# CS5260 Project: Unsupervised Data Augmentation for Object Detection

Zhang Yichen
School of Computing
National University of Singapore
<EMAIL_ADDRESS>
&Song Zeyang
School of Computing
National University of Singapore
<EMAIL_ADDRESS>
Li Wenbo
School of Computing
National University of Singapore
<EMAIL_ADDRESS>

###### Abstract

Data augmentation has always been an effective way to overcome overfitting
issue when the dataset is small. There are already lots of augmentation
operations such as horizontal flip, random crop or even Mixup. However, unlike
image classification task, we cannot simply perform these operations for
object detection task because of the lack of labeled bounding boxes
information for corresponding generated images. To address this challenge, we
propose a framework making use of Generative Adversarial Networks(GAN) to
perform unsupervised data augmentation. To be specific, based on the recently
supreme performance of YOLOv4, we propose a two-step pipeline that enables us
to generate an image where the object lies in certain position. In this way,
we can accomplish the goal that generating an image with bounding box label.
However, because of the lack of time and device, we are still working on its
experiments.

##  1 Introduction

Computer Vision is one of the most important areas in Artificial Intelligence.
It contains several basic tasks: image classification 2015arXiv151203385H ,
object detection 2013arXiv1311.2524G , instance segmentation
2017arXiv170306870H  and semantic segmentation 2014arXiv1412.7062C . Among
them, object detection has attracted widespread attention from researchers all
over the world. It solves two main problems in real world applications: "What
is the object?" and "Where is the object?". Current methods are simulating
human vision and cognition in specific application scenarios, such as
pedestrian detection, face detection and text detection. In recent years, with
the development of deep learning technology, significant breakthroughs have
been made in object detection, and it has also been applied to more fields
such as automatic driving, robot vision, and video surveillance.

Although object detection has developed rapidly and new algorithms have
constantly refreshing the scores of major rankings, one of its shortcomings is
its difficulty in obtaining labeled data compared to classification task
because of the tedious bounding box annotations process. A straight forward
solution to this is generating pseudo label from trained detection model with
confidence directly. However, it seems useless or even harmful since it’s
possible that the trained model will make mistake even if it gives out high
confidence value.

During the survey, we found that some methods start using Generative
Adversarial Network (GAN) for data augmentation recently, which achieve better
results than traditional methods. For example, DAGAN antoniou2018data  succeed
in generating augmented images for human face as shown in Figure 2, where only
the left top one is real image. And after training on those augmented images,
their model perform much better compared to those methods only using the real
images. Inspired by these methods, we raise an enlightening question about
this.

  * •

What if we can design a framework about GAN to generate augmented image
together with the bounding box labels?

To answer the above question, we develop a framework aiming to perform data
augmentation for object detection based on the success of YOLOv4 and GAN. To
be specific, we propose a two-step framework. In first step, our goal is to
train a GAN for normal data augmentation based on original dataset, i.e., be
able to generate real images. In second step, our goal is further extended to
generate images where object should lies in certain position. And we use
YOLOv4 to localize the position of object in the generated images and thus
instruct the training of generator. With this model, the only thing we need to
do is providing a position to the generator, then we can get the generated
images that meet the requirement. Further more, we can use this technique to
generate some hard samples, such as images with small or bias object. With
these hard samples, model may become more robust after training.

In summary, the contributions of this paper are three folds:

  * •

Propose a novel data augmentation framework to solve the essential problem of
lacking labeled sample in object detection.

  * •

Encode the position constraint into the generated images.

  * •

Provide comprehensive analysis and discussion on our results and experiment
process, which should be useful to other people who are also new to this area.

To our best knowledge, this is the first work that applying GAN to generate
the whole images directly based on a pure condition for object detection task.

##  2 Related Work

###  2.1 Object Detection

Most of the early object detection methods were constructed based on manual
features. Due to the lack of effective image representation at that time,
people could only design complex feature representations. A well-known method
is the Viola Jones detector 990517 . The VJ detector uses a sliding window to
view all possible window sizes and positions in an image. With the performance
of manual feature selection technology becoming saturated, object detection
has reached a stable development period after R. Girshick et al. proposed a
region with CNN features (RCNN) for object detection in 2014 7112511 . Since
then, object detection has developed at an unprecedented speed. In the era of
deep learning, object detection can be divided into two categories: "two-stage
detection" and "one-stage detection".

The former defines the detection frame as a "from coarse to fine" process, the
most typical methods are RCNN vandeSandeICCV2011 , Fast RCNN
2015arXiv150408083G , Faster RCNN NIPS2015_14bfa6bb , and later also proposed
a variety of improvements, including RFCN and Light head RCNN.

The latter defines it as "one step in place", and typical methods include
YOLO, SSD 2015arXiv151202325L , and RetinaNet 2017arXiv170802002L . We used
the YOLOv4 2020arXiv200410934B  framework this time, the backbone used is the
improved CSPDarknet 2019arXiv191111929W  and the new feature map usage method
BiFPN 2019arXiv191109070T  is used, both of which are robust methods in the
field of object detection.In Figure 2, We can see that even difficult sample
containing small objects, YOLOv4 can still provide excellent results. So we
can rely it on generating pseudo bounding box labels to some degree.

Figure 1: YOLOv4 Detection Result

Figure 2: Outcome of DAGAN

###  2.2 Generative Adversarial Network

GAN was first proposed by Ian Goodfellow in 2014. By introducing the concept
of adversarial learning, generator and discriminator trying their best to
compete with each other and gradually both become better during this process.
Because of its success in generating samples, GAN has been paid more and more
attention by researchers in recent years and has been applied to different
fields such as style transfer, portrait generation (StyleGAN
2018arXiv181204948K ), few-shot, one-shot, and zero-shot learning (AFHN
2020arXiv200313193L ).

GAN usually consists of a generator and a discriminator, where the generator
needs to generate fake samples to fool the discriminator while the
discriminator needs to have the ability to distinguish between real and fake
samples.

However, GAN has suffered from lots of problems such as model collapse,
difficulty in finding a balance between generator and discriminator, and the
lack of variety in generated samples. Many papers have tried to solve the
problem, but the effect is not satisfactory and the problem has not been
completely solved.

To solve the balance problem between generator and discriminator, which may
cause model collapse, Wasserstein GAN (WGAN) arjovsky2017wasserstein
introduce gradient penalty and Wasserstein distance, while the latter can also
act as the indicator for the training progress. It also has some other simple
modifications in GAN’s structure such as removing logarithm from the loss
function, removing the final Sigmoid layer of discriminator and outputs the
score of input samples rather than probability. All these simple changes
largely solve the above mentioned problems of GAN.

###  2.3 Data Augmentation

In object detection, traditional data augmentation methods include random
crop, changing picture size, changing picture tone, and random clip. Based on
the idea of traditional methods, this year there are also many outstanding
improvement methods such as mosic 2020arXiv200410934B , mixup
2017arXiv171009412Z , roimix 2019arXiv191103029L . The above-mentioned methods
all transform the original image to the augmented data. They can be well
explained, but most of them are not effective and cannot generate some unseen
images.

Besides these traditional methods, people also use deep learning methods, such
as GAN, for data augmentation. ASDN2017arXiv170403414W  use masks in the
feature layer to block some parts and generate images without some important
features. This method improves the network’s utilization of various features.
MTGANBai_2018_ECCV  use data augmentation for small objects detection, and
they use GAN to expand the resolution for small objects with low pixels. This
can improve the recognition accuracy of small objects, but both methods do not
expand the dataset. AugGanhuang2018auggan  converts images in daytime to other
time like night, which could improve generalization ability of object
detection model. However, this model requires semantic segmentation
information of input samples, which is hard to acquire in ordinary datasets.

##  3 Methodology

In this section, we will talk about our overall framework and its possible
useful application.

###  3.1 Overall Framework

Figure 3: Framework

Our model is proposed to generate visually plausible images with objects in
certain position. For this purpose, we design a 2-step training framework as
shown in Figure 3.

#### Step 1

In the first step, in order to generate indistinguishable fake images, we
first train a GAN which is able to translate noise to the target image domain.
In this part, we adopt Wasserstein GAN structure for the framework.

In this framework, there are two main components: generator and discriminator.
For generator, we used a network following the general architecture of
U-Netronneberger2015unet  as shown in Figure 4). Then for the discriminator, a
straight forward convolutional neural network is applied to predict the score
of input samples. Note that CNN model is able to get a probability for each
patch of input images. For example, if the size of input image is $256\times
256\times 3$, and the size of patch for discriminator is $8\times 8$, then the
output of discriminator will be $32\times 32\times 1$, representing the score
for each patch.

Figure 4: Unet structure used in our generator

For weight update of the GAN model, we used the similar way as WGAN, that is
minimizing the Wasserstein distance between real distance and generated
distance. In this model, a function $f_{w}$ with parameters $w$ are defined to
construct the loss:

$\mathcal{L}_{D}=E_{x~{}P_{r}}[f_{w}(x)]-E_{x~{}P_{g}}[f_{w}(x)]$ (1)

where the $P_{r}$ is the distribution of real samples and $P_{g}$ is the
distribution of generated samples, and $f_{w}$ is K-Lipschitz function. And
the generator loss $L_{G}$ and discriminitor loss $L_{D}$ are:

$\displaystyle\mathcal{L}_{G}$ $\displaystyle=-E_{x~{}P_{g}}[f_{w}(x)]$ (2)
$\displaystyle\mathcal{L}_{D}$
$\displaystyle=E_{x~{}P_{g}}[f_{w}(x)]-E_{x~{}P_{r}}[f_{w}(x)]$ (3)

By minimizing $L_{G}$ and $L_{D}$ iteratively, we then can get a model that
could generate visually plausible images.

#### Step 2

Then in the next step, the GAN is further trained to generate images with
target object conditioned on input position. A sequence of coordinate
$[x1,y1,x2,y2]$ is first converted to a mask, where the corresponding pixels
are set to one while the rest are zero. Then we fed it into an encoder and
concatenate its feature with noise as the input of generator.

For the generated images, if they get relatively high score from
discriminator, which means they are real enough and could be used for object
detection, thus we can use a trained detection model to get the position of
object from the generated images, compare it with the input bounding box data,
and calculate an object position loss to generator as the generation
constrain.

To achieve this, we applied a loss function similar to YOLO. It can be divided
into three parts: bounding box regression loss, confidence loss and
classification loss.

We applied CIoU loss zheng2020distance  as bounding box regression loss:

$\mathcal{L}_{CIoU}=1-IoU+\frac{\rho^{2}(B,B^{gt})}{c^{2}}+\alpha v$

where $IoU$ is the naive Intersection over Union(IoU) metrics, $B$ and
$B^{gt}$ denote the central points of detected bounding central points of
detected bounding box and target bounding box, $\rho(\cdot)$ is the Euclidean
distance and $x$ is the diagonal length of the smallest enclosing box covering
the two boxes. $v$ refers to width-height consistency ratio,

$v=\frac{4}{\pi^{2}}(arctan\frac{w^{gt}}{h^{gt}}-arctan\frac{w}{h})^{2}$

and $\alpha$ is a trade-off parameter defined as:

$\alpha=\frac{v}{(1-IoU)+v}$

The confidence loss is used to determine whether there are objects in the
prediction box, which is defined as below:

