QED Processes in Particle-In-Cell Simulations

tutorial
code
physics
Author

Shan Jin

Published

April 26, 2023

This post explains how to embed nonlinear Quantum Electrondynamics (QED) processes into the Particle-In-Cell code under local constant field approximation (LCFA).

量子电动力学过程在粒子模拟中的实现

本篇主要讲述如何将 QED 的散射截面，通过蒙特卡洛的方法嵌入到 PIC 模拟程序当中。不涉及散射截面的具体推导，且仅限于非线性 QED 的一阶过程。

散射截面

非线性 QED 有两个最基本的一阶过程：非线性 Compton 散射，和 非线性 Breit-Wheeler 过程。非线性 Compton 散射指的是一个高能电子 \(e^{-}\) 和 \(n\) 个激光光子 \(\gamma_l\) 发生散射，产生一个高能的 \(\gamma\) 光子，同时散射电子的能量降低。非线性 Breit-Wheeler 过程指的是一个高能的 \(\gamma\) 光子，和 \(n\) 个激光光子 \(\gamma_l\) 发生散射，产生一个正负电子对 (\(e^{-},e^{+}\)). 我们可以用公式表达这两个过程：

\[ e^{-} + n\gamma_l\to e^{-}+\gamma~~~~~~~~\gamma+n\gamma_l\to e^{-}+e^{+} \]

非线性效应体现在公式中的 \(n\), 当 \(n=1\) 时便是线性过程。

我们知道在粒子模拟程序中，电子的信息是以“粒子”这一种数据结构来储存的，而激光场 \(\gamma_l\) 则是通过在模拟盒内划分网格来储存信息的。那么 \(\gamma\) 光子应该采用哪种方式呢？一般飞秒激光的光子的能量在 1.2eV 左右，而 \(\gamma\) 光子 能量通常在 MeV 甚至 GeV 量级，它的波长要比激光波长小 \(10^6\sim 10^9\) 量级， 要在模拟网格上分辨 \(\gamma\) 光子几乎是不可能的。所以在 PIC 中 \(\gamma\) 光子以粒子的形式储存信息，它拥有粒子的位置、动量等信息，但没有类似于“振幅”的信息。

为了定量地描述这两个过程，需要先定义几个非线性参数。一个是我们最熟悉的无量纲化的激光振幅（等离子体文献中一般记为 \(a_0\), 在 QED 理论文献中一般记为 \(\xi\) ）： \[ a_0=\xi=\frac{eE}{m_e\omega_lc} \] 其中，\(E\) 是激光振幅，\(\omega_l\) 是激光的圆频率。当 \(\xi>1\) 时，多光子（非线性）效应开始变得显著，因而 \(\xi\) 被称为经典非线性参数。

另一个参数是量子非线性参数 \(\chi_e\) （QED模拟文献中会将其记为 \(\eta\)，也有少部分文献中记为\(\chi\)），它的定义是 \[ \chi_e=\left(e \hbar / m_e^3 c^4\right)\left|F_{\mu \nu} p^\nu\right| =\gamma\sqrt{(\vec{E}+\vec{v}\times\vec{B})^2-(\vec{E}\cdot\vec{v}/c)^2}/E_{cr} \] 其中，\(\gamma\) 指的电子的洛伦兹因子，注意与光子的记号区分，\(E_{cr}=m_e^2c^3/e\hbar\) 是 Schwinger 临界场，上式的物理意义是电子在其静止坐标系下感受到的激光场强与 Schwinger 临界场之比。 当 \(\chi_e>1\) 时 QED 效应开始变得显著。

\(\chi_e\) 描述的是非线性 Compton 散射的强度，与之类似的，描述非线性 Breit-Wheeler 过程时，我们采用另一个量子非线性参数 \(\chi_\gamma\) （QED模拟文献中会将其记为 \(\chi\)，也有少部分文献中记为\(\eta\)，注意与 \(\chi_e\) 的记法区分），它的定义是 \[ \chi_\gamma=\left(e \hbar^2 / 2 m_e^3 c^4\right)\left|F^{\mu \nu} k_\nu\right| =\frac{1}{2}\frac{\hbar\omega_\gamma}{m_ec^2}\sqrt{(\vec{E}+\frac{c\vec{k}}{\omega_\gamma}\times c\vec{B})^2-(\vec{E}\cdot\frac{c\vec{k}}{\omega_\gamma})^2}/E_{cr} \] 式中的 \(k_\nu\) 是激光光子的四矢量。

