# JefiPIC: A 3-D Full Electromagnetic Particle-in-Cell Simulator Based on
Jefimenko’s Equations on GPU

Jian-Nan Chen1,2, Jun-Jie Zhang1,2,∗, Xue-Ming Li3, Hai-Liang Qiao1,2 and
Yong-Tao Zhao4 1\. Northwest Institute of Nuclear Technology, Xi’an, 710024,
China  2\. State Key Laboratory of Intense Pulsed Radiation Simulation and
Effect, 710024, China  3\. Institute of Applied Physics and Computational
Mathematics, Beijing, 100094, China  4\. MOE Key Laboratory for Nonequilibrium
Synthesis and Modulation of Condensed Matter, School of Science, Xi’an
Jiaotong University, Xi’an, 710049, China<EMAIL_ADDRESS>
###### Abstract

This paper presents a novel 3-D full electromagnetic particle-in-cell (PIC)
code called JefiPIC, which uses Jefimenko’s equations as the electromagnetic
(EM) field solver through a full-space integration method. Leveraging the
power of state-of-the-art graphic processing units (GPUs), we have made the
challenging integral task of PIC simulations achievable. Our proposed code
offers several advantages by utilizing the integral method. Firstly, it offers
a natural solution for modeling non-neutral plasmas without the need for pre-
processing such as solving Poisson’s equation. Secondly, it eliminates the
requirement for designing elaborate boundary layers to absorb fields and
particles. Thirdly, it maintains the stability of the plasma simulation
regardless of the time step chosen. Lastly, it does not require strict charge-
conservation particle-to-grid apportionment techniques or electric field
divergence amendment algorithms, which are commonly used in finite-difference
time-domain (FDTD)-based PIC simulations. To validate the accuracy and
advantages of our code, we compared the evolutions of particles and fields in
different plasma systems simulated by three other codes. Our results
demonstrate that the combination of Jefimenko’s equations and the PIC method
can produce accurate particle distributions and EM fields in open-boundary
plasma systems. Additionally, our code is able to accomplish these
computations within an acceptable execution time. This study highlights the
effectiveness and efficiency of JefiPIC, showing its potential for advancing
plasma simulations.

Keywords: particle-in-cell, Jefimenko’s equation, plasma simulation, integral
method, GPU

Our code is published at:

1\. https://github.com/Juenjie/JefiPIC

2\. https://codeocean.com/capsule/5129018/tree

##  1 Introduction

###  1.1 Background

The particle-in-cell (PIC) method is a commonly used tool in plasma physics.
It utilizes macro-particles to describe charged particles in similar phase
space states and model the evolution of the particle distribution and
electromagnetic field. This method was initially proposed by Dawson in the
1960s [1] to study the Langmuir wave in 1-D electrostatic plasma. Later,
Langdon and Birdsall improved the PIC model by incorporating finite-size
particles [2] or particle clouds [3], which solved the issue of Coulomb
collision between particles. Marder [4] and Villasenor [5] addressed the
electric field divergence error in the current-driven method. Currently, PIC
is used across various fields, including simulating controlled/laser
thermonuclear fusion [6],[7], studying nuclear explosions [8]-[10] and space
physics effects [11], and designing vacuum electronic devices [12]-[14].

The PIC method offers an intuitive representation of charged particles, making
it easier for researchers to analyze simulative phenomena and data. As a
result, numerous commercial software and open-source PIC codes have emerged
for decades, including but not limited to MAGIC [15], XOOPIC [16], ICEPIC
[17], UNIPIC [18], and EPOCH [7],[19].

The PIC model comprises two primary components. The first part involves
updating the dynamics of EM fields according to Maxwell’s equations,

$\displaystyle\varepsilon\frac{\partial\mathbf{E}}{\partial t}$
$\displaystyle=$ $\displaystyle\nabla\times\frac{\mathbf{B}}{\mu}-\mathbf{J},$
(1) $\displaystyle\frac{\partial\mathbf{B}}{\partial t}$ $\displaystyle=$
$\displaystyle-\nabla\times\frac{\mathbf{E}}{\mu},$ (2)
$\displaystyle\nabla\cdot\mathbf{E}$ $\displaystyle=$
$\displaystyle\frac{\rho}{\varepsilon},$ (3)
$\displaystyle\nabla\cdot\mathbf{B}$ $\displaystyle=$ $\displaystyle 0,$ (4)

where $\mathbf{E}$ and $\mathbf{B}$ represent electric and magnetic fields,
$\varepsilon$ and $\mu$ donate the permittivity and permeability of the
medium, $\rho$ and $\mathbf{J}$ are the electric density and electric current
density.

The second part involves updating the positions and velocities of macro-
particles using the Newton-Lorentz force equations of motion,

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}(\gamma m\mathbf{v})$
$\displaystyle=$ $\displaystyle q(\mathbf{E}+\mathbf{v}\times\mathbf{B}),$ (5)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{r}$ $\displaystyle=$
$\displaystyle\mathbf{v},$ (6)

where $\gamma$, $m$, $q$, $\mathbf{r}$, and $\mathbf{v}$ represent the
relativistic factor, particle’s mass, charge, displacement, and velocity,
respectively.

After ensuring current continuity,

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot\mathbf{J}$
$\displaystyle=$ $\displaystyle 0,$ (7)

the electric field divergence Eq. (3) is implicitly obtained solely through
the first two curl Eqs. (1) and (2), while Eq. (4) holds at all time.

###  1.2 Motivation

The finite-difference time-domain (FDTD) method [20] is well-suited for
computing curl Eqs. (1) to (2) due to its simplicity and data sparsity. This
makes it more advantageous than other methods, such as the finite element
method (FEM) [21],[22] that necessitates the design of complex basic
functions. Thus, FDTD has been employed as the EM solver for several decades.
In the FDTD-based PIC framework, it is crucial that the finite-difference
implementation of the current continuity Eq. (7) should be consistent with
that of the electric field $\mathbf{E}$ and magnetic field $\mathbf{B}$ [5].
Furthermore, this current-driven method computes the electrostatic fields
indirectly by Eq. (7), rather than directly through the diverge Eq. (3) [23].
To ensure its accuracy, the current interpolation algorithm are meticulously
designed to uphold the charge-conservation law [24]-[27].

In this paper, we present JefiPIC, an alternative PIC method that utilizes
Jefimenko’s equations [28],[29], the general solution to Maxwell’s equations,
to compute EM fields in plasma simulations by the integral form. Jefimenko’s
equations are computationally time-consuming, particularly when used in
conjunction with the numerous particles updating in PIC methods. To enhance
the performance of this integral simulations, we employ the GPU-based package
JefiGPU [30], which directly calculate Jefimenko’s equation on the GPU with a
high-dimensional integration package ZMCintegral [31]. This approach allows us
to efficiently evolve plasma systems within an acceptable execution time.
Furthermore, we perform particle motion equations on GPU, making JefiPIC a
fully GPU-based plasma simulator. All variables are consistently defined in
GPU memory, ensuring swift execution of the simulations. Providing the initial
conditions of charged particles (neutral or not) and the background EM fields,
our code is supposed to automatically compute the particles’ space-time
distribution as well as the evolution of EM fields.

Compared to the traditional FDTD-based PIC simulation, our proposed JefiPIC
has the following advantages:

a) The numerical stability of the integral eliminates the need for the CFL
stability condition [32], a criterion in finite difference methods, to prevent
numerical divergence. This provides more flexibility and robustness in
selecting the computational time step.

b) Our method does not require the use of complex boundary layers to cut off
the propagation of EM fields [33],[34]. Besides, simply ’killing’ particle in
the boundary layers will not result in unintended charge deposition or errors
in the electrostatic components [35],[36]. Thus, our method is well-suited for
describing plasma systems with open boundaries.

c) Our method can handle non-neutral plasma without the need for additional
pretreatments, such as solving the Poisson equation [5], even for those new to
the field, making it easier to treat the initial conditions.

d) In JefiPIC, EM fields are precisely calculated by the utilization of both
charge and current densities, which strictly adheres to the charge
conservation law at every time step. Thus, instead of using strict charge
conservation algorithms [25] or implementing electric field divergence
amendments [4],[23], the simple linear interpolation method can be employed.

Apart from the above merits, our approach inevitably encounters some drawbacks
in its current version. For example, the large computational cost, and
inflextibility in dealing with the EM boundaries with matters, etc. We will
study these issues to make our code more practical in the future.

The paper is structured as follows. Section 2 outlines the detailed design
process of JefiPIC, including the Jefimenko’s equations, particle motion, the
interpolation method for particles and fields, and the implementation of
JefiPIC on GPU. In Section 3, we compare the results obtained from JefiPIC,
UNIPIC, EPOCH, and RBG-Maxwell [37],[38] through examining the evolutions of
three different plasma systems. We further analyze the advantages of JefiPIC
compared to other available codes. Finally, in Section 4, we present our
conclusions and discuss future research directions.

##  2 Algorithm and Implementation

As a PIC method, JefiPIC has a similar architecture to the traditional EM PIC,
shown in Figure 1. The method involves four processes: 1) computing EM fields
on each grid, 2) solving the field interpolation to calculate the forces on
the particles, 3) solving the particle interpolation to allocate the charges
on the grids, and 4) updating the particle motion. In this section, we will
introduce the details of JefiPIC in the sequence of Jefimenko’s equations,
particle motion, field and particle interpolation and the implementation on
GPU.

Figure 1: Process diagram for JefiPIC simulation. The overall computational
architecture is similar to traditional PIC methods, except for the use of a
different EM field solver and the negligible consideration of EM field and
particle boundary issues.

###  2.1 Jefimenko’s equations

In JefiPIC, the aim is to compute the EM fields over the entire computational
region by the Jefimenko’s equations with sources $\mathbf{J}$ and $\rho$,

$\displaystyle\mathbf{E}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi\varepsilon_{0}}\int\left[\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{3}}\rho(\mathbf{r}^{\prime},t_{\mathrm{r}})+\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{2}}\frac{1}{c}\frac{\partial\rho(\mathbf{r}^{\prime},t_{\mathrm{r}})}{\partial
t}\right.$ (8)
$\displaystyle\left.-\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|}\frac{1}{c^{2}}\frac{\partial\mathbf{J}(\mathbf{r}^{\prime},t_{\mathrm{r}})}{\partial
t}\right]\mathrm{d}^{3}\mathbf{r}^{\prime},$
$\displaystyle\mathbf{B}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle-\frac{1}{4\pi}\int\left[\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{3}}\times\mathbf{J}(\mathbf{r}^{\prime},t_{\mathrm{r}})\right.$
(9)
$\displaystyle\left.+\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|^{2}}\times\frac{1}{c}\frac{\partial\mathbf{J}(\mathbf{r}^{\prime},t_{\mathrm{r}})}{\partial
t}\right]\mathrm{d}^{3}\mathbf{r}^{\prime},$ $\displaystyle t_{\mathrm{r}}$
$\displaystyle=$ $\displaystyle t-|\mathbf{r}-\mathbf{r}^{\prime}|/c,$ (10)