$\mathcal{L}_{conf}=\sum_{i=0}^{S^{2}}\sum_{j=0}^{B}1_{ij}^{\mathrm{obj}}\left(C_{i}-\hat{C}_{i}\right)^{2}+\lambda_{\text{noobj
}}\sum_{i=0}^{S^{2}}\sum_{j=0}^{B}1_{ij}^{\text{noobj
}}\left(C_{i}-\hat{C}_{i}\right)^{2}$ (4)

The classification loss is used to determine whether the object in the
prediction frame is the predicted object type:

$\mathcal{L}_{clf}=\sum_{i=0}^{S^{2}}1_{i}^{\mathrm{obj}}\sum_{c\in\text{
classes }}\left(p_{i}(c)-\hat{p}_{i}(c)\right)^{2}$ (5)

Thus, the overall loss function is a linear combination between these three
and can be written as follow:

$\mathcal{L}=\alpha\mathcal{L}_{CIoU}+\beta\mathcal{L}_{conf}+\theta\mathcal{L}_{clf}$
(6)

###  3.2 Application

After training, we can feed the object position into this model to get
augmented image together with the bounding box information. Moreover, we can
generate some difficult training samples where object is small or object lies
in the corner.

With these generated images, we can construct an augmented dataset consists of
both original images and generated images from our generator and retrain the
object detection model on this augmented dataset.

##  4 Experiments

###  4.1 Experiment Setting

#### Dataset

We use VOC2007 as our dataset. It has totally 9k images where each image is
labeled with several bounding boxes. We resize those images to $64\times
64\times 3$ before feeding into the discriminator. Given the coordinate of
bounding box, we also generate mask of size $64\times 64\times 1$ to match the
size of input images and then feed into the generator together. We will
discuss the reason that we choose this setting in the Section 4.3.

#### GAN Model

We adopt Unet generator and N-Layer discriminator from Pix2Pix. We also adopt
gradient penalty as WGAN-GP because of its advantages in both performance and
stability. Generator and discriminator are trained iteratively. The goal of
training generator is to generate indistinguishable samples for discriminator
while the goal of training discriminator is to make sure it recognize real
samples correctly.

#### Detection Model

We adopt YOLOv4 as our detection model, which is both fast and accurate. As is
observed from our experiments, the YOLOv4 model pretrained on COCO dataset is
also suitable for detection on VOC2007 dataset since 20 classes of VOC2007 are
actually a subset of the 80 classes of COCO. So during training, we made use
of pretrained YOLOv4 model to help us calculate the object position loss.

###  4.2 Results Analysis

#### Step one

Results of the first step are shown in Figure 5, where the left image is real
and the right one is generated. As can be seen, our model can generate similar
image as real image. Although human being can recognize the object in
generated image, i.e., a horse, a rider and the prairie background, there are
still many noise-like pixels inside the generated image. These pixels may
largely affect the performance of detection model. This phenomenon might due
to the usage of WGAN, where the unsuitable magnitude of gradient penalty may
distract its weights. Then the generator may not be able to concentrate fully
on the competitive game itself, but also need to fulfill the requirement of
gradient penalty item.

Figure 5: Outcome of 1st-step model. Left: Real image from VOC dataset after
resizing into $64\times 64\times 3$. Right: Corresponding generated image from
our model.

#### Step two

The results of second step are shown in Figure 6. As shown in images, there is
a rectangle in each generated image. This is because in the input of
generator, we concatenate a mask to the noise in channel dimension. Without
telling what generator shall do with this mask, it just automatically
recognizes it as a shape it needs to display and thus there is an obvious
boundary around the rectangle in each image. It is an interesting phenomenon,
but it is actually not what we desired. In the future, one of our key
directions is finding a better way to help generator understand the purpose of
these input coordinates.

Figure 6: Temporary outcome of 2nd-step model.

###  4.3 Discussion

In this part, we will analyze our problems met during experiment.

Initially, we were planning to generate larger images with size $256\times
256\times 3$. However, because of the limit of both time and available GPUs,
after half month’s trial, we decided to decrease it to $64\times 64\times 3$.

It’s well-known that GAN training is unstable. To address this, we’ve tried
several methods to help stabilize the GAN training besides gradient penalty
after our failure at the very beginning: learning rate decay after certain
epochs and spectral normalization. Besides this, the largest problem is
definitely about finding the balance between generator and discriminator.

During training, we found that the discriminator outperforms generator most of
the time even if we use the architectures proved to have similar capacity
between generator and discriminator. The loss of such situation is shown in
Figure 7-left, where the generator cannot fool the discriminator after certain
epochs. To address this, we use multiple tricks from internet, including
adding Gaussian noise onto the input of discriminator, training generator with
more steps at each iteration or just directly increasing the number of
parameters of generator. However, it seems that this might relieve the problem
to some degree, but the balance is still too fragile to maintain in our
application. Most of the time, they compete with each other well at first, but
after maybe several hours, things suddenly go wrong and the model just falls
into the above situation that needs to stop and restart training.

On the contrary, if we train the generator more frequently than discriminator,
for example with a ratio of 3:1, then the loss curve will look like Figure
7-right. It seems like this figure shows a perfect competition process, but
actually, the generated images are even worse than the previous situation. The
reason is that in this situation, the ability of discriminator is too poor to
identify real images. A possible solution to these two situations is
increasing both the number of discriminator steps and generator steps per
iteration.

What’s more, compared to other application that use GAN for data augmentation,
it seems that their GAN training is not so difficult as ours. One possible
reason is that object detection dataset has more complex and distinct
background between different images, which might confuse the generator a lot.

Figure 7: Loss curve during our training. Red: score for real samples from
discriminator. Green and Blue: score for fake samples from discriminator.

###  4.4 Future Work

From our perspectives, this topic is interesting and creative. So we are going
to continue on this after this module if possible.

Our future work can be further extended in the following aspects:

  * •

Continue on finding a better way to tell generator the purpose of giving it
those coordinates.

  * •

Train a detection model based on augmented dataset to see if there are
performance improvement compared to the case when only using original dataset.

  * •

Further extend this framework, so that it can handle situation with multiple
input coordinates and generate multiple objects on the generated image at the
same time.

##  5 Conclusion

In this paper, we propose a novel framework to address one of the largest
shortages of object detection task: lack of labeled training data. With the
augmented dataset, we can further improve the performance of detection model
intuitively. However, because of the limit of time and available GPU, it’s a
pity that we cannot finish it during this module since we are all unfamiliar
with GAN before this project. And during this project, we entirely understand
why people always complain about the training of GAN. Anyway, it has been a
wonderful time playing with GAN structure. Thanks for your nice design of this
project that gives us an opportunity to work on this interesting area.

## References

  * [1] Antreas Antoniou, Amos Storkey, and Harrison Edwards.  Data augmentation generative adversarial networks, 2018.
  * [2] Martin Arjovsky, Soumith Chintala, and Léon Bottou.  Wasserstein gan, 2017.
  * [3] Yancheng Bai, Yongqiang Zhang, Mingli Ding, and Bernard Ghanem.  Sod-mtgan: Small object detection via multi-task generative adversarial network.  In Proceedings of the European Conference on Computer Vision (ECCV), September 2018.
  * [4] Alexey Bochkovskiy, Chien-Yao Wang, and Hong-Yuan Mark Liao.  YOLOv4: Optimal Speed and Accuracy of Object Detection.  arXiv e-prints, page arXiv:2004.10934, April 2020.
  * [5] Liang-Chieh Chen, George Papandreou, Iasonas Kokkinos, Kevin Murphy, and Alan L. Yuille.  Semantic Image Segmentation with Deep Convolutional Nets and Fully Connected CRFs.  arXiv e-prints, page arXiv:1412.7062, December 2014.
  * [6] Ross Girshick.  Fast R-CNN.  arXiv e-prints, page arXiv:1504.08083, April 2015.
  * [7] Ross Girshick, Jeff Donahue, Trevor Darrell, and Jitendra Malik.  Rich feature hierarchies for accurate object detection and semantic segmentation.  arXiv e-prints, page arXiv:1311.2524, November 2013.
  * [8] Ross Girshick, Jeff Donahue, Trevor Darrell, and Jitendra Malik.  Region-based convolutional networks for accurate object detection and segmentation.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(1):142–158, 2016.
  * [9] Kaiming He, Georgia Gkioxari, Piotr Dollár, and Ross Girshick.  Mask R-CNN.  arXiv e-prints, page arXiv:1703.06870, March 2017.
  * [10] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep Residual Learning for Image Recognition.  arXiv e-prints, page arXiv:1512.03385, December 2015.
  * [11] Sheng-Wei Huang, Che-Tsung Lin, Shu-Ping Chen, Yen-Yi Wu, Po-Hao Hsu, and Shang-Hong Lai.  Auggan: Cross domain adaptation with gan-based data augmentation.  In Proceedings of the European Conference on Computer Vision (ECCV), pages 718–731, 2018.
  * [12] Tero Karras, Samuli Laine, and Timo Aila.  A style-based generator architecture for generative adversarial networks, 2019.
  * [13] Kai Li, Yulun Zhang, Kunpeng Li, and Yun Fu.  Adversarial Feature Hallucination Networks for Few-Shot Learning.  arXiv e-prints, page arXiv:2003.13193, March 2020.
  * [14] Tsung-Yi Lin, Priya Goyal, Ross Girshick, Kaiming He, and Piotr Dollár.  Focal Loss for Dense Object Detection.  arXiv e-prints, page arXiv:1708.02002, August 2017.
  * [15] Wei-Hong Lin, Jia-Xing Zhong, Shan Liu, Thomas Li, and Ge Li.  RoIMix: Proposal-Fusion among Multiple Images for Underwater Object Detection.  arXiv e-prints, page arXiv:1911.03029, November 2019.
  * [16] Wei Liu, Dragomir Anguelov, Dumitru Erhan, Christian Szegedy, Scott Reed, Cheng-Yang Fu, and Alexander C. Berg.  SSD: Single Shot MultiBox Detector.  arXiv e-prints, page arXiv:1512.02325, December 2015.
  * [17] Shaoqing Ren, Kaiming He, Ross Girshick, and Jian Sun.  Faster r-cnn: Towards real-time object detection with region proposal networks.  In C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015.
  * [18] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation, 2015\.
  * [19] Mingxing Tan, Ruoming Pang, and Quoc V. Le.  EfficientDet: Scalable and Efficient Object Detection.  arXiv e-prints, page arXiv:1911.09070, November 2019.
  * [20] K. E. A. van de Sande, J. R. R. Uijlings, T. Gevers, and A. W. M. Smeulders.  Segmentation as selective search for object recognition.  In IEEE International Conference on Computer Vision, 2011.
  * [21] P. Viola and M. Jones.  Rapid object detection using a boosted cascade of simple features.  In Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. CVPR 2001, volume 1, pages I–I, 2001\.
  * [22] Chien-Yao Wang, Hong-Yuan Mark Liao, I-Hau Yeh, Yueh-Hua Wu, Ping-Yang Chen, and Jun-Wei Hsieh.  CSPNet: A New Backbone that can Enhance Learning Capability of CNN.  arXiv e-prints, page arXiv:1911.11929, November 2019.
  * [23] Xiaolong Wang, Abhinav Shrivastava, and Abhinav Gupta.  A-Fast-RCNN: Hard Positive Generation via Adversary for Object Detection.  arXiv e-prints, page arXiv:1704.03414, April 2017.
  * [24] Hongyi Zhang, Moustapha Cisse, Yann N. Dauphin, and David Lopez-Paz.  mixup: Beyond Empirical Risk Minimization.  arXiv e-prints, page arXiv:1710.09412, October 2017.
  * [25] Zhaohui Zheng, Ping Wang, Wei Liu, Jinze Li, Rongguang Ye, and Dongwei Ren.  Distance-iou loss: Faster and better learning for bounding box regression.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pages 12993–13000, 2020.


