量子隧穿效应模拟

10月7日，诺贝尔物理学奖公布，授予约翰・克拉克（John Clarke）、米歇尔・H・德沃雷（Michel H. Devoret）和约翰・M・马丁尼斯（John M. Martinis），以表彰他们“在电路中发现宏观量子力学隧穿效应与能量量子化现象”。

我在先前就已经使用Manim动画展现了量子隧穿效应，这里是理论的讲解和具体的模拟方法

量子隧穿效应

在量子力学中，粒子并不是简单地像经典物理那样沿着确定的轨迹运动，而是表现出波动性质。这种波动性质可以通过波函数来描述，波函数的平方 $|\psi(x,t)|^2$ 表示粒子在某个位置出现的概率密度。

其中，一维含时薛定谔方程 描述了粒子的时间演化过程：

$$i \frac{\partial \psi}{\partial t} = \left( -\frac{1}{2m} \frac{\partial^2}{\partial x^2} + V(x) \right) \psi(x,t)$$

这里，$\psi(x,t)$ 是波函数，$m$ 是粒子的质量，$V(x)$ 是势能。

量子隧穿：粒子穿越势垒

量子隧穿效应指的是，尽管粒子的能量 $E$ 小于某个势垒的高度 $V_0$，它依然有一定的概率穿越势垒。这是经典物理中无法解释的现象。

势垒可以被建模为一个有限的方势：

$$V(x) = \begin{cases} V_0 & \text{如果 } x \in [x_{\text{left}}, x_{\text{right}}] \\ 0 & \text{其他} \end{cases}$$

粒子在碰到势垒时，波函数会发生衰减，但并不是完全消失，存在一定的透过概率，称为隧穿概率。

波包的传播

量子波包的初始状态可以通过一个高斯包络与平面波的组合来描述：

$$\psi(x, 0) = A \cdot \exp \left( -\frac{(x - x_0)^2}{2 \sigma^2} \right) \cdot \exp(i k_0 x)$$

其中，$A$ 是振幅，$x_0$ 是波包中心，$\sigma$ 是波包的宽度，$k_0$ 是波数，表示粒子的动量。

分步傅里叶法

为了数值上求解波函数的时间演化，常用的数值方法之一是分步傅里叶法。其核心思想是将时间演化分为两个步骤：首先在实空间中考虑势能的影响，然后在动量空间中考虑动能的影响。具体公式如下：

$$\psi(t+\Delta t) \approx \exp\left(-i \frac{V(x)}{2} \Delta t\right) \cdot \mathcal{F}^{-1} \left[ \exp\left(-i \frac{k^2}{2m} \Delta t\right) \mathcal{F} \left( \exp\left(-i \frac{V(x)}{2} \Delta t\right) \psi(t) \right) \right]$$

其中，$\mathcal{F}$ 表示傅里叶变换，$k$ 是波数。

观察量子隧穿

在我们的模拟中，波包从左侧传入，并遇到一个有限高度的势垒，观察其在势垒前后的行为，尤其是隧穿效应的显现。通过数值方法模拟时间演化并绘制出波包的概率密度图 $|\psi(x,t)|^2$，可以清晰地看到当波包的能量低于势垒高度时，它依然有部分概率“穿透”势垒，从而展示量子力学中经典物理无法预测的现象。

模拟

采用Manim和Python脚本展示量子隧穿效应中粒子波是如何穿越障碍的

from manim import *
import numpy as np


# ===============================================================
# 1. 物理与数值：一维含时薛定谔方程 i dψ/dt = [-(1/2m)∂² + V(x)] ψ
#    单位化：ℏ = 1, m = 1
#    数值方法：分步傅里叶法（Split-Step Fourier, SSFM）
#    修改：模拟从无限远传播的波包
# ===============================================================