where $(\mathbf{r},t)$ represents the space-time point in the computational
region, while $\mathbf{r}=(x,y,z)$ and
$\mathbf{r}^{\prime}=(x^{\prime},y^{\prime},z^{\prime})$ separately refer to
the displacements of the EM fields and sources. For numerical computation, we
divide the computational region into the same structured cuboid grids. All the
EM quantities $\mathbf{J}$, $\rho$, $\mathbf{E}$ and $\mathbf{B}$ are defined
at the center of the grids. To perform integrations numerically, Eqs. (8) to
(11) are discretized as,

$\displaystyle\mathbf{E}(\mathbf{r}_{I,J,K},t^{n})$ $\displaystyle=$
$\displaystyle\frac{\mathrm{d}x^{\prime}\mathrm{d}y^{\prime}\mathrm{d}z^{\prime}}{4\pi\varepsilon_{0}}\sum_{i,j,k}\left[\frac{\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}}{|\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}|^{3}}\rho(\mathbf{r}_{i,j,k}^{\prime},t_{\mathrm{r}})\right.$
(11)
$\displaystyle+\frac{\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}}{|\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}|^{2}}\frac{1}{c}\frac{\rho(\mathbf{r}_{i,j,k}^{\prime},t_{\mathrm{r}})-\rho(\mathbf{r}_{i,j,k}^{\prime},t_{\mathrm{r}}-\mathrm{d}t)}{\mathrm{d}t}$
$\displaystyle\left.-\frac{1}{|\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}}\frac{1}{c^{2}}\frac{\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{\mathrm{r}})-\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{\mathrm{r}}-\mathrm{d}t)}{\mathrm{d}t}\right],$
$\displaystyle\mathbf{B}(\mathbf{r}_{I,J,K},t^{n})$ $\displaystyle=$
$\displaystyle-\frac{\mathrm{d}x^{\prime}\mathrm{d}y^{\prime}\mathrm{d}z^{\prime}}{4\pi}\sum_{i,j,k}\left[\frac{\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}}{|\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}|^{3}}\times\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{\mathrm{r}})\right.$
(12)
$\displaystyle\left.+\frac{\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}}{|\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}|^{2}}\times\frac{1}{c}\frac{\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{\mathrm{r}})-\mathbf{J}(\mathbf{r}_{i,j,k}^{\prime},t_{\mathrm{r}}-\mathrm{d}t)}{\mathrm{d}t}\right],$
$\displaystyle t_{r}$ $\displaystyle=$ $\displaystyle
t^{n}-|\mathbf{r}_{I,J,K}-\mathbf{r}_{i,j,k}^{\prime}|/c,$ (13)

where the subscripts $I$, $J$, and $K$ of displacement denote the grid index
of the unknown EM fields, while the subscripts $i$, $j$, and $k$ of
displacement denote the grid index of the source. The superscript n of time
refers to the n-th time step, and $\mathrm{d}x^{\prime}$,
$\mathrm{d}y^{\prime}$, and $\mathrm{d}z^{\prime}$ represent the size of a
grid.

###  2.2 Particle motion

The Newton–Lorentz force Eqs. (5) to (6) are discretized by the central leap-
frog difference method as,

$\displaystyle\frac{\gamma
m(\mathbf{v}^{n+1/2}-\mathbf{v}^{n-1/2})}{\mathrm{d}t}$ $\displaystyle=$
$\displaystyle
q\left(\mathbf{E}^{n}+\frac{\mathbf{v}^{n+1/2}+\mathbf{v}^{n-1/2}}{2}\times\mathbf{B}^{n}\right),$
(14) $\displaystyle\frac{\mathbf{r}^{n+1}-\mathbf{r}^{n}}{\mathrm{d}t}$
$\displaystyle=$ $\displaystyle\mathbf{v}^{n+1/2},$ (15)

This advancement of velocity involves an implicit equation. To figure out this
issue, the Boris particle pusher is employed to separate Eq. (15) into three
steps as,

$\displaystyle\mathbf{v}^{-}$ $\displaystyle=$
$\displaystyle\mathbf{v}^{n-1/2}+\frac{q\mathbf{E}^{n}}{m}\frac{\mathrm{d}t}{2},$
(16) $\displaystyle\frac{(\mathbf{v}^{+}+\mathbf{v}^{-})}{\mathrm{d}t}$
$\displaystyle=$ $\displaystyle\frac{q}{2\gamma
m}(\mathbf{v}^{+}+\mathbf{v}^{-})\times\mathbf{B}^{n},$ (17)
$\displaystyle\mathbf{v}^{n+1/2}$ $\displaystyle=$
$\displaystyle\mathbf{v}^{+}+\frac{q\mathbf{E}^{n}}{m}\frac{\mathrm{d}t}{2},$
(18)

where $\mathbf{v}^{-}$ and $\mathbf{v}^{+}$ donate the temporary variables to
update the velocity.

###  2.3 Charge, current and field interpolation method

Based on the discretization of Jefimenko’s equations, the EM fields, current
and charge densities all need to be allocated at the center of the grids,
which is different from that in FDTD. The traditional PIC methods usually
adopt two kinds of interpolation methods. One approach is the charge-
conservation method [39], which generate numerical noises in the solenoidal
part of current density [4],[23],[40]. The other method is the sample linear
interpolation method, which generates less noise [4], but violates the charge
conservation law, leading to inaccuracies in electric field. To compensate for
this loss in accuracy, electric field correction algorithms must be
implemented.

Because the integral-based PIC method calculates the fields directly from the
provided current and charge density, JefiPIC avoids the issue of charge
conservation. As a result, we can use the first-order linear interpolation
method for charge and field allocation in JefiPIC, combining the advantages of
both acceptable numerical noise and computational time.

We firstly exhibit how particle interpolation works. In a 3-D model, linear
interpolation is achieved through volume-weighted interpolation. For
simplicity, we use a 2-D model with the area-weighted interpolation method, as
shown in Figure 2. Assuming that there is one particle in the computational
region, labeled as a ’star’. The particle’s charge is distributed to the
nearest four center nodes indicated by the circles with indices $(i,j)$,
$(i+1,j)$, $(i,j+1)$ and$(i+1,j+1)$. The area-weighted interpolation method
specifies that the allocated charge on each node is proportional to the
weight, which is the ratio of the area of its thither rectangle to that of a
single grid. Accordingly, the weight can be computed by the corresponding
volume in 3-D models.

Figure 2: Area-weighted interpolation. The particle is indicated by a ’star’
and its charge is interpolated to the four neighboring grids according to the
ratio of the corresponding shaded area to a single grid area.

Since all the physical variables are defined at the center nodes of the grids,
we firstly need to find the indices of the grids which the particles belong
to,

$\displaystyle
i=\left[\frac{x-0.5\cdot\mathrm{d}x^{\prime}}{\mathrm{d}x^{\prime}}\right],$
$\displaystyle
j=\left[\frac{y-0.5\cdot\mathrm{d}y^{\prime}}{\mathrm{d}y^{\prime}}\right],$
$\displaystyle
k=\left[\frac{z-0.5\cdot\mathrm{d}z^{\prime}}{\mathrm{d}z^{\prime}}\right],$
(19)

where the square bracket means rounding down operation, and $x,$ $y$, and $z$
are the displacements in the three axes. Following the volume-weighted
interpolation method, a particle located at one grid $(i,j,k)$ should
interpolate its charge to the adjacent grid nodes as,

$Q_{i,j,k}=(1-\alpha)(1-\beta)(1-\chi)q,$

$Q_{i,j,k}=\alpha(1-\beta)(1-\chi)q,$

$Q_{i,j,k}=(1-\alpha)\beta(1-\chi)q,$

$Q_{i,j,k}=(1-\alpha)(1-\beta)\chi q,$

$Q_{i,j,k}=\alpha\beta(1-\chi)q,$

$Q_{i,j,k}=\alpha(1-\beta)\chi q,$

$Q_{i,j,k}=(1-\alpha)\beta\chi q,$

$Q_{i,j,k}=\alpha\beta\chi q,$ (20)

where

$\displaystyle\alpha=\frac{x-(i+0.5)\cdot\mathrm{d}x^{\prime}}{\mathrm{d}x^{\prime}},$
$\displaystyle\beta=\frac{y-(j+0.5)\cdot\mathrm{d}y^{\prime}}{\mathrm{d}y^{\prime}},$
$\displaystyle\chi=\frac{z-(k+0.5)\cdot\mathrm{d}z^{\prime}}{\mathrm{d}z^{\prime}},$
(21)

and $q$ is the charge of the moving particle. We can acquire the charge and
current density at node $(i,j,k)$ by,

$\displaystyle\rho_{i,j,k}$ $\displaystyle=$
$\displaystyle\frac{Q_{i,j,k}}{\mathrm{d}x^{\prime}\cdot\mathrm{d}y^{\prime}\cdot\mathrm{d}z^{\prime}},$
(22) $\displaystyle\left(J_{i,j,k,x},J_{i,j,k,y},J_{i,j,k,z}\right)$
$\displaystyle=$ $\displaystyle\rho_{i,j,k}\cdot(v_{x},v_{y},v_{z}),$ (23)

where $v_{x}$, $v_{y}$ and $v_{z}$ are the velocities of the particle. The
charge and current densities at the corresponding nodes can be obtained
according to Eqs. (20) to (23).

When a particle is located closer to a computational region boundary than half
a grid length, it is important to note that the adjacent nodes cannot
constitute a full grid. To address this issue, we utilize virtual grids (red
grids) to surround the computational region (black grids), as illustrated in
Figure 3. We can then implement the same charge interpolation process as
displayed in Figure 2, while excluding the computation of charge and current
densities on the virtual nodes.

Figure 3: Virtual grids in the charge allocation. The black grids and nodes
represent the real computational region, while the red dashed grids and nodes
are in the virtual one. The linked green dashed grids are used for regular 2-D
charge interpolation as in Figure 2. During this process, only the charge
densities on the black nodes will be computed and saved.

In this part, we will show the field interpolation. Once the EM fields are
established on all observation grids, it is required to compute the forces on
the particles using the field interpolation method. The field interpolation
process works in reverse to that of the charge interpolation method, whereby
the weighted EM fields from the adjacent eight grids are summed to calculate
the fields acting on a certain particle. The weight of electric field from one
grid on the particle is equal to the weight of the particles’ charge allocated
to that corresponding grid. Thus, the electric field on a particle located in
grid $(i,j,k)$ can be evaluated as,

$\displaystyle\mathbf{E}_{i,j,k}^{q}$ $\displaystyle=$
$\displaystyle(1-\alpha)(1-\beta)(1-\chi)\mathbf{E}_{i,j,k}+\alpha(1-\beta)(1-\chi)\mathbf{E}_{i+1,j,k}$
(24)
$\displaystyle+(1-\alpha)\beta(1-\chi)\mathbf{E}_{i,j+1,k}+(1-\alpha)(1-\beta)\chi\mathbf{E}_{i,j,k+1}$
$\displaystyle+\alpha\beta(1-\chi)\mathbf{E}_{i+1,j+1,k}+\alpha(1-\beta)\chi\mathbf{E}_{i+1,j,k+1}$
$\displaystyle+(1-\alpha)\beta\chi\mathbf{E}_{i,j+1,k+1}+\alpha\beta\chi\mathbf{E}_{i+1,j+1,k+1}.$