# Comparative Analysis of Molecular Dynamics and Method of Moments in Two-
Dimensional Concentric Circular Layers111This is the version of the article
before peer review or editing, as submitted by an author to Comparative
Analysis of Molecular Dynamics and Method of Moments in Two-Dimensional
Concentric Circular Layers IOP Publishing Ltd is not responsible for any
errors or omissions in this version of the manuscript or any version derived
from it. The Version of Record is available online at
https://doi.org/10.1088/1361-648X/ad5baf .

Robert Salazar ID<EMAIL_ADDRESS>Dirección de Ingeniería
Electrónica, Universidad ECCI, Bogotá, Colombia  Cristian Cobos
<EMAIL_ADDRESS>Dirección de Ingeniería Electrónica, Universidad
ECCI, Bogotá, Colombia  Diego Jaramillo ID<EMAIL_ADDRESS>Departamento de Ciencias, Universidad Antonio Nariño, Bogotá, Colombia  Camilo
Bayona-Roa ID<EMAIL_ADDRESS>Centro de Ingeniería Avanzada
Investigación y Desarrollo — CIAID, Bogotá, Colombia

###### Abstract

In this manuscript, we undertake an examination of a classical plasma deployed
on two finite co-planar surfaces: a circular region $\Omega_{in}$ into an
annular region $\Omega_{out}$ with a gap in between. It is studied both from
the point of view of statistical mechanics and the electrostatics of continua
media. We employ a dual perspective: the first one is by using Molecular
Dynamics (MD) simulations to find the system’s positional correlation
functions and velocity distributions. That by modeling the system as a
classical two-dimensional Coulomb plasma of point-like charged particles
$q_{1}$ and $q_{2}$ on the layers $\Omega_{in}$ and $\Omega_{out}$
respectively with no background density. The second one corresponds to a
finite surface electrode composed of planar metallic layers displayed on the
regions $\Omega_{in}$, $\Omega_{out}$ at constant voltages $V_{in}$, $V_{out}$
considering axial symmetry. The surface charge density is calculated by the
Method of Moments (MoM) under the electrostatic approximation. Point-like and
differential charges elements interact via a $1/r$ \- electric potential in
both cases. The thermodynamic averages of the number density, and electric
potential due to the plasma depend on the coupling and the charge ratio
$\xi=q_{1}/q_{2}$ once the geometry of the layers is fixed. On the other hand,
the fields due to the SE depend on the layer’s geometry and their voltage. In
the document, is defined a protocol to properly compare the systems. We show
that there are values of the coupling parameter, where the thermodynamic
averages computed via MD agree with the results of MoM for attractive $\xi=-1$
and repulsive layers $\xi=1$.

Keywords: molecular dynamics, method of moments, Coulomb systems, long-range
interaction.

##  1 Introduction

In this study, we examine a metallic configuration comprising distinct layers
denoted as $\Omega_{in}$ and $\Omega_{out}$, collectively forming a gapped
Surface Electrode (SE) devoid of background charge density. The inner layer is
a disk $\Omega_{in}=\left\\{(x,y):\sqrt{x^{2}+y^{2}}\leq
R_{1}\right\\}\subset\mathbb{R}^{2}$ of radius $R_{1}=R$ and the outer layer
is an annular region
$\Omega_{out}=\left\\{(x,y):R_{2}\leq\sqrt{x^{2}+y^{2}}\leq
R_{3}\right\\}\subset\mathbb{R}^{2}$. The system can be also modeled as two-
dimensional plasma where within $\Omega_{in}$, reside $N_{1}$ particles with
charge $q_{1}\in\mathbb{R}$, while $\Omega_{out}$ holds $N_{2}$ particles with
charge $q_{2}\in\mathbb{R}$ (see Fig. 1). The foregoing study draws
inspiration from [1], where the gapped surface electrode was characterized by
two-conductor flat regions with differing potentials and a distinct gap.
Electric field and surface charge density analysis in that system involved
resolving Laplace’s equation with specified conditions.

Fig. 1: The system. Two metallic sheets $\Omega_{in}$ and $\Omega_{out}$ with
no background density charge having $N_{1}$ particles of charge $q_{1}$ and
$N_{2}$ particles of charge $q_{2}$.

The $R_{1}\longrightarrow R_{2}$ (gapless) and $R_{3}\longrightarrow\infty$
(semi-infinite outer layer) limits of the gapped SE were studied in [2]. The
surface electrode was modeled as a two-dimensional classical Coulomb gas with
$+q$ and $-q$ charges on the inner and outer electrodes. Equilibrium states
and charge density were derived from Monte Carlo (MC) simulations. In that
study, interactions followed an inverse power law potential $1/r$, where the
coupling parameter $\Gamma$ has been found inversely proportional to the
system’s temperature and proportional to $q^{2}$.

Furthermore, authors of Ref. [1] explored the electric vector potential and
the Biot-Savart law analogy in electrostatics. These efforts set the zero
stage for the present study, aiming to further analyze the system in question
with the statistical mechanics approach. Subsequent sections detail our Monte
Carlo methodology, results, and implications within the previous
electrostatics descriptions of the system.

A strategy to study two-dimensional plasmas is by using Molecular Dynamics
(MD) simulations [3, 4, 5, 6, 7, 8, 9, 10, 11]. In the present study, we shall
perform simulations in the NVT-ensemble by using Nosé-Hoover and Langevin
thermostats. Our goal is to study the thermodynamic averages of number density
and the potential as a function of $\beta=k_{B}T$ for a fixed geometry of the
SE. We aim to know if a coulomb plasma with a finite number of particles can
be in agreement with the electrostatics of the SE.

The motivation behind this study stems from the electrostatic behavior of the
gapped surface electrode system comprising $\Omega_{in}$ and $\Omega_{out}$.
While prior research has shed light on this system’s electromagnetic fields,
our work extends the investigation into the realm of statistical mechanics. By
delving into the interplay between Molecular Dynamics simulations and the
Method of Moments technique, we aim to uncover the underlying principles
governing the system’s surface charge density and its correlation with the
Coulomb’s interaction between charges.

This document is structured as follows: We first delve into the Method of
Moments, focusing on the derivation of surface charge density for a single
circular layer and exploring the distinctive electrostatic behavior
originating from the interactions governed by a Coulomb potential. Then, we
describe the Molecular Dynamics techniques, together with the experimental
protocol that is followed to make comparisons. Subsequently, we present and
analyze results obtained through both the MoM and MD methods. The statistical
results are contrasted against analytic descriptions of SE. These findings are
discussed within the context of electrostatics and statistical mechanics,
highlighting their broader implications. Finally, we conclude by summarizing
our contributions.

##  2 Method of Moments

In this section, we delve into the main concepts of the Method of Moments and
its application to our system. From the point of view of the electrostatics of
continuum, (see Fig. 1-right) charge is continuously distributed over the SE
layers. The potential of the layers is known from the problem definition, but
the surface charge distribution $\sigma(\boldsymbol{r})$ is unknown. To
calculate the surface charge distribution on the layers we use the MoM. That
method is often employed to study classical electromagnetic radiation [12, 13,
14], but it is also useful to derive surface charge densities in electrostatic
systems [15, 16]. We start by examining the behavior of a single circular
layer and exploring the derivation of surface charge density. Furthermore, we
investigate the intricate interplay of interactions governed by the $1/r$
potential, shedding light on the unique electrostatic behavior of our system.

###  2.1 Single circular layer

First, we derive the surface charge density for a single circular layer
configuration. This first step serves as the foundational study in a single
component of the gapped SE system. The integral formula of the electrostatic
potential is given by

$\frac{1}{4\pi\epsilon_{o}}\int_{\Omega}\frac{\sigma(\boldsymbol{r}^{\prime})}{|\boldsymbol{r}-\boldsymbol{r}^{\prime}|}d^{2}\boldsymbol{r}=\Phi(\boldsymbol{r}\in\Omega)=V_{o},$

with $\epsilon_{o}$ the permittivity of vacuum and $V_{o}$ the voltage of the
circular layer $\Omega=\left\\{(x,y):\sqrt{x^{2}+y^{2}}\leq
R\right\\}\subset\mathbb{R}^{2}$ of radius $R$. Since the system has azimuthal
symmetry then potential electrostatic potential in the space is independent of
the angular $\phi$-coordinate $\Phi(\boldsymbol{r})=\Phi(u,z)$. Therefore, the
surface charge density is a function depending only of the radial
$u=\sqrt{x^{2}+y^{2}}$ coordinate $\sigma(\boldsymbol{r})=\sigma(u)$. Any
location on the plate is given by $\boldsymbol{r}(u,\phi)=u\hat{u}(\phi)$, and
the distance between two different points on the plate is

$|\boldsymbol{r}-\boldsymbol{r}^{\prime}|=\sqrt{u^{2}+(u^{\prime})^{2}-2uu^{\prime}\hat{u}(\phi)\cdot\hat{u}(\phi^{\prime})},$

with $\hat{u}(\phi)=(\cos\phi,\sin\phi,0)$ written in Cartesian coordinates.
Hence, the potential on the metallic layer can be calculated from

$\frac{1}{4\pi\epsilon_{o}}\int_{0}^{R}\int_{0}^{2\pi}\frac{\sigma(u)}{\sqrt{u^{2}+(u^{\prime})^{2}-2uu^{\prime}\cos(\phi-\phi^{\prime})}}u^{\prime}du^{\prime}d\phi^{\prime}=\Phi(\boldsymbol{r}\in\Omega)=V_{o}.$

We aim to obtain the surface charge density $\sigma$ from the previous
integral equation. The strategy is to expand $\sigma$ into a truncated series
of functions as follows,

$\sigma(u)=\sum_{n=0}^{N}\sigma_{n}f_{n}(u),$

with $\left\\{f_{n}(u)\right\\}_{n=1,\ldots,N}$ the set of basis functions
where $\sigma_{n}$ is the $n$-th coefficient for the surface charge density in
that basis. Then,

$\frac{1}{4\pi\epsilon_{o}}\sum_{n=0}^{N}\sigma_{n}\int_{0}^{R}\int_{-\phi}^{2\pi-\phi}\frac{f_{n}(u^{\prime})}{\sqrt{u^{2}+(u^{\prime})^{2}-2uu^{\prime}\cos(\beta)}}u^{\prime}du^{\prime}d\beta^{\prime}=\Phi\left((u,\phi,0)\in\Omega\right)=V_{o},$

where we have used the change of variable $\beta=\phi^{\prime}-\phi$. Note
that the choice of $\phi\in[0,2\pi]$ only ranges the surface of the plate,
then

$\frac{1}{4\pi\epsilon_{o}}\sum_{n=0}^{N}\sigma_{n}\int_{0}^{R}\left[\int_{0}^{2\pi}\frac{d\beta}{\sqrt{u^{2}+(u^{\prime})^{2}-2uu^{\prime}\cos(\beta)}}\right]f_{n}(u^{\prime})u^{\prime}du^{\prime}=\Phi\left((u,0,0)\in\Omega\right)=V_{o}.$