我们接下去的所有讨论都基于以下这两个假设 ：

  1. 静态场近似：通常被称为局域常数场近似 (LCFA)，它指的是当激光满足 \(\xi\gg1\) 时，光子的形成长度（即 QED 过程的空间尺度）近似为 \(\lambda_l/\xi\)， 远远小于激光的波长 \(\lambda_l\). 在此情况下 QED 过程可以认为是局域发生的，且在 QED 过程发生时，电子（或 \(\gamma\) 光子）感受到的电磁场为恒定交叉场（电场和磁场相互垂直、强度不变）。

  2. 弱场近似1：（实验室坐标系下）相互作用的电磁场要远小于 Schwinger 临界场。在此条件下，QED 截面只与 \(\chi_e\) 或 \(\chi_\gamma\) 相关。

在 弱场近似 下，非线性 Compton 散射的散射截面可以写为

\[ \frac{\mathrm{d}^2 W}{\mathrm{d} t \mathrm{d} \omega}=\frac{\alpha}{\sqrt{3} \pi \tau_c \gamma}\left\{\left[2+\frac{\omega^2}{(1-\omega)}\right] \mathrm{K}_{2 / 3}\left[\frac{2 \omega}{3 \chi_e(1-\omega)}\right] -\int_{2 \omega /[3 \chi_e(1-\omega)]}^{\infty} \mathrm{d} y \mathrm{~K}_{1 / 3}(y)\right\} \tag{1}\]

其中 \(\alpha=1/137\) 是精细结构常数，\(\tau_{c}=\hbar /\left(m c^{2}\right)\) 是康普顿时间，\(\gamma\) 是电子的洛伦兹因子，\(\omega\) 是辐射光子能量 (\(\epsilon_\gamma=\hbar\omega_\gamma\)) 与初态电子能量 (\(\epsilon=\gamma m_ec^2\)) 之比，即 \(\omega=\varepsilon_{\gamma} / \varepsilon\)，\(\mathrm{K}_{\nu}(x)\) 是第二类改良贝塞尔函数。上式的物理意义是电子在单位时间内辐射一个能量 为 \(\omega\epsilon\) 的光子的概率。

\(\omega\) 和 \(\chi_e,\chi_\gamma\) 之间有个重要的关系式 \[ \omega=\frac{2\chi_\gamma}{\chi_e} \] 利用此关系式和贝塞尔函数的递推关系 \[ \mathrm{~K}_{5 / 3}(x)=-\mathrm{~K}_{1 / 3}(x)-2\frac{\mathrm{d}\mathrm{~K}_{2 / 3}(x)}{\mathrm{d}x} \] 可以将 Equation 1 转换为文献(Ridgers et al. 2014) 中的式 (1). 通过艾里函数 \(\mathrm{Ai}(x)\) 和贝塞尔函数的关系 \[ \mathrm{Ai}(x)=\frac{\sqrt{x}}{\sqrt{3}\pi}\mathrm{~K}_{1 / 3}\left(\frac{2}{3}x^{3/2}\right) \]