def generate_spatial_grid(num_points: int, domain_half_length: float) -> tuple[np.ndarray, float]:
    """
    生成空间网格与网格间距。
    - num_points: 空间采样点数（建议为 2 的幂）
    - domain_half_length: 空间范围的一半 L，使得 x ∈ [-L, L]
    返回:
      x: 形如 [-L, ..., L) 的等距网格
      dx: 网格间距
    """
    x = np.linspace(-domain_half_length, domain_half_length, num_points, endpoint=False)
    dx = x[1] - x[0]
    return x, dx


def generate_wavenumbers(num_points: int, dx: float) -> np.ndarray:
    """
    生成与 FFT 对应的波数数组 k（与 numpy.fft.fftfreq 一致的排列）。
    """
    # 2π 因子放入动能传播子中统一处理，这里只生成频率再乘以 2π
    freqs = np.fft.fftfreq(num_points, d=dx)
    k = 2.0 * np.pi * freqs
    return k


def build_potential_barrier(x: np.ndarray, barrier_left: float, barrier_right: float, barrier_height: float) -> np.ndarray:
    """
    有限方势垒：
      V(x) = barrier_height, x ∈ [barrier_left, barrier_right]
              0,          其他
    """
    V = np.zeros_like(x)
    mask = (x >= barrier_left) & (x <= barrier_right)
    V[mask] = barrier_height
    return V


def build_initial_plane_wave(x: np.ndarray, k0: float, amplitude: float = 1.0) -> np.ndarray:
    """
    平面波 ψ(x,0) = A * exp(i k0 x)，不进行空间归一化。
    模拟从无限远传播的平面波。
    """
    psi = amplitude * np.exp(1j * k0 * x)
    return psi


def build_initial_wave_packet_from_left(x: np.ndarray, k0: float, amplitude: float = 1.0, 
                                       packet_center: float = -30.0, packet_width: float = 3.0) -> np.ndarray:
    """
    创建一个从左侧传播进来的波包，模拟平面波从无限远传播的效果。
    使用高斯包络来限制初始波包的位置，让它从左侧逐渐传播进来。
    """
    # 高斯包络
    gaussian_envelope = np.exp(-((x - packet_center) ** 2) / (2.0 * (packet_width ** 2)))
    # 平面波相位
    plane_wave_phase = np.exp(1j * k0 * x)
    # 组合
    psi = amplitude * gaussian_envelope * plane_wave_phase
    return psi


def build_initial_gaussian_packet(x: np.ndarray, center: float, width: float, k0: float) -> np.ndarray:
    """
    高斯波包 ψ(x,0) ~ exp(-(x-x0)^2/(4σ^2)) * exp(i k0 x)，并进行 L2 归一化。
    """
    gaussian_envelope = np.exp(-((x - center) ** 2) / (4.0 * (width ** 2)))
    plane_wave_phase = np.exp(1j * k0 * x)
    psi = gaussian_envelope * plane_wave_phase
    # L2 归一化
    norm = np.sqrt(np.trapz(np.abs(psi) ** 2, x))
    psi /= norm
    return psi