The angular integral can be evaluated straightforwardly since it is related to
the complete elliptic function of the first kind222This integral can be
simplified by introducing
$\zeta(u^{\prime},u):=\sqrt{(u^{\prime})^{2}+u^{2}-2u^{\prime}u}=|u-u^{\prime}|,\hskip
14.22636pt\mbox{and}\hskip
14.22636pt\xi:=\frac{4u^{\prime}u}{\zeta(u^{\prime},u)^{2}}=\frac{4u^{\prime}u}{|u-u^{\prime}|^{2}}.$
Changing the $\beta$ variable with $\beta=2t$ leads to
$\displaystyle(u^{\prime})^{2}+u^{2}-2u^{\prime}u\cos\beta$
$\displaystyle=(u^{\prime})^{2}+u^{2}-2u^{\prime}u+2u^{\prime}u(1-\cos\beta)$
$\displaystyle=((u^{\prime})^{2}+u^{2}-2u^{\prime}u)+4u^{\prime}u\sin^{2}t$
$\displaystyle=\zeta(u^{\prime},u)^{2}+\xi\zeta(u^{\prime},u)^{2}\sin^{2}t$
$\displaystyle=\zeta(u^{\prime},u)^{2}(1+\xi\sin^{2}t)=|u-u^{\prime}|^{2}(1+\xi\sin^{2}t).$
Hence, the integrand can be written as
$\int_{0}^{2\pi}\frac{d\beta}{\sqrt{u^{2}+u^{\prime
2}-2uu^{\prime}\cos(\beta)}}=\frac{4}{|u-u^{\prime}|}\int_{0}^{\pi/2}\frac{dt}{\sqrt{(1+\xi\sin^{2}t)}}=\frac{4}{|u-u^{\prime}|}K\left(-\xi\right).$
,

$\int_{0}^{2\pi}\frac{d\beta}{\sqrt{u^{2}+(u^{\prime})^{2}-2uu^{\prime}\cos(\beta)}}=\frac{4}{|u-u^{\prime}|}K\left(-\frac{4uu^{\prime}}{|u-u^{\prime}|^{2}}\right),$
(1)

where the elliptic function is given by

$K(\chi)=\int_{0}^{\pi/2}\frac{d\theta}{\sqrt{1-\chi\sin^{2}\theta}}.$

Fig. 2: Surface charge density of the singular circular plate at a constant
voltage according to MoM. (left) Profile of the charge density. (right)
Surface density charge for $\eta=1/2$ (transparent surface) and $\eta=3/2$
(solid surface).

The potential evaluated on the sheet is computed from

$\frac{1}{\pi\epsilon_{o}}\sum_{n=0}^{N}\sigma_{n}\int_{0}^{R}\frac{f_{n}(u^{\prime})}{|u-u^{\prime}|}K\left(-\frac{4uu^{\prime}}{|u-u^{\prime}|^{2}}\right)u^{\prime}du^{\prime}=\Phi\left((u,0,0)\in\Omega\right)=V_{o}.$

In this case, we can choose the value of $u$ in the interval $[0,R]$. We shall
define $\mathcal{N}+1$ equally spaced discrete values at the radial coordinate
$u_{0},\ldots,u_{\mathcal{N}}$ with $u_{n}=nR/\mathcal{N}$, to define
$\mathcal{N}$ intervals $U_{n}:=[u_{n},u_{n+1}]$ with centers located at
$u_{n}^{(c)}=(u_{n}+u_{n+1})/2$, then

$\frac{1}{\pi\epsilon_{o}}\sum_{n=0}^{N}\sigma_{n}\int_{0}^{R}\frac{f_{n}(u^{\prime})}{|u_{m}^{(c)}-u^{\prime}|}K\left(-\frac{4u_{m}^{(c)}u^{\prime}}{|u_{m}^{(c)}-u^{\prime}|^{2}}\right)u^{\prime}du^{\prime}=\Phi\left((u_{m}^{(c)},0,0)\in\Omega\right)=V_{o},$

which can be written as a linear set of algebraic equations that can be
inverted to find the $\sigma_{n}$ coefficients, as follows

$\sum_{n=0}^{\mathcal{N}}M_{mn}\sigma_{n}=\Phi_{m}\hskip
14.22636pt\therefore\hskip
14.22636pt\sigma_{m}=\sum_{n=0}^{\mathcal{N}}M_{mn}^{-1}\Phi_{n}=\sum_{n=0}^{\mathcal{N}}M_{mn}^{-1}\\{V_{o}\\}_{n},$

with

$M_{mn}=\frac{1}{\pi\epsilon_{o}}\int_{0}^{R}\frac{f_{n}(u^{\prime})}{|u_{m}^{(c)}-u^{\prime}|}K\left(-\frac{4u_{m}^{(c)}u^{\prime}}{|u_{m}^{(c)}-u^{\prime}|^{2}}\right)u^{\prime}du^{\prime}$

being a square matrix and
$\Phi_{m}=\Phi\left((u_{m}^{(c)},0,0)\in\Omega\right)$ a vector of constant
entries $V_{o}$.

There are several choice of basis functions
$\\{f_{u}(n)\\}_{n=0,\ldots,\mathcal{N}}$. One of the simplest choices is the
$C_{0}$ basis defined as

$f_{n}(u)=\begin{cases}1&u\in U_{n}=[u_{n},u_{n+1}],\\\
0&\mbox{otherwise},\end{cases}$

for which the $M_{mn}$ matrix becomes

$M_{mn}=\frac{1}{\pi\epsilon_{o}}\int_{U_{n}}\frac{u^{\prime}du^{\prime}}{|u_{m}^{(c)}-u^{\prime}|}K\left(-\frac{4u_{m}^{(c)}u^{\prime}}{|u_{m}^{(c)}-u^{\prime}|^{2}}\right).$

If the number of the interval divisions is large $\mathcal{N}>>1$ one can use
a mid-point approximation, leading to the simpler discrete term

$M_{mn}=\frac{1}{\pi\epsilon_{o}}\frac{u_{n}^{(c)}\Delta
u_{n}}{|u_{m}^{(c)}-u_{n}^{(c)}|}K\left(-\frac{4u_{m}^{(c)}u_{n}^{(c)}}{|u_{m}^{(c)}-u_{n}^{(c)}|^{2}}\right)\hskip
28.45274pt\mbox{for}\hskip 28.45274ptm\neq n,$

with $\Delta u_{n}=u_{n+1}-u_{n}$. However, it is often advisable to perform
numerical integration through sophisticated techniques to obtain accurate
values of $M_{mn}$. In this manuscript, these integrals are evaluated by using
the Global Adaptive method set by default in function NIntegrate of Wolfram
Mathematica 9.0 [17]. The MoM is a versatile method and it can be easily
generalized for other long-range potentials of the form $1/r^{\eta}$ with
$\eta\in\mathbb{R}^{+}$ (as demonstrated in Appendix A). For the case of a
single circular plate, the voltage $V_{o}$ plays a scaling role for the charge
density, thus one can set that voltage as one. In Fig. 2-left we show the
surface charge density profile for standard $1/r$-interaction potential (red-
solid line) for the circular plate. It occurs that surface charge density in
the electrostatic approximation diverges at the sheet border, that is
$\lim_{u\to R}\sigma(u)=\infty$. Numerically, we do not reach that limit since
we compute the vector $\sigma(u_{m}^{(c)})=\sigma_{m}$ and the largest value
$\sigma_{N}$ is located at $R-\Delta u_{\mathcal{N}}/2$ which approaches to
the border as $\mathcal{N}\longrightarrow\infty$, where
$\sigma_{\mathcal{N}\to\infty}\longrightarrow\infty$.

###  2.2 Two layers

The MoM for concentric circular layers (see Fig. 1) can be straightforwardly
generalized by following an analogous procedure to the one described for one
circular plate. A generalization of this procedure to the gapped SE is
presented in Appendix B. Summarizing, two basis functions
$\\{f_{n}(u)\\}_{n=1,\ldots,2\mathcal{N}}$ are employed to expand the surface
charge density: the first $n=1,\ldots,\mathcal{N}$ bases are employed for the
inner layer and the complementary $n=\mathcal{N}+1,\ldots,\mathcal{2N}$ are
used for the second layer. The surface charge density between these elements
is

$\sigma_{m}=V_{in}\sum_{n=1}^{\mathcal{N}}M_{m,n}^{-1}+V_{out}\sum_{n=\mathcal{N}+1}^{2\mathcal{N}}M_{m,n}^{-1}\hskip
28.45274ptm=1,\ldots,2\mathcal{N},$ (2)

with $M_{m,n}$ a $2\mathcal{N}\times 2\mathcal{N}$ matrix given by

$M_{mn}=\int_{u_{n-1}}^{u_{n}}I(u_{m}^{(c)},u^{\prime})u^{\prime}du^{\prime}\hskip
7.11317pt\mbox{{if}}\hskip 7.11317ptn=1,\ldots,\mathcal{N},\hskip
7.11317pt\mbox{{otherwise}}\hskip
7.11317pt\int_{u_{n}}^{u_{n+1}}I(u_{m}^{(c)},u^{\prime})u^{\prime}du^{\prime},$

and

$I(u,u^{\prime})=\frac{1}{\pi\epsilon_{o}}\frac{1}{|u-u^{\prime}|}K\left(-\frac{4uu^{\prime}}{|u-u^{\prime}|^{2}}\right).$

##  3 Molecular Dynamics

In this section, the SE system is modeled as $N_{1}$, $N_{2}$ distribution of
classical electric charges (with nominal values $q_{1}$ and $q_{2}$) confined
to planar geometries $\Omega_{1}$ and $\Omega_{2}$, respectively (see Fig. 1).
The coarse-grained method is used to model charged particles, which collide
elastically with the layer boundaries $\partial\Omega_{1}$,
$\partial\Omega_{2}$ and interacting with each other via a Coulomb-like
potential of the form

$\Phi(\boldsymbol{r}_{ij})=\frac{1}{4\pi\epsilon_{o}}\frac{q_{i}q_{j}}{|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|}\hskip
14.22636pt\mbox{if}\hskip
14.22636pt|\boldsymbol{r}_{i}-\boldsymbol{r}_{j}|>2b,\hskip
14.22636pt\mbox{otherwise}\hskip 14.22636pt\infty.$

We employ two different MD techniques. In the first one, we consider hard-disk
interactions with a radius of particles denoted by $b$, in addition to Coulomb
interactions. For that technique, we employ the Langevin thermostat. In the
second simulation, we eliminate hard-disk interactions and employ the Nose-
Hoover thermostat. In both cases, we implement a Verlet algorithm to
numerically integrate the equations of motion. Point-like charges remain on
the layers, ensuring that particles in $\Omega_{1}$ did not mix with particles
in $\Omega_{2}$.

For convenience, we introduce the parameter $\xi=q_{2}/q_{1}$, i.e., the ratio
between the two types of charges. We restrict that parameter in the range
$-1\leq\xi\leq 1$ fixing the charge of the inner layer $q_{1}=q$ and varying
the charge of the outer layer in the range $-q\leq q_{2}\leq q$.

We do not compensate the particles’ repulsion on the layers e.g. by
introducing a neutralizing background on each layer, since we look for charge
configurations that make layers equipotential surfaces, as it occurs in
electrostatics. In every simulation, we set the number of particles to
$N_{1}=N_{2}=N/2$ on the layers with $N\in 2\mathbb{N}$ the total number of
particles, thus being a globally neutral system when $\xi=-1$.

In MD simulations, we compute the number density $n(\boldsymbol{r})$ (i.e
number of particles per unit of area) that enables the computation of the
surface charge density, as follows,