\[ \mathrm{Ai}'(x)=-\frac{x}{\sqrt{3}\pi}\mathrm{~K}_{2 / 3}\left(\frac{2}{3}x^{3/2}\right) \] 还可以将 Equation 1 转换为文献 (Seipt and King 2020) 中的式 (42).

下图展示了 \(\chi_e=1,5,10\) 时，非线性 Compton 散射截面 (Equation 1) 随 \(\gamma\) 光子的相对能量 \(\omega\) 的变化趋势。可以看到 \(\chi_e\) 越大，反应概率越大。

点击查看代码 
import matplotlib.pyplot as plt
from scipy.integrate import quad
import numpy as np
from scipy.special import kv

def odWdw(w,eta):
    term1 = ( 2.0 + w**2/(1.0-w) ) * kv(2/3,2.0*w/3.0/eta/(1.0-w))
    term2 = -quad(lambda x: kv(1/3,x),2.0*w/3.0/eta/(1.0-w),np.inf,limit=1000,limlst=1000)[0]
    return term1+term2

OMEGA = np.linspace(0.01,0.99,200)
dWdtdw_eta1 = np.array([odWdw(w,1) for w in OMEGA])
dWdtdw_eta5 = np.array([odWdw(w,5) for w in OMEGA])
dWdtdw_eta10 = np.array([odWdw(w,10) for w in OMEGA])
plt.plot(OMEGA,dWdtdw_eta1,label=r'$\chi_e=1$')
plt.plot(OMEGA,dWdtdw_eta5,label=r'$\chi_e=5$')
plt.plot(OMEGA,dWdtdw_eta10,label=r'$\chi_e=10$')
plt.xlabel(r'$\omega$')
plt.ylabel(r'$d^2W/d\omega dt$')
plt.legend()
plt.yscale('log')
plt.ylim([1.0e-5,1.0e3])
plt.show()
Figure 1: Nonlinear Compton scattering sector

对 Equation 1 中的末态光子能量 \(\omega\) 进行积分，可以得到总的散射概率 \[ \frac{\mathrm{d}W}{\mathrm{d}t}=\int_0^1\frac{\mathrm{d}^2 W}{\mathrm{d} t \mathrm{d} \omega}\mathrm{d}\omega \tag{2}\]

对于非线性 Breit-Wheeler 过程，我们令 \(\omega\) 为末态负电子的能量和初态 \(\gamma\) 光子的能量之比，散射截面公式可以表示为 \[ \frac{\mathrm{d}^2 W}{\mathrm{d} t \mathrm{d} \omega}=\frac{\alpha m_ec^2}{\sqrt{3} \pi \tau_c \hbar\omega_\gamma} \left\{ \left[ \frac{1}{\omega(1-\omega)}-2 \right] \mathrm{K}_{2 / 3} \left[ \frac{2 }{3 \chi_\gamma\omega(1-\omega)} \right] +\int_{2/[3 \chi_\gamma\omega(1-\omega)]}^{\infty} \mathrm{d} y \mathrm{~K}_{1 / 3}(y) \right\} \tag{3}\] 同样的，利用艾里函数可以将 Equation 3 转换为文献 (Seipt and King 2020) 中的式 (61). 对 Equation 3 中的 \(\omega\) 积分之后，得到总的散射概率，形式和 Equation 2 相同。

蒙卡方法

粒子模拟程序中 QED 模块的蒙卡整体分为两部分，第一部分判断电子何时辐射（或者 \(\gamma\) 光子何时分裂），第二部分判断辐射光子的能量（或者分裂的正负电子对的能量）。下面以非线性 Compton 散射为例来讲述这个算法，非线性 Breit-Wheeler 过程可以类似地得到。

我们已知非线性 Compton 散射的概率为 \(\frac{\mathrm{d}W}{\mathrm{d}t}\)，由 Equation 2 得到，看似改式不含时间，实际上激光场中的电子能量是随时变化的，而 Equation 2 依赖于电子能量，因而 \(\frac{\mathrm{d}W}{\mathrm{d}t}\) 也是随时变化的。我们知道 \(\frac{\mathrm{d}W}{\mathrm{d}t}\) 的物理含义是电子在单位时间内辐射光子的概率，因此可以通过积分式来表示在一给定时间区间 \([0,t_0)\) 内电子发生辐射的总概率： \[ W(t_0)=\int_0^{t_0}\frac{\mathrm{d}W}{\mathrm{d}t}\mathrm{d}t \] 以此可以求得在时间区间 \([0,t_0)\) 内，电子在单位时间发生辐射的平均概率： \[ W(t_0)/t_0 \] 在选定 \(t_0\) 后，上式是不依赖于时间的。那么，我们可以求得在 \([0,t_0)\) 内，电子不发生任何辐射的概率为 \[ \lim_{\Delta t\to 0} \left( 1-\frac{W(t_0)\Delta t}{t_0} \right)^{\frac{t_0}{\Delta t}} =\lim_{\Delta t\to 0}\exp\left( \frac{t_0}{\Delta t}\ln \left( 1-\frac{W(t_0)\Delta t}{t_0} \right) \right)=\exp(-W(t_0)) \] 那么，在 \([0,t_0)\) 内发生辐射的概率为 \[ P(t_0)=1-e^{-W(t_0)} \tag{4}\] 蒙卡的实现是这样的：先产生一个0到1之间的随机数 \(P\)，通过 Equation 4 反解出辐射阈值概率 \(W(t_0)\)，再通过 \(W(t)=\int_0^{t}\frac{\mathrm{d}W}{\mathrm{d}t}\mathrm{d}t\) 在模拟中的每个时间步长计算 \(\frac{\mathrm{d}W}{\mathrm{d}t}\) 累加得到累计概率 \(W(t)\)，直到某个时间步长 \(W(t)\) 大于 \(W(t_0)\)，判断为电子在此时辐射出一个光子。此时还需要再进行一次蒙卡来确定辐射的光子的能量，利用截面公式 Equation 3 ，可以得到如下蒙卡过程：首先产生一个0到1之间的随机数 \(P_\omega\)，通过 \[ P_\omega = \frac{\int_0^\omega\frac{\mathrm{d}^2 W}{\mathrm{d} t \mathrm{d} \omega}\mathrm{d}\omega}{\int_0^1\frac{\mathrm{d}^2 W}{\mathrm{d} t \mathrm{d} \omega}\mathrm{d}\omega} \tag{5}\] 反解出 \(\omega\). 由 \(\epsilon_\gamma=\omega\epsilon\) 得到辐射光子能量。实际代码中一般会把 Equation 5 制成关于 \(\chi_e\) 和 \(\omega\) 的二维表格，在已知 \(P_\omega\) 和 \(\chi_e\) 的情况下通过线性插值得到 \(\omega\). 下面一段代码演示了这一过程，黑线是由 Equation 3 得到的解析结果，红线是蒙卡采样后的统计结果，可以看到蒙卡结果近乎完美地还原了散射截面。

点击查看代码 
import random
from math import floor

def find_idx_from_table(x_in, nx, x):
    x_value = x_in

    i1 = 0
    i2 = nx-1
    xdif1 = x[i1] - x_value
    xdif2 = x[i2] - x_value
    if(xdif1 * xdif2 < 0):
#        Use bisection to find the nearest cell
      while True:
        im = floor( (i1 + i2) / 2 )
        xdifm = x[im] - x_value
        if(xdif1 * xdifm < 0):
            i2 = im
        else:
            i1 = im
            xdif1 = xdifm
        if(i2 - i1 == 1):
            break
#        Interpolate in x
      fx = (x_value - x[i1]) / (x[i2] - x[i1])
    else:
#         if(warning and rank == 0):
#             print('*** WARNING ***')
#             print('Argument to "find_idx_from_table" outside the range of the table.')
#             print('Using truncated value. No more warnings will be issued.')
#             warning = False

        if(xdif1 >= 0):
            fx = 0.0
        else:
            fx = 1.0
    return [fx,i1,i2]

def find_value_from_table_1d(x_in, nx, x):
    fx,i1,i2 = find_idx_from_table(x_in, nx, x[:,0])
    return (1.0 - fx) * x[i1,1] + fx * x[i2,1]

def find_value_from_table_2d(x_in, p_value, nx,ny, x,y, values):
    value_interp = np.array([0.0 for i in range(3)])
    fx, ix1, ix2 = find_idx_from_table(x_in,nx,x)
    fy, iy1, iy2 = find_idx_from_table(p_value,ny,values[ix1,:])
    value_interp[0] = (1.0 - fy) * y[iy1] + fy * y[iy2]
    fy, iy1, iy2 = find_idx_from_table(p_value,ny,values[ix2,:])
    value_interp[1] = (1.0 - fy) * y[iy1] + fy * y[iy2]
    value_interp[2] = (1.0 - fx) * value_interp[0] + fx * value_interp[1]
    return value_interp[2]

drdw = np.loadtxt('./c0.table')
print(drdw.shape)
ETA = np.linspace(-3,2,100)
OMEGA = np.linspace(-4,-0.00001,500)
eta_log = 1.
p_values = np.array([random.random() for i in range(100000)])
w_check = np.array([0.0 for i in range(100000)])
for i in range(100000):
    w_check[i] = find_value_from_table_2d(eta_log, p_values[i], 100, 500, ETA, OMEGA, drdw)
H, edges = np.histogram(10**w_check,bins=np.linspace(0,1.0,100))

drdw_ana = np.array([odWdw(w,10**eta_log) for w in 10**OMEGA])

plt.plot(10**OMEGA,drdw_ana,'k-',label='analytical')
plt.plot(edges[:-1],H/100*np.pi,'r--',label='MC sampling')
plt.yscale('log')
plt.xlabel(r'$\omega$')
plt.ylabel(r'$d^2W/d\omega dt$')
plt.legend()
plt.show()
(100, 500)
Figure 2: Analytical results v.s. MC results

非线性 Breit-Wheeler 过程的操作与非线性 Compton 散射完全一致，只是换了散射截面，这里不再赘述。

Explore more:

References

Ridgers, C. P., J. G. Kirk, R. Duclous, et al. 2014. “Modelling Gamma-Ray Photon Emission and Pair Production in High-Intensity Laser–Matter Interactions.” Journal of Computational Physics 260: 273–85. https://doi.org/https://doi.org/10.1016/j.jcp.2013.12.007.
Seipt, D., and B. King. 2020. “Spin- and Polarization-Dependent Locally-Constant-Field-Approximation Rates for Nonlinear Compton and Breit-Wheeler Processes.” Phys. Rev. A 102 (November): 052805. https://doi.org/10.1103/PhysRevA.102.052805.

Footnotes

  1. 这里的所谓“弱场”是和 Schwinger 临界场（相当于 \(\xi_{cr}\sim 10^5\)）相比，现在（以及可预测的将来）实验室下的激光都满足“弱场”条件。同时，“弱场”不代表不会发生 QED 效应，因为在电子或 \(\gamma\) 光子静止坐标系下场强很容易达到 Schwinger 条件。即 \(\xi/\xi_{cr}\ll1\) 和 \(\chi_e>1\) 并不矛盾。↩︎