The dimensionless parameters in Eq. (24) correspond to those in Eq. (21), thus
allowing for the computation of the magnetic field $\mathbf{B}$ on the
particle through the same method. It is worth noting that the number of
adjacent nodes involved in driving the particle follows the same rule as
particle interpolation, and the EM fields outside the computational region are
not calculated. Thus, once the particles pass through the computational
boundary, the forces acting upon them are no longer considered, and they do
not return to the computational region. As we do not solve the charge
conservation equation, we can exclude their contribution to the source without
resulting in undesired charge deposition.

###  2.4 GPU Implementation

According to section 2.1 to 2.3, JefiPIC needs three major CUDA kernels, each
taking over of the Jefimenko’s equation, particle motion, and the charge and
field interpolation, respectively.

The first CUDA kernel is designed to perform the summations presented in Eqs.
(11) to (13) on GPU. Since the EM fields of a grid are produced by all
$\mathbf{J}$ and $\rho$ in the computational region with retarded time
$t_{r}$, we need to figure out the quantity of sources that should be saved.
Here is our solution. Firstly, we calculate the maximum length
($L_{\text{max}}$) across the computational region, while the value of
$L_{\text{max}}$ is typically determined by the diagonal length of the
computational region. We then divide $L_{\text{max}}$ by $c\mathrm{d}t$ to
obtain $N_{\mathrm{d}t}=\text{Celing}(L_{\text{max}}/(c\mathrm{d}t))$, which
represents the maximum number of time steps an EM field travelling across the
entire computational region. Consequently, the sources $\mathbf{J}$ and $\rho$
for at least the most recent $N_{\mathrm{d}t}$ time steps need to be stored to
guarantee the EM field computation. For instance, at time $t_{\mathrm{A}}$, we
need to save $\mathbf{J}$ and $\rho$ at all grids (the number of girds is M)
for $N_{\mathrm{d}t}$ time steps, i.e.,

$t_{\mathrm{A}}-\mathrm{d}t:[\rho(r_{1},t_{\mathrm{A}}-\mathrm{d}t),J(r_{1},t_{\mathrm{A}}-\mathrm{d}t)],\ldots,[\rho(r_{\mathrm{M}},t_{\mathrm{A}}-\mathrm{d}t),J(r_{\mathrm{M}},t_{\mathrm{A}}-\mathrm{d}t)],$

$t_{\mathrm{A}}-2\cdot\mathrm{d}t\text{:}[\rho(r_{1},t_{\mathrm{A}}-2\cdot\mathrm{d}t),J(r_{1},t_{\mathrm{A}}-2\cdot\mathrm{d}t)],\ldots,[\rho(r_{\mathrm{M}},t_{\mathrm{A}}-2\cdot\mathrm{d}t),J(r_{\mathrm{M}},t_{\mathrm{A}}-2\cdot\mathrm{d}t)],$

$\cdots\cdots$

$t\mathrm{{}_{A}}-N_{\mathrm{d}t}\cdot\mathrm{d}t:[\rho(r_{1},t\mathrm{{}_{A}}-N_{\mathrm{d}t}\cdot\mathrm{d}t),J(r_{1},t\mathrm{{}_{A}}-N_{\mathrm{d}t}\cdot\mathrm{d}t)],\ldots,[\rho(r_{\mathrm{M}},t\mathrm{{}_{A}}-N_{\mathrm{d}t}\cdot\mathrm{d}t),J(r_{\mathrm{M}},t\mathrm{{}_{A}}-N_{\mathrm{d}t}\cdot\mathrm{d}t)].$

For time steps $N$ that are ahead of the initial time
$t_{0}(t\mathrm{{}_{A}}-N\cdot\mathrm{d}t<t_{0})$, $\mathbf{J}$ and $\rho$ are
simply set to zero. Here $\mathbf{J}$ and $\rho$ are cached in the global GPU
memory. Obviously, this procedure requires a large amount of memory space on
GPU. In practice, we limit $N_{\mathrm{d}t}$ to 10000, i.e.,
$N_{\mathrm{d}t}=\text{min}(N_{\mathrm{d}t},10000)$, making the data storage
achievable for most GPU apparatus.

In the implementation of the integral, each thread in the CUDA kernel
calculates the EM fields corresponding to a certain grid. Thus, we loop
through each grid in the computational region. For convenience, we denote the
$i-th$ CUDA thread as $thread\text{\\_}i$, which calculates the EM fields on
the $i-th$ grid with index $grid\text{\\_}i$. To evaluate the EM fields of
$grid\text{\text{\\_}}i$, we loop through all the grids to compute the
Jefimenko’s equation through the local sources $\boldsymbol{\mathbf{J}}$ and
$\rho$ with retarded time tr (except for $grid\text{\\_}i$ itself to prevent
numerical divergence) and sum all the EM fields. This procedure is depicted in
Figure 4.

Figure 4: Process of the implementation of solving Jefimenko’s equations on
GPU. Supposing that the computational region contains 6 grids, and the process
involves two steps. In the first step, EM fields on one grid (for example
$grid\text{\\_}i$) are computed by looping through the other grids with
$\mathbf{J}$ and $\rho$. In the second step, the entire EM fields are computed
by looping through the rest grids repeating the operation in step 1.

The second CUDA kernel is to achieve the particle motion exhibited by Eqs.
(14) to (17). Here, the particle’s displacements for full time steps and
velocities for half time steps, i.e., $r^{n}$, $v^{n-1/2}$, are saved in the
global GPU memory. The information update of each particle is also handled by
an individual thread in the CUDA kernel. Hence, the maximum particle numbers
allowed in our code is restricted by the maximum number of threads. For one
NVIDIA A100 card, this number is about $10^{9}$, meaning that we cannot
calculate more than $10^{9}$ particles on one GPU card.

Since $\mathbf{J}$, $\rho$, $\mathbf{E}$ and $\mathbf{B}$ are defined at the
grid center node, we need to build the third CUDA kernel to calculate the
forces on the particles and the charge and current densities on all the grid
centers. To obtain $\mathbf{E}$ and $\mathbf{B}$, we first iterate through all
the particles in the computational region to compute $\mathbf{J}$ and $\rho$
for each grid. We then use the known $\mathbf{J}$ and $\rho$ with a retarded
time to calculate $\mathbf{E}$ and $\mathbf{B}$ by iterating through all the
grids. Here, $\mathbf{J}$, $\rho$, $\mathbf{E}$ and $\mathbf{B}$ are cached in
the global GPU memory and all intermediate variables like $\alpha$ and $\beta$
are set temporary local variables.

For plasma systems at different length scales, the numerical value of physical
quantities can vary significantly. Given that the numerical values in JefiPIC
must lie within the range of machine precision (float64), we must convert
these quantities from SI unit system to a new unit system. In this paper, we
use the Flexible Unit (FU) system (a brief introduction of FU is presented in
the appendix), with which our code can simulate plasma systems at different
length scales on GPU. By assigning proper values for constants such as vacuum
permittivity $\varepsilon_{0}$, reduced Planck constant $\hbar$, speed of
light $c$, and $\lambda$ (which relates the energy in SI and FU), we can limit
most of the numerical values to machine precision. Note that the use of FU
does not change the physics of the plasma system.

##  3 Computational Model and Result

To exhibit the comprehensive capabilities of JefiPIC, we performed a
comparative study of three plasma models against three other 3-D codes. We
simulated three non-neutral electron plasma models, including a) electrons
emitted from a point to verify the correctness of the particle and field cut-
off boundaries, b) electrons emitted with zero initial velocity to demonstrate
the ease to handle non-neutral plasma, and c) electrons emitted with random
initial velocity, showcasing the natural charge conservation feature. The
alternative codes we used include UNIPIC, a mature PIC code, EPOCH, an open-
source PIC code, and RGB-Maxwell, a plasma simulator. UNIPIC and EPOCH are
traditional FDTD-based PIC codes, while RGB-Maxwell solves the Boltzmann
equations. Here, we list some presets of these codes.

Particle boundary condition — The three PIC codes all use the cut-off boundary
to ’kill’ the particles that leave the computational region, and RGB-Maxwell
directly estimate the particle distribution function outside the computational
region.

Field boundary condition — Codes based on difference method utilize
convolutional perfect matched layers (CPML) as the field boundary, while those
based on integral method use a simple cut-off boundary.

Space grid and time step — All four codes divide the computational region with
the same spatial grid size of $\mathrm{d}x=\mathrm{d}y=\mathrm{d}z=10^{-5}$ m.
The total computational time is set to 1 ns. UNIPIC and EPOCH use the
difference methods, thus requiring the time step $\mathrm{d}t$ to at least
satisfy the CFL condition,

$\displaystyle c\cdot\mathrm{d}t$ $\displaystyle\leq$
$\displaystyle\frac{1}{\sqrt{\frac{1}{(\mathrm{d}x)^{2}}+\frac{1}{(\mathrm{d}y)^{2}}+\frac{1}{(\mathrm{d}z)^{2}}}},$
(25)

limiting $\mathrm{d}t$ to no more than around $1.92\times 10^{-14}$ s. JefiPIC
and RGB-Maxwell utilize Jefimenko’s equations instead, which are not limited
by numerical stability. Therefore, we set the time step as
$\mathrm{d}t=10^{-13}$ s, which is much larger than that used in difference
methods.

Computational region — For model 1, we set the computational region size as
$n_{x}\cdot n_{y}\cdot n_{z}=1\times 101\times 101$ in JefiPIC and RGB-
Maxwell. For model 2 and model 3, we set the computational region size as
$n_{x}\cdot n_{y}\cdot n_{z}=1\times 251\times 111$ in JefiPIC and RGB-
Maxwell, and $n_{x}\cdot n_{y}\cdot n_{z}=3\times 251\times 111$ in EPOCH and
UNIPIC.

Particle — In each simulation, we emit the total charge of $-5\times 10^{-14}$
C, with an electron count of $3.125\times 10^{5}$. To efficiently achieve the
computation, we use macro-particle to represent electrons in the similar
phrase space in PIC space. The ways to handle macro-particles are different in
different codes. JefiPIC implement individual tracking for all macro-particles
using CUDA kernels on an Nvidia A100 GPU card, which has a maximum limit of
around $10^{9}$ CUDA threads. Therefore, we are limited to running simulations
with a maximum of $10^{9}$ particles at once. To reduce the numerical noise,
JefiPIC simulation uses up to $3.125\times 10^{7}$ particles. In EPOCH, the
particle weight is dynamically adjusted by the code to ensure an approximate
count of 100 particles in each grid. In UNIPIC, a suitable particle weight is
set to balance numerical noise and execution time. On the other hand, RGB-
Maxwell used a distribution function on a six-dimensional phase space instead
of the concept of ’particles’. The initial velocity is represented in the
units of m/s in our simulation, and energy corresponding to the velocities
$v_{0}$ used in the following experiments equals to 10 eV.

###  3.1 Perfect Charge and Field Boundary Layer

The boundary of PIC contains two components: the field boundary and particle
boundary. In difference-based PIC, the accuracy of the EM fields is highly
sensitive to the choice of field boundary layers, and the simple method for
’killing’ particles in the boundary layer can introduce undesirable
electrostatic field errors. JefiPIC, on the other hand, effectively handles
the absorption of fields and particles at the computational boundary with a
simple cut-off boundary, aided by the integral operation. The evolution of
electrons emitted in all directions from the central grid $(i=1,j=50,k=50)$ is
simulated using JefiPIC and compared to RGB-Maxwell. The electrons are
initialized with a Gaussian distribution of velocities, with an average
velocity of $v_{0}$ and a standard deviation of 0.2$v_{0}$.