def split_step_fourier_time_evolution(
    psi0: np.ndarray,
    potential: np.ndarray,
    time_step_sim: float,
    num_sim_steps: int,
    dx: float,
    k_array: np.ndarray,
    store_stride: int,
) -> tuple[np.ndarray, np.ndarray]:
    """
    分步傅里叶法时间推进：
      ψ(t+dt) ≈ e^{-iV dt/2} F^{-1}[ e^{-i k^2 dt / (2m)} F[e^{-iV dt/2} ψ(t)] ]
    此处取 ℏ=1, m=1。

    - psi0: 初始 ψ(x)
    - potential: V(x)
    - time_step_sim: 数值仿真时间步长 dt（较小以稳定）
    - num_sim_steps: 总仿真步数 Nsteps
    - dx: 空间网格间距
    - k_array: 与 FFT 匹配的波数数组
    - store_stride: 每隔多少个仿真步保存一次数据

    返回:
      times_saved: shape (Nsaved,)
      psi_saved: shape (Nsaved, N), 复数 数组（可据此计算 |ψ|²）
    """
    num_points = psi0.size

    # 传播子（常量预计算）
    exp_halfV = np.exp(-1j * potential * (time_step_sim / 2.0))
    expT = np.exp(-1j * (k_array ** 2) * (time_step_sim / 2.0))  # 先保存半步，合并两次FFT可用同一数组
    # 实际上一次完整步长需要 expT^2；为了减少额外乘法，我们每个半步都乘一次 expT

    psi = psi0.astype(np.complex128, copy=True)

    # 预估保存次数
    num_saved = (num_sim_steps // store_stride) + 1
    psi_saved = np.zeros((num_saved, num_points), dtype=np.complex128)
    times_saved = np.zeros(num_saved, dtype=np.float64)

    save_index = 0
    psi_saved[save_index, :] = psi
    times_saved[save_index] = 0.0
    save_index += 1

    for step in range(1, num_sim_steps + 1):
        # V 半步
        psi *= exp_halfV
        # 动能半步（频域）
        psi_k = np.fft.fft(psi)
        psi_k *= expT
        psi = np.fft.ifft(psi_k)
        # 再做一次动能半步合成整步（T 半步 + T 半步），避免额外中间变量
        psi_k = np.fft.fft(psi)
        psi_k *= expT
        psi = np.fft.ifft(psi_k)
        # V 半步
        psi *= exp_halfV

        if (step % store_stride) == 0:
            psi_saved[save_index, :] = psi
            times_saved[save_index] = step * time_step_sim
            save_index += 1

    # 若末尾没有正好落在 stride 上，数组尾部仍为零；裁剪到实际保存长度
    psi_saved = psi_saved[:save_index, :]
    times_saved = times_saved[:save_index]
    return times_saved, psi_saved


# ===============================================================
# 2. Manim 场景：模拟从无限远传播的波包穿越势垒
#    - 扩大空间域以容纳更长的传播距离
#    - 调整初始条件使波包从更远的地方开始
#    - 优化动画参数以适应更长的传播时间
# ===============================================================


class QuantumBarrierInfinity(Scene):
    def construct(self) -> None:
        # ---------- 数值参数（大幅扩大空间域） ----------
        num_points = 4096            # 大幅增加空间采样点以支持更大的空间域
        L = 100.0                   # 大幅扩大空间半长，x ∈ [-L, L)
        x, dx = generate_spatial_grid(num_points, L)
        k = generate_wavenumbers(num_points, dx)

        # 势垒参数（有限高势垒）
        barrier_left = -7.0
        barrier_right = 7.0
        barrier_height = 2.0        
        V = build_potential_barrier(x, barrier_left, barrier_right, barrier_height)

        # 波包参数：从左侧传播进来，模拟平面波从无限远传播的效果
        k0 = 1.5                    # 动量（波数），入射动能
        amplitude = 0.5             # 波包振幅
        packet_center = -30.0       # 初始波包中心位置
        packet_width = 3.0          # 波包宽度
        psi0 = build_initial_wave_packet_from_left(x, k0, amplitude, packet_center, packet_width)

        # 时间参数（调整以适应波从左侧传播）
        t_max = 20.0                # 增加总时间，让波有足够时间从左侧传播到势垒
        dt_sim = 1e-3               # 调整数值仿真步长
        total_steps = int(t_max / dt_sim)

        # 存储步距：约 60 FPS，按 stride 对齐
        target_fps = 60.0
        dt_store = 1.0 / target_fps
        store_stride = max(1, int(round(dt_store / dt_sim)))

        # ---------- 预计算时间演化 ----------
        times, psi_saved = split_step_fourier_time_evolution(
            psi0=psi0,
            potential=V,
            time_step_sim=dt_sim,
            num_sim_steps=total_steps,
            dx=dx,
            k_array=k,
            store_stride=store_stride,
        )
        density_saved = (np.abs(psi_saved) ** 2).astype(np.float64)

        # 为纵轴选择一个合适的上限
        density_max = float(np.max(density_saved))
        y_max = max(0.3, min(0.8, density_max * 1.3))  # 调整范围以适应新的参数

        # ---------- 构建坐标轴与图元 ----------
        # 只显示与势垒相互作用的区域，而不是整个计算域
        display_range = 25.0  # 显示范围：x ∈ [-15, 15]
        axes = Axes(
            x_range=[-display_range, display_range, 3.0],   # 调整刻度间隔
            y_range=[0.0, y_max, y_max / 8.0],
            x_length=12.0,          # 增加坐标轴长度
            y_length=5.3,           # 降低动画高度
            axis_config={
                "include_numbers": True, 
                "font_size": 24,
                "decimal_number_config": {"num_decimal_places": 2}
            },
            tips=False,
        ).to_edge(DOWN, buff=0.8)
        
        # 手动调整纵坐标数字位置，让它们向左移动
        for number in axes.y_axis.numbers:
            number.shift(LEFT * 0.3)

        title = Tex("Quantum Tunneling through Finite Potential Barrier").to_edge(UP, buff=0.5)

        # 势垒：用半透明矩形表示，高度与势垒高度成比例
        # 将势垒高度映射到显示范围，使其更直观
        barrier_display_height = barrier_height * 0.3  # 直接使用缩放因子
        barrier_polygon = Polygon(
            axes.c2p(barrier_left, 0.0, 0.0),
            axes.c2p(barrier_right, 0.0, 0.0),
            axes.c2p(barrier_right, barrier_display_height, 0.0),
            axes.c2p(barrier_left, barrier_display_height, 0.0),
            color=YELLOW,
            fill_opacity=0.3,
            stroke_opacity=0.8,
        )

        # 添加势垒标签
        barrier_label = Tex("Barrier", font_size=24).next_to(barrier_polygon, UP, buff=0.1)

        # 概率密度曲线：基于预计算结果，使用 always_redraw 按时间索引重建
        time_tracker = ValueTracker(0.0)

        x_values = x

        def build_density_graph() -> Mobject:
            current_t = time_tracker.get_value()
            max_index = len(times) - 1
            # times 大致等间隔：使用 dt_store 对齐索引
            dt_store_effective = times[1] - times[0] if len(times) > 1 else 1.0 / 60.0
            index = int(np.clip(round(current_t / dt_store_effective), 0, max_index))
            y_values = density_saved[index]
            
            # 只显示显示范围内的数据
            mask = (x >= -display_range) & (x <= display_range)
            x_display = x[mask]
            y_display = y_values[mask]
            
            # 使用 line graph 连接离散点
            graph = axes.plot_line_graph(
                x_values=x_display,
                y_values=y_display,
                add_vertex_dots=False,
                line_color=BLUE,
                stroke_width=3,
            )
            return graph

        density_graph = always_redraw(build_density_graph)

        time_label = always_redraw(
            lambda: Tex(f"t = {time_tracker.get_value():.2f} s", font_size=32).next_to(axes, DOWN, buff=0.3)
        )

        params_text = VGroup(
            Tex(f"Initial Energy: {k0**2/2:.2f}", font_size=20),
            Tex(f"Potential Barrier Height: {barrier_display_height:.2f}", font_size=20),
        ).arrange(DOWN, aligned_edge=LEFT).to_edge(LEFT, buff=2.5)

        # ---------- 动画播放 ----------
        self.add(title, axes, barrier_polygon, barrier_label, density_graph, time_label, params_text)
        
        # 添加初始说明
        initial_text = Tex("Wave Propagate from left", font_size=28).to_edge(RIGHT, buff=1.5)
        self.add(initial_text)
        
        # 播放动画
        self.play(
            time_tracker.animate.set_value(float(times[-1])), 
            run_time=float(times[-1]), 
            rate_func=linear
        )
 
        self.remove(initial_text)
        self.wait(2)

本文只展示了最基本的粒子波穿越障碍的模拟，诺贝奖表彰给宏观的粒子隧穿效应的一个具体的应用，研究基础是“约瑟夫森结”，一定条件下，电流可以通过绝缘体流动。整体来说是旨在授予量子力学的奇异性质可以在宏观尺度上得到体现，为量子计算的应用做出贡献