For today’s recreational coding exercise, we will investigate plasma physics with particle-in-cell (PIC) simulations. We will create a simulation of two electron beams passing through each other in opposite directions. An interesting phenomenon occurs, called the Two-Stream Instability.

对于当今的娱乐性编码练习，我们将使用单元中粒子(PIC)模拟研究等离子体物理。 我们将模拟两个彼此相反的方向穿过的电子束。 发生一个有趣的现象，称为两流不稳定性。

You may find the accompanying Python code on github.

您可以在github上找到随附的Python代码。

Before we begin, below is a gif of what running our simulation looks like:

在开始之前，下面是模拟运行情况的图像：

单元中粒子(PIC)方法 (Particle-In-Cell (PIC) Method)

We will consider a one-dimensional system of electrons in an unmagnetized uniform medium of ions. The electrons will be described by N particles, indexed by i, each having a

我们将考虑在未磁化的离子均匀介质中的一维电子系统。 电子将由N个粒子(由i索引)描述，每个粒子都有一个

• position rᵢ

  位置rᵢ

• velocity vᵢ

  速度vᵢ

The electrons feel an acceleration aᵢ which is due to the electric field E at the location of the particle. The equations of motion for the electrons are given by:

电子感到加速度aᵢ ，这是由于粒子位置处的电场E引起的。 电子的运动方程式为：

What remains to be done is to calculate the electric field E, which depends on the number density n of the electrons. The electric field E is defined as the negative gradient of the electric potential ϕ, which is given by Poisson’s equation sourced by the density:

剩下要做的是计算电场E ，它取决于电子的数密度n 。 的电场E被定义为电势φ，它是由泊松方程由密度给定源的负梯度：

where n₀ is the average electron density (a constant).

式中n₀是平均电子密度(一个常数)。

The Particle-In-Cell (PIC) Method computes the density and subsequently the potential and electric field on a mesh. Hence it is called a particle-mesh method. The electric field is then extrapolated onto the location of the particles in order to advance them. That’s a long sequence of steps. Let’s dig in!

单元中粒子(PIC)方法可计算密度，然后计算网格上的电势和电场。 因此，它被称为粒子网格方法。 然后将电场外推到粒子的位置，以使粒子前进。 这是一个很长的步骤。 让我们开始吧！

密度计算 (Density Calculation)

The density is obtained by binning the particles onto the mesh gridpoints. Call the separation distance between mesh points Δx. A particle with position rᵢ will lie between a pair of coordinates xⱼ and x_{j+1} on the mesh. Then, the densities nⱼ and n_{j+1} at these two points gain a contribution:

通过将粒子合并到网格网格点上来获得密度。 称网格点之间的间隔距离Δx。 位置为rᵢ的粒子将位于网格上的一对坐标xⱼ和x_ {j + 1}之间。 然后，在这两个点上的密度nⱼ和n_ {j + 1}获得贡献：

After collecting contribution from all particles, the density has to be properly normalized so that the average density is the specified n₀.

在收集所有粒子的贡献后，必须适当地对密度进行归一化，以使平均密度为指定的n₀。

潜力计算 (Potential Calculation)

Knowing the density nⱼ on the mesh grid points, we now seek the potential ϕⱼ at these points. The second-derivative operator (called the ‘Laplacian’) in Poisson’s equation can be approximated on the mesh with a finite-different formula:

知道在网格点的密度nⱼ，我们现在寻求在这些点的潜力φⱼ。 泊松方程中的二阶导数算子(称为“拉普拉斯算子”)可以使用有限差分公式在网格上近似：

If all the values for nⱼ are arranged into a column vector, and similarly for ϕⱼ, then the equation can be written as a system of linear equations in matrix form:

如果所有用于nⱼ的值被设置成一个列向量，并且类似地对于φⱼ，则方程可以写为矩阵形式的线性方程系统：

The matrix represents the Laplacian operator and is sparse and tridiagonal. We have assumed periodic boundary conditions, i.e., the endpoints of the mesh wrap around to each other, which add the ‘1’s in the top-right and bottom-left corners of the matrix. This matrix equation can be efficiently numerically inverted and solved by using pre-built linear algebra packages such as the Python Scipy package’s “spsolve” function.

矩阵表示拉普拉斯算子，并且是稀疏的和三对角的。 我们假定了周期性边界条件，即网格的端点彼此环绕，这在矩阵的右上角和左下角加了“ 1”。 通过使用预构建的线性代数包(例如Python Scipy包的“ spsolve”函数)，可以有效地对这个矩阵方程进行数值反转和求解。