Figure 5: Initial particle distribution. The computational region was of size
$1\cdot\mathrm{d}x\times 101\cdot\mathrm{d}y\times 101\cdot\mathrm{d}z$ in
JefiPIC and RGB-Maxwell. The particles are placed in the shaded region
uniformly in the cuboid volume of $1\cdot\mathrm{d}x\times
1\cdot\mathrm{d}y\times 1\cdot\mathrm{d}z$ ranging from $50\cdot\mathrm{d}y$
to $51\cdot\mathrm{d}y$ in the $y-\mathrm{axis}$ and $50\cdot\mathrm{d}z$ to
$51\cdot\mathrm{d}z$ in the $z-\mathrm{axis}$. The velocities of the particles
obey the Gaussian distribution.

Figure 6 depicts the particle distribution of the two models at different time
instances. The top row displays the results from the particle model, while the
bottom row shows those from the distribution function. The figure demonstrates
that the two particle distributions are nearly identical, with most of the
particles forming a prominent light ring that spreads out gradually over time.
This outcome confirms the general results obtained from JefiPIC. Moreover,
since the integral method naturally filters out particles that are outside the
computational domain, no charge deposition occurs when the particles pass
through the computational boundary, and any electromagnetic force on those
escape particles would not be taken into account. Consequently, the motion of
the particles is not hindered, and the ring-like particle distribution is
sustained.

Figure 6: Particle distributions over time. Panel (a) shows the distribution
at 0.2 ns, followed by (b) at 0.4 ns, (c) at 0.6 ns, and (d) at 0.8 ns. The
two rows represent the particle distributions obtained from JefiPIC and RGB-
Maxwell, respectively.

In figure 7, the electric fields of the two models are plotted at grid
$(i=1,j=75,k=75)$, which is $25\cdot\mathrm{d}x$ and $25\cdot\mathrm{d}y$
distant from the particle emission point. As particles spread towards the
diagnostic point, the electric field initially becomes negative and then
transitions to positive before gradually decreasing towards zero as most of
the particles pass and move away from the diagnostic point. The electric
fields in both models display similar waveforms that are smooth with no
observable reflections, suggesting that the electric field does not interfere
with the particle distribution.

Figure 7: Comparison of y-direction electric field. The $y-\mathrm{direction}$
electric fields at point $(i=1,j=75,k=75)$ derived from JefiPIC and RGB-
Maxwell are compared. The black square and red circle represent the results
from RGB-Maxwell and JefiPIC.

Through figure 6 and figure 7, we can observe that despite utilizing the
first-order charge and field interpolation method, JefiPIC delivers results as
accurate as those obtained from the second-order accurate plasma simulator
RGB-Maxwell.

As a result, the integral-based PIC method JefiPIC is capable of naturally
cutting off fields and particles at the computational boundary without
reflecting the fields or causing charge deposition errors. This makes JefiPIC
a more effective approach compared to the difference-based PIC methods that
rely on boundary conditions. Thus, JefiPIC is particularly well-suited for
modeling open-boundary problems, such as space plasma [42], quark-gluon plasma
[38], or high-altitude nuclear explosions [43]-[45].

###  3.2 Non-neutral plasma condition

The integral-based PIC method has the advantage of its ability to handle non-
neutral plasma without extra pre-processing. Here, we compare the evolution of
a batch of electrons with zero velocity in JefiPIC, EPOCH, and RGB-Maxwell. To
better observe the evolution of particles, we apply a constant magnetic field
of $\mathbf{B}=10$ T along the $y$-direction to constrain the electrons’
transverse motion, considering velocity only in the $y$-direction. Initially,
these electrons are placed uniformly within a cuboid region of size
$1\cdot\mathrm{d}x\times 1\cdot\mathrm{d}y\times 101\cdot\mathrm{d}z$, ranging
from $0\cdot\mathrm{d}x$ to $1\cdot\mathrm{d}x$ ($1\cdot\mathrm{d}x$ to
$2\cdot\mathrm{d}x$ in EPOCH) in the $x$-axis, $124\cdot\mathrm{d}y$ to
$125\cdot\mathrm{d}y$ in the $y$-axis and $5\cdot\mathrm{d}z$ to
$106\cdot\mathrm{d}z$ in the $z$-axis, shown in Figure 8.

Figure 8: Initial particle distribution. The computational region was of size
$1\cdot\mathrm{d}x\times 251\cdot\mathrm{d}y\times 111\cdot\mathrm{d}z$ in
JefiPIC and RGB-Maxwell, and of size $3\cdot\mathrm{d}x\times
251\cdot\mathrm{d}y\times 111\cdot\mathrm{d}z$ in EPOCH. The particles are
placed in the shaded region uniformly in the cuboid volume of
$1\cdot\mathrm{d}x\times 1\cdot\mathrm{d}y\times 101\cdot\mathrm{d}z$ ranging
from $0\cdot\mathrm{d}x$ to $1\cdot\mathrm{d}x$ ($1\cdot\mathrm{d}x$ to
$2\cdot\mathrm{d}x$ in EPOCH) in the $x$-axis, $124\cdot\mathrm{d}y$ to
$125\cdot\mathrm{d}y$ in the $y$-axis and $5\cdot\mathrm{d}z$ to
$106\cdot\mathrm{d}z$ in the $z$-axis. The velocities of the particles are
zero.

Due to the mutual repulsion between particles, it is expected that they should
diffuse from the center to both sides along $y$-axis. As depicted in Figure 9,
particles simulated by JefiPIC and RGB-Maxwell consistently spread out as
expected. However, EPOCH fails to demonstrate a reasonable distribution when
pre-processing tasks, such as Poisson’s equation, are not executed correctly.
This outcome can be explained by the fact that JefiPIC and RGB-Maxwell utilize
Jefimenko’s equations, which can model the EM field as long as charge and
current densities are supplied. Therefore, even in the absence of initial
current, the electrostatic fields calculated from the charge density can still
drive the charges apart. In contrast, EPOCH is a current-based PIC code. To
achieve the same result, we must perform the pre-processing to offer the
equivalent electrostatic fields, such as solving the Poisson’s equation.
However, this essential operation will lead to additional computation time.

Figure 9: Particle distributions over time. The three rows represent the
particle distributions obtained from JefiPIC, RGB-Maxwell, and EPOCH,
respectively. Panel (a) shows the distribution at 0.0 ns, followed by (b) at
0.4 ns, (c) at 0.6 ns, and (d) at 1.0 ns.

To further validate the accuracy of JefiPIC, we compare its particle number
density distribution with the analytical solution [46],

$\displaystyle\psi$ $\displaystyle\sim$
$\displaystyle\frac{1}{\tau\left[y^{2}/2\tau^{2}+\psi_{0}\right]},$ (26)

where $\psi$,$\,\tau$, and $y$ represent the normalized number density, time
and distance. As illustrated in Figure 10, the particle distributions produced
by JefiPIC closely match the analytical solution.

Figure 10: Simulated and analytical particle number density distribution. The
figure displays the simulated particle number density distribution obtained
through JefiPIC and the analytical solution from Reference [46] at $t=0.4$ ns,
0.6 ns, and 1.0 ns. The data are extracted from the grids ranging from [0·d$x$
to 1·d$x$, 125·d$y$ to 225·d$y$, 54·d$z$ to 55·d$z$]. The gray hollow square,
pink hollow triangle and blue hollow rhombus represent the results by JefiPIC,
while the black, red, and dark blue lines represent the analytical results.

###  3.3 Obeying Charge Conservation Law

In this subsection, we choose another FDTD-based PIC code, UNIPIC, to further
demonstrate the advantage of the integral-based PIC codes when dealing with
the charge-conservation problem. In this model, only electrons are emitted
from the shaded rectangle face, shown in Figure 11. The electrons’ initial
velocities in the $y$-direction are randomly from 0 to $v_{0}$. In line with
section 3.2, we apply a constant magnetic field of $\mathbf{B}$ = 10 T in the
$y$-direction again to constrain the electrons’ transverse motion. Thus, only
velocities in the $y$-direction are considered.

Figure 11: Initial electron distribution in the simulation. The total
computational region was of size 1·d$x$×251·d$y$×111·d$z$ in JefiPIC and
3·d$x$×251·d$y$×111·d$z$ in UNIPIC. The particles were placed uniformly on the
shaded rectangle face of $y=9\cdot\mathrm{d}y$ and 1·d$x$×101·d$z$, with
0·d$x$ to 1·d$x$ (1·d$x$ to 2·d$x$ in UNIPIC) in the $x$-axis and 5·d$z$ to
106·d$z$ in the $z$-axis. Particle velocities obeyed a random distribution
from 0 to $v_{0}$.

Since UNIPIC emits electrons from the metal emission surface, the produced
positive charges according to the charge-conservation law is shielded [27].
Therefore, we can regard the model in UNIPIC as a pure electron plasma.
However, UNIPIC employs the simple linear interpolation method rather than the
charge-conservation method, and thus, it uses the Langdon-Marder correction
algorithm from time to time.

Figure 12 displays the particle distribution obtained from JefiPIC and UNIPIC
at two different time points. A comparison of the two distributions reveals
their overall similarity, validating our code to a certain extent. However,
particles simulated by JefiPIC exhibits a larger $y$-direction displacement
compared to UNIPIC. Table 1 lists the maximum displacement and velocity of
particles under corresponding time steps in Figure 12. It is shown that the
particles simulated by JefiPIC attained larger velocities than those simulated
by UNIPIC, leading to more extensive movements. We attribute the difference to
dissimilar Coulomb forces exhibited by the particles in each simulation. This
suggests that since UNIPIC does not strictly adhere to charge conservation at
every time step, it is not capable of capturing the accurate electric fields
all the time. In contrast, JefiPIC does not face such issues and can simulate
particles with greater accuracy.

Figure 12: Comparison of particle distributions at t = 0.4 ns and t = 0.6 ns in (a) JefiPIC and (b) UNIPIC simulations. Initial velocities were applied only in the $y$-direction, randomly ranging from 0 to $v_{0}$. This illustrates the differences in Coulomb forces between the two simulations, highlighting the more accuracy of JefiPIC in simulating particle movements. Table 1: Comparison of the maximum particle displacement and velocity between JefiPIC and UNIPIC Time/ns | Maximum Displacement/m | Maximum Velocity/(m·s-1)
---|---|---
JefiPIC | UNIPIC | JefiPIC | UNIPIC
0.4 | $1.01\times 10^{-3}$ | $8.49\times 10^{-4}$ | $2.47\times 10^{6}$ | $2.26\times 10^{6}$
0.6 | $1.51\times 10^{-3}$ | $1.30\times 10^{-3}$ | $2.53\times 10^{6}$ | $2.30\times 10^{6}$