$\sigma(u)=\begin{cases}n(u)q&0\leq u\leq R_{1},\\\ n(u)\xi q&R_{2}\leq u\leq
R_{3},\\\ 0&\mbox{otherwise}.\end{cases}$ (3)

The number density is computed by counting the particles in

$\boldsymbol{ann}(u,\Delta
u):=\left\\{(u^{\prime},\phi)\in\mathbb{R}^{2}:u<|(u^{\prime},\phi)|<u+\Delta
u\right\\}$

e.g. concentric annular regions of thickness $\Delta u$ and dividing by the
area of those regions.

##  4 Numerical experiment protocol

In order to perform proper comparisons between the results of the approaches
in the study, we follow an experimental design. This is necessary because the
starting parameters of MoM and MD are different, even for the common geometry
defined by $R_{1}$, $R_{2}$ and $R_{3}$. Regarding the MoM, one has to set the
layers’ potential $V_{in}$ and $V_{out}$. On the other hand, MD requires
establishing the number of particles of each layer $N_{1}$, $N_{2}$, the
charge ratio $\xi$, and the temperature $T$, since the system is in the
canonical ensemble. Of course, the starting parameters of MoM and MD are
related by the electrostatic energy

$U_{e}=\int_{\Omega_{in}\cup\Omega_{out}}\Phi(\boldsymbol{r})\sigma(\boldsymbol{r})d^{2}\boldsymbol{r}=V_{in}\int_{\Omega_{in}}\sigma(\boldsymbol{r})d^{2}\boldsymbol{r}+V_{out}\int_{\Omega_{out}}\sigma(\boldsymbol{r})d^{2}\boldsymbol{r}=V_{in}Q_{in}+V_{out}Q_{out},$

with $Q_{in}=N_{1}q_{1}$ and $Q_{out}=N_{2}q_{2}$ the total charge of each
layer. This is the energy stored, or conversely, the the amount of energy
spent to charge both layers. That energy is always positive since

$U_{e}=\frac{1}{2}\epsilon_{o}\int_{\mathbb{R}^{3}}\boldsymbol{E}^{2}(\boldsymbol{r})d^{3}\boldsymbol{r}\geq
0,$

with $\boldsymbol{E}(\boldsymbol{r})$ the electric field on the
$\mathbb{R}^{3}$-space due to the surface charge distribution. It implies that

$V_{in}N_{1}q_{1}+V_{out}N_{2}q_{2}\geq 0,$

thus the potential employed in the MoM and the total charge on the layer set
in MD are conditioned by the electrostatic energy. In order to avoid the use
of scaling factors in the computation of surface charge density comparisons we
proceed as follows:

  * •

Step 1. We perform a MD simulation of the system by setting the number of
particles $N=2N_{1}=2N_{2}$, the charge ratio $\xi$, and the thermostat target
temperature.

  * •

Step 2. Once the system is in equilibrium, the particles’ density $\langle
n(u)\rangle$ is computed from the thermodynamic average of configurations. The
surface charge density is calculated using Eq. (3).

  * •

Step 3. The average potential of those several configurations is computed from

$\langle\Phi(u,z)\rangle=\frac{1}{4\pi\epsilon_{o}}\sum_{n}\langle\sigma(u_{n})\rangle
u_{n}\delta
u_{n}\frac{4}{\sqrt{(u-u_{n})^{2}+z^{2}}}K\left(-\frac{4uu_{n}}{(u-u_{n})^{2}+z^{2}}\right),$
(4)

with $u_{n}$’s radial positions on the layers. Brackets $\langle\cdot\rangle$
indicate thermodynamic average on equilibrated configurations.

  * •

Step 4. The potential of the layers is estimated as follows

$V_{in}=\frac{1}{R_{1}}\int_{0}^{R_{1}}\langle\Phi(u,\delta)\rangle du\hskip
7.11317pt\mbox{and}\hskip
7.11317ptV_{out}=\frac{1}{R_{3}-R_{2}}\int_{R_{2}}^{R_{3}}\langle\Phi(u,\delta)\rangle
du,$

with small value of $\delta\approx 0$. These integrals are numerically
approximated.

  * •

Step 5. The MoM is fed with the potential values $V_{in}$ and $V_{out}$ from
the previous step.

##  5 Results

In this section, we present the comparisons between MD and MoM simulations. We
perform MD simulations in the NVT-ensemble by using Nosé-Hoover and Langevin
Thermostats [18, 19], by placing $N_{1}=500$ particles in the inner layer and
$N_{2}=500$ in the outer one, varying the charge ratio $\xi$ and the parameter
$\beta=k_{B}T$. Typical runs of MD simulations consist of $5\times 10^{6}$
steps with a time step of $\Delta t=0.001\tau$. The number of MD steps to
reach equilibrium are in the order of $1.5\times 10^{6}$.

For each set of parameters e.g. $\xi$, $\beta$ and fixed geometry, we feed the
MoM with the thermodynamic average scalar potential near the layers calculated
via MD. The MoM is run with 50 basis functions for each layer.

In Subsection 5.1, the coupling effect on the system behavior is discussed. In
general, a finite N-system which is far from the thermodynamic limit cannot be
described by the classical and continuum electrostatics in the
$\beta\to\infty$ (weak coupling) and $\beta\to 0$ (strong coupling) limits.
However, it is possible to find finite values of $\beta$ where the system
behaves according to continuum electrostatics. Subsection 5.2 is devoted to
checking the MD simulations and the MoM calculations of surface density with
the exact electrostatic result in the SE gapless limit. The last subsection
compares MD and MoM results in the SE gapped case.

###  5.1 The coupling effect over the finite SE system

In general, the MD thermodynamic averages and the electrostatic predictions
may diverge far from the thermodynamic limit $N\to\infty$, depending on the
coupling between particles. For finite-size systems, the coupling between
particles can change drastically the particle density distribution on the
sheets, and the behavior of the system cannot be simply described by standard
electrostatics of continuous media. To exemplify this situation, the behavior
of the finite system in a strong coupling regime is considered. Fig. 3 shows
the results of a MD simulation with a Nosé-Hoover thermostat at
$\beta=k_{B}T=0.03$. The parameters of the system are set as
$N_{1}=N_{2}=500$, $R_{1}=4$, $R_{2}=5R_{1}/4$, and $R_{3}=2R_{1}$, with
opposite change $\xi=-1$. At this value of $\beta$, the electric interaction
between particles is large enough in comparison with their kinetic energy,
such that it freezes the system.

Technically, a non-crystalline solid is obtained at the $\beta\rightarrow 0$
limit via MD. At $\beta=0.03$ the speed distribution $f(v)$ found through MD
does not match the speed distribution of the classical two-dimensional ideal
system:

$f(v)=\frac{mv}{\beta}\exp(-\frac{1}{2}\beta mv^{2})\hskip
14.22636pt\mbox{Maxwell speed distribution of a 2D ideal gas}$ (5)

at the same $\beta$, since the histogram shows that there are fixed particles
(Fig. 3-b) due to the coupling.

(a) Integrated charge via MD and MoM

(b) Speed distribution

(c) Potential at $z=z_{1}=R_{1}/100$.

(d) Averaged Potential

Fig. 3: Strong coupling behaviour for $\xi=-1$ and $\beta=0.03$ MD results.
The MD results corresponds to the Nosé-Hoover thermostat. The values of $z$
for (d) are $z_{1}=R_{1}/100$, $z_{2}=R_{1}/50$, $z_{3}=R_{1}/2$,
$z_{4}=R_{1}$,$z_{5}=2R_{1}$, and $z_{6}=10R_{1}$.

At the strong coupling regime $\beta\rightarrow 0$, the particles at the outer
layer are arranged in concentric rings, as shown in the inset plot of Fig.
3-d. On the other hand, the particles in the inner layer do not generate ring
arrangements at $\beta=0.03$.

Fig. 3-c shows the normalized integrated charge $q/Q$ on each layer given by

$q(u)=2\pi\begin{cases}\int_{0}^{u}\sigma(u^{\prime})du^{\prime}&\mbox{inner
sheet,}\\\ \int_{R_{2}}^{u}\sigma(u^{\prime})du^{\prime}&\mbox{outer
sheet,}\end{cases}$ (6)

with $Q=q(R_{1})$ and $Q=q(R_{3})$ the total charge of each layer.

One can observe that the integrated charge via MD and MoM are in agreement in
the inner layer, but they are not for the outer one as it is shown in Fig.
3-a. The plot of $q/Q$ at the outer layer via MD has a stairs-like profile.
Those jumps on $q(u)$ for $u\in[R_{2},R_{3}]$ are due the circular
arrangements of the particles at the outer layer. These circular arrangements
of particles also affect the thermodynamic average potential near the outer
layer, as it is shown in Fig. 3-c.

Fig. 4: MD average potential near the outer layer at the plane
$z=z_{1}=R_{1}/100$. Molecular dynamics simulations with $\beta=0.03,0.06,$
and $0.09$ (left to right).

The thermodynamic average potential $\langle\Phi(u,z)\rangle$ computed via MD
is shown in Fig. 3-d. At $z=0.004$ and $z=0.04$, the average potential shows
peaks corresponding to the discrete nature of the system in the strong
coupling regime. Those peaks practically disappear at $z=0.08$ where the curve
of $\langle\Phi(u,z)\rangle$ tends to a flat region $[0,R_{1}]$ corresponding
to the inner layer. This is still consistent with electrostatics since the
inner layer is metallic and it corresponds to an equipotential surface. The
average potential evaluated near the outer layer $[R_{2},R_{3}]$ also exhibits
peaks for $z=0.004$ and $z=0.04$. The number of peaks on the potential profile
Fig. 3-c & d corresponds to the number of rings of the typical configurations
near the ground state (one of them plotted in the inset image of Fig. 3-d) and
the stair-like profile of $q/Q$ at the outer layer. The behavior of the system
at the strong coupling regime is therefore not described by the standard
electrostatics due to the discrete effects. As the $\beta$ parameter is
increased, the circular arrangements of particles on the outer (typical of the
ground state of the finite system) disappear due to the decoupling among
particles. The progressing growth of the $\beta$ parameter also affects the
thermodynamic average of the potential evaluated near the outer layer, as it
is shown in Fig. 4. In that figure, one can observe that $<\phi(z,u)>$ tends
to be flat as it is evaluated at $z=z_{1}=R_{1}/100$ near the outer metallic
layer. This is according to electrostatics since metallic layers are
equipotential surfaces. We have to remark that $<\phi(z,u)>$ starts to be an
equipotential function near both metallic layers if particles move enough to
demonstrate a velocity distribution (consistent with the Maxwell-Boltzmann
distribution) such as the one in Eq. 5 (the one of the ideal gas), but
sufficiently electrically coupled to maintain the electrostatic features.

Fig. 5: MD velocity distribution for different couplings.

In Fig. 5, the MD velocity distribution is shown for $\beta=0.25$ and
$\beta=2.0$, and charge ratios of $\xi=1$ (repulsive layers) and $\xi=-1$
(attractive layers). In both cases the MD speed distribution $f(v)$ computed
with the Nosé-Hoover thermostat is according to the Maxwell speed distribution
given by Eq. (5). However, only the value $\beta=0.25$ ensures that moving
particles are sufficiently coupled to generate collectively a flat enough
potential $<\phi(z,u)>$ on the layers, as it is shown in Fig. 6 for attractive
and repulsive layers. On the other hand, if the parameter is set as
$\beta=2.0$, then the thermodynamic average of the potential does not generate
equipotentials near the layers (see Fig. 6).