To start, we will construct the matrix for now with some code:

首先，我们将使用一些代码来构建矩阵：

# Construct matrix Lmtx to computer Laplacian (2nd derivative)
# on mesh size Nx, separation dx
e = np.ones(Nx)
diags = np.array([-1,0,1])
vals  = np.vstack((e,-2*e,e))
Lmtx = sp.spdiags(vals, diags, Nx, Nx);
Lmtx = sp.lil_matrix(Lmtx) # tranform mtx type to modify entries
Lmtx[0,Nx-1] = 1
Lmtx[Nx-1,0] = 1
Lmtx /= dx**2
Lmtx = sp.csr_matrix(Lmtx) # transform mtx type

电场计算(Electric Field Calculation)

Given that we have solved for ϕⱼ on the mesh, we can next solve for the electric field Eⱼ by discretizing the first-derivative operator (called the ‘gradient’):

鉴于我们解决了为φⱼ网格上，我们可以下一个求解由离散一阶导数算子(称为“梯度”)的电场Eⱼ：

Again we can write out the system of equations in matrix form:

同样，我们可以用矩阵形式写出方程组：

The system is simpler this time. There is no need to invert a matrix. The electric field is directly obtained by applying the matrix multiplication.

这次系统更简单。 无需反转矩阵。 通过应用矩阵乘法直接获得电场。

The code to construct the matrix for the gradient operator is shown below:

下面显示了构造用于梯度算子的矩阵的代码：

# Construct matrix G to computer Gradient  (1st derivative)
# on mesh size Nx, separation dx
e = np.ones(Nx)
diags = np.array([-1,1])
vals  = np.vstack((-e,e))
Gmtx = sp.spdiags(vals, diags, Nx, Nx);
Gmtx = sp.lil_matrix(Gmtx) # tranform mtx type to modify entries
Gmtx[0,Nx-1] = -1
Gmtx[Nx-1,0] = 1
Gmtx /= (2*dx)
Gmtx = sp.csr_matrix(Gmtx) # transform mtx type

加速度计算(Acceleration Calculation)

So far we have gone from particles to mesh values of density to potential to electric field. To complete the loop, we need the electric field at the location of the particles. Fortunately, this can now be done with interpolation from the 2 nearest gridpoints xⱼ and x_{j+1}. For particle i, the electric field Eᵢ is:

到目前为止，我们已经从粒子到密度的网格值，再到电场的电势。 为了完成循环，我们需要在粒子位置处施加电场。 幸运的是，现在可以通过从2个最近的网格点xⱼ和x_ {j + 1}进行插值来完成此操作。 对于粒子i，电场Eᵢ为：

using similar weights as when be binned to get density.

使用与合并时类似的权重以获取密度。

Putting all the steps above together, we can write a function to calculate the acceleration on each particle, given the particle distribution:

将以上所有步骤放在一起，我们可以编写一个函数来计算给定粒子分布的每个粒子的加速度：

def getAcc( pos, Nx, boxsize, n0, Gmtx, Lmtx ):"""Calculate the acceleration on each particle due to electric fieldpos      is an Nx1 matrix of particle positionsNx       is the number of mesh cellsboxsize  is the domain [0,boxsize]n0       is the electron number densityGmtx     is an Nx x Nx matrix for calculating the gradient on the gridLmtx     is an Nx x Nx matrix for calculating the laplacian on the grida        is an Nx1 matrix of accelerations"""# Calculate Electron Number Density on the Mesh by # placing particles into the 2 nearest bins (j & j+1, with proper weights)# and normalizingN          = pos.shape[0]dx         = boxsize / Nxj          = np.floor(pos/dx).astype(int)jp1        = j+1weight_j   = ( jp1*dx - pos  )/dxweight_jp1 = ( pos    - j*dx )/dxjp1        = np.mod(jp1, Nx)   # periodic BCn  = np.bincount(j[:,0],   weights=weight_j[:,0],   minlength=Nx);n += np.bincount(jp1[:,0], weights=weight_jp1[:,0], minlength=Nx);n *= n0 * boxsize / N / dx # Solve Poisson's Equation: laplacian(phi) = n-n0phi_grid = spsolve(Lmtx, n-n0)# Apply Derivative to get the Electric fieldE_grid = - Gmtx @ phi_grid# Interpolate grid value onto particle locationsE = weight_j * E_grid[j] + weight_jp1 * E_grid[jp1]a = -Ereturn a

时间整合(Time Integration)