Figure 13 exhibits the $y$-direction electric fields at point
$(i=1,j=50,k=55)$ in JefiPIC and at point $(i=2,j=50,k=55)$ in UNIPIC for
detail. The electric fields of the two simulations differ in three aspects: a)
JefiPIC’s field experiences a little earlier peak time, b) the pulse width of
JefiPIC’s field is larger, and c) JefiPIC’s field peak partly exceeds that of
the UNIPIC. Therefore, the fields calculated via integral equations
effectively accelerate particles over a longer duration, leading to a wider
particle distribution and an earlier peak time of electric fields. Besides,
JefiPIC shows a significantly smoother electric field profile in comparison to
UNIPIC. Thus, in the aspect of charge conservation, the integral method
provides a natural advantage in terms of accuracy and simplicity.

Figure 13: Comparison of $y$-direction electric field. The $y$-direction
electric fields at point $(i=0,j=50,k=55)$ derived from JefiPIC and at point
$(i=1,j=50,k=55)$ from UNIPIC are compared. The black square and red circle
represent the results from JefiPIC and UNIPIC.

We also compared the computational performance of the two codes. UNIPIC was
run on a workstation with Intel Xeon Processor E7450 CPU across 200 threads,
and it took over 6 hours to complete the calculations. In contrast, using a
larger time step, JefiPIC completed the computations within only half an hour
on a single GPU A100 card, achieving over 90% time-saving benefits.

###  3.4 Choice of Macro-Particle

We conduct test on JefiPIC’s time consumption with varying macro-particle
numbers to provide readers with additional reference. It can be concluded from
Figure 14 that when the particle number is less than one percent of the CUDA
kernel limitation ($\sim$just under $10^{7}$), the computational time keeps
relatively short due to the high degree of parallelism on GPU and a small
number of data transfer between CPU and GPU. However, once the particle number
exceeds $10^{7}$, the time cost increase significantly. Note that the particle
number in JefiPIC (i.e., the occupied kernel number) must not exceed the CUDA
limitation ($\leavevmode\nobreak\ 10^{9}$), or the computation will fail.
Fortunately, $10^{7}\sim 10^{8}$ macro-particles are typically sufficient for
most simulations within the controlled noises.

Figure 14: Time cost of JefiPIC with varying particle numbers. The maximum
CUDA threads limits the maximum number of particles to $10^{9}$.

##  4 Conclusion

This paper purposes JefiPIC, a powerful plasma simulation package that uses a
3-D particle-in-cell method in conjunction with Jefimenko’s equations to
accurately model plasma systems. Despite its complexity, we have successfully
implemented JefiPIC on GPU for practical calculations. JefiPIC is more user-
friendly and easier to get started with, as it doesn’t require many
complicated amendment algorithms. This makes it an alternative option for
newcomers who want to dive into the field of particle-in-cell simulations.

Our comparative study between UNIPIC, EPOCH, and RGB-Maxwell has revealed that
JefiPIC has certain advantages in the following aspects. First, JefiPIC’s use
of integral operation enables accurately employing the linear particle-to-grid
apportionment technique, which simplifies the numerical scheme and brings less
noise. Secondly, JefiPIC naturally cuts off charge and field at the
computational boundary, preventing field reflection or charge deposition.
Third, JefiPIC does not require extra pre-processing, making it a superior
choice for calculating non-neutral plasma cases. Fourth, the integral
equations used in JefiPIC free the simulation from the CFL condition, enabling
larger time steps and helping to conserve computational resources. Finally,
the results obtained from JefiPIC are comparable to those obtained using
second-order accurate Boltzmann equations. Overall, our study demonstrates
that JefiPIC is a promising and time-efficient method for simulating plasma
physics problems with open boundary conditions.

However, to increase JefiPIC’s versatility, the next version will incorporate
interactions between electrons and neutral molecules or charged ions. Also,
measures will be taken to improve the efficiency of the integral operation or
to accomplish the parallel computing. We believe that these improvements will
extend the applicability of JefiPIC further and help advance our understanding
of plasma behavior.

We express our gratitude to Prof. Jian-Guo Wang, Associate Researcher Zai-Gao
Chen, and Assistant Researcher Ze-Ping Ren from NINT for their valuable
contributions and insightful discussions on the PIC simulations. The work is
supported by the National Key Research and Development Program of China under
Grant No. 2020YFA0709800 and the National Natural Science Foundation of China
(NSFC) under Grants No. 12105227.

##  Appendix A The unit converion table

In Flexible Unit (FU), all physical quantities have the dimension of energy E,
shown in Table 2. One needs to choose proper values for $\lambda$, $c$,
$\hbar$, $\varepsilon_{0}$ so that all numerical quantities are in the range
of float64 on GPU, where $c$, $\hbar$, $\varepsilon_{0}$ denote the speed of
light, the reduced Planck constant, and the vacuum permittivity. $\lambda$ is
constant relating the energy quantity between SI and FU.

Table 2: The transform of the physical quantities from SI Unit to Flexible Unit Physical quantity | SI Unit | Flexible Unit
---|---|---
Magnetic field | T | $5.01398\times 10^{-36}/(\lambda^{2}\hbar^{3/2}\varepsilon_{0}^{1/2}c^{5/2})\mathrm{E}^{2}$
Length | m | $3.16304\times 10^{25}\hbar c\lambda\mathrm{E}^{-1}$
Time | s | $9.48253\times 10^{33}\hbar\lambda\mathrm{E}^{-1}$
Electric charge | C | $1.89032\times 10^{18}(\hbar\varepsilon_{0}c)1/2$
Momentum | $\mathrm{kg\cdot m\cdot s^{-1}}$ | $2.99792\times 10^{8}/(c\lambda)E$
Energy | J | $1/\lambda E$
Mass | kg | $8.89752\times 10^{8}/(c^{2}\lambda)E$
Electric current | A | $1.99347\times 10^{-16}(\varepsilon_{0}c/\hbar)^{1/2}/\lambda E$
Electric field | $\mathrm{V\cdot m^{-1}}$ | $1.67249\times 10^{-44}/(\lambda^{2}\hbar^{3/2}\varepsilon_{0}^{1/2}c^{3/2})\mathrm{E}^{2}$
Force | $\mathrm{kg\cdot m\cdot s^{-2}}$ | $3.16153\times 10^{-26}/(\hbar\varepsilon_{0}c)\mathrm{E}^{2}$
Unit charge | 1 | $0.30286(\hbar\varepsilon_{0}c)^{1/2}$

## ——————

## References

  * [1] J. Dawson, One-dimensional plasma model, The Physics of Fluids, 5 (1962), 445-459.
  * [2] A. B. Langdon and C. K. Birdsall, Theory of plasma simulation using finite-size particles, The Physics of Fluids, 13 (1970) 2115-2122.
  * [3] C. K. Birdsall and D. Fuss, Clods-in-clouds, clouds-in-cells physics for many-body plasma simulation, Journal of Computational Physics, 3 (1969) 494-511.
  * [4] B. Marder, A method for incorporating Gauss’ law into electromagnetic PIC codes, Journal of Computational Physics, 68 (1987) 48-55.
  * [5] J. Villasenor, O. Buneman, Rigorous charge conservation for local electromagnetic field solvers, Computer Physics Communications, 69 (1992) 306-316.
  * [6] B. Wang, et al., Modern gyrokinetic particle-in-cell simulation of fusion plasmas on top supercomputers, The International Journal of High Performance Computing Applications, 33 (2019) 168-188.
  * [7] T. D. Arber, et al., Contemporary particle-in-cell approach to laser-plasma modelling, Plasma Physics and Controlled Fusion, 57 Article. 113001 (2015).
  * [8] J. Chen and J Zhang, Impact of geometrical parameters on SGEMP responses in cylinder model, Nuclear Engineering and Technology, 54 (2022) 3415-3421.
  * [9] H. Zhang, et al. Particle-in-cell simulations of low-pressure air plasma generated by pulsed x rays, Journal of Applied Physics, 130 Article. 173303 (2021).
  * [10] Z. Xu, et al., A code verification for the cavity SGEMP simulation code LASER-SGEMP, IEEE Transactions on Nuclear Science, 68 (2021) 1251-1257.
  * [11] Q. Lu, et al., Particle-in-cell simulations of whistler waves excited by an electron К distribution in space plasma, Journal of Geophysical Research Space Physics, 115 Article. A02213 (2010).
  * [12] J. Chen, et al., 2D planar PIC simulation of space charge limited current with geometrical parameters, varying temporal-profile and initial velocities, IEEE ACCESS, 10 (2022) 28499-28508.
  * [13] J. Wang, A megawatt-level surface wave oscillator in Y-band with larger oversized structure driven by annular relativistic electron beam, Scientific Reports, 8 Article. 6978 (2018).
  * [14] D. R. Welch, Fast hybrid particle-in-cell technique for pulsed-power accelerators, Physical review accelerators and beams, 23 Article. 110401 (2020).
  * [15] B. Goplen, et al., User-configurable MAGIC for electromagnetic PIC calculations, Computer Physics Communications, 87 (1995) 54-86.
  * [16] J. P. Verboncoeur, et al., An object-oriented electromagnetic PIC code, Computer Physics Communications, 87 (1995) 199-211.
  * [17] J. D. Blahovec, et al., 3-D ICEPIC simulations of the relativistic klystron oscillator, IEEE Transactions on Plasma Science, 28 (2000) 821-829.
  * [18] J. Wang, et al., UNIPIC code for simulations of high power microwave devices, Phys. Plasmas, 16 Article. 033108 (2009).
  * [19] C. P. Ridgers, et al., Modelling gamma-ray photon emission and pair production in high-intensity laser–matter interactions, Journal of Computational Physics, 260 (2014) 273-285.
  * [20] K. Yee, Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media, IEEE Transactions on Antennas and Propagation, 1 (1966) 302-307.
  * [21] D-Y. Na, Y. A. Omelchenko, H. Moon, B. V. Borges, F. L. Teixeira, Axisymmetric Charge-Conservative Electromagnetic Particle Simulation Algorithm on Unstructured Grids: Application to Microwave Vacuum Electronic Devices, Journal of Computational Physics, 346 (2017) pp. 295-317.
  * [22] D-Y. Na, H. Moon, Y. A. Omelchenko, and F. L. Teixeira, Relativistic extension of a charge-conservative finite element solver for time-dependent Maxwell-Vlasov equations, Physics of Plasmas, 25(1) Article. 013109 (2018).
  * [23] A. B. Langdon, On enforcing Gauss’ law in electromagnetic particle-in-cell codes, Computer Physics Communications 70 (1992) 447-450.
  * [24] J. P. Verboncoeur, Symmetric Spline Weighting for Charge and Current Density in Particle Simulation, Journal of Computational Physics, 174 (2001) 421-427.
  * [25] T. Zh. Esirkepov, Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor, Computer Physics Communications, 135 (2001) 144-153.
  * [26] J. W. Eastwood, The virtual particle electromagnetic particle-mesh method, Computer Physics Communications, 64 (1991) 252-266.
  * [27] Y. Wang, et al., Charge conserving emission technique for three-dimensional conform particle-in-cell simulations, High Power Laser and Particle Beams, 28 (2016) Article. 033020.
  * [28] O. D. Jefimenko, Electricity and magnetism: an introduction to the theory of electric and magnetic fields, Electret Scientific Co, Star City, W. Va, (1989).
  * [29] D. Griffiths, Introduction to electrodynamics, Prentice Hall, Upper Saddle River, N.J, (1999).
  * [30] J. Zhang, et al., JefiGPU: Jefimenko’s equations on GPU, Computer Physics Communications, 276 Article. 108328 (2022).
  * [31] H. Wu, et al., ZMCintegral: A package for multi-dimensional Monte Carlo integration on multi-GPUs, Computer Physics Communications, 248 Article. 106962 (2020).
  * [32] R. Courant, K. Friedrichs and H. Lewy, On the partial difference equations of mathematical physics, IBM Journal of Research and Development, 11 (1967) 215-234.
  * [33] J. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics, 114 (1994) 185-200.
  * [34] G. Mur, Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations, IEEE Transactions on Electromagnetic Compatibility, EMC-23(4) (1981) 377-382.
  * [35] M. F Pasik, D. B Seidel, and R. W Lemke, A Modified Perfectly Matched Layer Implementation for Use in Electromagnetic PIC Codes, Journal of Computational Physics, 148 (1999) 125-132.
  * [36] G. Jost, T.M. Tran, K. Appert, and S. Wüthrich, Effect of electron-cyclotron instabilities on gyrotron beam quality, Computer Physics Communications, 100 (1997) 47-55.
  * [37] J. Zhang, et al., Towards a full solution of the relativistic Boltzmann equation for quark-gluon matter on GPUs, Physical Review D, 102 Article. 074011 (2020).
  * [38] J. Zhang, et al., Charge-dependent directed flows in heavy-ion collisions by Boltzmann-Maxwell equations, Physical Review Research, 4 (2022) Article. 033138.
  * [39] O. Buneman and W. B. Pardo, Relativistic Plasmas, Benjamin, New York, (1968).
  * [40] B. Goplin, et al., Users Manual for MAGIC, Mission Research Corporation Report No. MRC/WDC-R-068, Alexandria, VA, (1983).
  * [41] C. K. Birdsall and A. B. Landon, Plasma Physics via Computer Simulation, Mc-Graw-hill, New York, (1985).
  * [42] G. G. Howes, Laboratory space physics: Investigating the physics of space plasmas in the laboratory, 25 Article. 055501 (2018).
  * [43] G. Peng, et al., Two typical collective behaviors of the heavy ions expanding in cold plasma with ambient magnetic field, Physics of Fluids, 33 (2021) Article. 076602.
  * [44] S. H. Brecht and D. J. Larson, The physics of ion decoupling in magnetized plasma expansions, Journal of Geophysical Research: Space Physics, 116 (2011) Article. A11310.
  * [45] J. Chen, et al., Simulation of SGEMP using particle-in-cell method based on conformal technique, IEEE Transactions on Nuclear Science, 66 (2019) 820-826.
  * [46] K.E. Lonngren, Expansion of an electron cloud, Physical Letters A, 59 (1976) 285-286.