Fig. 6: MD potential profiles as a function of $\beta$. (left) Attractive
layers $\xi=-1$ and (right) repulsive layers $\xi=1$.

The $\beta\rightarrow\infty$ limit case corresponds to a plasma, where kinetic
energy greatly exceeds the electric interaction and the surface charge density
on the layers must be uniform. A close study of the $\beta\rightarrow\infty$
limit for the single circular layer starts from the surface charge density,

$<\sigma(\boldsymbol{r})>=<n(\boldsymbol{r})>q=\frac{1}{Z_{N,\beta}}\frac{q}{(N-1)!}\prod_{n=2}^{N}\int_{0}^{2\pi}d\phi_{n}\int_{0}^{R}u_{n}du_{n}\prod_{1\leq
i<j\leq
N}\exp\left(-\frac{1}{4\pi\epsilon_{o}}\frac{q_{i}q_{j}/\beta}{\sqrt{u_{i}^{2}+u_{j}^{2}-2u_{i}u_{j}\cos(\phi_{ij})}}\right),$

with $q$ the charge of the particles (assumed to be identical), $R$ the radius
of the layer, $\phi_{ij}=\phi_{i}-\phi_{j}$, and

$Z_{N,\beta}=\frac{1}{N!}\prod_{n=1}^{N}\int_{0}^{2\pi}d\phi_{n}\int_{0}^{R}u_{n}du_{n}\prod_{1\leq
i<j\leq
N}\exp\left(-\frac{1}{4\pi\epsilon_{o}}\frac{q_{i}q_{j}/\beta}{\sqrt{u_{i}^{2}+u_{j}^{2}-2u_{i}u_{j}\cos(\phi_{ij})}}\right)$

the canonical partition function. In the $\beta\rightarrow\infty$ limit, the
partition function simplifies to a two-dimensional ideal gas and

$\lim_{\beta\to\infty}<\sigma(\boldsymbol{r})>=\frac{Nq}{\pi R^{2}}$

is distributed uniformly on the layer. The same occurs for the two-layer
system. Then, the potential can be developed from

$\displaystyle\lim_{\beta\to\infty}\langle\Phi(\boldsymbol{r})\rangle$
$\displaystyle=\frac{1}{4\pi\epsilon_{o}}\int_{\Omega_{in}\cup\Omega_{out}}\lim_{\beta\to\infty}\frac{\langle\sigma(\boldsymbol{r}^{\prime})\rangle}{|\boldsymbol{r}-\boldsymbol{r}^{\prime}|}d^{2}\boldsymbol{r}$
$\displaystyle=\frac{1}{4\pi\epsilon_{o}}\int_{u^{\prime}\in[0,R_{1}]\cup[R_{2},R_{3}]}\lim_{\beta\to\infty}\langle\sigma(u^{\prime})\rangle
u^{\prime}du^{\prime}\int_{0}^{2\pi}\frac{d\phi^{\prime}}{\sqrt{u^{2}+u^{\prime
2}-2uu^{\prime}\cos(\phi-\phi^{\prime})+z^{2}}}$
$\displaystyle=\frac{1}{4\pi\epsilon_{o}}\int_{u^{\prime}\in[0,R_{1}]\cup[R_{2},R_{3}]}\lim_{\beta\to\infty}\langle\sigma(u^{\prime})\rangle
u^{\prime}du^{\prime}\frac{4}{\sqrt{(u-u^{\prime})^{2}+z^{2}}}K\left(-\frac{4uu^{\prime}}{(u-u^{\prime})^{2}+z^{2}}\right).$

Since the surface charge density at each layer is uniform at the
$\beta\to\infty$ limit, then

$\Phi^{(unif)}(u,z)=\lim_{\beta\to\infty}\langle\Phi(\boldsymbol{r})\rangle=\frac{1}{\pi\epsilon_{o}}\sum_{k=1,2}\frac{N_{k}q_{k}}{A_{k}}\int_{U_{k}}\frac{u^{\prime}du^{\prime}}{\sqrt{(u-u^{\prime})^{2}+z^{2}}}K\left(-\frac{4uu^{\prime}}{(u-u^{\prime})^{2}+z^{2}}\right),$
(7)

with $U_{1}=[0,R_{1}]$, $U_{1}=[R_{2},R_{3}]$, $A_{1}=\pi R_{1}^{2}$,
$A_{2}=\pi(R_{3}^{2}-R_{2}^{2})$ the area of the layers. A plot of Eq (7) is
shown in Fig. 6. We can observe that MD simulations for $\beta=2.0$ (far from
the thermodynamic limit) and Eq. (7) for $\beta\to\infty$ do not generate
equipotential surfaces since the surface charge distribution tends to be
uniform, as occurs with the ideal gas.

###  5.2 The gapless limit

The exact density charge in the gapless limit is given by

$\sigma^{Gapless}(u)=\lim_{R_{3}\to\infty}\lim_{R_{2}\to
R_{1}}\sigma(u)=\frac{2\epsilon_{o}(V_{in}-V_{out})}{\pi}\left\\{\frac{1}{R_{1}-u}E\left[\frac{4R_{1}u}{(R_{1}+u)^{2}}\right]+\frac{1}{R_{1}+u}K\left[\frac{4R_{1}u}{(R_{1}+u)^{2}}\right]\right\\},$
(8)

with $K(z)$ and $E(z)$ the complete elliptic integrals of the first and second
kind respectively. Eq. (8) can be obtained by using an electric vector
potential approach described in [20, 21]. This limit is achieved in MD when we
set $N_{1}=N_{2}$ and $\xi=-1$ e.g. the same number of particles in both
layers with opposite signs.

Fig. 7: Surface density charge in the gapless limit. The inner sheet is at
$V_{in}=6.2043V$ and the outer one is grounded $V_{out}=-0.1028V$. The exact
solution corresponds to Eq. (8). MD simulations are performed with $N=1000$
particles at $\beta=0.25$. The inset figure corresponds to a equilibrated
configuration of the particles in both layers via MD simulations.

One can use Eq. (8) as a benchmark for the current problem in the gapless
limit. Fig. 7 compares the density profiles due to MoM and the exact solution
for $R_{1}=0.999R_{2}$ and $R_{3}=3R_{2}$. Note that, formally, the Eq. (8)
works in the $R_{3}\rightarrow\infty$ limit. Nonetheless, the MoM gives an
accurate approximation for a finite value of $R_{3}$, since particles of
different charges are strongly coupled near the gap as $R_{1}\rightarrow
R_{2}$ and practically no particle reaches the boundary at $R_{3}$. Similarly,
the MD results for the density profile calculated with the Langevin and Nose-
Hoover thermostats are also in agreement with the exact solution, as shown in
Fig. 7. Note that Eq. (8) strictly works for $\xi=-1$, this is, opposite
charges in each layer where the system is globally neutral. In this scenario,
charges on the metallic layers are highly concentrated near the
$R=(R_{1}+R_{2})/2$ circle due to electric attraction (see the inset image in
Fig. 7). Other scenarios concerning equal charge as $\xi=1$ cannot be modeled
in the gapless limit since particles on the outer layer will escape to the
infinite.

###  5.3 The gapped case

This subsection presents the results for the finite SE system with a gap. The
parameters which define the geometry were set as follows: $R_{1}=4$,
$R_{2}=5R_{1}/4$ and $R_{3}=2R_{1}$.

Fig. 8: Surface density charge for $\xi=-1$ and $\beta=0.25$. Symbols
correspond to MD simulations and the solid line represents the MoM solution.

The plot in Fig. 8 shows the density profile on both sheets for $\xi=-1$, this
is, inner particles having positive charge $q_{1}=q$ and outer ones carrying
opposite charge $q_{1}=-q$ thus the system being globally neutral. The solid
line in Fig. 8 corresponds to the MoM approximation calculated according to
Eq. (2). MoM solution assumes that the density charge is continuous and it
differs from the one in Eq. (8) (valid in the gapless limit) since MoM
predicts a charge concentration in the $R_{3}$ boundary. One observes that
numerical solutions provided by MD simulations with Langevin and Nosé-Hoover
thermostats are in agreement with the MoM solution. The parameter $\beta$ in
the MD simulation is set as 0.3 for both thermostats. The inset semi log-plot
in Fig. 8 is the integrated charge (see Eq. (9)), this is, the total charge
stored in a circular region of radius $u$ in each layer scaled by the total
charge $Q$ of the layer.

Fig. 9: Scalar electric potential for $\xi=-1$ and $\beta=0.25$. The inset
plot is the potential in the $\mathbb{R}^{3}$-space close to the $xy$-plane at
$z=z_{1}$. Metallic layers are represented as black surfaces. The symbols are
the MD simulations with the Langevin thermostat and solid lines correspond to
MoM solution. The set of heights chosen for this plot are $z_{1}=R_{1}/100$,
$z_{2}=R_{1}/50$, $z_{3}=R_{1}/2$, $z_{4}=R_{1}$,$z_{5}=2R_{1}$, and
$z_{6}=10R_{1}$.

Fig. 9 shows the scalar electric potential profiles due to the charged layers
at different heights $z_{1}<z_{2}<\ldots<z_{6}$. Methodologically we start by
running a MD simulation with a fixed geometry, number of particles, and
$\beta$ parameter. Later the surface charge density is calculated by
thermodynamic averages of the number density of each layer. We build the
scalar potential profile very close to the $xy$-plane by using Eq. (4), which
enables us to compute the inner $V_{in}$ and outer $V_{in}$ voltages. The
voltages computed via MD simulations are used as input parameters in Eq. (2)
to compute surface charge density with the MoM. That procedure not only
ensures that both numerical methods generate profiles on the same scale but
also the same profile distributions since a random choice of voltages in the
MoM may develop profiles corresponding to several ratios of charges.

Fig. 10: Scalar electric potential for $\xi=1$ and $\beta=0.25$. The inset
plot is the potential in the $\mathbb{R}^{3}$-space close to the $xy$-plane at
$z=z_{1}$. Metallic layers are represented as black surfaces. The symbols are
the MD simulations with the Langevin thermostat and solid lines correspond to
MoM solution.

Fig. 10 shows the scalar electric potential for the $\xi=1$ case, this is, the
same charge on both layers. The MD simulations are performed with a thousand
of identically charged particles, $N_{1}=500$ charges in the inner layer and
$N_{2}=500$ in the outer one. In general, the results of the DM potential and
the MoM are also consistent for the case of repulsive layers. However, unlike
the attractive case, there is a slight discrepancy in the potential over the
outer plate near the edge $R_{3}$ as can be observed in Fig. 10, where
$z=z_{1}$. This occurs because the particles on the inner plate, having the
same charge as the particles on the outer plate, repel the latter in such a
way that these charges less frequently visit the boundary $R=R_{3}$. This
results in the statistics near this boundary not being as good and the
potential appearing less constant in this region of the conductor.

The electric field due to configurations for $\xi=-1$ (opposite charge) and
$\xi=1$ (equal charge) are shown in Fig. 11. The electric field is obtained
from computing the numerical gradient of Eq. (4). That is,

$E_{u}(u,z)=-\frac{\partial}{\partial u}\langle\Phi(u,z)\rangle=-\lim_{\Delta
u\to 0}\frac{\langle\Phi(u+\Delta u,z)\rangle-\langle\Phi(u,z)\rangle}{\Delta
u}$

and

