[1] Iqraa Ehsan

[1]Department of Computer Science, Technical University of Kaiserslautern,
Kaiserslautern, 67663, Germany

2]Mindgarage, Technical University of Kaiserslautern, Kaiserslautern, 67663,
Germany

3]Computer Vision, German Research Institute for Artificial Intelligence
(DFKI), Kaiserslautern, 67663, Germany

#  End-to-End Semi-Supervised approach with Modulated Object Queries for Table
Detection in Documents

<EMAIL_ADDRESS>Tahira Shehzadi  Didier Stricker  Muhammad Zeshan Afzal  *
[  [

###### Abstract

Table detection, a pivotal task in document analysis, aims to recognize and
locate tables within document images precisely. Although deep learning has
shown remarkable progress in this realm, it typically requires an extensive
dataset of labeled data for proficient training. Current CNN-based semi-
supervised table detection approaches use the anchor generation process and
Non-Maximum Suppression (NMS) in their detection process, limiting training
efficiency. Meanwhile, transformer-based semi-supervised techniques adopted a
one-to-one match strategy that provides noisy pseudo-labels, limiting overall
efficiency. This study presents an innovative transformer-based semi-
supervised table detector. It improves the quality of pseudo-labels through a
novel matching strategy combining one-to-one and one-to-many assignment
techniques. This approach significantly enhances training efficiency during
the early stages, ensuring superior pseudo-labels for further training. Our
semi-supervised approach is comprehensively evaluated on benchmark datasets,
including PubLayNet, ICADR-19, and TableBank. It achieves new state-of-the-art
results, with a mAP of 95.7% and 97.9% on TableBank (word) and PubLaynet,
respectively. The results clearly show the superiority of our semi-supervised
approach, surpassing all existing state-of-the-art methods by substantial
margins. This research represents a significant advancement in semi-supervised
table detection methods, offering a more efficient and accurate solution for
practical document analysis tasks.

###### keywords:

Table Detection, Document Analysis, Semi-Supervised Learning, Detection
Transformer

##  1 Introduction

A visual summary is crucial for many applications, such as document summaries
and identifying graphical components in the visualization pipeline. Therefore,
a crucial stage in the summary and analysis of the document will be localizing
and identifying graphical elements. Tables are one of the integral elements in
document analysis, adeptly consolidating crucial information into a format
that maximizes spatial efficiency. Precise text extraction, especially within
intricate table layouts, remains a recurring challenge. Manually extracting
table objects has become impractical due to the growing number of documents.
Automated methodologies [mori1992historical, nguyen2021survey] offer
dependable, impactful, and streamlined resolutions for manual tasks. OCR
[singh2012survey] is of paramount importance in the field of document
analysis, particularly in the context of tables. It plays a crucial role in
addressing the complexities associated with handling structured data.
Conventional OCR[bourbakis1991ocr] systems often falter when confronted with
the structured format of tables, leading to misinterpretations, data loss, and
formatting disparities. Early attempts [mori1992historical] to address these
challenges relied on external metadata embedded within documents. However,
these methods proved inadequate for adapting to novel table configurations,
such as borderless tables.

Recent advancements in deep learning methodologies have transformed table
detection [bhatt2018optical, prasad2020cascadetabnet, shekar2021optical,
shehzadi2d], notably through the integration of object detection concepts from
computer vision. As technology advances, new challenges have emerged within
the domain of table detection. A notable challenge is that training deep
learning models necessitates extensive annotation of a substantial amount of
data, a task that demands substantial resources. While supervised approaches
[xue2022language, amrhein2018supervised] excel on standardized benchmarks,
their application in industrial contexts is according to the availability of
domain-specific annotated datasets. The process of manually generating or
preprocessing labeled data is not only time-consuming but also introduces a
risk of inaccuracies. As a result, there’s a growing demand for networks that
can effectively function with limited labeled datasets.

Semi-supervised learning presents an invaluable solution precisely in
situations where labeled data is limited. These techniques utilize both
labeled and unlabeled data, subsequently providing pseudo-labels for the
unlabeled instances. The goal of semi-supervised techniques is to improve
model performance and generalization abilities by utilizing the advantages of
both types of data. By doing so, it harnesses the inherent patterns and
structures in the unlabeled data, enhancing the model’s learning capacity and
potentially improving its accuracy, all while reducing the dependency on
extensive labeled datasets. Earlier, Semi-Supervised Object Detection (SSOD)
[softTeacher56, app98] relies heavily on object detectors containing anchor
generation process and non-maximum suppression (NMS) [hosang2017learning] for
post-processing. Later on, DETR [detr34] has been adapted for semi-supervised
settings, capitalizing on its inherent design which eliminates the need for
NMS and efficiently utilizes both labeled and unlabeled data. The integration
of semi-supervised learning in DETR presents challenges. DETR’s dependence on
a fixed set of object queries, established during training, may lead to
inefficiencies in scenes with varying object numbers, potentially causing
missed detections or false positives. Furthermore, the one-to-one assignment
can exacerbate errors when pseudo-labels are inaccurate. For instance, if a
pseudo-label, though not entirely precise, is incorporated into the training
process, the one-to-one assignment scheme may compel a predicted bounding box
to closely align with this inaccurate pseudo-label.

In response to the challenges addressed, we introduce a novel Semi-supervised
framework based on the Teacher-Student model. Our innovative matching module
integrates two training stages, utilizing both one-to-many [fang2023feataug]
and one-to-one [li2023one, harris1986teaching] assignment strategies. The
former ensures NMS-free end-to-end detection [hosang2017learning] advantages,
while the latter generates a higher-quality candidate proposal set,
contributing to a more efficiently optimized detector. In our semi-supervised
learning approach, the student and teacher modules are equipped to handle both
one-to-many [fang2023feataug] and one-to-one [li2023one, harris1986teaching]
queries. This unified processing entails directing object queries from the
weakly augmented images to the teacher module and those from strongly
augmented images to the student module. Furthermore, in the one-to-many
matching scheme, we implement augmented ground truth, resulting in the
duplication of instances for high-quality pseudo labels.

In summary, this paper presents significant contributions in the following:

  * •

We introduce a novel end-to-end DETR-based semi-supervised approach for table
detection. The approach utilizes a dual query assignment mechanism that
integrates the one-to-many matching strategy with the one-to-one matching
strategy in the training phase.

  * •

We introduce a query filtering module in a one-to-many matching strategy to
provide high-quality pseudo-labels. To the best of our knowledge, this
research is the first to focus on object queries in a semi-supervised setting
for table detection tasks.

  * •

We conduct a comprehensive evaluation using diverse datasets, including
PubLayNet [PubLayNet3], ICDAR-19 [icdar19], and TableBank [tablebank8]. Our
approach surpasses both CNN-based and transformer-based semi-supervised
approaches due to its superior efficiency and ability to generate high-quality
pseudo-labels without using Non-maximum Suppression (NMS), resulting in a
comparable performance with improved accuracy.

The subsequent sections of the paper are organized in the following manner.
Section 2 provides an overview of prior research conducted in the field of
detecting table objects inside document images. Section 3 explains the
approach, while Section 4 provides a brief discussion of the dataset. We
examine the experimental outcomes in Section 5. The ablation study is covered
in Section 6. Finally, Section 7 provides a comprehensive overview of the
experimental work conducted and offers insights into potential future research
topics.

Figure 1: Comparing our approach, which utilizes the TableBank-Both and
Publaynet datasets, with the prior method, both supervised and semi-
supervised, we conducted experiments with a 30$\%$ data split. The semi-
supervised deformable transformer is referred to as the Def. Semi
[shehzadi2023towards]. For an extensive summary of the results, please refer
to Table 6.

##  2 Related Work

The identification and extraction of tables in document image analysis is a
crucial and significant task. Numerous studies have introduced various methods
for identifying tables with diverse structures in document images. In the
past, a significant number of these methodologies either depended on pre-
established rules or required further meta-data input to tackle table
detection tasks [tupextract3, DPmatch4]. In recent times, researchers have
increasingly utilized statistical methodologies [tsr87] and included deep
learning techniques to augment the flexibility and adaptability of table
detection systems [DeepDeSRT3, CasTab45, Hyb65]. This section provides an in-
depth overview of these methodologies and offers insights into the strategies
used for semi-supervised object detection based on Convolutional Neural
Networks (CNNs).

###  2.1 Rule-based Strategies

Itonori et al. [tsk5] pioneered table detection in document images by
representing tables as text blocks with predefined rules. Subsequent studies
by the same authors introduced a method considering both horizontal and
vertical lines [strctRTB3]. Pyreddy et al. [TINTIN67] used tailored heuristics
for tabular section detection, while Pivk et al. [PIV67] transformed HTML
formatted tables into logical forms. Hu et al. [DRR6] proposed a method based
on identifying white regions and vertically connected elements. Rule-based
methods, effective for specific formats, are discussed in [RecogTable5,
extractTab9, TSsurvey32]. While effective, rule-based strategies may lack
universal applicability, highlighting the need for more adaptable systems in
the challenge of detecting tables in document images.

###  2.2 Learning-Based Methodologies

Cesarini et al.’s [trainTD5] supervised learning approach is intended to
reliably detect and categorize table elements within document images. It
entails converting an image of a document into an MXY tree structure, then
recognizing and categorizing individual components as tables according to the
presence of vertical and horizontal lines. To detect tables, researchers
utilize Hidden Markov Models [RichMM3, FTabA4] and a Support Vector Machine
(SVM) classifier with traditional heuristics [lineln6]. While machine learning
outperforms rule-based approaches, it depends on additional information like
ruling lines. In contrast, deep learning-based methodologies show higher
precision and effectiveness. These techniques can be categorized into object
identification, semantic segmentation, and bottom-up approaches.

Approaches based on semantic segmentation. Table detection is framed as a
segmentation problem in techniques like [Xi17, He761, Ik36, Paliw9],
leveraging existing semantic segmentation networks to create pixel-level
segmentation masks. These masks are then consolidated for the ultimate table
detection results, showcasing superior performance over conventional methods
on benchmark datasets [icdar19, PubLayNet3, tablebank8, pubtables5]. In Yang
et al.’s work [Xi17], a fully convolutional network (FCN) [FCNseg4] integrates
linguistic and visual features, enhancing segmentation for tables and other
page elements. In order to generate segmentation masks for text and table
sections as well as their outlines, He et al. [He761] suggested a multi-scale
FCN that was further improved for the final table elements.

Inductive Approach The approaches employed in this study consider identifying
tables as a problem of graph labeling, whereby the nodes of the graph
represent various page components and the edges indicate the relationships
between them. In their study, Li et al. [Li64] employed a conventional layout
analysis approach to extract line regions from a given document. These line
areas were further categorized into four distinct classes, including text,
figure, formula, and table. This classification process was accomplished
through the utilization of two CNN-CRF networks. The researchers made further
predictions on the cluster that would correlate to pairs of line areas. In
contrast, Holecek et al. [Martin66] addressed the analysis of text sections by
seeing them as individual nodes within a graph structure, enabling the
identification and understanding of the overall document architecture. Graph-
neural networks were later utilized for node-edge categorization. It is
important to acknowledge that these methodologies are contingent upon certain
assumptions, such as the incorporation of text line boxes as supplementary
input.

Approaches Based on Object Detection. The identification and localization of
tables inside document images have often been addressed as a task of object
detection, initially employing R-CNN methods like those by Yi et al. [Yi77]
and Hao et al. [Hao789], relying on heuristic rules. Advancements introduced
efficient single-stage [retinaNet68], [yolos6] and two-stage [fast15],
[faster23], [mask86], [cascadercnn8] detectors for document elements. Some
studies [Azka62, Ayan29, arif48] applied image transformations for
improvement, while Siddiqui et al. [Sidd32] enhanced Faster R-CNN with
deformable-convolution and RoI-Pooling for geometrical modifications. The
researchers Agarwal et al. [Agarwal52] employed a composite network using
deformable convolution to augment the performance of Cascade R-CNN. In
contrast, the presented semi-supervised technique, treating detection as a set
prediction task [bridging_sheh], eliminates stages like proposal generation
and post-processing (e.g., NMS [hosang2017learning]), simplifying and
streamlining the detection process for increased efficiency.

###  2.3 Object Detection with Semi-supervised

In the domain of object detection’s semi-supervised learning, there are two
primary categories: consistency-based methods [consemi6, propsemi8] and
pseudo-label generation-based methods [omnisup8, rethink98, semimask1,
simplesemi76]. The approach we employ can be classified as pertaining to the
pseudo-label. In previous studies [omnisup8, rethink98], pseudo-labels for
unlabeled data were generated by combining prediction results using various
data augmentation techniques. Another method [selfsup6] trained a SelectiveNet
to create pseudo-labels, involving transferring bounding boxes from
unannotated to annotated data and assessing localization consistency. However,
this complex detection process manipulates images. STAC’s technique
[simplesemi76] uses robust augmentation for pseudo-labels and less intense
augmentation for model training. Our semi-supervised table identification uses
a deformable transformer in an end-to-end framework [Deformable54], employing
a multi-level training strategy to eliminate anchor generation and post-
processing like NMS [hosang2017learning]. Our method revolutionizes semi-
supervised learning by eliminating the need for NMS, a traditional post-
processing step for removing overlapping bounding boxes. Simultaneously, it
enhances pseudo-label precision, making them more effective in semi-supervised
scenarios, leading to improved model performance and reduced annotation
efforts.

##  3 Methods

###  3.1 PRELIMINARY

The DETR model [detr34], which is commonly used for applications such as
object detection, consists of an encoder, a transformer decoder, and
predictors for object class and box position[detr34, att75]. The process
initiates by passing an input image through the backbone and transformer
encoder [att75], resulting in the extraction of image features. The decoder
subsequently takes the object queries and image features, the input consists
of a matrix $Q$ which is written as $\left\\{q_{1},q_{2},...,q_{N}\right\\}$.
The DETR model utilizes task-specific prediction heads on the altered
embeddings of object queries ($\tilde{Q}$) at every transformer decoder layer.
These prediction heads are then followed by the predictors, resulting in an
output implied as $Y=\left\\{y_{1},y_{2},...,y_{N}\right\\}$.

The decoder utilizes self-attention mechanisms on individual object queries,
facilitating interactions between queries to gather information pertinent to
duplicate detection. Then, it utilizes cross-attention between queries and
image features to extract relevant data from the image features for object
detection. Lastly, it incorporates a feed-forward network that independently
analyzes the queries to improve the recognition of entities.

During the training process of the DETR model [detr34], the predictions are
structured assets and do not possess a direct connection with ground-truth
objects. The DETR employs a one-to-one strategy [harris1986teaching],
establishing a bipartite matching between the predictions and the ground-truth
objects:

$\left(y_{\sigma\left(1\right)},\bar{y}_{1}\right),\left(y_{\sigma\left(2\right)},\bar{y}_{2}\right),...,\left(y_{\sigma\left(N\right)},\bar{y}_{N}\right)$
(1)

In this case, the optimal permutation of $N$ indices is denoted by
$\sigma\left(\cdot\right)$, whereas the ground truth is denoted by
$\bar{Y}=\left\\{\bar{y}_{1},\bar{y}_{2},...,\bar{y}_{N}\right\\}$. Hence, the
expression of the loss function is formulated as:

$\mathcal{L}=\sum_{n=1}^{N}\mathscr{l}\left(y_{\sigma\left(n\right)},\bar{y}_{n}\right)$
(2)

where the classification loss and the box regression loss between the
predicted value $y$ and the ground truth object $\bar{y}$. The objective of
optimization using one-to-one assignment is to evaluate the accuracy of
predictions by prioritizing one prediction for each ground-truth object while
simultaneously reducing the prominence of duplicate predictions. This scoring
process leverages the comparison of a single prediction against several
others, and the relevant information is provided by the decoder’s self-
attention over queries.

In this end-to-end process, each ground-truth is associated with the
corresponding prediction having the least matching cost, and the relevant
supervision is applied accordingly. The utilization of a one-to-one matching
mechanism in this process imposes limitations, resulting in suboptimal
efficiency in training. One-to-one matching is a critical aspect of DETR,
pivotal in enabling its end-to-end capability, thereby obviating the need for
manual NMS [hosang2017learning] to eliminate duplicate detections. The
versatility of DETR relies on its end-to-end characteristic, which has been
extended to encompass a broader spectrum of vision tasks. However, it’s
important to note that only a restricted number of queries are designated as
positive samples, and the training effectiveness of these positives is
significantly diminished by the one-to-one matching.

Figure 2: The decoder utilizes both one-to-one and one-to-many assignment
methods, passing queries through multiple transformer decoder layers to
generate high-quality pseudo-labels, thereby improving object detection
precision. In one-to-one matching, the initial query set establishes a direct
link between predictions and true objects. In one-to-many matching, a second
set of queries results in multiple sets of predictions, achieved by
replicating ground truth information.

###  3.2 Dual Assignment Strategies for Pseudo-Label Generation

The one-to-one query matching approach [harris1986teaching] has its set of
challenges. Tables, given their varied sizes, structures, and content
densities, may not always align perfectly with a rigid query mechanism. The
inherent constraint of one-to-one matching can lead to suboptimal table
identifications, especially when the number of tables in a document doesn’t
match the designated query count. As a result, certain tables might remain
undetected, or there may be wasteful queries without corresponding tables.
When combined with pseudo-label generation in a semi-supervised setup, this
constraint might amplify inaccuracies. If the initial model makes a suboptimal
match due to the one-to-one limitation, the generated pseudo-label for an
unlabeled table might be erroneous. This can subsequently misguide the model
during further training, potentially hampering its effectiveness in detecting
tables across diverse documents.

To address the limitations of one-to-one matching, our approach combines both
one-to-one and one-to-many matching strategies as shown in Figure 2. This
approach combines the strengths of both techniques to improve the overall
efficiency of the object detection model. Within the matching strategies, two
query sets are upheld:
$Q=\left\\{q_{1},q_{2},\textperiodcentered\textperiodcentered\textperiodcentered,q_{n}\right\\}$
and
$\hat{Q}=\left\\{\hat{q_{1}},\hat{q_{2}},\textperiodcentered\textperiodcentered\textperiodcentered,\hat{q_{T}}\right\\}$.
Subsequently, predictions based on $Q$ or $\hat{Q}$ undergo either one-to-one
or one-to-many matching techniques.

In one-to-one matching, the initial query set, $Q$, traverses through
transformer decoder layers [att75], generating predictions for each layer’s
output. Following this, bipartite matching is conducted between the pairs of
predictions and ground truth across every layer. The loss is then determined
using the formula:

$\vspace{-0.8em}\mathcal{L}_{o2o}=\sum_{l=1}^{L}\mathcal{L}_{H}\left({y}_{p},y_{g}\right)$
(3)

Here, ${y}_{p}$ signifies the outcomes generated by the transformer decoder
layer at position l, while $y_{g}$ denotes the ground-truth. The optimal
assignment is efficiently calculated through the application of the Hungarian
algorithm. One-to-one matching strategies are applied during both the model’s
training and testing phases.

In the one-to-many matching, the second set of queries, $\hat{Q}$, go through
the same $L$ transformer decoder layers, resulting in $L$ sets of predictions.
To enable one-to-many matching [fang2023feataug], the ground-truth was
replicated $K$ times, leading to an extended target
$\hat{y_{g}}=\left\\{y_{g}^{1},y_{g}^{2},\textperiodcentered\textperiodcentered\textperiodcentered,y_{g}^{K}\right\\}$.
Following this, bipartite matching was conducted between the predictions and
the augmented ground truth pairs across each layer. The loss is then
determined using the formula:

$\vspace{-0.8em}\mathcal{L}_{o2m}=\sum_{l=1}^{L}\mathcal{L}_{H}\left(\hat{y}_{p},\hat{y}_{g}\right)$
(4)

Here, $\hat{y}_{p}$ signifies the outcomes generated by the transformer
decoder layer at position l. The Hungarian algorithm is employed for the
proficient computation of the best assignment. The one-to-many matching
strategy is exclusively applied during the model training phase.

By integrating these matching strategies, the model gains the flexibility to
adapt to diverse object configurations. It can accurately handle both
scenarios where objects are distinctly separate and situations where they
exhibit complex interactions. Our approach represents a sophisticated solution
that balances the advantages of one-to-one [harris1986teaching] and one-to-
many matching [fang2023feataug], resulting in a more robust and adaptable
object detection system. This capability allows the model to perform very well
in various practical situations, therefore rendering it a valuable asset in
the field of computer vision applications.

###  3.3 Semi-Supervised Learning

In this section, we introduce an enhanced approach for semi-supervised object
detection built upon the DETR framework [detr34]. For a detailed visual
representation, please refer to Figure 3. This innovative method smoothly
incorporates a transformer architecture while also elaborating on a thorough
pseudo-labeling approach. Our approach embodies a comprehensive learning
paradigm that strategically leverages both annotated and unannotated data to
enhance object detection performance. During training, we have access to a
labeled image set $D_{l}=\left\\{x^{l}_{i},y^{l}_{i}\right\\}^{N_{l}}_{i=1}$
and an unlabeled image set $D_{u}=\left\\{x^{u}_{i}\right\\}^{N_{u}}_{i=1}$,
denoted as $N_{l}$ and $N_{u}$ respectively, indicating the number of images
with labels and those without labels. For the labeled images $x_{l}$, the
annotations $y_{l}$ offer information about bounding box coordinates and
object categories.

Figure 3: Our Semi-supervised approach: It focuses on combining one-to-one and
one-to-many matching strategies in a semi-supervised setting and incorporates
two unique categories of data for training: labeled data and unlabeled data.
The framework consists of two modules: one for students and one for teachers.
The teacher module exclusively handles unlabeled images, employing weak
augmentation prior to input. The student module is designed to handle both
labeled and unlabeled images, applying strong augmentation techniques
specifically to the unlabeled images. During the training process, the student
module utilizes an EMA technique to consistently update the teacher module.

The methodology comprises two fundamental components: the student module and
the teacher module. The student module handles both labeled and unlabeled
images, primarily applying intensive augmentation to the unlabeled data and
employing a combination of strong and weak augmentation techniques for the
labeled data. It undergoes training utilizing detection losses derived from
both labeled and unlabeled data, guided by pseudo-boxes. Pseudo-boxes with
class labels and their associated bounding boxes are attached to the unlabeled
data. By comparison, the teacher’s module only works with unlabeled, weakly
augmented images. In contrast, the teacher module exclusively handles weakly
augmented, unlabeled images. It generates pseudo-labels, which are
subsequently employed for supervised training by the student module. Weak
augmentation for the teacher module enhances the precision of the generated
pseudo-labels, whereas strong augmentation for the student module fosters a
more intricate learning process. Additionally, the student module processes a
limited set of labeled images, subjecting them to both strong and weak
augmentation. Both modules initiate the training process with randomized
parameters. During the training process, the student module integrates pseudo-
labels supplied by the teacher, prioritizing those that demonstrate high
confidence scores. The student’s parameters are updated through back-
propagation, while an EMA [zhao2023does] technique consistently updates the
teacher module. The training process unfolds in two distinct phases. In the
initial phase, the student module is trained using data with labels and
simultaneously generates pseudo-labels through the teacher module. Following
this, both modules are combined to produce the final predictions.

Our primary contribution focuses on combining one-to-one and one-to-many
matching strategies in a semi-supervised setting, a pivotal advancement in
addressing the core challenges of DETR-based Semi-Supervised Object Detection
(SSOD) [zhang2023semi]. This technique generates multiple positive proposals
for each pseudo-label, refining high-quality suggestions and expediting
convergence. Post-training, refined pseudo-labels are acquired, reducing
unnecessary predictions. This approach provides adaptability to various object
configurations, handling both separate and complex interactions. It achieves a
balance between one-to-one [harris1986teaching] and one-to-many matching
[fang2023feataug], resulting in a more robust and versatile object detection
system.

In the Pseudo-Labeling Framework, we adopt a conventional method from SSOD to
assign pseudo-labels to unlabeled data using the teacher module’s output.
Object detectors typically produce a confidence score vector $c_{k}$ for each
bounding box $b_{k}$. By applying a threshold to these scores, pseudo-labels
can be created. In a basic pseudo-labeling filter, if the confidence value
$c_{k}\in\left[0,1\right]^{G_{i}}$ for the ground-truth class $g_{k}$ is not
higher than the predicted value, as the pseudo-label, the prediction with the
highest value is utilized.

In our approach, we conceptualize the assignment of pseudo-labels as a task
involving both one-to-one and one-to-many matching, where the predictions from
the teacher module are aligned with the generated semi-labels. Achieving a
one-to-one assignment involves implementing the Hungarian algorithm based on
the student’s predictions and the pseudo-labels generated by the teacher.

$\vspace{-0.8em}{\rho}_{o2o}=\underset{\sigma\in\zeta_{N}}{argmin}\sum_{i=1}^{N}\mathcal{L}_{H}\left({y}_{i}^{T},{y}_{\sigma\left(i\right)}^{S}\right)$
(5)

Here, $\zeta_{N}$ represents the collection of permutations involving $N$
elements. Additionally,
$\mathcal{L}_{H}\left({y}_{i}^{T},{y}_{\sigma\left(i\right)}^{S}\right)$
stands for the cost incurred during the matching process between the pseudo
labels ${y}^{T}$ and the prediction made by the student identified by the
index $\sigma\left(i\right)$.

In the early stages of self-supervised object detection (SSOD) training, it is
typical to encounter incorrect and unreliable pseudo-labels that are generated
by the teacher. The utilization of the one-to-one assignment approach
[harris1986teaching] has the potential for generating proposals that are both
limited and of low quality. To improve the effectiveness of semi-supervised
learning, we suggest shifting from the one-to-one assignment approach to the
one-to-many assignment method. [fang2023feataug].

${\rho}_{o2m}=\left\\{\underset{\sigma_{i}\in\chi}{argmin}\sum_{i=1}^{N}\mathcal{L}_{H}\left({y}_{i}^{T},{y}_{\sigma_{i}\left(k\right)}^{S}\right)\right\\}_{i=1}^{\left|{y}^{T}\right|}$
(6)

where $\chi$ is used to represent the assignment of a subset of proposals to
each pseudo box ${y}_{i}^{T}$. In the one-to-one matching branch, the initial
query set passes through transformer decoder layers [att75], generating
predictions for each layer’s output. In the one-to-many matching, a second set
of queries traverses the same transformer decoder layers [att75], resulting in
sets of predictions. To facilitate one-to-many matching [fang2023feataug], the
ground truth is replicated K times, leading to an extended target. Following
this, bipartite matching is carried out between the predictions and the
augmented ground truth pairs across each layer.

Furthermore, two losses are concurrently integrated throughout the training
process to enhance training speed in parallel. Additionally, they incorporate
masked multi-head self-attention to prevent interference, incurring no
additional training cost during experimentation. The process of choosing
correspondences between the predictions of the teacher module and the semi-
labels, achieved by utilizing a threshold, parallels the heuristic criteria
for matching proposals or anchors by utilizing CNNs in object detectors that
capture actual objects [bappy2016cnn].

A significant contribution is our matching strategy, effectively addressing
core challenges in DETR-based SSOD [zhang2023semi]. This leads to high-quality
pseudo-labels, improving object detection accuracy. Furthermore, transitioning
from one-to-one [harris1986teaching] to one-to-many [fang2023feataug]
assignment enhances the efficiency of semi-supervised learning. Through
careful integration of loss functions and attention mechanisms, our approach
optimizes training speed without incurring additional costs. This refined
methodology assures a more robust and adaptable object detection system,
proficient in handling diverse object configurations.

##  4 Experimental Setup

###  4.1 Dataset

TableBank: TableBank [tablebank8] is a significant dataset in document
analysis, being the second-largest dedicated to table recognition. It
comprises 417,000 document images, meticulously annotated using the arXiv
database. These images include tables from LaTeX (253,817), Word (163,417),
and a combined set, totaling 417,234 images.

In our experiments, we exclusively used the TableBank dataset due to its
comprehensive coverage of diverse document formats and table structures. This
deliberate choice allowed us to engage with a dataset of substantial depth,
poised to drive advancements in table detection within document analysis.

Figure 4: Comparative Performance Analysis of our approach using the
Tablebank, Publaynet, and Icdar datasets. We employed a ResNet-50 backbone and
conducted experiments with three distinct data splits: 10$\%$, 30$\%$, and
50$\%$. The metric mAP denotes the mean Average Precision within the IoU
threshold range of 50$\%$ to 95$\%$.

PubLayNet: The PubLayNet dataset [PubLayNet3] is a substantial and publicly
available resource, containing 335,703 training images, along with 11,240 and
11,405 images in the validation and test sets, respectively. It is richly
annotated with polygonal segmentation and bounding boxes, encompassing various
elements like figures, list titles, tables, and textual content from diverse
research papers. Evaluation of this dataset employs the coco analytical
technique for a rigorous and standardized assessment. In a specific
experiment, the focus was on analyzing and utilizing 102,514 out of the total
86,460 table annotations in PubLayNet. This targeted approach provided a
detailed examination and validation of table-related elements, contributing to
a deeper understanding of document layout and structure.

ICDAR-19: The 2019 ICDAR competition marked a significant milestone with the
introduction of the Table Detection and Recognition (cTDaR) track [icdar19].
This track introduced two meticulously curated datasets, one contemporary and
one historical, dedicated to serving as benchmarks for table detection in the
competition’s TRACK A. To facilitate a direct comparison with established
state-of-the-art methods, we present our results on these datasets using
Intersection over Union (IoU) levels ranging from 0.5 to 0.9. This evaluation
method provides crucial insights into the model’s accuracy in table detection
tasks, providing a detailed understanding of its performance across different
degrees of overlap between predicted and ground truth bounding boxes.

###  4.2 Evaluation Criteria

We employ distinct assessment measures to evaluate the efficacy of our semi-
supervised table identification methodology. This section provides an overview
of the evaluation measures utilized, specifically precision, recall, and
F1-score. Precision [pr61] is a metric that measures the ratio of accurately
detected instances (True Positives) to instances that were mistakenly
identified as positive (False Positives) concerning the True Positives.
According to [pr61], the recall metric measures the proportion of correctly
recovered True Positives to the total of True Positives and False Negatives.
By computing their harmonic mean [pr61], the F1-score is a metric that
measures how well Precision and Recall are balanced. Furthermore, the IoU is
computed by dividing the intersection of the ground-truth box $B_{g}$ and the
expected bounding box $B_{p}$ by their union.

$IoU=\frac{area(B_{g}\cap B_{p})}{area(B_{g}\cup B_{p})}$ (7)

The calculation of IoU is utilized to ascertain the classification of a
detected table object as either a true positive or a false positive. The
computation of the average precision (AP) involves the utilization of a
precision-recall (PR) curve, which is consistent with the evaluation
methodology employed in MS COCO [coco14]. This involves computing the area
under the PR curve, as expressed by the subsequent formula:

$\vspace{-0.8em}AP=\sum_{k=1}^{N}(Re_{k+1}-Re_{k})P_{intr}(Re_{k+1})$ (8)

Figure 5: Comparative Visualization of Predictions of our Approach. It
highlights that the inclusion of both one-to-one and one-to-many queries
enhances the model’s accuracy. The ground truth is depicted by the blue boxes,
while the green boxes display the results obtained through our approach.

Here, Re1, Re2, . . . , $Re_{k}$ denote the recall parameters. The mean
average precision (mAP) serves as a common metric for assessing detection
method performance. It is determined by averaging the individual average
precisions across all classes within a dataset. However, it’s crucial to
recognize that alterations in the accuracy of individual classes, resulting
from class mapping, can impact the mean Average Precision (mAP). This is a
recognized limitation of this metric. For our evaluation, we establish IoU
threshold values at 0.5 and 0.95. The calculation of mAP is as follows:

$\vspace{-0.8em}mAP=\frac{1}{S}\sum_{s=1}^{S}AP_{s}$ (9)

Where S denotes the overall number of classes.

###  4.3 Implementation Details

We employ the one-to-one and one-to-many matching with a ResNet-50 [resnet45]
backbone, pre-trained on the ImageNet [NId6] dataset, as our detection model
for assessing the effectiveness of the semi-supervised approach. Training is
conducted on PubLayNet, ICDAR-19, and all three partitions of the TableBank
dataset, utilizing 10$\%$, 30$\%$, and 50$\%$ of labeled data, with the
remainder designated as unlabeled. The pseudo-labeling threshold is set at
0.7. Training spans 150 epochs, with a learning rate reduction by a factor of
0.1 at the 140th epoch. Strong augmentations including horizontal flips,
resizing, patch removal, cropping, grayscale, and Gaussian blur are applied,
alongside weak augmentation through horizontal flipping. The optimal
configuration for the number of queries to the decoder input for one-to-one
matching (N) is set at 30. Additionally, 400 queries for one-to-many matching
(T) are utilized, and the ground truth K is replicated six times, all
contributing to the optimal outcome. Unless otherwise specified, model
evaluations employ mAP (AP50:95) metrics. The weights $\alpha_{1}$ and
$\alpha_{2}$ for balancing the classification loss ($L_{cls}$) and regression
loss ($L_{box}$) are set at 2 and 5, respectively. To expedite training, input
image dimensions are set to 600 pixels, with the standard size of 800 pixels
retained for comparative purposes with other methodologies. Figure 5
illustrates the enhanced table detection capabilities in documents achieved
through our approach as opposed to alternative methods.

Dataset | Label-splits | mAP | AP50 | AP75 | ARL
---|---|---|---|---|---
TableBank-word | 10$\%$ | 84.1 | 87.1 | 85.5 | 96.2
30$\%$ | 95.7 | 96.9 | 96.2 | 98.1
50$\%$ | 95.5 | 98.5 | 97.2 | 98.7
TableBank-latex | 10$\%$ | 78.0 | 87.5 | 85.5 | 90.9
30$\%$ | 88.0 | 94.3 | 93.1 | 94.1
50$\%$ | 90.9 | 97.2 | 95.9 | 95.9
TableBank-both | 10$\%$ | 89.4 | 94.9 | 93.6 | 95.8
30$\%$ | 93.3 | 96.8 | 95.8 | 97.1
50$\%$ | 95.1 | 98.1 | 97.1 | 97.5
Table 1: Evaluation of our semi-supervised approach across TableBank-word, TableBank-latex, and TableBank-both data divisions with varying proportions of labeled data. In this context, mAP signifies the mean Average Precision within the IoU range of (50:95), AP50 stands for the Average Precision at an IoU threshold of 0.5, and AP75 represents the Average Precision at an IoU threshold of 0.75. ARL, on the other hand, stands for the Average Recall for larger objects. Method | Technique | Detector | 10$\%$ | 30$\%$ | 50$\%$
---|---|---|---|---|---
S. Ren [ren2015faster]. | Sup | Faster R-CNN | 80.1 | 80.6 | 83.3
X. Zhu [Deformable54]. | Sup | Def. DETR | 80.8 | 82.6 | 86.9
Sohn, K. [simplesemi76] | Semi-sup | Faster R-CNN | 82.4 | 83.8 | 87.1
Liu, Y. [liu2021unbiased] | Semi-sup | Faster R-CNN | 83.9 | 86.4 | 88.5
Tang, Y. [tang2021humble] | Semi-sup | Faster R-CNN | 83.4 | 86.2 | 87.9
Xu, M. [softTeacher56] | Semi-sup | Faster R-CNN | 83.6 | 86.8 | 89.6
Shehzadi, T.[shehzadi2023towards] | Semi-sup | Def. DETR | 84.2 | 86.8 | 91.8
Our | Semi-sup | Def. DETR | 94.9 | 96.8 | 98.1
|  |  |  |  |
Table 2: Comparison of the performance of previous supervised and semi-
supervised methodologies. The ResNet-50 backbone was used to train the
supervised models, Deformable-DETR(Def. DETR) and Faster R-CNN, using 10$\%$,
30$\%$, and 50$\%$ of the TableBank dataset. The same splits and datasets are
utilized for semi-supervised assessment. It’s important to note that all
results are given in terms of mAP (0.5:0.95).

##  5 Results and Discussion

###  5.1 TableBank

This section presents the experimental results that include all partitions of
the TableBank dataset [tablebank8], taking into account different proportions
of labeled data. We perform a comprehensive analysis, contrasting our
transformer-based semi-supervised approach [att75] with prior methodologies
encompassing both supervised and semi-supervised techniques based on deep
learning. Table 1 displays the results obtained from our approach, which
leverages the hybrid deformable transformer [Deformable54] across the
TableBank-word, TableBank-latex, and TableBank-both data partitions. These
experiments encompass 10$\%$, 30$\%$, and 50$\%$ labeled data, with the
remaining portion designated as unlabeled data. For a comparative analysis
with other datasets employed, please refer to Figure 4. Notably, the
TableBank-word data split demonstrates the highest AP50 score of 98.5$\%$. For
a detailed assessment of the precision in comparison with other splits, please
consult and 6(a). Following this, TableBank-latex achieves 97.2$\%$ at 50$\%$
labeled data. For a detailed evaluation of the precision in comparison with
other splits, please consult Figure 6(b). Lastly, TableBank-both achieves
98.1$\%$ at 50$\%$ labeled data, for a comparative analysis between precision,
please refer to Figure 6(c).

In this study, we compare our semi-supervised approach, which incorporates the
deformable transformer [Deformable54]. The supervised model is trained on
several subsets of labeled data, specifically 10$\%$, 30$\%$, and 50$\%$,
inside the TableBank-both data split. The results are shown in Table 2.
Significantly, our proposed methodology, which utilizes an attention
mechanism, demonstrates comparable outcomes in terms of quality. This
eliminates the need for generating proposals and performing additional post-
processing procedures like NMS[hosang2017learning].

Figure 6: Precision-Based Analysis of our Model’s Performance. This assessment
covers four distinct datasets: (a) Tablebank-Word, (b) Tablebank-Latex, (c)
Tablebank-Both and (d) Publaynet. Utilizing ResNet-50 as the foundational
architecture, experiments were conducted across three different data splits:
10$\%$, 30$\%$, and 50$\%$. The metric mAP, denoting mean Average Precision
within the IoU threshold range of 50$\%$ to 95$\%$, was employed in the
evaluation.

###  5.2 PubLayNet

In this section, we critically examine the experimental findings pertaining to
the dataset for the table class in PubLayNet[PubLayNet3], considering varying
proportions of labeled data. Furthermore, we perform a comprehensive
comparative analysis between our transformer-based [att75] semi-supervised
approach and earlier supervised and semi-supervised methodologies grounded in
deep learning.

Dataset | Label-percent | mAP | AP50 | AP75
---|---|---|---|---
publaynet | 10$\%$ | 97.4 | 98.5 | 98.3
30$\%$ | 97.9 | 98.7 | 98.5
50$\%$ | 98.6 | 98.9 | 99.3
Table 3: Evaluation Results for the Table Class in the PubLayNet Dataset. In this context, mAP signifies the mean Average Precision within the IoU threshold range of (50:95), AP50 corresponds to the Average Precision at an IoU threshold of 0.5, and AP75 signifies the Average Precision at the IoU threshold of 0.75. ARL denotes the Average Recall specifically for larger objects. Method | Technique | Detector | 10$\%$ | 30$\%$ | 50$\%$
---|---|---|---|---|---
S. Ren [ren2015faster] | Sup | Faster R-CNN | 93.6 | 95.6 | 95.9
X. Zhu [Deformable54] | Sup | Def. DETR | 93.9 | 96.2 | 97.1
Sohn, K. [simplesemi76] | Semi-sup | Faster R-CNN | 95.8 | 96.9 | 97.8
Liu, Y. [liu2021unbiased] | Semi-sup | Faster R-CNN | 96.1 | 97.4 | 98.1
Tang, Y.[tang2021humble] | Semi-sup | Faster R-CNN | 96.7 | 97.9 | 98.0
Xu, M. [softTeacher56] | Semi-sup | Faster R-CNN | 96.5 | 98.1 | 98.5
Shehzadi, T.[shehzadi2023towards] | Semi-sup | Def. DETR | 98.5 | 98.8 | 98.8
Our | Semi-sup | Def. DETR | 98.5 | 98.7 | 98.9
|  |  |  |  |
Table 4: The performance evaluation contrasts earlier supervised and semi-
supervised methodologies. Deformable-DETR (Def. DETR) and Faster R-CNN were
exclusively trained on 10$\%$, 30$\%$, and 50$\%$ of table data, whereas semi-
supervised networks utilized the same proportions of the PubLayNet dataset,
designating them as labeled, while the remainder was considered unlabeled. The
reported results are uniformly presented in terms of AP50 at an IoU threshold
of 0.5.

The outcomes of our semi-supervised approach for the table class in PubLayNet
are presented in Table 3, employing different fractions of labeled data, while
the remainder serves as unannotated data. Notably, the AP50 metrics are
recorded at 98.5$\%$, 98.7$\%$, and 98.9$\%$ for 10$\%$, 30$\%$, and 50$\%$
labeled data, respectively. For a detailed comparison of results with other
datasets, please consult Figure 4, and for visual representations of the
precision, refer to Figure 6(d).

Table 4 presents a comprehensive evaluation of both supervised and semi-
supervised networks utilizing the ResNet-50 architecture, with specific
emphasis on the table class from the PubLayNet dataset [PubLayNet3]. We
additionally compare a supervised deformable-DETR model, which was trained on
varying proportions of labeled data (10$\%$, 30$\%$, and 50$\%$) from the
PubLayNet table class, with our semi-supervised approach employing the
deformable transformer [Deformable54]. Notably, our methodology produces
comparable outcomes, obviating the necessity for proposal generation and post-
processing procedures like NMS.

###  5.3 ICDAR-19

The efficacy of our table-detecting system is assessed on the Modern Track. A
portion of the dataset was used for table detection in the cTDaR competition
at ICDAR 2019 [icdar19]. Table 5 displays the outcomes of our semi-supervised
methodology used on the ICDAR 2019 dataset, utilizing 10$\%$, 30$\%$, and
50$\%$ of annotated data, while the remaining portion is assigned as
unannotated data.

Dataset | Label-percent | mAP | AP50 | AP75 | ARL
---|---|---|---|---|---
icdar | 10$\%$ | 60.2 | 64.7 | 61.4 | 94.8
30$\%$ | 64.8 | 69.0 | 65.6 | 93.7
50$\%$ | 65.4 | 68.6 | 67.1 | 96.4
Table 5: Comprehensive Evaluation Results for the ICDAR-19 Dataset. Here, mAP
refers to the mean Average Precision within the IoU threshold span of (50:95),
AP50 denotes the Average Precision at an IoU threshold of 0.5, and AP75
represents the Average Precision at an IoU threshold of 0.75.

In addition, we conducted a comparative analysis of our semi-supervised method
using 30$\%$ labeled data against earlier techniques in table detection and
document analysis on TableBank and PubLayNet datasets (refer to Table 6).
Notably, we found that even with only 30$\%$ labeled data, our approach
yielded results comparable to those obtained by previous supervised methods.

Method | Technique | Labels | TableBank | PubLayNet
---|---|---|---|---
CDeC-Net [Agarwal52] | Sup | 100% | 96.5 | 97.8
EDD [zhong2020image] | Sup | 100% | 86.0 | 89.9
TableFormer [nassar2022tableformer] | Sup | 100% | 89.6 | 96.7
CasTabDetectoRS [CasTab45] | Sup | 100% | 95.3 | -
VSR [vsr45] | Sup | 100% | - | 95.9
Shehzadi, T.[shehzadi2023towards] | Semi-sup | 30% | 86.8 | 90.3
Our | Semi-sup | 30% | 96.8 | 98.7
Table 6: We evaluated the performance of our semi-supervised network in
document analysis by comparing it to earlier supervised and semi-supervised
methods. Our approach leverages 30$\%$ of the data with labels, while the
remaining 70$\%$ is unlabeled. The results are presented in terms of mean
average precision (mAP) across a range from 0.5 to 0.95.

##  6 Ablation Study

In this segment, we evaluate the pivotal elements of the design. It’s worth
noting that, unless explicitly stated otherwise, all ablation studies utilize
a ResNet-50 backbone along with a dataset subset consisting of 10$\%$ labeled
images sourced from PubLayNet.

Comparing the utilization of one-to-one and one-to-many matching: We analyze
the predictions based solely on one-to-one or one-to-many matching in Table 7.
One-to-one matching, which associates each ground truth object with the most
similar predicted object, doesn’t cause issues during testing but may lead to
incomplete learning during training. On the other hand, one-to-many matching
allows multiple predicted objects to correspond to a single ground truth
object, enhancing training but potentially causing overlapping predictions
during testing. To balance these factors, our strategy involves using both
one-to-one and one-to-many matching during training and exclusively one-to-one
matching during testing to ensure accurate predictions without NMS problems.

One-to-One | One-to-Many | NMS-free | mAP | AP50 | AP75 | ARL
---|---|---|---|---|---|---
✓ | ✗ | ✓ | 96.5 | 98.1 | 97.1 | 98.9
✗ | ✓ | ✗ | 95.8 | 97.0 | 96.7 | 98.6
✓ | ✓ | ✓ | 96.9 | 98.5 | 97.5 | 98.9
Table 7:  Comparative Analysis of Performance among One-to-One, One-to-Many
Training Schemes, and NMS methods. The best results are highlighted in bold.

The Impact of the Quantity of Learnable Queries: The investigation examines
the effects of manipulating the quantity of queries supplied to the decoder.
Optimal performance is achieved when employing 30 one-to-one queries and 400
one-to-many queries. Deviating from these values leads to a decline in
performance. Table 9 provides a thorough examination of the results derived
from various quantities of object queries. Opting for a lower value might
result in the model’s incapacity to precisely detect certain objects, thereby
adversely affecting its overall performance. Conversely, selecting a large
value may lead to overfitting, causing the model to misclassify certain
regions as objects.

o2o | o2m | mAP | AP50 | AP75 | ARL
---|---|---|---|---|---
10 | 500 | 96.1 | 97.7 | 96.9 | 98.8
30 | 400 | 96.9 | 98.5 | 97.5 | 98.9
300 | 1000 | 96.4 | 97.8 | 97.0 | 99.1
Table 8: Evaluation of Model Performance Based on Varying Learnable Query
Quantities. Optimal results across different settings are highlighted in bold
for clarity and emphasis.

T | mAP | AP50 | AP75
---|---|---|---
0.5 | 96.9 | 98.3 | 97.2
0.6 | 96.2 | 97.6 | 96.0
0.7 | 97.1 | 98.5 | 97.5
0.8 | 96.8 | 98.3 | 97.2
Table 9: Analysis of Pseudo-Labeling Confidence Thresholds and Their Impact on
Performance. The table delineates outcomes across different thresholds, with
the top-performing results accentuated in bold.

Pseudo-Labeling Confidence Threshold: This threshold is a critical parameter
that influences the quality and quantity of pseudo-labels assigned to data
points during semi-supervised learning. Increasing the threshold reduces the
number of examples that meet the criteria, but those that do tend to have
higher quality. Conversely, lowering the threshold leads to more cases meeting
the criteria, but this also raises the likelihood of incorrect results. Table
9 presents a comprehensive summary of the effects observed across several
threshold values, spanning from 0.5 to 0.8. According to the results of our
investigation, it has been decided that the most suitable threshold value is
0.7.

##  7 Conclusion

In conclusion, our study addresses the challenges inherent in DETR-based
object detectors within semi-supervised object detection. We recognized the
limitations of the conventional one-to-one assignment approach, particularly
its inefficiency in handling imprecise pseudo-labels. To overcome these
issues, we introduced a semi-supervised, transformer-based approach for semi-
supervised object detection that operates end-to-end. Our method integrates
one-to-one and one-to-many matching strategies, capitalizing on the strengths
of both techniques.

We refine the pseudo-class labels and pseudo-bounding boxes by employing two
distinctive modules, student and teacher. Through the iterative update process
facilitated by an Exponential Moving-Average (EMA) function, these modules
enhanced one another, resulting in accurate classification and bounding box
predictions. The experimental findings underscored the effectiveness of our
approach, surpassing the performance of supervised models across various
labeling ratios on TableBank and the PubLayNet training data.

Future research will investigate how different proportions of annotated data
impact overall performance. Additionally, we aim to develop models capable of
operating effectively with limited amounts of labeled data. Furthermore, we
aim to extend the semi-supervised learning framework based on transformers to
tackle the task of table structure recognition. This research represents
significant progress in advancing more robust and efficient methodologies for
semi-supervised object detection within document analysis tasks.


Further author information: (Send correspondence to A.L.)
A.L.: E-mail<EMAIL_ADDRESS>Telephone: +91 9031 953805

# A numerical integration scheme for vectorised phase-space of one-dimensional
collision-free, electrostatic systems

Allen Lobo  Department of Physics, Sikkim Manipal Institute of Technology,
Sikkim Manipal University, Sikkim-737136, India.  Vinod Kumar Sayal
Department of Physics, Sikkim Manipal Institute of Technology, Sikkim Manipal
University, Sikkim-737136, India.

###### Abstract

The kinetic analyses of many-particle soft matter often employ many simulation
studies of various physical phenomena which supplement the experimental
limitations or compliment the theoretical findings of the study. Such
simulations are generally conducted by the numerical integration techniques of
the governing equations. In the typical case of collisionless electrostatic
systems such as electrostatic plasmas, the Vlasov-Poisson (VP) equation system
governs the dynamical evolution of the particle phase-space. The one-
dimensional position-velocity (1D-1V) particle phase-space, on the other hand,
is known to exhibit direct analogy with ordinary two-dimensional fluids,
wherein the Vlasov equation resembles the fluid continuity equation of an in-
compressible fluid. On the basis of this fluid-analogy, we present, in this
work, a new numerical integration scheme which treats the 1D-1V phase-space as
a two-dimensional fluid vector space. We then perform and present analyses of
numerical accuracy of this scheme and compare its speed and accuracy with the
well-known finite splitting scheme, which is a standardised technique for the
numerical Vlasov-Poisson integration. Finally, we show some simulation results
of the 1D collisionless electrostatic plasma which highlight the higher speed
and accuracy of the new scheme. This work presents a fast and sufficiently
accurate numerical integration technique of the VP system which can be
directly employed in various simulation studies of many particle systems,
including plasmas.

###### keywords:

Computational Plasma Physics, Numerical Techniques, Numerical Vlasov
Integration, ANI scheme

##  1 INTRODUCTION

The kinetic theory of matter is well-known to be a meticulous approach to
understand various aspects of the dynamics of many-particle systems. This
approach employs a statistical distribution of the mechanical states of
various particles in the system in order to explain the macroscopic phenomena
which occur in it. Plasmas, which are high-energy, partially or completely
ionised states, are not unknown to the kinetic approach – their dynamics have
been assiduously explored using the plasma kinetic theories[1, 2, 3, 4, 5]
which includes certain phenomena that can not be addressed by the magneto-
hydrodynamics approach. In the plasma kinetic theory, the dynamical evolution
of the system is exhibited by the well-known Boltzmann equation. In the hot
plasma case, collisions and inter-particle interactions can be ignored and the
equation modifies into the Vlasov equation[6],

$\frac{\partial f}{\partial t}+\vec{v}\cdot\frac{\partial
f}{\partial\vec{r}}+\vec{a}\cdot\frac{\partial f}{\partial\vec{v}}=0.$ (1)

In the above equation, $f(\vec{r},\vec{v},t)$ represents the probability
distribution function of the plasma species in its position-velocity phase-
space. The terms $\vec{r}$, $\vec{v}$ and $\vec{a}$ represent position,
velocity and acceleration vectors and $t$ represents time, respectively. The
Vlasov equation (1) describes the dynamical evolution of the system in terms
of its phase-space distribution function. This evolution is regulated by the
presence of field interactions in the system, which collectively contribute to
the acceleration term $\vec{a}$. In the case of purely electrostatic
interactions, this acceleration term can be represented by the contribution of
the electric field $\vec{E}$, using the Poisson equation,

$\vec{a}=\frac{q}{m}\vec{E},\quad\vec{\nabla}_{\vec{r}}\cdot\vec{E}=-\nabla_{\vec{r}}^{2}\phi=\frac{1}{\varepsilon_{0}}\sum_{i=1}^{N}q_{i}\int
f_{i}(\vec{r},\vec{v})dv.$ (2)

In the above equation, $\varepsilon_{0}$ is the spatial electrostatic
permittivity, $\phi$ is the electrostatic potential, $q_{i}$ is the charge of
the plasma species and $m$ is the mass of the species. Equations (1) and (2)
collectively describe the Vlasov-Poisson (VP) system, which governs the
dynamical evolution of the phase-space kinetics of the plasma species. The
kinetic approach encompasses various linear and nonlinear phenomena which
occur in plasmas, including wave propagation, instabilities and structure
formations. The Vlasov equation (1) exhibits remarkable similarity, in its
form, to the continuity equation of an ordinary in-compressible fluid. This
similarity has been explored by some authors [7, 8, 9] and it was recently
shown[10] that this analogy extends beyond the Vlasov equation, that is, a
1D-1V phase-space and a two-dimensional conventional fluid surface exhibit
similarities in terms of the evolution dynamics, flow-like behaviour and
deformations due to external fields.

Kinetic analyses of the plasma phenomena are often accompanied by numerical
simulations of the system[11], which, apart from visualising and showcasing
the studied phenomena, also exhibit the dynamics of the field waves which may
not be observable in experimental studies, and also help in validating the
study and its findings. The simulation is generally conducted by numerically
integrating the VP equation system. Various numerical schemes have been
developed and analysed in this regard [12, 13, 14, 15, 16, 17, 18]. One of the
well-known benchmark numerical schemes, developed by Cheng and Knorr
(1976)[19], known popularly as the finite splitting (FS) scheme, has been used
in various numerical investigations of nonlinear processes in plasma
physics[15, 18]. The FS scheme employs a leap-frog integration scheme of the
Vlasov-equation along position and velocity spaces, each at different, half-
time steps. These steps, each marked by asterisks $(^{*})$, are as follows:

$f(x,v_{x},t)\rightarrow f^{*}(x-v_{x}\Delta t/2,t)\rightarrow
f^{**}(x,v_{x}-\Delta tq\vec{E}/m,t)\rightarrow f^{***}(x-v_{x}\Delta
t/2,t)=f(x,v_{x},t+\Delta t).$ (3)

In the above numerical equation, $\Delta t$ represents the time-step value and
$m$ represents the mass of the plasma species. This 3-step integration scheme
of the Vlasov equation employs various interpolation schemes for shifting the
distribution function along respective position and velocity axes. In the
above equation (3), each arrow $(\rightarrow)$ indicates this interpolating
shift of the phase-space distribution function (df). This means that for the
evolution of the phase-space by one time-step, a three-step algorithm should
be implemented, which may prove to be an overburden on computation systems,
when dealing with huge systems. This is due to the dependency of the accuracy
of this scheme on the grid-spacing[20]. Therefore, large-sized systems with
high sensitivity to numerical accuracy might result in longer simulation times
and slower data generation. This limitation serves as the chief motivation of
this work.

In this article, we present a new numerical algorithm based on vectorisation
of the plasma phase-space. By vectorisation, it is connoted that the phase-
space df is observed as a three-dimensional vector field, with independent
evolving components along each axis of the dynamical phase-space: position
$(x)$, velocity $(v_{x})$ and time $(t)$. We then expand the evolution of the
phase-space df separately along each axis and present an agile, numerical
integration scheme which evolves the phase-space by one-time step in one
numerical step, which significantly reduces the computational time of the
evolution equation when compared to the FS scheme. We then study and present
the stability of our code and compare its speed and accuracy with the FS
scheme for some well-known plasma phenomena. This article is planned as
follows: in section 2, the three-dimensional outlook of the time-evolving
dynamical phase-space is discussed and the phase-space is described as a
three-dimensional vector space. Phase-space evolution is then exhibited as a
flow of the phase-space hydrodynamic behavior along the three dimensions
independently. In section 3, using this idea of independent flow-like
evolution of the phase-space, a new numerical scheme is developed, and its
stability and accuracy is discussed. In section 4, the new numerical scheme is
employed to simulate various plasma phenomena, and its numerical accuracy and
speed are tested against the FS scheme. Finally, section 5 is devoted to the
concluding remarks. This article presents a faster and accurate numerical
scheme which can be used in computational set-ups to observe and study various
linear and nonlinear plasma phenomena with much less computational burden than
the existing techniques.

##  2 Evolving phase-space as a three-dimensional vector field

The 1D-1V particle phase-space has often been compared[9, 7] to a two-
dimensional fluid surface, by observing the well-known analogy between the
Vlasov equation and the continuity equation. While the latter represents the
nature of flow and conservation of fluid density, the former states that a
particle in its phase-space moves only along a trajectory where-in its
probability density remains a constant, that is, the total time derivative of
the phase-space df is zero, specifically in the collisionless case. This
analogy, however, lies beyond the Vlasov-Continuity equations, as shown by
Lobo and Sayal[10]. It is shown that the particle phase-space also possesses
its own vorticity field, and experiences a statistical pressure due to non-
uniform probability distributions, and can also form vortices[21, 22, 10] and
other coherent structures. Kinetic hydrodynamic theories can then be developed
in order to study these fluid-analogous structures in the phase-space[10].

The dynamical evolution of the phase-space is described by the time-dependence
of the phase-space df $f(x,v_{x},t)$. While this time-dependence can be both
extrinsic or intrinsically through the space and velocity coordinates, it is
well-known that the evolution of the phase-space df can be described by the
motion of the particle along a phase-space trajectory which conserves this
density, that is,

$\frac{df(x,v_{x},t)}{dt}=0,\quad\text{for collisionless case.}$ (4)

At a certain time, an infinitesimal phase-space volume element $(\delta
V=\tau\delta x\delta v_{x})$ can be assigned a phase-space coordinate set
$(x,\tau v_{x})$, where $\tau$ is a characteristic time of the system. For
plasmas, it is the inverse of the plasma frequency
$\omega_{p}=\sqrt{n_{0}q^{2}/\varepsilon_{0}m}$, where $n_{0}$ is the
unperturbed spatial particle density. Therefore, the spatio-temporal location
of the phase-space volume element becomes $(ct,x,\tau v_{x})$, where $c$ is a
characteristic speed, at which information travels in time, and is closely
related to the Courant-Friedrichs-Lewy (CFL) number [20]. Usually, $c$ is is
taken to be equal to $1$. The temporal-phase-space evolution is then
observable in a three-dimensional phase-space, representing a one-dimensional
system. As an example of the same, we present the case of a simple harmonic
oscillator with unit angular frequency –

$\frac{dv_{x}}{dt}=-x,\qquad\frac{dx}{dt}=v_{x}.$ (5)

Fig, 1 shows the phase-space trajectory of the harmonic oscillator, for a
total of eight periods. At any instant, the mechanical state of the particle
is represented by the coordinate set $(x,\tau v_{x})$.

Figure 1: Phase-space trajectories of a simple harmonic oscillator. A two-
dimensional $(x,\tau v_{x})$ phase-space portrait representing motion of the
particle along the phase-space trajectory. At any given time, the state of the
particle lies in the phase-space trajectory (indicated by the red arrow).
However, the exact position of the particle in its phase-space can not be
found at a given time in this representation.

Fig. 2 showcases the same phase-space, evolving with time in three-dimensions.
While in both cases, the information relayed is the same, the three-
dimensional representation of the phase-space is a better representation of
the dynamically evolving system.

Figure 2: Phase-space trajectories of a simple harmonic oscillator, in a
three-dimensional representation $(ct,x,\tau v_{x})$. The mechanical state of
the particle at a specific time can be easily located in this representation
(red, continuous arrow). The red, dotted arrow represents the unidirectional
flow of time, from the origin $(ct=0,x=0,\tau v_{x}=0)$.

Having described the phase-space in a three-dimensional set-up of position,
velocity and time, the state of a phase-pace volume element $(\delta V)$ can
be described in a three coordinate set $(ct,x,\tau v_{x})$. Therefore, at any
instant, the position vector $\bm{\Gamma}$ of the $\delta V$ element can be
described as:

$\displaystyle\bm{\Gamma}=ct\hat{t}+x\hat{x}+\tau v_{x}\hat{v}_{x}.$ (6)
Differentiating w.r.t. time, we get the velocity field $\bm{u}$ of the volume
element $\delta V$,
$\displaystyle\bm{u}=c\hat{t}+v_{x}\hat{x}+\frac{q}{m}\vec{E}\hat{v}_{x}.$ (7)

The phase-space velocity field $\bm{u}$ defined in equation (7) describes the
evolution of the particle phase-space in terms of a flow-like behaviour of the
phase-space df in phase-space and time. This hydrodynamic behaviour of the
phase-space in a time-independent case was introduced earlier by Berk, Nielson
and Roberts[7] and later by Lobo and Sayal[10]. While the velocity components
along $\hat{x}$ and $\hat{v}_{x}$ directions are known from the system and
system-specific interacting fields[10], the phase-space velocity component
along the time axis $(\hat{t})$ describes the rate of evolution of the
dynamical system, as follows:

For the system evolves for a time $\Delta t$, the phase-space must also evolve
accordingly,

$\displaystyle\allowdisplaybreaks t\xrightarrow{t+\Delta t}t+\Delta t,$ (8)
$\displaystyle x\xrightarrow{t+\Delta t}x+\hat{x}\cdot\bm{u}\Delta t=x+\Delta
x\bigr{|}_{t}^{t+\Delta t},$ (9) $\displaystyle\tau v_{x}\xrightarrow{t+\Delta
t}\tau v_{x}+\tau\hat{v}_{x}\cdot\bm{u}\Delta t=\tau v_{x}+\tau\Delta
v_{x}\bigr{|}_{t}^{t+\Delta t},$ (10) $\displaystyle\therefore f(ct,x,\tau
v_{x})\xrightarrow{t+\Delta t}f(ct+c\Delta t,x+\Delta x\bigr{|}_{t}^{t+\Delta
t},\tau v_{x}+\tau\Delta v_{x}\bigr{|}_{t}^{t+\Delta t}).$ (11) Therefore, $c$
is the rate at which the phase-space fluid evolves along the time axis
$\hat{t}$. If $c$ is less than $1$, it would imply that the phase-space
evolves slower than the information propagating along-with the evolving time.
This means that there would exist an asymmetry in the development of the
phase-space df itself and its coordinates $x$ and $\tau v_{x}$ which develop
along-with time. On the other hand, if $c$ is greater than $1$, it implies
that the speed at which the information associated with the system’s evolution
travels faster than the system itself. In a physical system, due to the idea
of causality, this condition becomes meaningless. However, this condition may
be allowed in a numerically simulated system following an implicit time-
integration scheme. Therefore, $\displaystyle c\geq 1.$ (12)

The restriction imposed on the value of $c$ in equation (12) is also in
agreement to the CFL condition[20], and $c$ itself is therefore identical to
the Courant number.

Having described the dynamically evolving phase-space as a three-dimensional
vector space, the Vlasov equation can be rewritten in a vector form as the
divergence of the phase-space flux density $\bm{\mathcal{J}}$, which can be
defined as the phase-space velocity field $\bm{u}$ times the phase-space df
$f(ct,x,\tau v_{x})$.

$\displaystyle\bm{\nabla}=\frac{1}{c}\frac{\partial}{\partial
t}\hat{t}+\frac{\partial}{\partial
x}\hat{x}+\frac{1}{\tau}\frac{\partial}{\partial v_{x}}\hat{v}_{x}.$ (13)
$\displaystyle\bm{\mathcal{J}}=f\bm{u}=cf\hat{t}+v_{x}f\hat{x}-\frac{q\tau}{m}\frac{\partial\phi}{\partial
x}f\hat{v}_{x},$ (14) where, the electric field $\vec{E}$ is substituted by
the electric potential gradient $(-\partial\phi/\partial x)$. The Vlasov
equation then takes the form,
$\displaystyle\bm{\nabla}\cdot\bm{\mathcal{J}}=0.$ (15)

Equation (15) presents the Vlasov equation in a non-divergent, differential
form. It is well-known from conventional vector calculus that equation (15)
can also be represented in an integral form as follows:

$\oiint\bm{\mathcal{J}}\cdot d\bm{S}=0.$ (16)

In the above equation, $d\bm{S}$ represents surface elements in the three
dimensional dynamical phase-space, and is equal to

$d\bm{S}=\tau dxdv_{x}\hat{t}+c\tau dtdv_{x}\hat{x}+cdtdx\hat{v}_{x}.$ (17)

The integral form of the Vlasov equation (16) presents the well-known[10]
solenoidal nature of the phase-space flux density field. It states that as the
system evolves in time, there is no net flux of the phase-space fluid across a
closed dynamical phase-space volume $\int_{\delta t}\oiint\delta V\cdot
c\delta t$.

Therefore. the dynamical evolution of the particle phase-space can be
conveniently observed as a kinetic-hydrodynamic flow along three independent
coordinates, of time, length and speed, in a three-dimensional vector space.
Figs. 3 and 4 portray the temporal evolution of the dynamical phase-space in
terms of the phase-space df and the magnitude of the flux density field
$\bm{\mathcal{J}}$. In both the cases, formation of a phase-space vortex from
two-stream instability has been simulated using Vlasov kinetic simulation
employing the FS scheme. The hydrodynamic nature of the phase-space evolution
shown in the figures portray the phase-space evolution as a fluid flowing
along the positive time axis, deforming and evolving in accordance with the
physics of interacting electric field and the plasma[16].

Figure 3: Evolution of the plasma dynamic phase-space. Phase-space df portrait
shown during the evolution of an initial two-stream plasma into a phase-space
vortex. Color gradient representing magnitude of the phase-space density
$f(ct,x,\tau v_{x}$ and its evolution with time.

The vectorisation of the dynamical phase-space and accordingly the Vlasov
equation presented in this section dispenses that the evolution of the phase-
space df occurs independently along each direction of the dynamical phase-
space. This orthogonal evolution of the phase-space implies that the evolution
of the phase-space df (and the flux density $\bm{\mathcal{J}}$) can be
simulated by individually evolving it along each direction, instead of
shifting the evolution along two axis-sets, as has been done classically[19,
18]. Using this outcome, we devise and present a new numerical scheme for the
dynamical evolution of the plasma phase-space in the next section.

Figure 4: Evolution of the plasma dynamic phase-space flux density magnitude
$|\bm{\mathcal{J}}|$ from an initial two-stream plasma into a phase-space
vortex, with magnitudes represented bu the color gradient.

##  3 A Numerical Integration scheme for the Vlasov Equation and its stability
analysis

As discussed above, the vectorised phase-space makes it convenient to
numerically evolve the phase-space df along individual coordinates
independently of each other. This is highlighted by the vector form of the
Vlasov equation (15). Expanding the differential to present the classical form
of the VP system,

$\displaystyle\allowdisplaybreaks\bm{\nabla}\cdot\bm{\mathcal{J}}=\frac{\partial
f}{\partial t}+v_{x}\frac{\partial f}{\partial
x}-\frac{q}{m}\frac{\partial\phi}{\partial x}\frac{\partial f}{\partial
v_{x}}=0.$ (18) In the above equation, the terms represent independent
gradients along each axis of the dynamical phase-space. Discretizing the
equation, we get – $\displaystyle\frac{\Delta f}{\Delta t}\Bigr{|}_{\Delta
x=0,\Delta v_{x}=0}+v_{x}\frac{\Delta f}{\Delta x}\Bigr{|}_{\Delta t=0,\Delta
v_{x}=0}-\frac{q}{m}\frac{\partial\phi}{\partial x}\frac{\Delta f}{\Delta
v_{x}}\Bigr{|}_{\Delta x=0,\Delta t=0}=0.$ (19) Evolving the dynamical phase-
space coordinates in time $\Delta t$ as shown in equations (8), (9) and (10),
we get – $\displaystyle\frac{\Delta f}{\Delta t}\Delta t+(v_{x}\Delta
t)\frac{\Delta f}{\Delta x}+\left(-\frac{q}{m}\frac{\partial\phi}{\partial
x}\Delta t\right)\frac{\Delta f}{\Delta v_{x}}=\Delta f\Bigr{|}_{t}^{t+\Delta
t}+\frac{\Delta f}{\Delta x}\Delta x\Bigr{|}_{t}^{t+\Delta t}+\frac{\Delta
f}{\Delta v_{x}}\Delta v_{x}\Bigr{|}_{t}^{t+\Delta t}=0.$ (20)
$\displaystyle\Rightarrow\Delta f\Bigr{|}_{t}^{t+\Delta t}+\Delta
f\Bigr{|}_{x}^{x+\Delta x}+\Delta f\Bigr{|}_{v_{x}}^{v_{x}+\Delta v_{x}}=0.$
(21)

Equation (21) presents the complete evolution of the phase-space df in a time
interval $\Delta t$ in terms of a set of flow-like progressions along each
axis of the dynamical phase-space. This evolution, as shown in equation (21),
occurs independently along each axis, as discussed above. It can also be seen
that the equation represents a time-explicit integration scheme with one
unknown. Using indices $n$, $j$ and $k$ for descretized steps of time
$(c\Delta t)$, position $(\Delta x)$ and velocity $(\tau\Delta v_{x})$,
equation (21) obtains the form –

$f^{n+1}_{j,k}=f^{n}-\Delta_{x}f^{n}-\Delta_{v_{x}}f^{n}.$ (22)

Equation (22) presents a numerical algorithm to evolve the phase-space df in
time. While the temporal evolution of the phase-space df itself is presented
in terms of its flow along the position and velocity axes, the evolution of
the position-velocity coordinates themselves can be determined directly using
equations (9) and (10). The phase-space density at these points can be
determined by numerical interpolation. The complete algorithm can be described
in the following step, in correspondence to the FS scheme[19] shown in
equation (3):

$\displaystyle f(x,v_{x},t+\Delta t)=f(x,v_{x},t)+[f(x\rightarrow
x-v_{x}\Delta t,v_{x},t)-f(x,v_{x},t)]+[f(x,v_{x}-qE(x)\Delta
t/m,t)-f(x,v_{x},t)],$ (23) $\displaystyle\Rightarrow f(x,v_{x},t+\Delta
t)=f(x\rightarrow x-v_{x}\Delta t,v_{x},t)+f(x,v_{x}-qE(x)\Delta
t/m,t)-f(x,v_{x},t).$ (24)

In the above equations (19 \- 24), $c$ and $\tau$ are normalised to be $1$. As
can be seen, the numerical scheme described by equation (24) presents a one-
step algorithm with two interpolation shifts, which can be performed
simultaneously and are independent of each other. It can be seen that the
numerical accuracy of the scheme depends on the accuracy of the chosen
interpolation technique, and does not decrease beyond it since consecutive
interpolations are not performed on the phase-space df. if $\Delta\epsilon$ is
the numerical error produced in the phase-space df due to each interpolation,
the FF scheme compounds this error by the the successive interpolations, as
shown in equation (3).

$f^{n}\rightarrow f^{*}(=f+\Delta\epsilon)\rightarrow
f^{**}(=f+\Delta\epsilon+\Delta\epsilon)\rightarrow
f^{***}(=f+\Delta\epsilon+\Delta\epsilon+\Delta\epsilon)=f^{n+1}.$ (25)

Therefore the numerical error at one time-step becomes $3\Delta\epsilon$.
However, this compounding collective error does not occur in the ANI scheme,
which does not utilise the successive interpolation algorithm.

$f^{n+1}=f^{n}_{j^{\prime},k}+f^{n}_{j,k^{\prime}}(=f+\Delta\epsilon)-f^{n}_{j,k}.$
(26)

This agile numerical integration (ANI) scheme is therefore an accuracy
preserving, single-step semi-Lagrangian scheme, with comparatively more
accuracy, even though the order of the accuracy remains same. The stability of
the scheme can be analysed using the well-known Von-Neumann technique[23, 24].
Representing the numerical error at a phase-space grid $(j,k)$ as
$\epsilon_{j,k}^{n}$ for the $n^{\text{th}}$ time-step, and expressing it in
terms of a Fourier series, we get

$\displaystyle\epsilon(t,x,v_{x})=\sum_{\sigma,\eta}A_{\sigma,\eta}(t)\exp(ik_{\sigma}x)\cdot\exp(ih_{\eta}v_{x}).$
(27) $\displaystyle\text{Here,}\quad
k_{\sigma}=\frac{\pi\sigma}{L}\quad\text{and}\quad\sigma=-\frac{L}{\Delta
x},-\frac{L}{\Delta x}+1,\quad.\quad.\quad 0\quad.\quad.\quad\frac{L}{\Delta
x}-1,\frac{L}{\Delta x}.$ (28) $\displaystyle\text{Similarly},\quad
h_{\eta}=\frac{\pi\eta}{2v_{0}}\quad\text{and}\quad\eta=-\frac{2v_{0}}{\Delta
v_{x}},-\frac{2v_{0}}{\Delta v_{x}}+1,\quad.\quad.\quad
0\quad.\quad.\quad\frac{2v_{0}}{\Delta v_{x}}-1,\frac{2v_{0}}{\Delta v_{x}}.$
(29) Here, $L$ and $v_{0}$ represent the length and velocity amplitude of the
system, and $A_{\sigma,\eta}$ is the amplitude of the error function.
Therefore, we get –
$\displaystyle\epsilon_{j,k}^{n+1}=A_{\sigma,\eta}(t+\Delta
t)\exp(ik_{\sigma}x)\cdot\exp(ih_{\eta}v_{x}),\quad\text{where
}t^{n+1}=t^{n}+\Delta t.$ (30)
$\displaystyle\epsilon_{j^{\prime},k}^{n}=A_{\sigma,\eta}(t)\exp(ik_{\sigma}(x-v_{x}\Delta
t))\cdot\exp(ih_{\eta}v_{x}),\quad\text{where
}x_{j^{\prime}}=x_{j}-v_{x}\Delta t.$ (31)
$\displaystyle\epsilon_{j,k^{\prime}}^{n}=A_{\sigma,\eta}(t)\exp(ik_{\sigma}x)\cdot\exp(ih_{\eta}(v_{x}-q\vec{E}(x_{j})\Delta
t/m)),\quad\text{where
}v_{x_{k^{\prime}}}=v_{x_{k}}-\frac{q}{m}\vec{E}(x_{j})\Delta t.$ (32)

Inserting the above equations into the numerical scheme presented in equation
(24), we get –

$A_{\sigma,\eta}(t+\Delta
t)\exp(ik_{\sigma}x)\cdot\exp(ih_{\eta}v_{x})=A_{\sigma,\eta}(t)\exp(ik_{\sigma}(x-v_{x}\Delta
t))\cdot\exp(ih_{\eta}v_{x})\\\
+A_{\sigma,\eta}(t)\exp(ik_{\sigma}x)\cdot\exp(ih_{\eta}(v_{x}-q\vec{E}(x)\Delta
t/m))-A_{\sigma,\eta}(t)\exp(ik_{\sigma}x)\cdot\exp(ih_{\eta}v_{x}).$ (33)

Equation (33) simplifies into –

$\frac{A_{\sigma,\eta}(t+\Delta
t)}{A_{\sigma,\eta}(t)}=\exp({-ik_{\sigma}v_{x}\Delta
t})+\exp(-ih_{\eta}q\vec{E}(x)\Delta t/m)-1.$ (34)

The ratio derived in equation (33) describes the growth rate of numerical
error at each time-integration and must be less than or equal to 1 for the
numerical scheme to be stable. Expanding equation (33) be introducing
$\theta_{1}=k_{\sigma}v_{x}\Delta t$ and $\theta_{2}=h_{\eta}q\vec{E}(x)\Delta
t/m$,

$\displaystyle\exp(-i\theta_{1})+\exp(-i\theta_{2})-1\leq
1\Rightarrow\exp(-i\theta_{1})+\exp(-i\theta_{2})\leq 2.$ (35) Using the
identity $e^{-i\theta}=\cos(\theta)-i\sin(\theta)$, we get –
$\displaystyle\cos\theta_{1}+\cos\theta_{2}-i\left(\sin\theta_{1}+\sin\theta_{2}\right)\leq
2\Rightarrow\cos\left(\frac{\theta_{1}+\theta_{2}}{2}\right)\cos\left(\frac{\theta_{1}-\theta_{2}}{2}\right)-i\sin\left(\frac{\theta_{1}+\theta_{2}}{2}\right)\cos\left(\frac{\theta_{1}-\theta_{2}}{2}\right)\leq
1.$ (36) Taking $(\theta_{1}+\theta_{2})/{2}=\Theta_{A}$ and
$(\theta_{1}-\theta_{2})/{2}=\Theta_{B}$ and finding the squared-modulus on
both sides, we get –
$\displaystyle\cos^{2}\Theta_{A}\cos^{2}\Theta_{B}+\sin^{2}\Theta_{A}\cos^{2}\Theta_{B}=\cos^{2}\Theta_{B}(\cos^{2}\Theta_{A}+\sin^{2}\Theta_{A})$
(37) $\displaystyle=\cos^{2}\Theta_{B}\leq 1,\quad\text{which is true for all
values of }\Theta_{B}.$ (38)

Equation (38) therefore presents the unconditional stability of the numerical
algorithm presented in this work in equation (24). In this section, a semi-
Lagrangian, one-step ANI scheme for the time-integration of the Vlasov
Equation has been introduced and its stability is discussed. As stated
earlier, the numerical accuracy of this scheme is dependent on the
interpolation technique. When coupled with the Poisson equation, the numerical
accuracy then also depends on the choice of the Poisson solver. In the next
section, we employ this ANI scheme in order to simulate some well-known
phenomena in plasma physics, and compare the results obtained and the speed of
the simulation with the FS scheme.

##  4 Numerical Results and Discussions

In this part of our work, we employ the one-step, ANI scheme presented in
section 3 equation (24) to numerically simulate some well-known collisionless
plasma phenomena. This includes linear wave propagation and its Landau
damping, and the plasma echo phenomena. These results are compared with the
simulation results of the FS scheme[19]. We test the numerical accuracy and
speed of our reduced scheme and present the results. In order to perform the
numerical simulation, we employ the cubic spline interpolation in our scheme,
which is well-known to exhibit higher-order accuracy up to the $4^{th}$ order.
We employ a phase-space grid of size $M\times N$, where $M$ represents the
number of grid points in the position $(x)$ space and $N$ represents the
number of grid points in the velocity $(\tau v_{x})$ space. We take $M=512$
and $N=512$, therefore dealing with a total of $2,62,144$ grid points. This
ensures a small grid-size, leading to significantly reduced numerical
dissipation. We specifically deal with the case of electron plasma waves, and
simulate the electron phase-space. We normalise the position with electron
Debye length $\lambda_{De}=\sqrt{\varepsilon_{0}K_{B}T_{e}/n_{0}e^{2}}$, with
$K_{B}$ and $T$ being the Boltzmann’s constant and electron temperature,
respectively. The electron charge is represented by $e$. Similarly, the
particle velocity is normalised with electron thermal velocity
$v_{Te}=\sqrt{2K_{B}T_{e}/m}$ and time with inverse plasma frequency
$\omega_{pe}=\sqrt{n_{0}e^{2}/\varepsilon_{0}m}$. Singly charged ions are
assumed to remain stationary and form a uniform, neutralising background. The
resultant VP system is as follows:

$\displaystyle\frac{\partial f}{\partial t}+v_{x}\frac{\partial f}{\partial
x}-E(x)\frac{\partial f}{\partial v_{x}}=0,$ (39) $\displaystyle\frac{\partial
E}{\partial x}=-\frac{\partial^{2}\phi}{\partial
x^{2}}=1-\int_{-\infty}^{\infty}fdv_{x}.$ (40)

At each time-step, the Vlasov equation is numerically integrated and the
Poisson equation is solved after integrating the modified phase-space df. We
employ the inverse fast-Fourier technique in order to solve the Poisson
equation, and use periodic boundary conditions. The following CFL conditions
are employed for additional stability at each time-integration step[20] –

$\displaystyle dt=\frac{\Delta x}{v_{0}}\quad\text{or}\quad dt=\frac{\Delta
v_{x}}{E_{0}},\quad\text{whichever is smaller}.$ (41)

Here, $v_{0}$ is the maximum velocity and $E_{0}$ is the absolute electric
field amplitude. We start with the numerical simulation of linear Landau
damping[2] of an electron plasma wave, which is a collisionless decay in the
wave energy occurring due to transfer of the wave energy to the resonant
particles. For the same, the initial distribution of the electrons in phase-
space is perturbed by adding an electron distribution as follows:

$f(x,v_{x},t=0)=\frac{1}{\sqrt{\pi}}\exp(-v_{x}^{2})\cdot(1+A\cos(\kappa x)),$
(42)

where $A$ is the perturbation amplitude and $\kappa$ is the wave number. We
use a small amplitude of $A=0.01$ in normalised units in order to restrict the
damping in the linear domain. We choose wave number $\kappa$ as $0.5$. The
values of observed wave angular frequency $\omega$ and damping factor
$\gamma$, along-with the time of the simulation in both finite-splitting
scheme case and the ANI scheme are reported in table 1. The observation has
been shown in fig. 5. For the simulation, the same computation set-up has been
used for both numerical schemes and the time of the complete code-run has been
calculated digitally, accurate to the third decimal. Three runs of the codes
have been used and the average of the simulation times of the three trials has
been reported for each case.

Table 1: Data of numerical simulation of linear Landau damping of electron plasma wave. Comparison of numerical simulation schemes in terms of accuracy and speed for wave angular frequency $\omega$ and damping factor $\gamma$ are reported for $\kappa=0.5$. | $\kappa$ |  $\omega$ (error) |  $\gamma$ (error) | simulation time (seconds)
---|---|---|---|---
Theoretical | 0.5 | 1.4156 | -0.153359 | –
FS scheme | 0.5 | 1.4159 (0.02119%) | -0.152356 (-0.653%) | 862.601
ANI scheme | 0.5 | 1.4158 (0.01412%) | -0.153264 (-0.061%) | 469.441
Figure 5: Linear Landau damping of electron Langmuir wave. (Wave energy vs.
time plot. Graph representing damping simulation using the ANI scheme (black,
continuous curve) and the finite-splitting scheme (green, dashed curve).
Damping factor calculated by the slope of the red line.

From the results presented in table 1, it is clear that the ANI scheme
developed herein reduces the computational time significantly, to almost half
when compared with the FS scheme. The accuracy of the two schemes are
comparable to each other. For the nonlinear Landau damping case, the wave
perturbation amplitude $A$ is increased to $0.5$ to shift the interaction into
nonlinearity. This is the transfer of energy initially from the wave to the
slower resonant particles, which then accelerate beyond the wave phase-speed.
This is followed by the transfer of energy from these higher speed particles
back to the wave, thus presenting an oscillation of the wave energy and
formation of phase-space df troughs centered at the wave phase-speed. Using
$\kappa=0.5$, we simulate the electron phase-space. Table 2 and fig. 6
presents the results. The wave energy can be seen to decay initially with a
decay factor $\gamma=-0.575377$ and then grow with a growth factor
$\gamma=0.180057$, to slowly reach a saturation. This is a typical behaviour
of nonlinear collisionless damping and has been reported in many previous
works [19, 3]. These results are in agreement with the data reported by Gazdag
[25] and later by Cheng and Knorr[19], with errors of $\leq 0.236\%$. Fig. 7
shows the formation of phase-space holes due to the nonlinear trapping of
electrons in the nonlinear Landau damping phenomena. These are regions of
reduced particle densities in the phase-space, travelling in either directions
with the wave phase-speeds.

Table 2: Data of numerical simulation of nonlinear Landau damping of electron plasma wave, reported for $\kappa=0.5$ and $A=0.5$. Decay and growth rates of the electrostatic wave amplitude reported, along-with the simulation times of the finite-splitting and ANI scheme. $\kappa$ |  $\gamma$ (decay) |  $\gamma$ (growth) | simulation time FS (seconds) | simulation time ANI (seconds)
---|---|---|---|---
0.5 | -0.287689 | 0.0900289 | 783.298 | 404.373
Figure 6: Nonlinear Landau damping of Langmuir waves, with amplitude $A=0.5$
and wave-number $\kappa=0.5$. Wave damping rate and growth rate shown.
Simulation performed using both ANI (red, dashed curve) and FS (green,
continuous curve) schemes.

We next present the numerical simulation of the temporal plasma echo
phenomena[26, 27] using the ANI scheme. Plasma echo is well-known to be a
highly nonlinear phenomena which is quite difficult to simulate numerically
due to its dependence on the numerical stability[28]. This wave-wave
interaction phenomena occurs due to the persistence of the electron velocity
df oscillations even after the Landau damping of the waves’ macroscopic
properties. Thus, two completely damped waves exhibit resonant interactions
resulting in higher order oscillations after some measurable time intervals.

Figure 7: Phase-space density portrait exhibiting presence of electron phase-
space holes due to nonlinear Landau damping. Simulation performed using the
ANI scheme with $A=\kappa=0.5$. Colour gradient representing magnitude of the
phase-space df.

This time interval depends on the wave-nature of the two waves. For initial
perturbations of the form represented in equation (42), presented in the
plasma after a time interval of $\mathcal{T}$, having wave-numbers
$\kappa_{1}$ and $\kappa_{2}$, harmonic modes $\alpha_{1}$ and $\alpha_{2}$
and amplitudes $A_{1}$ and $A_{2}$ respectively, the echo occurs at a time
interval $t_{echo}$ after the first perturbation, such that

$t_{echo}=\frac{\alpha_{2}\kappa_{2}}{\alpha_{2}\kappa_{2}-\alpha_{1}\kappa_{1}}\mathcal{T}.$
(43)

For the simulation of the plasma temporal echo, the velocity grid spacing of
the phase-space should be small enough for recurrence effect to occur much
later than the echo time. The recurrence time $T_{R}$ is equal to
$2\pi/k\Delta v_{x}$. The plasma echo phenomena is also demonstrated to be
dependent on the nature of the damping of the waves[28]. In case of nonlinear
damping, due to particle trapping, the initial oscillations produced by the
perturbations undergo phase-mixing, resulting in loss of initial information
of the waves. This causes plasma echoes to disappear as the nonlinearity
increases in the damping nature. Hence, small amplitude perturbations must be
used.

For the simulation of temporal plasma echo, we introduce two perturbations
similar to the form represented by equation (42). We use a small amplitude of
$A_{1}=A_{2}=0.005$ in order to restrict the damping to the linear case. We
use the waves of wave numbers $\kappa_{1}=0.5$ and $\kappa_{2}=1.0$ and
introduce them at a time interval $\mathcal{T}$ of $50.0\omega_{pe}^{-1}$, at
their natural (first) harmonic cases. From equation (43), the predicted echo
occurrence time $t_{echo}=100.0\omega_{pe}^{-1}$. Fig. 8 presents the
numerical simulation of the temporal plasma echo using the ANI scheme. It can
be seen that in the numerical simulation, the plasma temporal echo occurs at
time $t_{echo}=100.918\omega_{pe}^{-1}$ (represented by the peak of the wave
energy).

Figure 8: A temporal plasma echo. Wave amplitudes $A_{1}=A_{2}=0.005$ and
$\kappa_{2}=2\kappa_{1}=1.0$ (both in $1^{st}$ harmonics). Echo occurring at
time $t_{echo}=100.918\omega_{pe}^{-1}$, as against the theoretical prediction
of $t_{echo}=100.918\omega_{pe}^{-1}$ (error $=0.918\%$.

The accuracy of the ANI scheme is evident with an error in the echo time of
$0.918\%$, as compared to the FS scheme which presents an error of $1.732\%$.
The simulation time using ANI scheme is found to be $2023.415$ sec., as
compared to the FS scheme ($3940.127$ sec). The later peaks formed in the
simulation, as visible fig. 8 occur due to the numerical recurrence
phenomena[28].

In this section, the accuracy, speed and applicability of the ANI scheme has
been established in order to numerically simulate various linear and nonlinear
plasma phenomena using the kinetic Vlasov simulation. In the next section, we
conclude our work by briefly discussing the work presented in this article.

##  5 Conclusion

In this article, the phase-space of collisionless electrostatic one-
dimensional plasmas has been presented and discussed as a three-dimensional
vector space by observing its evolution in time. The phase-space evolution is
then showcased as a hydrodynamic flow described by the phase-space velocity
field, defined in equation (7) and the flux density vector field, defined in
equation (14). This flow-like evolution presents an intuitive, though
unconventional outlook of the phase-space dynamics of the system. Using this
approach, the formation of coherent kinetic structures and phase-space
contortions can be analysed as fluid-analogous behaviour.

The vectorisation of the dynamic phase-space also permits the development of a
numerical scheme for the Vlasov equation integration. This new agile numerical
integration (ANI) scheme, as developed in this work, is an accurate and stable
numerical scheme for the Vlasov equation and reduces the computation burden to
almost half, when compared to the finite splitting scheme[19], which is a
benchmark numerical technique for the kinetic Vlasov simulation. It has been
shown in this work that this new ANI scheme reduces the computational times to
almost half, when compared with the finite splitting scheme. It has high
accuracy which is comparable to the finite splitting scheme, and is capable to
simulate both linear and nonlinear phenomena occurring in the plasma.

The ANI scheme can be employed in various numerical simulations of the plasma
kinetic theory to study collisionless plasma phenomena. Due to its reduced
computation load, it can be utilised for faster simulations in less-powered
computation set-ups. In this work, we have employed the cubic spline
interpolation technique in the numerical simulations, which has a high
accuracy. However, other interpolation techniques can be employed to further
enhance the speed or accuracy of the numerical simulation. The scheme exhibits
unconditional stability, unlike the conservative Eulerian scheme which depends
on the spatio-temporal grid sizes. Hence, the numerical scheme presents a fast
and accurate numerical integration technique of the Vlasov integration without
compromising on the stability of the simulation.

## References

  * [1] Van Kampen, N. G., “On the theory of stationary waves in plasmas,” Physica 21(6-10), 949–963 (1955).
  * [2] Landau, L. D., “On the vibrations of the electronic plasma,” J. Phys. (USSR) 10, 25–34 (1946).
  * [3] Manfredi, G., “Long-time behavior of nonlinear landau damping,” Physical Review Letters 79, 2815–2818 (10 1997).
  * [4] Brodin, G., “Nonlinear Landau Damping,” Physical Review Letters 78, 1263 (2 1997).
  * [5] Sayal, V. K. and Sharma, S. R., “Electron kinetic effects on an ion-acoustic rarefactive soliton in a two-electron temperature plasma,” Physics Letters A 151(1-2) (1990).
  * [6] Vlasov, A. A., “On vibration properties of electron gas,” Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki 8(3), 291 (1938).
  * [7] Berk, H. L., Nielsen, C. E., and Roberts, K. V., “Phase Space Hydrodynamics of Equivalent Nonlinear Systems: Experimental and Computational Observations,” The Physics of Fluids 13, 980–995 (4 1970).
  * [8] Schamel, H., “Electron holes, ion holes and double layers: Electrostatic phase space structures in theory and experiment,” Physics reports 140(3), 161–191 (1986).
  * [9] Schamel, H., “Cnoidal electron hole propagation: Trapping, the forgotten nonlinearity in plasma and fluid dynamics,” Physics of Plasmas 19(2), 20501 (2012).
  * [10] Lobo, A. and Sayal, V. K., “Theory of electron and ion holes as vortices in the phase-space of collision-less plasmas - NASA/ADS,” (2023).
  * [11] Forslund, D. W., “Fundamentals of plasma simulation,” Space Science Reviews 42, 3–16 (10 1985).
  * [12] Watanabe, T. H., Sugama, H., and Sato, T., “A Nondissipative Simulation Method for the Drift Kinetic Equation,” Journal of the Physical Society of Japan 70, 3565–3576 (11 2001).
  * [13] Aminmansoor, F. and Abbasi, H., “Hybrid (Vlasov-Fluid) simulation of ion-acoustic soliton chain formation and validity of Korteweg de-Vries model,” Physics of Plasmas 22(8) (2015).
  * [14] Turikov, V. A., “Computer simulation of the formation of langmuir solitons and holes in a cylindrical magnetized plasma column,” tech. rep., Risoe national lab. (1978).
  * [15] Filbet, F., Sonnendrücker, E., and Bertrand, P., “Conservative numerical schemes for the Vlasov equation,” Journal of Computational Physics 172, 166–187 (9 2001).
  * [16] Morse, R. L. and Nielson, C. W., “One-, two-, and three-dimensional numerical simulation of two-beam plasmas,” Physical Review Letters 23(19), 1087 (1969).
  * [17] Colombi, S. and Touma, J., “Vlasov-Poisson: The waterbag method revisited,” Communications in Nonlinear Science and Numerical Simulation 13(1 SPEC. ISS.) (2008).
  * [18] Filbet, F. and Sonnendrücker, E., “Comparison of Eulerian Vlasov solvers,” Computer Physics Communications 150, 247–266 (2 2003).
  * [19] Cheng, C. Z. and Knorr, G., “The integration of the vlasov equation in configuration space,” Journal of Computational Physics 22, 330–351 (11 1976).
  * [20] Courant, R., Friedrichs, K., and Lewy, H., “Über die partiellen Differenzengleichungen der mathematischen Physik,” Mathematische Annalen 100, 32–74 (12 1928).
  * [21] Guio, P., Børve, S., Daldorff, L. K. S., Lynov, J. P., Michelsen, P., Pécseli, H. L., Rasmussen, J. J., Saeki, K., and Trulsen, J., “Nonlinear Processes in Geophysics Phase space vortices in collisionless plasmas,” Nonlinear Processes in Geophysics 10, 75–86 (2003).
  * [22] Hutchinson, I. H., “Electron holes in phase space: What they are and why they matter,” Physics of Plasmas 24, 055601 (5 2017).
  * [23] Crank, J. and Nicolson, P., “A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type,” Advances in Computational Mathematics 6(1), 207–226 (1996).
  * [24] Charney, J. G., Fjörtoft, R., and Neumann, J. v., “Numerical Integration of the Barotropic Vorticity Equation,” Tellus A: Dynamic Meteorology and Oceanography 2, 237–254 (1 1950).
  * [25] Gazdag, J., “Numerical solution of the Vlasov equation with the accurate space derivative method,” Journal of Computational Physics 19(1) (1975).
  * [26] Gould, R. W., O’Neil, T. M., and Malmberg, J. H., “Plasma Wave Echo,” Physical Review Letters 19, 219–222 (7 1967).
  * [27] Malmberg, J. H., Wharton, C. B., Gould, R. W., and O’Neil, T. M., “Plasma wave echo experiment,” Physical Review Letters 20(3) (1968).
  * [28] Hou, Y. W., Ma, Z. W., and Yu, M. Y., “The plasma wave echo revisited,” Physics of Plasmas 18, 012108 (1 2011).