11institutetext:  University of Copenhagen, Department of Computer Science,
Denmark,
<EMAIL_ADDRESS>22institutetext: Technical University of
Denmark, DTU Compute, Denmark,
<EMAIL_ADDRESS>
$\dagger$ Authors contributed equally

# Is segmentation uncertainty useful?

Steffen Czolbe †1†1 Kasra Arnavaz †1†1 Oswin Krause 11 Aasa Feragen 22

###### Abstract

Probabilistic image segmentation encodes varying prediction confidence and
inherent ambiguity in the segmentation problem. While different probabilistic
segmentation models are designed to capture different aspects of segmentation
uncertainty and ambiguity, these modelling differences are rarely discussed in
the context of applications of uncertainty. We consider two common use cases
of segmentation uncertainty, namely assessment of segmentation quality and
active learning. We consider four established strategies for probabilistic
segmentation, discuss their modelling capabilities, and investigate their
performance in these two tasks. We find that for all models and both tasks,
returned uncertainty correlates positively with segmentation error, but does
not prove to be useful for active learning.

###### Keywords:

Image segmentation Uncertainty quantification Active learning.

##  1 Introduction

00footnotetext: Code available at
github.com/SteffenCzolbe/probabilistic_segmentation

Image segmentation – the task of delineating objects in images – is one of the
most crucial tasks in image analysis. As image acquisition methods can
introduce noise, and experts disagree on ground truth segmentations in
ambiguous cases, predicting a single segmentation mask can give a false
impression of certainty. Uncertainty estimates inferred from the segmentation
model can give some insight into the confidence of any particular segmentation
mask, and highlight areas of likely segmentation error to the practitioner. It
adds transparency to the segmentation algorithm and communicates this
uncertainty to the user. This is particularly important in medical imaging,
where segmentation is often used to understand and treat disease.
Consequently, quantification of segmentation uncertainty has become a popular
topic in biomedical imaging [6, 11].

Training segmentation networks requires large amounts of annotated data, which
are costly and cumbersome to attain. Active learning aims to save the
annotator’s time by employing an optimal data gathering strategy. Some active
learning methods use uncertainty estimates to select the next sample to
annotate [10, 7, 9]. While several potential such data gathering strategies
exist [13, 16], a consistent solution remains to be found [8].

While several methods have been proposed to quantify segmentation uncertainty
[4, 6, 11], it is rarely discussed what this uncertainty represents, whether
it matches the user’s interpretation, and if it can be used to formulate a
data-gathering strategy. We compare the performance of several well-known
probabilistic segmentation algorithms, assessing the quality and use cases of
their uncertainty estimates. We consider two segmentation scenarios: An
unambiguous one, where annotators agree on one underlying true segmentation,
and an ambiguous one, where a set of annotators provide potentially strongly
different segmentation maps, introducing variability in the ground truth
annotation.

We investigate the degree to which the inferred uncertainty correlates with
segmentation error, as this is how reported segmentation uncertainty would
typically be interpreted by practitioners. We find that uncertainty estimates
of the models coincide with likely segmentation errors and strongly correlate
with the uncertainty of a set of expert annotators. Surprisingly, the model
architecture used does not have a strong influence on the quality of
estimates, with even a deterministic U-Net [12] giving good pixel-level
uncertainty estimates.

Second, we study the potential for uncertainty estimates to be used for
selecting samples for annotation in active learning. Reducing the cost of data
annotation is of utmost importance in biomedical imaging, where data
availability is fast-growing, while annotation availability is not. We find
that there are many pitfalls to an uncertainty-based data selection strategy.
In our experiment with multiple annotators, the images with the highest model
uncertainty were precisely those images where the annotators were also
uncertain. Labeling these ambiguous images by a group of expert annotators
yielded conflicting ground truth annotations, providing little certain
evidence for the model to learn from.

Figure 1: Segmentation uncertainty is often interpreted as probable
segmentation error, as seen near the lesion boundary in the first two
examples. In the third example, however, model bias leads to a very certain,
yet incorrect segmentation.

##  2 Modelling segmentation uncertainty

Image segmentation seeks to estimate a well-defined binary111For simplicity,
we consider binary segmentation; the generalization to multi-class
segmentation is straightforward. segmentation $g\colon\Omega\to\\{0,1\\}$ for
a discrete image domain $\Omega$. Typically, a predictive model
$h(\mathbf{x},\mathbf{w})$ with parameters $\mathbf{w}$, such as a neural
network, is fitted to binary annotation data $a\colon\Omega\to\\{0,1\\}$ by
minimizing a loss $\mathcal{L}(a,h(\mathbf{x},\mathbf{w}))$. Here,
$\mathbf{x}\in\mathbb{R}^{\Omega}$ is the image, and
$\mathbf{y}=h(\mathbf{x},\mathbf{w})$ defines an image of pixel-wise
segmentation probabilities, such as the un-thresholded softmax output of a
segmentation network $h$.

Typically, the annotation is assumed to be error-free, that is $a=g$, and
predictors are typically trained on a single annotation per image. We assume
that the trained neural network $h(\mathbf{x},\mathbf{w})$ satisfies

$h(\mathbf{x},\mathbf{w})=g(\mathbf{x})+b+err\enspace,$

where $b$ and $err$ denote bias and segmentation error. Segmentation
uncertainty is often interpreted as correlating with this error, although this
is primarily realistic for small bias. Such segmentation tasks are called
_unambiguous_ ; we consider a running example of skin lesion segmentation from
dermoscopic images [3, 15], where the lesion boundary is clearly visible in
the image (Fig. 1).

Recent work has considered _ambiguous_ segmentation tasks [6, 11], where there
is no accessible “ground truth” segmentation, either because the data is not
sufficient to estimate the segmentation, or because there is subjective
disagreement. Examples include lesions in medical imaging, where the boundary
can be fuzzy due to gradual infiltration of tissue, or where experts disagree
on whether a tissue region is abnormal or not.

In such tasks, we make no assumption on the underlying segmentation $g$ or the
errors $err$, but regard the observed annotations as samples from an unknown
“ground truth” distribution $p(a|\mathbf{x})$ over annotations $a$ conditioned
on the image $\mathbf{x}$. The goal of segmentation is to estimate the
distribution $p(a|\mathbf{x})$, or its proxy distribution
$p(\mathbf{y}|\mathbf{x})$ over pixel-wise class probabilities
$\mathbf{y}\colon\Omega\to[0,1]$, as accurately as possible for a given image
$\mathbf{x}$. If successful, such a model can sample coherent, realistic
segmentations from the distribution, and estimate their variance and
significance. As a running example of an ambiguous segmentation task, we
consider lung lesions [1, 2, 6]. For such tasks, predictors are typically
trained on multiple annotators, who may disagree both on the segmentation
boundary and on whether there is even an object to segment.

From the uncertainty modelling viewpoint, these two segmentation scenarios are
rather different. Below, we discuss differences in uncertainty modelling for
the two scenarios and four well-known uncertainty quantification methods.

##  3 Probabilistic Segmentation Networks

A probabilistic segmentation model seeks to model the distribution
$p(\mathbf{y}|\mathbf{x})$ over segmentations given an input image
$\mathbf{x}$. Here, our annotated dataset $(\mathbf{X},\mathbf{A})$ consists
of the set $\mathbf{X}$ of $N$ images
$\big{\\{}\mathbf{x}_{n}~{}|~{}n=1,...,N\big{\\}}$, and $L$ annotations are
available per image, so that $\mathbf{A}=\big{\\{}a_{n}^{(l)}\sim
p(\mathbf{y}|\mathbf{x}_{n})~{}|~{}(n,l)=(1,1),...,(N,L)\big{\\}}$.

Taking a Bayesian view, we seek the distribution

$p(\mathbf{y}|\mathbf{x},\mathbf{X},\mathbf{A})=\int
p(\mathbf{y}|\mathbf{x},\mathbf{w})p(\mathbf{w}|\mathbf{X},\mathbf{A},h)\,d\mathbf{w}\enspace,$
(1)

over segmentations $\mathbf{y}$ given image $\mathbf{x}$ and data
$(\mathbf{X},\mathbf{A})$, which can be obtained by marginalization with
respect to the weights $\mathbf{w}$ of the model $h$.

In most deep learning applications, our prior belief over the model $h$,
denoted $p(h)$, is modelled by a Dirac delta distribution indicating a single
architecture with no uncertainty. In the context of uncertain segmentation
models, however, we would like to model uncertainty in the parameters
$\mathbf{w}$. Denoting our prior belief over the parameters $\mathbf{w}$ by
$p(\mathbf{w}|h)$, Bayes’ theorem gives