$E_{z}(u,z)=-\frac{\partial}{\partial z}\langle\Phi(u,z)\rangle=-\lim_{\Delta
z\to 0}\frac{\langle\Phi(u,z+\Delta z)\rangle-\langle\Phi(u,z)\rangle}{\Delta
z}$

by the forward finite differences of the scalar potential thermodynamic
average. MD thermodynamic averages of the electric field agree with the
electrostatic prediction of the streamlines of $\boldsymbol{E}$. Thus we can
observe agreement between the discrete MD and continuous MoM approaches at
$\beta=0.25$ since thermodynamic averages of the electric potential near the
layer tend to a constant value (see Figs. 9 and Figs. 10), and consequently,
the electric field streamlines become orthogonal to the layers (see Fig. 11).

Fig. 11: Electric field from the thermodynamic average $<\sigma>$ at
$\beta=0.25$ computed via MD simulations. The field is evaluated on the plane
$y=0$. (left) Layers having opposite charges and (right) the same charge.
Metallic layers are represented as solid lines on the $x$-axis. The blue and
red color corresponds to the positive and negative charge respectively.

##  6 Conclusions

In this manuscript, we conducted molecular dynamics simulations of a classical
plasma in two dimensions, incorporating Coulomb interactions within a circular
geometry devoid of background density. The system’s electrostatic analog is a
circular surface electrode in the context of finite metallic layers held at
fixed potentials. Our methodology emphasized the application of molecular
dynamics simulations and the Method of Moments to solve the system, alongside
developing a protocol for their comparative analysis. This rigorous approach
revealed the MD method’s efficacy in correctly describing the system under
specific coupling regimes, highlighting the nuanced behavior of the two-
dimensional plasma as opposed to the simpler SE system.

The complexity of the plasma system extends beyond the scope of continuum
media electrostatics, especially when considering scenarios far from the
thermodynamic limit. MD simulations near the ground state showed finite size
effects leading to circular arrangements of particles, which significantly
influence the thermodynamic average of potential near the metallic layers.
Notably, a comparison of the integrated charge on the layers between MD and
MoM at the strong coupling regime ($\beta\rightarrow 0$) showed a staircase-
like behavior, absent in electrostatics, for which the integrated charge
profile is smooth. Moreover, at the limit $\beta\rightarrow\infty$, the system
generates uniform surface charge distributions, resulting in potential
profiles on the layers that defy the equipotential surfaces inherent in
electrostatics of metallic layers.

Our numerical analysis determined that achieving the equipotential
characteristic of metallic layers from the Coulomb plasma requires increasing
the value of $\beta$ from the ground state. This discovery facilitated the
identification of a $\beta$ range where the thermodynamic average of surface
charge density via MD closely resembles that predicted by the MoM for the SE
system. Within this range, particularly at $\beta=0.25$, finite-size effects
of the Coulomb plasma, even with a thousand particles, become negligible in
replicating the SE’s electrostatics. This condition was further corroborated
in the gapless limit, where MD simulations, employing Nosé-Hoover and Langevin
thermostats, the MoM, and the analytical benchmark provided by Eq. (8),
demonstrated remarkable congruence.

While the gaped case lacks an exact solution, implementing the MoM allowed for
the approximation of surface charge density on the metallic layers. The MoM’s
outcomes, particularly at $\beta=0.25$ and considering attractive layers
($\xi=-1$), aligned well with those derived from MD simulations. This
comprehensive study not only underscores the critical role of methodological
adherence in analyzing complex systems but also illuminates the intricate
behavior of two-dimensional plasmas, thereby contributing valuable insights
into the field of electrostatics and plasma physics.

## Acknowledgments

This work was partially supported by Vicerrectoría de investigación,
Universidad ECCI through project IN-08-28. Camilo Bayona expresses gratitude
to the Centro de Ingeniería Avanzada, Investigación y Desarrollo - CIAID for
the technical support provided for the present study. Cristian Cobos thanks
the financial support of Centro de investigación e innovación en tecnología y
ciencia CEINTECCI project CEINTECCI-CPJI-2023-ID-14. Diego Jaramillo thanks
Universidad Antonio Nariño for financial support through project No 2023202.

## References

  * [1] R. Salazar, C. Bayona-Roa, and D. Jaramillo, “Electric vector potential and the biot-savart like law in electrostatics,” Revista de la Academia Colombiana de Ciencias Exactas, Físicas y Naturales, vol. 46, no. 181, pp. 939–946, 2022.  doi: 10.18257/raccefyn.1671
  * [2] R. Salazar, C. Bayona, and G. Tellez, “Monte carlo simulations of two-component coulomb gases applied in surface electrodes,” Journal of Physics: Condensed Matter, 2022.  doi: 10.1088/1361-648X/ac4aa8
  * [3] G. Dimonte and J. Daligault, “Molecular-dynamics simulations of electron-ion temperature relaxation in a classical coulomb plasma,” Physical review letters, vol. 101, no. 13, p. 135001, 2008.
  * [4] B. Zelener, B. Zelener, E. Manykin, S. Y. Bronin, A. Bobrov, and D. Khikhlukha, “Self-diffusion and conductivity in an ultracold strongly coupled plasma: Calculation by the method of molecular dynamics,” in Journal of Physics: Conference Series, vol. 946, p. 012126, IOP Publishing, 2018.
  * [5] Y. V. Arkhipov, A. Askaruly, A. Davletov, D. Y. Dubovtsev, Z. Donkó, P. Hartmann, I. Korolov, L. Conde, and I. Tkachenko, “Direct determination of dynamic properties of coulomb and yukawa classical one-component plasmas,” Physical Review Letters, vol. 119, no. 4, p. 045001, 2017.
  * [6] C. Durniak, D. Samsonov, N. P. Oxtoby, J. F. Ralph, and S. Zhdanov, “Molecular-dynamics simulations of dynamic phenomena in complex plasmas,” IEEE transactions on plasma science, vol. 38, no. 9, pp. 2412–2417, 2010\.
  * [7] G. Zwicknagel, C. Toepffer, and P.-G. Reinhard, “Molecular dynamic simulations of ions in electron plasmas at strong coupling,” Hyperfine Interactions, vol. 99, no. 1, pp. 285–291, 1996.
  * [8] A. Calisti and B. Talin, “Classical molecular dynamics model for coupled two-component plasmas–ionization balance and time considerations,” Contributions to Plasma Physics, vol. 51, no. 6, pp. 524–528, 2011.
  * [9] B. Scheiner and S. D. Baalrud, “Testing thermal conductivity models with equilibrium molecular dynamics simulations of the one-component plasma,” Physical Review E, vol. 100, no. 4, p. 043206, 2019.
  * [10] J. Mithen, J. Daligault, and G. Gregori, “Molecular dynamics simulations for the shear viscosity of the one-component plasma,” Contributions to Plasma Physics, vol. 52, no. 1, pp. 58–61, 2012.
  * [11] A. Calisti, S. Ferri, and B. Talin, “Classical molecular dynamics model for coupled two component plasmas,” High Energy Density Physics, vol. 5, no. 4, pp. 307–311, 2009.
  * [12] W. C. Gibson, The method of moments in electromagnetics.  Chapman and Hall/CRC, 2021.  doi: 10.1201/9780429355509
  * [13] L.-W. Li, M.-S. Yeo, and M.-S. Leong, “Method of moments analysis of em fields in a multilayered spheroid radiated by a thin circular loop antenna,” IEEE transactions on antennas and propagation, vol. 52, no. 9, pp. 2391–2402, 2004.  doi: 10.1109/TAP.2004.834018
  * [14] L.-W. Li, C.-P. Lim, and M.-S. Leong, “Method-of-moments analysis of electrically large circular-loop antennas: Nonuniform currents,” IEE Proceedings-Microwaves, Antennas and Propagation, vol. 146, no. 6, pp. 416–420, 1999.  doi: 10.1049/ip-map:19990784
  * [15] R. F. Harrington, “The method of moments in electromagnetics,” Journal of Electromagnetic waves and Applications, vol. 1, no. 3, pp. 181–200, 1987\.  doi: 10.1163/156939387X00018
  * [16] C. A. Balanis, Antenna theory: analysis and design.  John wiley & sons, 2016.
  * [17] M. Wolfram, “Version 9.0,” Champaign, IL, 2012.
  * [18] G. Kumar, R. R. Mishra, and A. Verma, “Introduction to molecular dynamics simulations,” in Forcefields for Atomistic-Scale Simulations: Materials and Applications, pp. 1–19, Springer, 2022.  doi: 10.1007/978-981-19-3092-8
  * [19] B. Leimkuhler and C. Matthews, “Molecular dynamics,” Interdisciplinary applied mathematics, vol. 39, p. 443, 2015.  doi: 10.1007/978-3-319-16375-8
  * [20] R. Salazar, C. Bayona, and G. Téllez, “Electric vector potential formulation in electrostatics: Analytical treatment of the gaped surface electrode,” Eur. Phys. J. Plus, vol. 135, no. 878, 2020.  doi: https://doi.org/10.1140/epjp/s13360-020-00864-0
  * [21] R. Salazar, C. Bayona, and J. Chaves, “Electrostatic field of angular-dependent surface electrodes,” Eur. Phys. J. Plus, vol. 135, no. 93, 2019.  doi: 10.1140/epjp/s13360-019-00090-3

##  Appendix A The $1/r^{\eta}$ interaction

Let us consider the electrostatic problem with the following particle-particle
interaction potential $\Phi_{\eta}(r)=q\nu_{\eta}(r)$, with

$\nu_{\eta}(\boldsymbol{r},\boldsymbol{r^{\prime}})=\nu_{\eta}(|\boldsymbol{r}-\boldsymbol{r^{\prime}}|)=\frac{L^{\eta-1}}{4\pi\epsilon_{o}}\frac{1}{|\boldsymbol{r}-\boldsymbol{r}^{\prime}|^{\eta}}$

a scalar function, $L$ a characteristic parameter with units of length, and
$\eta$ a real positive number. The surface charge density is found from the
following integral expression

$\int_{\Omega}\sigma(\boldsymbol{r}^{\prime})\nu(\boldsymbol{r},\boldsymbol{r^{\prime}})d^{2}\boldsymbol{r}=\Phi(\boldsymbol{r}\in\Omega)=V_{o},$

which can be written as

$\sum_{n=0}^{N}\sigma_{n}\int_{0}^{R}\mathcal{I}_{\eta}(u_{m}^{(c)},u^{\prime})f_{n}(u^{\prime})u^{\prime}du^{\prime}=\Phi\left((u_{m}^{(c)},0,0)\in\Omega\right)=V_{o}.$

by introducing

$\mathcal{I}_{\eta}(u,u^{\prime})=\frac{L^{\eta-1}}{4\pi\epsilon_{o}}\int_{0}^{2\pi}\nu_{\eta}(\sqrt{u^{2}+(u^{\prime})^{2}-2uu^{\prime}\cos(\beta)})d\beta.$

Fig. 12: Total charge.

The angular integral in the previous expression can be written in terms of
special functions, by following an analogous procedure than the one in Eq.
(1). The angular integral is

$\int_{0}^{2\pi}\frac{d\beta}{(\sqrt{u^{2}+(u^{\prime})^{2}-2uu^{\prime}\cos(\beta)})^{\eta}}=\frac{2\pi}{|u-u^{\prime}|^{\eta}}{}_{2}F_{1}\left(\frac{1}{2},\frac{\eta}{2},1;-\frac{4uu^{\prime}}{|u-u^{\prime}|^{2}}\right)$

where