We are in the home stretch now! We will loop over timesteps and evolve the electrons using a leap-frog scheme (‘kick-drift-kick’). For each timestep Δt, each particle receives a half-step ‘kick’:

我们现在在家中！ 我们将在时间步长上循环，并使用跳越方案('kick-drift-kick')释放电子。 对于每个时间步长Δt ，每个粒子都接收一个半步“踢”：

followed by a full-step ‘drift’:

然后是整步“漂移”：

followed by another half-step ‘kick’.

接下来是半步“踢”。

The evolution is performed in the code using a For-loop and our function for the acceleration from earlier:

使用For循环在代码中执行演化，而我们的函数则用于从早期开始进行加速：

# Simulation Main Loop
for i in range(Nt):# (1/2) kickvel += acc * dt/2.0# drift (and apply periodic boundary conditions)pos += vel * dtpos = np.mod(pos, boxsize)# update accelerationsacc = getAcc( pos, Nx, boxsize, n0, Gmtx, Lmtx )# (1/2) kickvel += acc * dt/2.0# update timet += dt

In case this looks familiar, it might be because this time-stepping algorithm is a common one for particle-based simulations, and we have also used it in our tutorial on gravitational N-body simulations.

如果看起来很熟悉，可能是因为这种时步算法是基于粒子的模拟中的常用算法，并且我们还在引力N体模拟的教程中也使用了它。

初始条件 (Initial Conditions)

The only thing left to do is specify the initial conditions (positions and velocities at time t=0) of the electrons in the simulation. In the code, we implement the example that we have a uniform beam traveling to the right, and another one superimposed, traveling left. Particle positions are drawn from a random uniform distribution. Velocities are drawn from a Gaussian centered about the beam velocity (one positive, one negative). Small initial perturbations in the form of a cosine function are added to the velocity distribution. Feel free to check out the code and modify or the initial conditions or explore your own.

剩下要做的就是指定模拟中电子的初始条件(在时间t = 0时的位置和速度)。 在代码中，我们实现了一个示例，我们有一个均匀的光束向右传播，另一个光束向左传播。 粒子位置是从随机均匀分布中得出的。 速度是从以光束速度为中心的高斯中得出的(一个正，一个负)。 余弦函数形式的小的初始扰动被添加到速度分布中。 随时检查代码并进行修改或初始条件，或者自行探索。

Running the code allows you to visualize the simulation in real time and will yield the figure:

运行代码可以使您实时可视化仿真，并产生图：

We plot the configuration of electrons in what’s called phase-space, where the x-axis is the particle’s position and the y-axis is the velocity. This allows us to visualize all the variables in the system directly. Particles initially belonging to the right-moving beam have a blue color, while particles in the left-moving beam are red. The two beams interact with each other and go unstable in mesmerizing fashion.

我们在所谓的相空间中绘制电子的构图，其中x轴是粒子的位置，y轴是速度。 这使我们可以直接可视化系统中的所有变量。 最初属于向右移动光束的粒子为蓝色，而向左移动的光束中的粒子为红色。 两束光相互影响，并以令人着迷的方式变得不稳定。

So there you have it. We’ve created a plasma PIC simulation from scratch. I have here assumed the simplest scenario where electrons are non-relativistic and ions form a static background, and we’ve already learned an interesting instability can occur from counter-streaming plasma beams. Astro- and plasma physicists use these types of simulations to study all sorts of physical systems, from magnetic reconnection in the Sun to the design of tokamaks. In tokamak design, for example, the goal is to figure out a way to confine plasma to produce energy, without triggering disruptive plasma instabilities. It’s a very hard problem, but some day the technology may be an efficient way to produce clean energy.

所以你有它。 我们从头开始创建了等离子PIC仿真。 我在这里假设了最简单的情况，其中电子是非相对论的，并且离子形成了静态背景，而且我们已经了解到，逆流等离子体束会产生有趣的不稳定性。 天体物理学家和等离子物理学家使用这些类型的模拟来研究各种物理系统，从太阳的磁重联到托卡马克的设计。 例如，在托卡马克设计中，目标是找到一种方法来限制等离子体产生能量，而不会触发破坏性的等离子体不稳定性。 这是一个非常棘手的问题，但是总有一天该技术可能是生产清洁能源的有效方法。

演示地址

Download the Python code on github for our PIC algorithm to visualize the simulation in real time and play around with the initial conditions. Enjoy!

在github上为我们的PIC算法下载Python代码，以实时可视化仿真并处理初始条件。 请享用！