$p(\mathbf{w}|\mathbf{X},\mathbf{A},h)=\frac{p(\mathbf{w}|h)p(\mathbf{A}|\mathbf{X},\mathbf{w},h)}{p(\mathbf{A}|\mathbf{X},h)},$
(2)

where the likelihood update function is given by

$p(\mathbf{A}|\mathbf{X},\mathbf{w},h)=\exp{\left(\sum_{n=1}^{N}\sum_{l=1}^{L}\mathbf{A}^{(l)}\log\left(h(\mathbf{x}_{n},\mathbf{w})\right)+(1-\mathbf{A}^{(l)})\log\left(1-h(\mathbf{x}_{n},\mathbf{w})\right)\right)}$

and normalizing constant

$p(\mathbf{A}|\mathbf{X},h)=\int
p(\mathbf{w}|h)p(\mathbf{A}|\mathbf{X},\mathbf{w},h)\,d\mathbf{w}\enspace.$

This integral is generally intractable, making it impossible to obtain the
proper posterior (2). Below, we discuss how empirical approximations
$\hat{p}(\mathbf{y}|\mathbf{x},\mathbf{X},\mathbf{A})$ to the distribution
$p(\mathbf{y}|\mathbf{x},\mathbf{X},\mathbf{A})$ found in (1) are performed in
four common segmentation models. Note that both $p$ and $\hat{p}$ can be
degenerate, depending on the number of annotations available and models used.

U-Net$h(\mathbf{x},\hat{\mathbf{w}})$Ensemble$h(\mathbf{x},\mathbf{w}^{(1)})$…$h(\mathbf{x},\mathbf{w}^{(M)})$$h(\mathbf{x},\mathbf{w}^{(r)})$$\mathbf{w}^{(r)}\sim
p(\mathbf{w};\theta)$MC-
Dropout$h(\mathbf{x},\mathbf{z}^{(i)},\hat{\mathbf{w}})$$\mathbf{z}^{(i)}\sim
p_{\text{prior}}(\mathbf{z}|\mathbf{x})$Prob. U-Net Figure 2: Schematic
overview (adapted from [6]) of the evaluated models. Blue: residual blocks .
Orange: Dropout layers essential to the networks’ functionality.

#### U-Net with softmax output.

The well established U-Net [12] architecture with a softmax output layer
yields class-likelihood estimates. As the model is deterministic,
$p(h|\mathbf{X},\mathbf{A})$ is degenerate. Parameters are selected by a
maximum a posteriori (MAP) estimate i.e.
$p(\mathbf{w}|\mathbf{X},\mathbf{A},h)\approx\delta(\mathbf{w}-\hat{\mathbf{w}})$
in which $\hat{\mathbf{w}}=\text{argmax
}p(\mathbf{w}|\mathbf{X},\mathbf{A},h)$. The model output (1) is approximated
by the degenerate distribution
$\hat{p}(\mathbf{y}|\mathbf{x},\mathbf{X},\mathbf{A})\approx
p(\mathbf{y}|\mathbf{x},\hat{\mathbf{w}})$. The softmax output layer predicts
a pixel-wise class probability distribution
$p(\mathbf{y}_{(i,j)}|\mathbf{x},\mathbf{X},\mathbf{A})$. As no co-variance or
dependencies between pixel-wise estimates are available, segmentation masks
sampled from the pixel-wise probability distributions are often noisy [11]. An
alternative approach followed by our implementation is the thresholding of
pixel-wise probability values, which leads to a single, coherent segmentation
map.

#### Ensemble methods

combine multiple models to obtain better predictive performance than that
obtained by the constituent models alone, while also allowing the sampling of
distinct segmentation maps from the ensemble. We combine $M$ U-Net models
$h(\mathbf{x},\mathbf{w}^{(m)})$ where, if labels from multiple annotators are
available, each constituent model is trained on a disjoint label set
$\mathbf{A}^{(m)}$. When trained on datasets with a single label, all
constituent models are trained on the same data and their differences stem
from randomized initialization and training. Treating the models as samples,
we obtain an empirical distribution approximating (1) by drawing from the
constituent models at random.

#### Monte-Carlo Dropout

[4] is a Bayesian approximation technique based on dropout, where samples from
the posterior over dropout weights give a better approximation of the true
posterior than a MAP estimation. Given a selected model $h$, one can
approximate (1) as
$\hat{p}(\mathbf{y}|\mathbf{x},\mathbf{X},\mathbf{A})\approx
1/R\sum_{r=1}^{R}p(\mathbf{y}|\mathbf{x},\mathbf{w}^{(r)})$ when
$\mathbf{w}^{(r)}\sim p(\mathbf{w}|\mathbf{X},\mathbf{A},h)$. Since
$p(\mathbf{w}|\mathbf{X},\mathbf{A},h)$ is intractable, it is approximated [4]
by a variational distribution $p(\theta)$ as $\theta_{i}=w_{i}\cdot
z_{i},z_{i}\sim\text{Bernoulli}(p_{i})$, where $p_{i}$ is the probability of
keeping the weight $w_{i}$ in a standard dropout scheme.

#### The Probabilistic U-Net

[6] fuses the output of a deterministic U-Net with latent samples from a
conditional variational auto-encoder modelling the variation over multiple
annotators. Test-time segmentations are formed by sampling a latent
$\mathbf{z}$, which is propagated with the image through the U-Net.
Predictions are made as
$\hat{p}(\mathbf{y}|\mathbf{x},\mathbf{X},\mathbf{A})\approx
p(\mathbf{y}|\mathbf{x},\mathbf{z}^{(i)},\hat{\mathbf{w}})$, with
$\mathbf{z}^{(i)}\sim p_{\text{prior}}(\mathbf{z}|\mathbf{x})$.

Skin LesionsLung CancerImageU-NetEnsembleMC-DropoutProb. U-Net Figure 3:
Segmentation Uncertainty. Samples from the test set of the two datasets.
Images in row one, model uncertainty (entropy) heat-maps in rows 2-5. Outline
of mean ground truth annotations in Blue, mean model predictions in Orange.

##  4 Experiments

###  4.1 Data

Practical applications of uncertainty in segmentation tasks differ both in the
type of ambiguity, and the availability of expert annotations. We select two
representative datasets for our evaluation.

The ISIC18 dataset consists of skin lesion images with a single annotation
available [3, 15], and is used as an example of unambiguous image
segmentation. We rescale the images to $256\times 256$ pixels and split the
dataset into $1500$ samples for the train-set and $547$ each for the
validation and test sets.

The LIDC-IDRI lung cancer dataset [1, 2] contains 1018 lung CT scans from 1010
patients. For each scan, 4 radiologists (out of 12) annotated abnormal
lesions. Anonymized annotations were shown to the other annotators, who were
allowed to adjust their own masks. Significant disagreement remains between
the annotators: Among the extracted patches where at least one annotator
marked a lesion, an average of $50\%$ of the annotations are blank. We pre-
processed the images as in [6], resampled to $0.5$mm resolution, and cropped
the CT-slices with lesions present to $128\times 128$ pixels. The dataset is
split patient-wise into three groups, $722$ for the training-set and $144$
each for the validation and test sets.

###  4.2 Model tuning and training

To allow for a fair evaluation, we use the same U-Net backbone of four encoder
and decoder blocks for all models. Each block contains an up-/down-sampling
layer, three convolution layers, and a residual skip-connection. The ensemble
consists of four identical U-Nets. The latent-space encoders of the
probabilistic U-Net are similar to the encoding branch of the U-Nets, and we
choose a six-dimensional latent space size, following the original paper’s
recommendation.

All models were trained with binary cross-entropy. The probabilistic U-Net has
an additional $\beta$-weighted KL-divergence loss to align the prior and
posterior distributions, as per [6]. The optimization algorithm was Adam, with
a learning rate of $10^{-4}$ for most models, except the probabilistic U-Net
and MC-Dropout models on the skin lesion dataset, where a lower learning rate
of $10^{-5}$ gave better results. We utilized early stopping to prevent over-
fitting, and define the stopping criteria as 10 epochs without improvement of
the validation loss, 100 epochs for models trained with the reduced learning
rate. For the MC-Dropout and probabilistic U-Net models we performed a hyper-
parameter search over the dropout probability $p$ and the loss function
weighting factor $\beta$, selecting the configuration with the lowest
generalized energy distance on the validation set. We arrived at $p=0.5$,
$\beta=0.0005$.

(a) Skin Lesion Dataset

(b) Lung Cancer Dataset

Figure 4: Pixelwise uncertainty by prediction correctness (True Positive,
False Positive, False Negative, True Negative). The scatter plot shows
individual pixels, with the median circled. For the lung cancer dataset, we
discarded pixels with annotator disagreement.

Figure 5: Pixel-wise model uncertainty on the lung cancer dataset, grouped by
agreement of expert annotations. Experts agree: $H(p)=0$, somewhat agree
$0<H(p)<1$, disagree $H(p)=1$.
Figure 6: Generalized Energy Distance of models on the lung cancer dataset,
approximation by 1 to 16 samples, median highlighted. Lower distances are
better.

###  4.3 Uncertainty Estimation

For all models, our uncertainty estimates are based on non-thresholded pixel-
wise predictions. For the U-Net, we take the final softmax predictions; for
the remaining models we average across $16$ non-thresholded samples. We
quantify the pixel-wise uncertainty of the model by the entropy

$\mathrm{H}(p(\mathbf{y}_{(i,j)}|\mathbf{x},\mathbf{X},\mathbf{A}))=\sum_{c\in
C}p(\mathbf{y}_{(i,j)}=c|\mathbf{x},\mathbf{X},\mathbf{A})\log_{2}\frac{1}{p(\mathbf{y}_{(i,j)}=c|\mathbf{x},\mathbf{X},\mathbf{A})}$

with $p(\mathbf{y}_{(i,j)}=c|\mathbf{x})$ as the pixel-wise probability to
predict class $c\in C$. We plot the resulting uncertainty map for random
images $\mathbf{x}$ from both datasets in Fig. 3. For visual reference, we
overlay the mean expert annotation in Blue, and the mean model prediction in
Orange. Darker shades indicate higher uncertainty.

We quantitatively assess the quality of uncertainty estimates by examining
their relation to segmentation error in Fig. 4. On both datasets, models are
more certain when they are correct (true positive, true negative) compared to
when they are incorrect (false positive, false negative). A repeated measure
correlation test finds a significant $(\alpha=0.01)$ correlation between
segmentation error and model uncertainty on both datasets, for all methods.
The relation holds, but is less strong, for MC-dropout on the skin dataset,
which retains high uncertainty even when it is correct. On the lung cancer
dataset, all models have high uncertainty on true positive predictions. This
might be caused by the imbalance of the dataset, where the positive class is
strongly outweighed by the background and annotators often disagree. We tried
training the models with a class-occurrence weighted loss function, which did
produce true positive predictions with higher certainty but suffered an
overall higher segmentation error.

We assess the correlation of model uncertainty with the uncertainty of the
annotators on the lung cancer dataset in Fig. 5. For all models, this
correlation is significant $(\alpha=0.01)$ . The median model uncertainty is
very low ($<0.1$) when all annotators agree, but high ($>0.7$) when they
disagree. There is a minor difference in model uncertainty between partial
agreement (annotators split 3 – 1) and full disagreement (annotators split 2 –
2).