${}_{2}F_{1}\left(a,b,c;\xi\right)=\frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_{0}^{1}\frac{t^{b-1}(1-t)^{c-b-1}}{(1-t\xi)^{a}}dt$

is Gauss hypergeometric function, and $\Gamma(z)$ is the gamma function. Once
more, the problem can be reduced to a linear set of algebraic equations

$\sum_{n=0}^{N}M^{(\eta)}_{mn}\sigma_{n}=\Phi_{m},$

with

$M^{(\eta)}_{mn}=\int_{U_{n}}\mathcal{I}_{\eta}(u_{m}^{(c)},u^{\prime})f_{n}(u^{\prime})u^{\prime}du^{\prime}=\frac{L^{\eta-1}}{4\pi\epsilon_{o}}\int_{U_{n}}\left[\frac{2\pi}{|u_{m}^{(c)}-u^{\prime}|^{\eta}}{}_{2}F_{1}\left(\frac{1}{2},\frac{\eta}{2},1;-\frac{4u_{m}^{(c)}u^{\prime}}{|u_{m}^{(c)}-u^{\prime}|^{2}}\right)\right]f_{n}(u^{\prime})u^{\prime}du^{\prime},$

being $U_{n}\in[0,R]$ the interval where the $n$-th basis function is
different from zero. The surface charge density is shown in Fig. 2 for several
$\eta$ values.

Another quantity of interest is the charge inside a circular region, that can
be computed from

$q(u)=2\pi\int_{0}^{u}\sigma(u^{\prime})du^{\prime}=2\pi\sum_{n=0}^{N}\sigma_{n}\int_{0}^{u}f_{n}(u^{\prime})u^{\prime}du^{\prime},$

where the total charge is $Q=q(R)$.

##  Appendix B The Method of Moments for two circular layers

Let us generalize the MoM approach for two concentric layers on the plane. By
studying the interplay between charges, we aim to elucidate the intricate
correlation between the surface charge density and the underlying Coulomb
interactions.

We proceed by writing the electrostatic integral formula for the scalar
potential

$\frac{1}{4\pi\epsilon_{o}}\int_{\Omega_{in}\cup\Omega_{out}}\frac{\sigma(\boldsymbol{r}^{\prime})}{|\boldsymbol{r}-\boldsymbol{r}^{\prime}|}d^{2}\boldsymbol{r}=\Phi(\boldsymbol{r}\in\Omega_{in}\cup\Omega_{out})=\begin{cases}V_{in}&\boldsymbol{r}\in\Omega_{in}\\\
V_{out}&\boldsymbol{r}\in\Omega_{out}\end{cases}$

with $V_{in}$ and $V_{out}$ the voltage of the inner and outer layers. One can
use the angular symmetry of the system to write

$\int_{u^{\prime}\in[0,R_{1}]\cup[R_{2},R_{3}]}\sigma(u^{\prime})I(u,u^{\prime})u^{\prime}du^{\prime}=\Phi(u,0,0)=\begin{cases}V_{in}&u\in[0,R_{1}]\\\
V_{out}&u\in[R_{2},R_{3}]\end{cases}$

with

$I(u,u^{\prime})=\frac{1}{4\pi\epsilon_{o}}\int_{0}^{2\pi}\frac{d\beta}{\sqrt{u^{2}+u^{\prime
2}-2uu^{\prime}\cos(\beta)}}=\frac{1}{4\pi\epsilon_{o}}\frac{4}{|u-u^{\prime}|}K\left(-\frac{4uu^{\prime}}{|u-u^{\prime}|^{2}}\right).$

The potential profiles at different heights of the system are shown in Fig.
13. We are interested in the surface density charge $\sigma$ on both layers.

Fig. 13: Illustrative potential profiles at different heights. Both layers are
at the same voltage $V_{in}=V_{out}=1V$ and their cross sections are
represented with black rectangles.

We define two basis: one having $\mathcal{N}$-basis functions on the inner
layer $\left\\{f_{n}^{(in)}(u)\right\\}_{n=1,\ldots,\mathcal{N}}$ and another
having $\mathcal{N}$-basis functions on the outer layer
$\left\\{f_{n}^{(out)}(u)\right\\}_{n=\mathcal{N}+1,\ldots,2\mathcal{N}}$.
Thus, the surface charge density can be computed from

$\sigma(u)=\sum_{n=1}^{\mathcal{N}}\sigma_{n}f_{n}^{(in)}(u)+\sum_{n=\mathcal{N}+1}^{2\mathcal{N}}\sigma_{n}f_{n}^{(out)}(u)=\sum_{n=1}^{2\mathcal{N}}\sigma_{n}f_{n}(u),$

with $f_{n}(u)=f_{n}^{(in)}(u)\hskip 7.11317pt\mbox{{if}}\hskip 7.11317pt1\leq
n\leq\mathcal{N},\hskip 7.11317pt\mbox{{otherwise}}\hskip
7.11317ptf_{n}^{(out)}(u)$. For simplicity, we shall take piece-wise functions
as bases

$f_{n}^{in}(u)=\begin{cases}1&u\in[u_{n-1},u_{n}]:=U_{n}^{(in)}\\\ 0&u\notin
U_{n}^{(in)}\end{cases}\hskip 14.22636pt\mbox{and}\hskip
14.22636ptf_{n}^{out}(u)=\begin{cases}1&u\in[u_{n},u_{n+1}]:=U_{n}^{(out)}\\\
0&u\notin U_{n}^{(out)}\end{cases}$

where node points are

$u_{n}=\begin{cases}u_{n}^{(in)}=n\frac{R_{1}}{\mathcal{N}}&n=0,\ldots,\mathcal{N},\\\
u_{n}^{(out)}=R_{2}+[n-(N+1)]\frac{R_{3}-R_{2}}{\mathcal{N}}&n=\mathcal{N}+1,\ldots,2\mathcal{N}+1.\end{cases}$

The basis functions are represented schematically in Fig. 14. Note that even
when the bases are split into two groups of functions, simultaneously solving
the problem for both layers $\Omega_{1}$ and $\Omega_{2}$ is required.

Fig. 14: Illustrative discrete bases. (left) Layers and nodes
$\\{u_{n}\\}_{n=1,\ldots,2\mathcal{N}+1}$ in purple color. (right) Basis
functions. Red points are located at the element’s center
$\\{u_{m}^{(c)}\\}_{m=1,\ldots,2\mathcal{N}}$.

Fig. 15: Illustrative electric potential on the $z=$ constant plane. Potential
$\Phi$ evaluated at $z=5R$, $z=0.3R$ and $z=0.025R$ (left to right). Both
layers are set to the same voltage $V_{in}=V_{out}=1V$.

Therefore,

$\sum_{n=1}^{2\mathcal{N}}\sigma_{n}\int_{u^{\prime}\in[0,R_{1}]\cup[R_{2},R_{3}]}f_{n}(u^{\prime})I(u,u^{\prime})u^{\prime}du^{\prime}=\Phi(u,0,0)=\begin{cases}V_{in}&u\in[0,R_{1}]\\\
V_{out}&u\in[R_{2},R_{3}]\end{cases}$

Since the radial distance $u$ is not constrained, we may evaluate the previous
equation at the centers
$\left\\{u_{m}^{(c)}\right\\}_{m=1,\ldots,2\mathcal{N}}$ of the elements
$U_{m}^{(in)}$ and $U_{m}^{(out)}$, hence

$\displaystyle\sum_{n=1}^{2\mathcal{N}}\left\\{\int_{u^{\prime}\in[0,R_{1}]\cup[R_{2},R_{3}]}f_{n}(u^{\prime})I(u_{m}^{(c)},u^{\prime})u^{\prime}du^{\prime}\right\\}\sigma_{n}$
$\displaystyle=\Phi(u_{m}^{(c)},0,0)=\begin{cases}V_{in}&u_{m}^{(c)}\in[0,R_{1}]\\\
V_{out}&u_{m}^{(c)}\in[R_{2},R_{3}]\end{cases}$
$\displaystyle\sum_{n=1}^{2\mathcal{N}}M_{mn}\sigma_{n}$
$\displaystyle=\Phi_{m}$

with the $M_{mn}$ matrix given by

$M_{mn}=\begin{cases}\int_{u_{n-1}}^{u_{n}}I(u_{m}^{(c)},u^{\prime})u^{\prime}du^{\prime}&n=1,\ldots,\mathcal{N}\\\
\int_{u_{n}}^{u_{n+1}}I(u_{m}^{(c)},u^{\prime})u^{\prime}du^{\prime}&n=N+1,\ldots,2\mathcal{N}.\\\
\end{cases}$

Fig. 16: Illustrative surface density charge (left) and integrated charge
profiles (right). Both sheets are at the same voltage $V_{in}=V_{out}=1V$ and
their cross sections are represented with black rectangles.

Fig. 16-left shows the profile of the surface charge density $\sigma$ when
both sheets are at same potential $V_{in}=V_{out}=1V$ with $R_{1}=R$,
$R_{2}=3R$ and $R_{3}=4R$. The potential in the space due to the layers at the
same voltage is shown in Fig. 15 and the integrated charge $q(u)$ of each
layer is shown Fig. 13-left. The last quantity is computed by

$q(u)=2\pi\begin{cases}\int_{0}^{u}\sigma(u^{\prime})du^{\prime}&\mbox{inner
sheet}\\\ \int_{R_{2}}^{u}\sigma(u^{\prime})du^{\prime}&\mbox{outer
sheet}\end{cases}$ (9)

with $Q$ the total charge of each layer, this is $Q_{in}=q(R_{1})$ and
$Q_{out}=q(R_{3})$ to normalize the plot. The ratio between total charge is
$Q_{out}/Q_{in}=9.06972$. The scalar potential in the $\mathbb{R}^{3}$ space
is computed from

$\displaystyle\Phi(\boldsymbol{r})$
$\displaystyle=\frac{1}{4\pi\epsilon_{o}}\int_{\Omega_{in}\cup\Omega_{out}}\frac{\sigma(\boldsymbol{r}^{\prime})}{|\boldsymbol{r}-\boldsymbol{r}^{\prime}|}d^{2}\boldsymbol{r}$
$\displaystyle=\frac{1}{4\pi\epsilon_{o}}\int_{u^{\prime}\in[0,R_{1}]\cup[R_{2},R_{3}]}\sigma(u^{\prime})u^{\prime}du^{\prime}\int_{0}^{2\pi}\frac{d\phi^{\prime}}{\sqrt{u^{2}+u^{\prime
2}-2uu^{\prime}\cos(\phi-\phi^{\prime})+z^{2}}}$
$\displaystyle=\frac{1}{4\pi\epsilon_{o}}\int_{u^{\prime}\in[0,R_{1}]\cup[R_{2},R_{3}]}\sigma(u^{\prime})u^{\prime}du^{\prime}\frac{4}{\sqrt{(u-u^{\prime})^{2}+z^{2}}}K\left(-\frac{4uu^{\prime}}{(u-u^{\prime})^{2}+z^{2}}\right).$

Here we may use a mid-point approximation, such that the potential takes the
form:

$\Phi(u,z)=\frac{1}{4\pi\epsilon_{o}}\sum_{n=1}^{2N}\sigma(u_{n}^{(c)})u_{n}^{(c)}\delta
u_{n}\frac{4}{\sqrt{(u-u_{n}^{(c)})^{2}+z^{2}}}K\left(-\frac{4uu_{n}^{(c)}}{(u-u_{n}^{(c)})^{2}+z^{2}}\right).$