ImageAnnotationImageAnnotationsSamplesSamplesU-
Net$0.00$$0.00$Ensemble$0.10$$0.27$MC-Dropout$0.26$$0.39$Prob.
U-Net$0.07$$0.43$ Figure 7: Samples from the probabilistic models. First row:
Image and ground truth annotations from the skin dataset (left) and lung
nodule dataset (right). Following rows: samples
$\mathbf{y}\sim\hat{p}(\mathbf{y}|\mathbf{x},\mathbf{X},\mathbf{A})$ drawn
from the various models. Sample diversity over the entire dataset shown next
to the model name.

###  4.4 Sampling Segmentation Masks

Fig. 7 shows segmentation masks $\mathbf{y}$ sampled from the trained models
$\hat{p}(\mathbf{y}|\mathbf{x})$. The U-Net model is fully deterministic and
does not offer any variation in samples. The sample diversity of the ensemble
is limited by the number of constituent models (four in our experiment). The
MC-Dropout and probabilistic U-Net allow fully random sampling and achieve a
visually higher diversity. On the skin lesion dataset, where only one export
annotation per image is available, models still produce diverse predictions.
On the lung cancer dataset, samples from the MC-Dropout and probabilistic
U-Net represent the annotator distribution well.

We measure the distance between the model distribution
$\hat{p}(\mathbf{y}|\mathbf{x},\mathbf{X},\mathbf{A})$ and the annotator
distribution $p(\mathbf{y}|\mathbf{x})$ with the Generalized Energy Distance
[6, 11, 14]. The distance measure is calculated as

$D^{2}_{GED}(p,\hat{p})=2\mathbb{E}_{y\sim
p,\hat{y}\sim\hat{p}}\left[d(y,\hat{y})\right]-\mathbb{E}_{y,y^{\prime}\sim
p}\left[d(y,y^{\prime})\right]-\mathbb{E}_{\hat{y},\hat{y}^{\prime}\sim\hat{p}}\left[d(\hat{y},\hat{y}^{\prime})\right]\enspace.$
(3)

We use $1-\text{IoU}(\cdot,\cdot)$ as the distance $d$. A low $D^{2}_{GED}$
indicates similar distributions of segmentations. We approximate the metric by
drawing up to 16 samples from both distributions, and sample with replacement.
The results are shown in Fig. 6. We observe that the annotator distribution is
best approximated by the probabilistic U-Net, with MC-dropout and Ensemble
closely behind; these pairwise ranks are significant $(\alpha=0.01)$ with
left-tailed t-tests. A deterministic U-Net architecture is not able to
reproduce the output distribution. Our results are consistent with [6],
verifying our implementation. Following [11], we use the last term of (3) to
assess the diversity of samples drawn from the model and note them in Fig. 7.
They reinforce the qualitative observations of sample diversity.

###  4.5 Uncertainty estimates for active learning

Instead of training the models with all available data
$\\{\mathbf{X},\mathbf{A}\\}$, we now start with a small random subset
$\\{\mathbf{X}_{0},\mathbf{A}_{0}\\}$. We train the model with this subset at
iteration $t=0$, and then add a set of $k$ images from
$\\{\mathbf{X},\mathbf{A}\\}$ to form
$\\{\mathbf{X}_{t+1},\mathbf{A}_{t+1}\\}$. Samples are selected based on the
sum of pixel-wise entropies [7]. We repeat for $T$ iterations, benchmarking
against a random sample selection strategy.

Figure 8: Learning curves for the four algorithms on both datasets. Note that
the Probabilisitic U-net only applies to the ambiguous segmentation.

For both skin lesion and lung cancer datasets, we start with a training size
of 50 images, add $k=25$ images at each iteration, and repeat $T=10$ times.
The models are trained for 5000 gradient updates with a batch size of 16 and
32 for the respective datasets. Since annotations are costly and to speed up
computations, no validation-loss based early stopping is used. The
experimental setup has been picked to ensure meaningful model uncertainties
for the data selection policy and to ensure convergence within each active
learning iteration.

The learning curves in Fig. 8 show that random-based sampling leads to a
faster reduction in test loss over the uncertainty-based sampling strategy for
both datasets. We further investigated the samples selected by the
uncertainty-based strategy by looking at the images which caused a large
increase in the test error. One such image is shown in Fig. 9.

ImageAnnotator 1Annotator 2Annotator 3Annotator 4 Figure 9: An example of the
ambiguous samples frequently selected for inclusion into the training set
under the uncertainty-based data gathering strategy. This unseen sample was
selected when 150 annotations were revealed. The group of expert annotators
provided disagreeing segmentation masks, confirming the model uncertainty but
providing little additional information to learn from.

##  5 Discussion & Conclusion

Our results in Fig. 4 show that there is a clear relation between uncertainty
estimates and segmentation error. The examples in Fig. 3 further highlight
that areas of high uncertainty are not merely distributed around class
boundaries, but also encompass areas with ambiguous labels. Fig. 5 shows that
the uncertainty estimates obtained from the model are a good representation of
the uncertainty of a group of expert annotators. We conclude that pixel-wise
model uncertainty estimates give the practitioner a good indication of
possible errors in the presented segmentation mask, allowing those predictions
to be examined with care.

The learning curves in Fig. 8 show that estimated uncertainty is not generally
useful for selecting active learning samples, for any model or dataset. Our
results depend on using the sum of pixel-wise entropies as a per-image
entropy, which is correct for the softmax model, but only an approximation for
the other models. This might impact our results. For the Lung Cancer dataset,
all models estimate high uncertainty for the positive class, and the active
learner thus selects images with a large foreground, skewing the proportion of
classes represented in the training set. Furthermore, the selected images
often have high annotator disagreement, illustrated in Fig. 9. If the active
learner prefers sampling ambiguous images, it will be presented with
inconsistent labels leading to harder learning conditions and poor
generalisation. This may stem from an incorrect active learning assumption
that annotations are noise-free and unambiguous, which is often not true. In
conclusion, for a fixed budget of annotated images, we find no advantage in
uncertainty-based active learning.

We observed similar behaviour of pixel-wise uncertainty estimates across all
four segmentation models. The models differ in their ability to generate a
distribution of distinct and coherent segmentation masks, with only the MC-
dropout and probabilistic U-Net offering near unlimited diversity (see Fig.
7). But these models are harder to implement, more resource-intensive to
train, and require hyperparameter tuning. The choice of model is ultimately
application dependent, but our experiments show that even a simple U-net is
competitive for the common task of assessing segmentation error. This agrees
with [5], which compared uncertainty quantification models for unambiguous
segmentation.

Our division of segmentation tasks into ambiguous and unambiguous considers it
as ”unambiguous” when a fundamentally ambiguous segmentation task is covered
by a single annotator – or potentially several annotators, but with only one
annotator per image, as for the Skin Lesion dataset. Even if the underlying
task _is_ ambiguous, the models considered in this paper inherently assume
that it is _not_ , as there is no mechanism to detect annotator variance when
every image is only annotated once. More fundamental modelling of segmentation
ambiguity and uncertainty thus remains a highly relevant open problem.

To conclude – is segmentation uncertainty useful? We find that uncertainty,
even in the simplest models, reliably gives practitioners an indication of
areas of an image that might be ambiguous, or wrongly segmented. Using
uncertainty estimates to reduce the annotation load has proven challenging,
with no significant advantage over a random strategy.

#### Acknowledgements.

Our data was extracted from the “ISIC 2018: Skin Lesion Analysis Towards
Melanoma Detection” grand challenge datasets [3, 15]. The authors acknowledge
the National Cancer Institute and the Foundation for the National Institutes
of Health, and their critical role in the creation of the free publicly
available LIDC/IDRI Database used here. This work was funded in part by the
Novo Nordisk Foundation (grants no. NNF20OC0062606 and NNF17OC0028360) and the
Lundbeck Foundation (grant no. R218-2016-883).

## References

  * [1] Armato III, S.G., McLennan, G., Bidaut, L., McNitt-Gray, M.F., Meyer, C.R., Reeves, A.P., Zhao, B., Aberle, D.R., Henschke, C.I., Hoffman, E.A., et al.: The lung image database consortium (LIDC) and image database resource initiative (IDRI): a completed reference database of lung nodules on CT scans. Medical physics 38(2), 915–931 (2011)
  * [2] Clark, K., et al.: The cancer imaging archive (TCIA): maintaining and operating a public information repository. Journal of digital imaging 26(6), 1045–1057 (2013)
  * [3] Codella, N.C., Gutman, D., Celebi, M.E., Helba, B., Marchetti, M.A., Dusza, S.W., Kalloo, A., Liopyris, K., Mishra, N., Kittler, H., et al.: Skin lesion analysis toward melanoma detection: A challenge at the 2017 international symposium on biomedical imaging (ISBI), hosted by the international skin imaging collaboration (ISIC). In: 2018 IEEE 15th International Symposium on Biomedical Imaging. pp. 168–172
  * [4] Gal, Y., Ghahramani, Z.: Dropout as a Bayesian approximation: Representing model uncertainty in deep learning. In: international conference on machine learning. pp. 1050–1059 (2016)
  * [5] Jungo, A., Reyes, M.: Assessing reliability and challenges of uncertainty estimations for medical image segmentation. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 48–56. Springer (2019)
  * [6] Kohl, S., Romera-Paredes, B., Meyer, C., De Fauw, J., Ledsam, J.R., Maier-Hein, K., Eslami, S., Jimenez Rezende, D., Ronneberger, O.: A probabilistic u-net for segmentation of ambiguous images. Advances in neural information processing systems 31, 6965–6975 (2018)
  * [7] Lewis, D.D., Catlett, J.: Heterogeneous uncertainty sampling for supervised learning. In: Machine learning proceedings 1994, pp. 148–156. Elsevier (1994)
  * [8] Loog, M., Yang, Y.: An empirical investigation into the inconsistency of sequential active learning. In: 2016 23rd international conference on pattern recognition (ICPR). pp. 210–215. IEEE (2016)
  * [9] MacKay, D.J.: The evidence framework applied to classification networks. Neural computation 4(5), 720–736 (1992)
  * [10] MacKay, D.J.: Information-based objective functions for active data selection. Neural computation 4(4), 590–604 (1992)
  * [11] Monteiro, M., Le Folgoc, L., Coelho de Castro, D., Pawlowski, N., Marques, B., Kamnitsas, K., van der Wilk, M., Glocker, B.: Stochastic segmentation networks: Modelling spatially correlated aleatoric uncertainty. In: Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M.F., Lin, H. (eds.) Advances in Neural Information Processing Systems. vol. 33, pp. 12756–12767 (2020)
  * [12] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [13] Settles, B.: Active learning literature survey. Tech. rep., University of Wisconsin-Madison Department of Computer Sciences (2009)
  * [14] Székely, G.J., Rizzo, M.L.: Energy statistics: A class of statistics based on distances. Journal of statistical planning and inference 143(8), 1249–1272 (2013)
  * [15] Tschandl, P., Rosendahl, C., Kittler, H.: The HAM10000 dataset, a large collection of multi-source dermatoscopic images of common pigmented skin lesions. Scientific data 5, 180161 (2018)
  * [16] Yang, Y., Loog, M.: A benchmark and comparison of active learning for logistic regression. Pattern Recognition 83, 401–415 (2018)