# Converting Epics/Stories into Pseudocode using Transformers

Gaurav Kolhatkar  SCTR’s Pune Institute of Computer Technology
Pune, India
<EMAIL_ADDRESS>Akshit Madan  Nidhi Kowtal  SCTR’s Pune Institute of
Computer Technology
Pune, India
<EMAIL_ADDRESS>SCTR’s Pune Institute of Computer Technology
Pune, India
<EMAIL_ADDRESS>Satyajit Roy  Sheetal Sonawane  SCTR’s Pune Institute
of Computer Technology
Pune, India
<EMAIL_ADDRESS>Associate Professor
SCTR’s Pune Institute of Computer Technology
Pune, India
<EMAIL_ADDRESS>
###### Abstract

The conversion of user epics or stories into their appropriate representation
in pseudocode or code is a time-consuming task, which can take up a large
portion of the time in an industrial project. With this research paper, we aim
to present a methodology to generate pseudocode from a given agile user story
of small functionalities so as to reduce the overall time spent on the
industrial project. Pseudocode is a programming language agnostic
representation of the steps involved in a computer program, which can be
easily converted into any programming language. Leveraging the potential of
Natural Language Processing, we want to simplify the development process in
organizations that use the Agile Model of Software Development. We present a
methodology to convert a problem described in the English language into
pseudocode. This methodology divides the Text to Pseudocode conversion task
into two stages or subtasks, each of which is treated like an individual
machine translation task. Stage 1 is Text to Code Conversion and Stage 2 is
Code to Pseudocode Conversion. We find that the CodeT5 model gives the best
results in terms of BLEU score when trained separately on the two subtasks
mentioned above. BLEU score is a metric that is used to measure the similarity
between a machine-translated text and a set of reference translations.

###### Index Terms:

Text to code generation, Code to Pseudocode generation, Transformers

##  I Introduction

Efficiency of work is of the highest importance in modern organizations and
businesses. A majority of the workplaces today use the Agile Model for
software development. Agile is a software development approach based on
iterative development, wherein tasks are divided into smaller iterations or
sprints. In Agile project management tools such as Jira are used to document
the user requirements in the form of epics or user stories. Developers need to
understand these requirements and write code for the same. However, a
significant amount of development time and efforts can be saved by automating
the process of code/pseudo code generation, especially for simple or
repetitive problems that have been solved before. The motivation of our
research paper is to simplify the work of developers so that they can focus on
more complex tasks and in the process, to optimize the software development
lifecycle.

Jira is a software application used for issue tracking and project management.
It is widely used by agile development teams to track bugs, stories, epics,
and other tasks. Epics are large bodies of work that can be broken down into a
number of smaller tasks (called stories). Stories, also called “user stories,”
are short requirements or requests written from the perspective of an end
user. Our aim is to convert epics/stories to pseudo code.

Despite the advantages of the Agile Model, developing software may still be a
difficult and drawn-out process, especially when it comes to translating user
requirements into functional code. Developers often translate user epics or
stories manually into code as part of this process, which can take a lot of
time and work.

Our study intends to investigate the potential of utilising machine learning
methods and natural language processing to automate the process of generating
code and pseudocode from user stories in order to address this difficulty. By
doing this, we hope to streamline developers’ tasks, enhance the software
development lifecycle, and boost the effectiveness of the entire industrial
project.

Recent advancements in the field of natural language processing have made it
possible to automate a variety of formerly manual operations. Recent
developments in deep learning, in particular, have made it possible to create
sophisticated natural language models that can extract context and meaning
from text input. By utilising these models, we can quickly and efficiently
create code or pseudocode from user stories, relieving the pressure on
developers.

Our study will examine a variety of currently used methodologies, including as
rule-based systems, statistical models, and deep learning techniques, for
producing code or pseudocode from user stories. We will assess the benefits
and drawbacks of each methodology and suggest a novel strategy that makes the
most of their advantages.

Overall, our research paper aims to make software development more efficient
and effective by automating the process of code/pseudocode generation from
user stories. By doing so, we hope to free up developers to focus on more
complex tasks, reduce the risk of errors, and ultimately deliver better
software products to users.

##  II Literature Survey

Saeki, Horai, and Enomoto’s groundbreaking work in 1989 [15] marked the
beginning of the path to close the gap between natural language and
programming languages. With the aid of this fundamental discovery, software
development could be sped up by converting natural language specifications
into code. This work served as a conceptual starting point but also provided a
framework for further investigation.

Cohn’s important book ”User Stories Applied” from 2004 [13] became a
cornerstone as the software development industry turned towards more agile
approaches. User stories are critical to Agile development, and Cohn
underlined the value of quickly turning them into code. User stories have
since taken on a key role in Agile workflows, signaling a dramatic turn in the
field.

The emergence of automated code generation tools accelerated the shift from
user stories written in plain language to programming code. Notably, Pseudogen
was presented by Oda et al. in 2015 [19]. Even though it wasn’t just concerned
with text-to-pseudocode conversion, this work was a significant advancement in
automating the conversion of code from descriptions written in plain language.

With the introduction of deep learning and neural networks, the conversion of
text to pseudocode reached a turning point. In 2018, Devlin et al. unveiled
BERT [1], a revolutionary approach for interpreting natural language. Code
generation underwent a revolution thanks to BERT’s ability to record context,
which made it possible to comprehend programming-related language more
sophisticatedly.

By pre-training on data flow, GraphCodeBERT [2] expanded BERT’s coding
capabilities in 2021. This ground-breaking method, which made a substantial
advancement in the discipline, used deep bidirectional transformers to improve
code representation and comprehension. By addressing the need for models that
could comprehend both the code and the data flow within it, generated
pseudocode was of higher quality.

A unified approach for program understanding and citation generation was
proposed by Ahmad et al. in 2021 [3], signaling the convergence of various
research strands. By providing a thorough understanding of program structures,
this work aimed to improve code generation. It emphasized the possibility of
improving software development through an understanding of the entire
development cycle, from user stories to code generation.

The value of standardized datasets became clear as the field developed. By
serving as benchmarks for text-to-code and code-to-pseudocode conversion
methods, datasets like the MBPP Dataset [26] and the ASE15 Django Dataset [27]
significantly contributed to the advancement of technology. Researchers have
been able to rigorously assess and compare their models thanks to these
datasets.

Models like Code2Text [8] and DeepPseudo [10] have become more popular
recently. These models improve the precision and applicability of text-to-
pseudocode conversion by combining pre-training, deep learning, and code
feature extraction. They offer reliable solutions for this task and are at the
cutting edge of technology. These models are made to comprehend both the
structure of the programming context and the code, allowing for more precise
pseudocode generation.

##  III Methodology

Figure 1: Stages of the proposed approach

###  III-A Stages of our proposed approach

Our research paper aims to generate pseudocode from a given English language
prompt. In order to achieve this task, we have divided the workflow into two
stages, namely Text to Code Conversion and Code to Pseudocode Conversion. We
have fine-tuned CodeT5 Model to get the required output. The performance of
the model was measured in the form of BLEU score.

####  III-A1 Text to Code Conversion

This stage converts the initial Natural Language text input into Python code.
Python was chosen as the intermediate representation language due to its
common use and availability of datasets on Python code. Text-to-code
generation produces program code in a programming language from a natural
language description of a program as its input. An encoder (such as a
transformer-based design) creates a set of hidden states from the input text
after preprocessing it (for example, tokenization and part-of-speech tagging).
The output program code is produced by a decoder (such as an LSTM-based
architecture) using these concealed states as input. The final result is then
returned after postprocessing (for example, to fix syntax mistakes).

  * •

Input Text (natural language description):
This refers to the initial text input that is to be converted into code or
pseudocode. This could be a sentence, a paragraph, or a longer piece of text
that describes a computational task or problem to be solved.

  * •

Preprocessing Layers (tokenization, POS tagging):
This involves preparing the input text for further processing by breaking it
down into individual words or symbols (tokenization) and identifying the part
of speech of each token (POS tagging). Tokenization is the process of dividing
text into meaningful units, or tokens, while POS tagging is the process of
labeling each token with its corresponding part of speech, such as noun, verb,
adjective, etc.

  * •

Encoder (transformer-based architecture):
The encoder takes in the preprocessed text and processes it using a
transformer-based architecture, which learns to map the input text to a
numerical representation that captures its meaning. The transformer-based
architecture is a type of neural network architecture that has achieved state-
of-the-art performance in many natural language processing tasks.

  * •

Decoder (transformer-based architecture):
The decoder takes in the encoded representation of the input text and
generates the corresponding code or pseudocode output. The decoder also uses a
transformer-based architecture and is trained to generate code or pseudocode
that is consistent with the input text.

  * •

Postprocessing Layers (code syntax correction):
This stage involves checking the generated code for syntax errors and
correcting them to ensure that it adheres to the syntax rules of the target
programming language. This is an important step to ensure that the generated
code is syntactically correct and can be executed without errors.

  * •

Output Code (programming language code):
This refers to the final output of the text-to-code generation process, which
is a block of code written in a programming language. The programming language
used for the output code depends on the task or problem being solved.

####  III-A2 Code to Pseudocode Generation

The Python code generated in stage one is converted to pseudocode form. In
code-to-pseudocode generation, a programming language’s programme code serves
as the input, and a simplified, high-level description of the code is produced
as the output in everyday English. The input code is first preprocessed (for
example, by parsing the syntax), and after that it is sent via an encoder (for
example, a transformer-based architecture) to produce a set of hidden states.
The decoder (also a transformer-based architecture) uses these hidden states
as input and produces the output pseudocode. The final result is then
postprocessed (for example, to make the text simpler) and returned as the
output pseudocode.
Both stages of the proposed approach are analogous to a language translation
task. We use an encoder-decoder transformer model to first convert the English
text to Python code, and then the Python code to pseudocode.

  * •

Preprocessing Layers (syntax parsing):
This stage involves preparing the input code for further processing by
breaking it down into its component parts and identifying the relationships
between those parts (syntax parsing). Syntax parsing is the process of
analyzing the structure of the code and identifying its constituent parts,
such as variables, functions, and control structures.

  * •

Encoder (transformer-based architecture):
The encoder takes in the parsed code and processes it using a transformer-
based architecture to learn a numerical representation that captures its
meaning. This encoder is similar to the one used in the text-to-code
generation process, but it is trained on code rather than natural language
text.

  * •

Decoder (transformer-based architecture):
The decoder takes in the encoded representation of the input code and
generates the corresponding pseudocode output. This decoder is similar to the
one used in the text-to-code generation process, but it is trained to generate
pseudocode instead of code.

  * •

Postprocessing Layers (text simplification):
This stage involves simplifying the generated pseudocode to make it more
easily understandable to humans. This may involve removing redundant or
ambiguous information, simplifying complex expressions, or rephrasing complex
statements in simpler terms.

  * •

Output Pseudocode (natural language description):
This refers to the final output of the code-to-pseudocode generation process,
which is a simplified natural language description of the input code. The
output pseudocode is designed to be more easily understandable to humans than
the original code, and may involve rephrasing complex statements in simpler
terms, removing redundant information, and simplifying complex expressions.
The pseudocode may be used as a higher-level description of the code, making
it easier to understand and maintain.

TABLE I: Datasets used Stages |  Dataset |  Programming Language |  Number of Samples in the Dataset
---|---|---|---
Text to Code |  MBPP |  Python |  974
Code to Pseudocode |  Django |  Python |  16000

##### As shown in Table, our research utilized two datasets for the Text to
Code and Code to Pseudocode stages. The first dataset, MBPP, consisted of 974
samples and was used for the Text to Code stage. The second dataset, Django,
consisted of 16,000 samples and was used for the Code to Pseudocode stage.
This table provides important information about the datasets used in our
research, including their size and the programming language they were written
in.

###  III-B Transformers

Transformers are attention-based models that don’t use the typical recurrent
layers found in encoder-decoder designs, but rather use multi-headed self-
attention. Word embeddings from the input sequence are sent to the first
encoder.

The data is then transformed and transmitted to the next encoder. The final
encoder in the encoder-stack sends its output to every decoder in the stack of
decoders. For translation tasks, the Transformer may be trained significantly
more quickly than designs based on recurrent or convolutional layers.

The encoder and decoder layers in the Transformer architecture
LABEL:fig:t5_architecture each employ multi-headed self-attention processes as
opposed to conventional recurrent or convolutional layers. In order to produce
a sequence of context embeddings, the encoder takes in a sequence of input
embeddings and processes it through several layers. One token at a time, the
decoder creates an output sequence using the context embeddings and a series
of target embeddings. The model is tuned during training to reduce the
discrepancy between the goal sequence and the anticipated output sequence. The
Transformer is ideally suited for jobs that involve comprehending long-range
dependencies because of its self-attention mechanism, which enables it to pay
attention to various input and output sequences based on their relevance to
the present prediction.

###  III-C CodeT5 Model

Modern neural language model CodeT5 was created with the goal of producing
excellent source code from natural language inquiries. It is founded on the T5
architecture, which uses a framework for encoders and decoders based on
transformers. The model is pre-trained using extensive corpora of natural
language and code, which enables it to accurately capture the nuanced
relationships between normal language and code. An encoder that handles
natural language inputs and a decoder that produces outputs in the form of
code make up the architecture of CodeT5. Each token’s contextual
representation is produced by the encoder, a multi-layer transformer that
analyses the incoming text. The decoder, which transforms the encoder’s output
into a series of code tokens, is also a multi-layer device. CodeT5 also
employs a novel copy mechanism that allows it to directly copy tokens from the
input text to the output code, improving the model’s ability to handle rare or
out-of-vocabulary words. We train CodeT5 for 40 epochs for the Text to Code
conversion task, and for 5 epochs for the Code to Pseudocode conversion task.

###  III-D Rule Based Approach

This method makes use of a Python script to convert python code given to the
script as input in the form of a .py file. The Python code given as input is
converted into pseudocode with the help of a fixed set of rules.
There are three types of rules:

  * •

Basic conversion rules

  * •

Prefix conversion rules

  * •

Advanced conversion rules

The Python code is scanned line by line and pseudocode is generated for every
corresponding set of code. The output of the script is the pseudocode
generated in the form of a .txt file

###  III-E Disadvantages of Rule-based approach

The proposed approach is more robust, flexible and dynamically adaptive as
compared to the rule-based approach. The scope of the rule-based approach is
limited as compared to the extensive scope of our proposed approach. On human
evaluation and making use of evaluation metrics like the BLEU score, the
proposed approach performs better than the rule-based approach. The rule-based
approach fails to generate appropriate pseudocode if keywords beyond the scope
of the fixed set of conversion rules exist in the code. However, such a case
can be appropriately handled by the proposed approach that makes use of
transformers. In summary, the proposed approach is more efficient and accurate
as compared to the rule-based approach.

##  IV Performance Analysis

We have used BLEU score to analyse the performance of our Machine Learning
Models.

###  IV-A BLEU Score

A popular technique for assessing the quality of machine-generated
translations against a set of reference translations created by humans is
called BLEU (Bilingual Evaluation Understudy) [b28]. Its foundation is the
comparison of the n-gram overlap between the machine translation and the
reference translations.
The precision score, which assesses how many n-grams in the machine-generated
translation also present in the reference translations, is first computed for
each n-gram up to a specific maximum length (usually 4), before the BLEU score
is determined. In order to discourage short and inadequate translations, the
precision score is then adjusted by a brevity penalty factor that considers
how long the machine-generated translation is in comparison to the reference
translations. The final step is to integrate the updated accuracy scores into
a geometric mean to create the overall BLEU score.

TABLE II: BLEU Score Stages |  Dataset Used |  BLEU Score Value
---|---|---
Text to Code |  MBPP |  0.4
Code to Pseudocode |  Django |  0.74

The BLEU score values obtained in our study, as shown in Table, were 0.4 for
Text to Code and 0.74 for Code to Pseudocode. The BLEU score for the first
stage of the project denotes understandable to good translations, which is
evident from the syntactical correctness of the generated Python code.
However, there may be logical errors in the Python code due to insufficient
training data. The BLEU score for the second stage points to high quality
translations, with the generated pseudocode exhibiting a high level of logical
and structural correctness. This clearly shows that, given sufficient training
data, this methodology promises to give excellent results for the Text to
Pseudocode conversion task.

##  V Conclusion

This paper examines the enormous potential for programming work automation
offered by natural language processing. We propose a two-stage methodology to
convert English language user stories into pseudocode. The advantage of
pseudocode is that it allows easy conversion into any programming language of
the developer’s choice. The two stages of this approach are text to code
conversion and code to pseudocode conversion. Each of these stages is treated
as a language translation task. We use the CodeT5 model for this task, getting
a BLEU score of 0.4 for Stage 1 and 0.74 for Stage 2. Our proposed system
simplifies the software development process in organizations.
In the future, we plan to curate a larger dataset for English text to Python
code conversion. We can obtain higher accuracy and better generalisation to
new examples with a larger and more varied dataset. This would necessitate
considerable data gathering, but it might open the door to more efficient and
useful text to code conversion methods. Additionally, a pertinent text to
pseudocode dataset might be produced to help simplify the conversion
architecture.

## References

  * [1] Devlin, J., Chang, M.-W., Lee, K., Toutanova, K.: BERT: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805 (2018)
  * [2] GraphCodeBERT: Pre-training Code Representations with Data Flow ICLR 2021 · Daya Guo, Shuo Ren, Shuai Lu, Zhangyin Feng, Duyu Tang, Shujie Liu, Long Zhou, Nan Duan, Alexey Svyatkovskiy, Shengyu Fu, Michele Tufano, Shao Kun Deng, Colin Clement, Dawn Drain, Neel Sundaresan, Jian Yin, Daxin Jiang, Ming Zhou
  * [3] Unified Pre-training for Program Understanding and Generation Wasi Uddin Ahmad, Saikat Chakraborty, Baishakhi Ray, Kai-Wei Cha (2021)
  * [4] Code Generation from Natural Language with Less Prior Knowledge and More Monolingual Data Sajad Norouzi, Keyi Tang, Yanshuai Cao (2021)
  * [5] Retrieval Augmented Code Generation and Summarization Md Rizwan Parvez, Wasi Uddin Ahmad, Saikat Chakraborty, Baishakhi Ray, Kai-Wei Chang (2021)
  * [6] Fine-grained Pseudo-code Generation Method via Code Feature Extraction and Transformer Guang Yang, Yanlin Zhou, Xiang Chen, Chi Yu (2021)
  * [7] Learning to Generate Pseudo-code from Source Code using Statistical Machine Translation Yusuke Oda, Hiroyuki Fudaba, Graham Neubig, Hideaki Hata, Sakriani Sakti, Tomoki Toda, and Satoshi Nakamura (2021)
  * [8] Code2Text: Dual Attention Syntax Annotation Networks for Structure-Aware Code Translation Authors : Yun Xiong, Shaofeng Xu, Keyao Rong, Xinyue Liu, Xiangnan Kong, Shanshan Li, Philip Yu, Yangyong Zhu
  * [9] Home Database Systems for Advanced Applications Conference paper Code2Text: Dual Attention Syntax Annotation Networks for Structure-Aware Code Translation Yun Xiong, Shaofeng Xu, Keyao Rong, Xinyue Liu, Xiangnan Kong, Shanshan Li, Philip Yu & Yangyong Zhu
  * [10] DeepPseudo: Deep Pseudo-code Generation via Transformer and Code Feature Extraction Guang Yang, Xiang Chen, +1 author Chi Yu (2021)
  * [11] Dense Passage Retrieval for Open-Domain Question Answering Vladimir Karpukhin, Barlas Oğuz, Sewon Min, Patrick Lewis, Ledell Wu, Sergey Edunov, Danqi Chen, Wen-tau (2020)
  * [12] Shah, V.: The implication of agile traditional method as a practice in pharmaceutical industry (2017)
  * [13] Cohn, M.: User Stories Applied: For Agile Software Development. Addison- Wesley Professional (2004)
  * [14] Iyer, S., Konstas, I., Cheung, A., Zettlemoyer, L.: Mapping language to code in programmatic context. arXiv preprint arXiv:1808.09588 (2018)
  * [15] Saeki, M., Horai, H., Enomoto, H.: Software development process from natural language specification. In: Proceedings of the 11th International Conference on Software Engineering, pp. 6473 (1989)
  * [16] Gilson, F., Galster, M., Georis, F.: Extracting quality attributes from user stories for early architecture decision making. In: 2019 IEEE International Conference on Software Architecture Companion (ICSA-C), pp. 129136 (2019). IEEE
  * [17] Robeer, M., Lucassen, G., Van DerWerf, J.M.E., Dalpiaz, F., Brinkkemper, S.: Automated extraction of conceptual models from user stories via nlp. In: 2016 IEEE 24th International Requirements Engineering Conference (RE), pp. 196–205 (2016). IEEE
  * [18] Husain, H.,Wu, H.-H., Gazit, T., Allamanis, M., Brockschmidt, M.: Code- SearchNet challenge: Evaluating the state of semantic code search. arXiv preprint arXiv:1909.09436 (2019)
  * [19] Fudaba, H., Oda, Y., Akabe, K., Neubig, G., Hata, H., Sakti, S., Toda, T., Nakamura, S.: Pseudogen: A tool to automatically generate pseudo-code from source code. In: 2015 30th IEEE/ACM International Conference on Automated Software Engineering (ASE), pp. 824{829 (2015). IEEE
  * [20] Alhefdhi, A., Dam, H.K., Hata, H., Ghose, A.: Generating pseudo-code from source code using deep learning. In: 2018 25th Australasian Software Engineering Conference (ASWEC), pp. 21{25 (2018). IEEE
  * [21] Rai, S., Gupta, A.: Generation of pseudo code from the python source code using rule-based machine translation. arXiv preprint arXiv:1906.06117 (2019)
  * [22] Sutskever, I., Vinyals, O., Le, Q.V.: Sequence to sequence learning with neural networks. In: Advances in Neural Information Processing Systems, pp. 3104{3112 (2014)
  * [23] Devlin, J., Chang, M.-W., Lee, K., Toutanova, K.: BERT: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805 (2018)
  * [24] Liu, Y., Ott, M., Goyal, N., Du, J., Joshi, M., Chen, D., Levy, O., Lewis, M., Zettlemoyer, L., Stoyanov, V.: RoBERTa: A robustly optimized BERT pretraining approach. arXiv preprint arXiv:1907.11692 (2019)
  * [25] Dataset for text to code generation
  * [26] Dataset for code to pseudo-code generation
  * [27] BLEU: a Method for Automatic Evaluation of Machine Translation Kishore Papineni, Salim Roukos, Todd Ward, and Wei-Jing Zh


# Improved uniform error bounds on time-splitting methods for the long-time
dynamics of the Dirac equation with small potentials

Weizhu Bao  Department of Mathematics, National University of Singapore,
Singapore, 119076, Singapore
<EMAIL_ADDRESS>Yue Feng  Department of Mathematics, National University
of Singapore, Singapore, 119076, Singapore
<EMAIL_ADDRESS>Jia Yin  Applied Mathematics and Computational Research
(AMCR) Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
USA
<EMAIL_ADDRESS>
((Day Month Year); (Day Month Year))

###### Abstract

We establish improved uniform error bounds on time-splitting methods for the
long-time dynamics of the Dirac equation with small electromagnetic potentials
characterized by a dimensionless parameter $\varepsilon\in(0,1]$ representing
the amplitude of the potentials. We begin with a semi-discritization of the
Dirac equation in time by a time-splitting method, and then followed by a
full-discretization in space by the Fourier pseudospectral method. Employing
the unitary flow property of the second-order time-splitting method for the
Dirac equation, we prove uniform error bounds at $C(t)\tau^{2}$ and
$C(t)(h^{m}+\tau^{2})$ for the semi-discretization and full-discretization,
respectively, for any time $t\in[0,T_{\varepsilon}]$ with
$T_{\varepsilon}=T/\varepsilon$ for $T>0$, which are uniformly for
$\varepsilon\in(0,1]$, where $\tau$ is the time step, $h$ is the mesh size,
$m\geq 2$ depends on the regularity of the solution, and
$C(t)=C_{0}+C_{1}\varepsilon t\leq C_{0}+C_{1}T$ grows at most linearly with
respect to $t$ with $C_{0}\geq 0$ and $C_{1}>0$ two constants independent of
$t$, $h$, $\tau$ and $\varepsilon$. Then by adopting the regularity
compensation oscillation (RCO) technique which controls the high frequency
modes by the regularity of the solution and low frequency modes by phase
cancellation and energy method, we establish improved uniform error bounds at
$O(\varepsilon\tau^{2})$ and $O(h^{m}+\varepsilon\tau^{2})$ for the semi-
discretization and full-discretization, respectively, up to the long-time
$T_{\varepsilon}$. Numerical results are reported to confirm our error bounds
and to demonstrate that they are sharp. Comparisons on the accuracy of
different time discretizations for the Dirac equation are also provided.

###### keywords:

Dirac equation; long-time dynamics; time-splitting methods; improved uniform
error bound; regularity compensation oscillation.

(xxxxxxxxxx)

AMS Subject Classification: 35Q41, 65M70, 65N35, 81Q05

##  1 Introduction

In this paper, we consider the Dirac equation in one or two dimensions (1D or
2D), which can be represented in the two-component form with wave function
$\Phi:=\Phi(t,{\mathbf{x}})=(\phi_{1}(t,{\mathbf{x}}),\phi_{2}(t,{\mathbf{x}}))^{T}\in\mathbb{C}^{2}$
as[13, 14, 33]

$i\partial_{t}\Phi=\left(-i\sum_{j=1}^{d}\sigma_{j}\partial_{j}+\sigma_{3}\right)\Phi+\varepsilon\left(V(\mathbf{x})I_{2}-\sum_{j=1}^{d}A_{j}(\mathbf{x})\sigma_{j}\right)\Phi,\quad\mathbf{x}\in\Omega,$
(1)

where $\Omega\subset\mathbb{R}^{d}$ $(d=1,2)$ is a bounded domain equipped
with periodic boundary condition. Here, $i=\sqrt{-1}$, $t\geq 0$ is time,
$\mathbf{x}=(x_{1},\ldots,x_{d})^{T}\in\mathbb{R}^{d}$,
$\partial_{j}=\frac{\partial}{\partial x_{j}}(j=1,\ldots,d)$,
$\varepsilon\in(0,1]$ is a dimensionless parameter, $V(\mathbf{x})$ and
$A_{j}({\mathbf{x}})$ are the given real-valued time-independent electric and
magnetic potentials, respectively. $I_{2}$ is the $2\times 2$ identity matrix,
and $\sigma_{1},\sigma_{2},\sigma_{3}$ are the Pauli matrices defined as

$\sigma_{1}=\begin{pmatrix}0&\ \ 1\\\ 1&\ \
0\end{pmatrix},\quad\sigma_{2}=\begin{pmatrix}0&\ \ -i\\\ i&\ \
0\end{pmatrix},\quad\sigma_{3}=\begin{pmatrix}1&\ \ 0\\\ 0&\ \
-1\end{pmatrix}.$ (2)

In order to study the dynamics of the Dirac equation (1), the initial
condition is taken as

$\Phi(t=0,\mathbf{x})=\Phi_{0}(\mathbf{x}),\quad{\bf x}\in\overline{\Omega}.$
(3)

For the Dirac equation (1) with $\varepsilon=1$, i.e., the classical regime,
there are comprehensive analytical and numerical results in the literatures.
Along the analytical front, for the existence and multiplicity of bound states
and/or standing wave solutions, we refer to Refs. ES,GGT,Gross and references
therein. In the numerical aspect, different numerical methods have been
proposed and analyzed, such as the finite difference time domain (FDTD)
methods[4, 28], exponential wave integrator Fourier pseudospectral (EWI-FP)
method[4, 5], and time-splitting Fourier pseudospectral (TSFP) method[7, 8],
etc. For more details related to the numerical schemes, we refer to Refs.
BSG,FLB,Gosse,GSX,MK and references therein.

However, to our knowledge, much less research has been done for the long-time
dynamics when $\varepsilon\in(0,1]$. With a small $\varepsilon$, we can study
the long-time dynamics of the Dirac equation (1) for $t\in[0,T_{\varepsilon}]$
with $T_{\varepsilon}:=T/\varepsilon$ for $T>0$, while either $\varepsilon$ or
$T$ can be fixed (or changed) in the analysis. When $\varepsilon\in(0,1]$ is
fixed, e.g. $\varepsilon=1$, we would like to investigate how the errors of
different numerical schemes perform for $t\in[0,T_{\varepsilon}]$ with
increasing $T$. On the other hand, when $T$ is fixed, we are interested in
analyzing the dependency of errors for different numerical methods on the
parameter $\varepsilon$ for $t\in[0,T_{\varepsilon}]$. In our recent works, we
examined the long-time error bounds for finite difference time domain (FDTD)
methods and finite difference Fourier pseudospectral (FDFP) methods, where the
finite difference was applied to discretize the Dirac equation (1) in time and
different discretizations were taken in space. We rigorously proved that the
FDTD methods share an error bound at
$O(\frac{h^{2}}{\varepsilon}+\frac{\tau^{2}}{\varepsilon})$ up to the long-
time $T_{\varepsilon}$ with $h$ the spatial mesh size and $\tau$ the temporal
step size, while the FDFP methods exhibit an error bound at
$O(h^{m}+\frac{\tau^{2}}{\varepsilon})$ up to the long-time[18]. The uniform
spatial error bound with respect to $\varepsilon$ for the FDFP methods
indicate that these methods have better spatial resolution than the FDTD
methods to solve the Dirac equation with small potentials in the long-time
regime. Nevertheless, both types of methods face severe numerical burdens when
$\varepsilon\to 0^{+}$ due to the $\varepsilon$-dependence of the temporal
error. To deal with this problem, we adopted the exponential wave integrator
(EWI) method for time discretization, and obtained a uniform error bound at
$O(h^{m}+\tau^{2})$ up to the long-time[19].

To further improve the temporal error for $\varepsilon\in(0,1]$, we consider
the time-splitting methods, which have been widely used to numerically solve
dispersive partial differential equations (PDEs)[1, 21, 26, 29]. It has been
proved that the temporal error bounds grow linearly when the time-splitting
methods are applied to the Maxwell’s equations[11, 12]. In addition, the
unitary flow property of the time-splitting methods contributes to the time-
dependent uniform error bounds for the Schrödinger equation and the improved
uniform error bounds for the Schrödinger/nonlinear Schrödinger equation with
the help of the regularity compensation oscillation (RCO) technique[2, 3].
This improves the analysis in the previous studies and obtains surprising
error estimates for the time-splitting methods[2, 3]. The aim of this paper is
to establish the improved uniform error bounds on time-splitting methods for
the long-time dynamics of the Dirac equation with small potentials. First, we
prove a uniform error bound where the error constant grows linearly with
respect to the time $t$. Based on the error bound, for a given accuracy
$\delta_{0}$ and time step $\tau$, the second order time-splitting (Strang
splitting) could be applied to simulate the dynamics of the Dirac equation up
to $O(\delta_{0}/\tau^{2})$ when $\varepsilon=1$, i.e., the smaller the time
step $\tau$, the longer the dynamics that can be calculated. Then by employing
the RCO technique[2, 3], we establish improved uniform error bounds at
$O(\varepsilon\tau^{2}+\tau_{0}^{m})$ and
$O(h^{m}+\varepsilon\tau^{2}+\tau_{0}^{m})$ for the semi-discretization and
full-discretization, respectively, for the Dirac equation with
$O(\varepsilon)$-potentials up to the long-time at $O(1/\varepsilon)$. In the
error bounds, $m\geq 2$ depends on the regularity of the exact solution and
$\tau_{0}\in(0,1)$ is a fixed chosen parameter. When the solution is smooth,
i.e. $m\to\infty$, then the error bounds collapse to $O(\varepsilon\tau^{2})$
and $O(h^{m}+\varepsilon\tau^{2})$ for the semi-discretization and full-
discretization, respectively.

The rest of this paper is organized as follows. In Section 2, the second-order
time-splitting methods including the semi-discretization and full-
discretization for the long-time dynamics of the Dirac equation with small
potentials are presented. In Section 3, uniform error bounds for the time-
splitting methods are established up to the long-time at $O(1/\varepsilon)$
and the errors are shown to grow linearly with respect to the time $t$. In
Section 4, we prove the improved uniform error bounds rigorously by adopting
the RCO technique. Extensive numerical results are reported in Section 5.
Finally, some conclusions are drawn in Section 6. Throughout this paper, we
adopt the notation $A\lesssim B$ to represent that there exists a generic
constant $C>0$, which is independent of the mesh size $h$ and time step $\tau$
as well as the parameter $\varepsilon$ such that $|A|\leq CB$.

##  2 The time-splitting methods

For simplicity, in the following sections, we focus on the Dirac equation (1)
in 1D, i.e., $d=1$ in (1), for the numerical methods and corresponding
analysis. The methods and results can be easily generalized to (1) in 2D,
i.e., $d=2$, and to the four-component Dirac equation given in Refs. BCJT and
BCJY.

The Dirac equation (1) in 1D on the bounded computational domain
$\Omega=(a,b)$ with periodic boundary condition collapses to

$\displaystyle
i\partial_{t}\Phi=(-i\sigma_{1}\partial_{x}+\sigma_{3})\Phi+\varepsilon(V(x)I_{2}-A_{1}(x)\sigma_{1})\Phi,\qquad
x\in\Omega,\quad t>0,$ (4) $\displaystyle\Phi(t,a)=\Phi(t,b),\quad t\geq
0;\qquad\Phi(0,x)=\Phi_{0}(x),\quad x\in\overline{\Omega},$ (5)

where $\Phi:=\Phi(t,x)$ and $\Phi_{0}(a)=\Phi_{0}(b)$.

###  2.1 The semi-discretization

Define the operators

$\mathbf{T}:=-\sigma_{1}\partial_{x}-i\sigma_{3},\quad\mathbf{V}:=-i(V(x)I_{2}-A_{1}(x)\sigma_{1}),$
(6)

then (4) can be expressed by

$\partial_{t}\Phi(t,x)=(\mathbf{T}+\varepsilon\mathbf{V})\Phi(t,x),\
x\in\Omega,\ t>0.$ (7)

Take a time step size $\tau>0$ and denote the time grids as $t_{n}=n\tau$ for
$n=0,1,\ldots$. As both $\mathbf{T}$ and $\mathbf{V}$ are time-independent,
the solution to (7) with (5) can be propagated through

$\Phi(t_{n+1},x)=e^{\tau(\mathbf{T}+\varepsilon\mathbf{V})}\Phi(t_{n},x),\quad
n=0,1,\ldots.$ (8)

To approximate the operator $e^{\tau(\mathbf{T}+\varepsilon\mathbf{V})}$, we
apply the second-order time-splitting (Strang splitting)[32], which gives

$e^{\tau(\mathbf{T}+\varepsilon\mathbf{V})}\approx
e^{\frac{\tau}{2}\mathbf{T}}e^{\varepsilon\tau\mathbf{V}}e^{\frac{\tau}{2}\mathbf{T}}.$
(9)

Therefore, the semi-discretization of the Dirac equation (4) via Strang
splitting can be expressed as

$\Phi^{[n+1]}(x)=\mathcal{S}_{\tau}(\Phi^{[n]}(x)):=e^{\frac{\tau}{2}\mathbf{T}}e^{\varepsilon\tau\mathbf{V}}e^{\frac{\tau}{2}\mathbf{T}}\Phi^{[n]}(x),\quad
n=0,1,\ldots,$ (10)

where we take the initial condition $\Phi^{[0]}(x):=\Phi_{0}(x)$ for
$x\in\overline{\Omega}$. Here, $\Phi^{[n]}(x)$ is the approximation of
$\Phi(t_{n},x)$.

###  2.2 The full-discretization

By noticing the definition of $\mathbf{V}$ in (6), it is easy to derive that

$e^{\varepsilon\tau\mathbf{V}}\Phi(t,x)=e^{-i\varepsilon\tau\left(V(x)I_{2}-A_{1}(x)\sigma_{1}\right)}\Phi(t,x),\quad
x\in\overline{\Omega},\quad t>0.$ (11)

On the other hand, to get $e^{\tau\mathbf{T}}\Phi(t,x)$, we can discretize (7)
in space by the Fourier spectral method, and then it is possible to integrate
the operator analytically in the phase space. We take $M+1$ uniformly sampled
grid points in $\overline{\Omega}$ with $M$ a positive even integer

$x_{j}=a+jh,\quad h=\frac{b-a}{M},\quad j=0,1,\ldots,M,$ (12)

and we denote the sets $X_{M}$, $Y_{M}$, $Z_{M}$ as

$\displaystyle
X_{M}=\left\\{U=(U_{0},U_{1},\ldots,U_{M})^{T}~{}|~{}U_{j}\in\mathbb{C}^{2},\,j=0,1,\ldots,M,\,U_{0}=U_{M}\right\\},$
$\displaystyle Y_{M}=Z_{M}\times Z_{M},\quad
Z_{M}={\rm{span}}\\{\phi_{l}(x)=e^{i\mu_{l}(x-a)},\,l\in\mathcal{T}_{M}\\},$

where the index set
$\mathcal{T}_{M}=\\{l~{}|~{}l=-M/2,-M/2+1,\ldots,M/2-1\\}$, and $\mu_{l}=2\pi
l/(b-a)$ for $l\in\mathcal{T}_{M}$. The projection operator
$P_{M}:\left(L^{2}(\Omega)\right)^{2}\rightarrow Y_{M}$ is defined as

$\left(P_{M}U\right)(x):=\sum_{l\in\mathcal{T}_{M}}\widehat{U}_{l}e^{i\mu_{l}(x-a)},\quad
U(x)\in\left(L^{2}(\Omega)\right)^{2},$

where

$\widehat{U}_{l}=\frac{1}{b-a}\int_{a}^{b}U(x)e^{-i\mu_{l}(x-a)}dx,\quad
l\in\mathcal{T}_{M};$ (13)

and by taking $\left(C_{\rm
per}(\overline{\Omega})\right)^{2}=\left\\{U\in\left(C(\overline{\Omega})\right)^{2}~{}|~{}U(a)=U(b)\right\\}$,
the interpolation operator $I_{M}:\left(C_{\rm
per}(\overline{\Omega})\right)^{2}\rightarrow Y_{M}$ or
$I_{M}:X_{M}\rightarrow Y_{M}$ is defined as

$\left(I_{M}U\right)(x):=\sum_{l\in\mathcal{T}_{M}}\widetilde{U}_{l}e^{i\mu_{l}(x-a)},\qquad
U(x)\in\left(C_{\rm per}(\overline{\Omega})\right)^{2}\quad{\rm{or}}\quad U\in
X_{M},$

where

$\widetilde{U}_{l}=\frac{1}{M}\sum_{j=0}^{M-1}U_{j}e^{-2ijl\pi/M},\quad
l\in\mathcal{T}_{M}.$ (14)

Here we take $U_{j}=U(x_{j})$ if $U$ is a function.

Let $\Phi_{j}^{n}$ be the numerical approximation of $\Phi(t_{n},x_{j})$ and
denote
$\Phi^{n}=\left(\Phi_{0}^{n},\Phi_{1}^{n},\ldots,\Phi_{M}^{n}\right)^{T}\in
X_{M}$ as the solution vector at $t=t_{n}$. Take the initial value
$\Phi_{j}^{0}=\Phi_{0}(x_{j})$ for $j=0,\ldots,M$, then the time-splitting
Fourier pseudospectral (TSFP) method for discretizating the Dirac equation (4)
is given as

$\begin{split}&\Phi^{(1)}_{j}=\sum_{l\in\mathcal{T}_{M}}e^{-i\frac{\tau\Gamma_{l}}{2}}\widetilde{(\Phi^{n})}_{l}\;e^{i\mu_{l}(x_{j}-a)}=\sum_{l\in\mathcal{T}_{M}}Q_{l}\;e^{-i\frac{\tau
D_{l}}{2}}(Q_{l})^{T}\widetilde{(\Phi^{n})}_{l}\;e^{\frac{2ijl\pi}{M}},\\\
&\Phi^{(2)}_{j}=e^{-i\varepsilon\tau
G(x_{j})}\Phi^{(1)}_{j}=Pe^{-i\varepsilon\Lambda_{j}}P^{T}\Phi^{(1)}_{j},\\\
&\Phi^{n+1}_{j}=\sum_{l\in\mathcal{T}_{M}}e^{-i\frac{\tau\Gamma_{l}}{2}}\widetilde{(\Phi^{(2)})}_{l}\;e^{i\mu_{l}(x_{j}-a)}=\sum_{l\in\mathcal{T}_{M}}Q_{l}\;e^{-i\frac{\tau
D_{l}}{2}}(Q_{l})^{T}\widetilde{(\Phi^{(2)})}_{l}\;e^{\frac{2ijl\pi}{M}},\end{split}$
(15)

for $n=0,1,\ldots$, where
$\Gamma_{l}=\mu_{l}\sigma_{1}+\sigma_{3}=Q_{l}D_{l}(Q_{l})^{T}$ with
$\delta_{l}=\sqrt{1+\mu^{2}_{l}}$, $(Q_{l})^{T}$ is the transpose of $Q_{l}$
and

$\Gamma_{l}=\begin{pmatrix}1&\ \mu_{l}\\\ \mu_{l}&\ -1\end{pmatrix},\quad
Q_{l}=\begin{pmatrix}\frac{1+\delta_{l}}{\sqrt{2\delta_{l}(1+\delta_{l})}}&\
-\frac{\mu_{l}}{\sqrt{2\delta_{l}(1+\delta_{l})}}\\\
\frac{\mu_{l}}{\sqrt{2\delta_{l}(1+\delta_{l})}}&\
\frac{1+\delta_{l}}{\sqrt{2\delta_{l}(1+\delta_{l})}}\end{pmatrix},\quad
D_{l}=\begin{pmatrix}\delta_{l}&\ 0\\\ 0&\ -\delta_{l}\end{pmatrix},$ (16)

and $G(x_{j})=V(x_{j})I_{2}-A_{1}(x_{j})\sigma_{1}=P\Lambda_{j}P^{T}$ with
$\Lambda_{j}=\textrm{diag}(\Lambda_{j,+},\Lambda_{j,-})$ and
$\Lambda_{j,\pm}=V(x_{j})\pm A_{1}(x_{j})$, $P=I_{2}$ if $A_{1}(x_{j})=0$ and
otherwise

$P=\begin{pmatrix}\frac{1}{\sqrt{2}}\ &\frac{1}{\sqrt{2}}\\\
-\frac{1}{\sqrt{2}}\ &\frac{1}{\sqrt{2}}\end{pmatrix}.$

######  Remark 2.1.

If the electromagnetic potential $V$ and/or $A_{1}$ are time-dependent, the
time-ordering technique[35] should be applied when we implement (10), as the
operator $\mathbf{T}+\varepsilon\mathbf{V}$ is no longer commutable for
different time coordinates $t_{1}\neq t_{2}$.

##  3 Uniform error bounds

In this section, we prove the uniform error bounds for the second-order time-
splitting method in propagating the Dirac equation with small potentials in
the long-time regime up to $T_{\varepsilon}=T/\varepsilon$ for any given
$T>0$. We will start with the results for the semi-discretized scheme, and
then extend it to the full-discretization.

###  3.1 For semi-discretization

Suppose there exists a positive integer $m\geq 2$, such that for the
potentials, we have

$\textrm{(A)}\qquad\qquad V(x)\in{W^{m,\infty}_{\rm per}}(\Omega),\quad
A_{1}(x)\in{W^{m,\infty}_{\rm per}}(\Omega),$

where $W^{m,\infty}_{\rm per}(\Omega):=\\{u|u\in
W^{m,\infty}(\Omega),\,\partial^{l}_{x}u(a)=\partial^{l}_{x}u(b),\,l=0,\,\ldots,\,m-1\\}$.
In addition, we assume that the exact solution $\Phi:=\Phi(t,x)$ of the Dirac
equation (4) up to the long-time $T_{\varepsilon}=T/\varepsilon$ satisfies

$\textrm{(B)}\qquad\qquad\|\Phi\|_{L^{\infty}([0,T_{\varepsilon}];(H_{\rm
per}^{m}(\Omega))^{2})}\lesssim 1,$

where similarly, $H^{m}_{\rm per}(\Omega):=\\{u|u\in
H^{m}(\Omega),\,\partial^{l}_{x}u(a)=\partial^{l}_{x}u(b),\,l=0,\,\ldots,\,m-1\\}$,
with the equivalent $H^{m}$-norm on $H^{m}_{\rm per}(\Omega)$ given as

$\left\|\phi\right\|_{H^{m}}=\left(\sum\limits_{l\in\mathbb{Z}}(1+\mu_{l}^{2})^{m}|\widehat{\phi}_{l}|^{2}\right)^{1/2}.$
(17)

Taking into account the assumptions (A) and (B), we can easily derive

$\|\partial_{t}\Phi\|_{L^{\infty}([0,T_{\varepsilon}];(H_{\rm
per}^{m-1}(\Omega))^{2})}\lesssim 1,$ (18)

from the Dirac equation (1). By taking the semi-discretized second-order time-
splitting given by (10) for the Dirac equation (7) with the operators
$\mathbf{T}$ and $\mathbf{V}$ defined in (6), we have the following error
estimate.

######  Theorem 3.1.

Let $\Phi^{[n]}(x)$ be the numerical approximation obtained from the semi-
discretized second-order time-splitting (10) for the Dirac equation (7), then
under the assumptions (A) and (B), for any $0<\varepsilon\leq 1$, we have the
uniform error estimate

$\|\Phi(t_{n},x)-\Phi^{[n]}(x)\|_{L^{2}}\leq C_{1}\varepsilon
t_{n}\tau^{2}\leq C_{1}T\tau^{2},\quad 0\leq n\leq\frac{T/\varepsilon}{\tau},$
(19)

where $C_{1}$ is a positive constant independent of $h,\tau,n$ and
$\varepsilon$.

######  Proof 3.2.

We notice that $\mathbf{T}$ generates a unitary group in $(H^{k}_{\rm
per}(\Omega))^{2}$ $(k\geq 0)$. Denote the exact solution flow
$\Phi(t_{n})\to\Phi(t_{n+1})$ as

$\Phi(t_{n+1})=\mathcal{S}_{e,\tau}(\Phi(t_{n})),\quad 0\leq
n\leq\frac{T/\varepsilon}{\tau},$ (20)

where we take $\Phi(t_{n}):=\Phi(t_{n},x)$ for simplicity.

In order to prove the convergence, we adopt the approach via formal Lie
calculus introduced in Ref. Lubich and split the proof into the following two
steps.

Step 1 (Bounds for local truncation error). We begin with the local truncation
error, i.e., to estimate the error generated by one time step computed via
(10). By using Taylor expansion for $e^{\varepsilon\tau\mathbf{V}}$, we have

$\mathcal{S}_{\tau}(\Phi_{0})=e^{\tau{\mathbf{T}}}\Phi_{0}+\varepsilon\tau
e^{\frac{\tau{\mathbf{T}}}{2}}{\mathbf{V}}e^{\frac{\tau{\mathbf{T}}}{2}}\Phi_{0}+\varepsilon^{2}\tau^{2}\int^{1}_{0}(1-\theta)e^{\frac{\tau{\mathbf{T}}}{2}}e^{\varepsilon\theta\tau{\mathbf{V}}}{\mathbf{V}}^{2}e^{\frac{\tau{\mathbf{T}}}{2}}\Phi_{0}d\theta.$

By Duhamel’s principle, we can write

$\displaystyle\mathcal{S}_{e,\tau}(\Phi_{0})=$ $\displaystyle\
e^{\tau{\mathbf{T}}}\Phi_{0}+\varepsilon\int^{\tau}_{0}e^{(\tau-s)\mathbf{T}}{\mathbf{V}}e^{s\mathbf{T}}\Phi_{0}ds$
$\displaystyle+\varepsilon^{2}\int^{\tau}_{0}\int^{s}_{0}e^{(\tau-s)\mathbf{T}}{\mathbf{V}}e^{(s-w)\mathbf{T}}{\mathbf{V}}\Phi(w)dwds.$

Denote

$Y(s)=e^{(\tau-s)\mathbf{T}}{\mathbf{V}}e^{s\mathbf{T}}\Phi_{0},\quad
B(s,w)=e^{(\tau-s)\mathbf{T}}{\mathbf{V}}e^{(s-w)\mathbf{T}}{\mathbf{V}}e^{w\mathbf{T}}\Phi_{0},$
(21)

then the local truncation error can be written as

$\displaystyle\mathcal{S}_{\tau}(\Phi_{0})-\mathcal{S}_{e,\tau}(\Phi_{0})=$
$\displaystyle\ \varepsilon\tau
Y\left(\frac{\tau}{2}\right)-\varepsilon\int^{\tau}_{0}Y(s)ds+\frac{\varepsilon^{2}\tau^{2}}{2}B\left(\frac{\tau}{2},\frac{\tau}{2}\right)$
$\displaystyle-\varepsilon^{2}\int^{\tau}_{0}\int^{s}_{0}B(s,w)dsdw+\varepsilon^{2}R_{1}+\varepsilon^{2}R_{2},$

with

$\begin{split}&R_{1}=\tau^{2}\int^{1}_{0}(1-\theta)e^{\frac{\tau{\mathbf{T}}}{2}}e^{\varepsilon\theta\tau{\mathbf{V}}}{\mathbf{V}}^{2}e^{\frac{\tau{\mathbf{T}}}{2}}\Phi_{0}d\theta-\frac{\tau^{2}}{2}B\left(\frac{\tau}{2},\frac{\tau}{2}\right),\\\
&R_{2}=-\int^{\tau}_{0}\int^{s}_{0}e^{(\tau-s)\mathbf{T}}{\mathbf{V}}e^{(s-w)\mathbf{T}}{\mathbf{V}}\Phi(w)-B(s,w)dwds.\end{split}$

It is easy to check that

$\displaystyle\|R_{1}\|_{L^{2}}$
$\displaystyle\lesssim\tau^{2}\max_{\theta\in(0,1)}\\{\|\partial_{\theta\theta}((1-\theta)e^{\frac{\tau{\mathbf{T}}}{2}}e^{\varepsilon\theta\tau{\mathbf{V}}}{\mathbf{V}}^{2}e^{\frac{\tau{\mathbf{T}}}{2}}\Phi_{0})\|_{L^{2}}\\}$
$\displaystyle\lesssim\varepsilon\tau^{3}\left(\|{\mathbf{V}}^{3}\Phi_{0}\|_{L^{2}}+\varepsilon\tau\|{\mathbf{V}}^{4}\Phi_{0}\|_{L^{2}}\right)\lesssim\varepsilon\tau^{3}.$

In view of the properties of $\mathbf{T}$ and $B(s,w)$, the quadrature rule
implies

$\displaystyle\|R_{2}\|_{L^{2}}$
$\displaystyle\lesssim\tau^{3}\max_{s,w\in(0,\tau)}\\{\|e^{(\tau-s){\mathbf{T}}}{\mathbf{V}}e^{(s-w)\mathbf{T}}{\mathbf{V}}\partial_{w}\Phi(w)\|_{L^{2}}\\}$
$\displaystyle\lesssim\tau^{3}\|\partial_{s}\Phi(\cdot)\|_{L^{\infty}([0,\tau];({L^{2}})^{2})}\lesssim\tau^{3},$

and

$\left\|\frac{\tau^{2}}{2}B\left(\frac{\tau}{2},\frac{\tau}{2}\right)-\int^{\tau}_{0}\int^{s}_{0}B(s,w)dwds\right\|_{L^{2}}\lesssim\tau^{3}\max_{0\leq
w\leq
s\leq\tau}(\|\partial_{s}B\|_{L^{2}}+\|\partial_{w}B\|_{L^{2}})\lesssim\tau^{3}.$

Finally, we estimate the last term, which contains the major part of the local
error

$\mathcal{F}(\Phi_{0}):=\varepsilon\tau{Y}(\frac{\tau}{2})-\varepsilon\int^{\tau}_{0}Y(s)ds=\varepsilon\tau^{3}\int^{1}_{0}\text{ker}(\theta){Y}^{\prime\prime}(\theta\tau)d\theta,$
(22)

where $\text{ker}(\theta)$ is the Peano kernel for midpoint rule. In addition,
we have

$Y^{\prime\prime}(s)=e^{(\tau-s)\mathbf{T}}[{\mathbf{T}},[\mathbf{T},{\mathbf{V}}]]e^{s\mathbf{T}}\Phi_{0}.$

For the double commutator $[{\mathbf{T}},[{\mathbf{T}},{\mathbf{V}}]]$, we
have

$\left\|[{\mathbf{T}},[\mathbf{T},{\mathbf{V}}]]\Psi\right\|_{L^{2}}\lesssim\left(\|V(\cdot)\|_{W^{2,\infty}}+\|A_{1}(\cdot)\|_{W^{2,\infty}}\right)\|\Psi\|_{H^{1}}.$

Combining all the results above, we find the one step local error as

$\mathcal{S}_{\tau}(\Phi_{0})-\mathcal{S}_{e,\tau}(\Phi_{0})=\mathcal{F}(\Phi_{0})+\mathcal{R}^{0},$
(23)

where $\|\mathcal{R}^{0}\|_{L^{2}}\lesssim\varepsilon^{2}\tau^{3}$ and

$\mathcal{F}(\Phi_{0})=\varepsilon\tau^{3}\int^{1}_{0}\rm{ker}(\theta)e^{\tau(1-\theta){\mathbf{T}}}[{\mathbf{T}},[{\mathbf{T}},{\mathbf{V}}]]e^{\tau\theta{\mathbf{T}}}\Phi_{0}d\theta.$
(24)

Define the local truncation error at $t_{n}$ for $0\leq
n\leq\frac{T/\varepsilon}{\tau}-1$ as

$\mathcal{E}^{n}(x)=\mathcal{S}_{\tau}(\Phi(t_{n},x))-\mathcal{S}_{e,\tau}(\Phi(t_{n},x)),\quad
a\leq x\leq b,$ (25)

then from the above computation, we can get

$\mathcal{E}^{n}(x)=\mathcal{F}(\Phi(t_{n}))+\mathcal{R}^{n},\quad a\leq x\leq
b,\quad 0\leq n\leq\frac{T/\varepsilon}{\tau}-1,$ (26)

where for $0\leq n\leq\frac{T/\varepsilon}{\tau}-1$,

$\|\mathcal{F}(\Phi(t_{n}))\|_{L^{2}}\lesssim\varepsilon\tau^{3}\|\Phi(t_{n})\|_{H^{1}},\quad\|\mathcal{R}^{n}\|_{L^{2}}\lesssim\varepsilon^{2}\tau^{3}.$
(27)

Step 2 (Bounds for the global error). We are going to prove the error bound
(19). Denote ${\bf e}^{[n]}(x)=\Phi^{[n]}-\Phi(t_{n})$, then $\|{\bf
e}^{[0]}(x)\|_{L^{2}}=0$ by definition. For $0\leq
n\leq\frac{T/\varepsilon}{\tau}-1$, we have

${\bf
e}^{[n+1]}=\mathcal{S}_{\tau}(\Phi^{[n]})-\mathcal{S}_{\tau}(\Phi(t_{n}))+\mathcal{S}_{\tau}(\Phi(t_{n}))-\mathcal{S}_{e,\tau}(\Phi(t_{n})).$
(28)

By the error bound (27) for the local truncation error, we obtain for $0\leq
n\leq\frac{T/\varepsilon}{\tau}-1$,

$\left\|{\bf e}^{[n+1]}\right\|_{L^{2}}\leq\left\|{\bf
e}^{[n]}\right\|_{L^{2}}+C_{1}\varepsilon\tau^{3},$ (29)

where $C_{1}>0$ is a constant and independent of $h,n,\tau$ and $\varepsilon$.
Then it is straightforward to derive

$\left\|{\bf e}^{[n+1]}\right\|_{L^{2}}\leq\left\|{\bf
e}^{[0]}\right\|_{L^{2}}+C_{1}(n+1)\varepsilon\tau^{3}\leq C_{1}\varepsilon
t_{n+1}\tau^{2},$ (30)

which completes the proof of Theorem 3.1.

######  Remark 3.3.

According to Theorem 3.1, the uniform error bound on the time-splitting method
for the Dirac equation grows linearly with respect to $t$. In fact, given an
accuracy bound $\delta_{0}>0$, the time (for simplicity, assume
$\varepsilon=1$ here) for the second-order splitting method to violate the
accuracy requirement $\delta_{0}$ is $O(\delta_{0}/\tau^{2})$. The results can
be extended to other time-splitting methods. For the first-order Lie-Trotter
splitting and fourth-order compact splitting or partitioned Runge-Kutta
splitting (PRK4)[9, 34], the times are at $O(\delta_{0}/\tau)$ and
$O(\delta_{0}/\tau^{4})$, respectively. In other words, higher order splitting
method performs much better in the long-time simulation not only regarding the
higher accuracy but also longer simulation time to produce accurate solutions.
In addition, extensions to 2D/3D are straightforward.

###  3.2 For full-discretization

For the full-discretization given in (15) by the second-order time-splitting
method for the Dirac equation (7), we could further derive the following
uniform error estimate.

######  Theorem 3.4.

Let $\Phi^{n}$ be the approximation obtained from the TSFP (15) for the Dirac
equation (7). Under the assumptions (A) and (B), for any $0<\varepsilon\leq
1$, we have

$\|\Phi(t_{n},x)-I_{M}\Phi^{n}\|_{L^{2}}\leq
C(t_{n})\left(h^{m}+\tau^{2}\right),\qquad 0\leq
n\leq\frac{T/\varepsilon}{\tau},$ (31)

where $C(t)=C_{0}+C_{1}\varepsilon t\leq C_{0}+C_{1}T$ with $C_{0}$ and
$C_{1}$ two positive constants independent of $h,\tau,n$ and $\varepsilon$.

_Proof._ Noticing that

$I_{M}\Phi^{n}-\Phi(t_{n})=I_{M}\Phi^{n}-P_{M}\Phi(t_{n})+P_{M}\Phi(t_{n})-\Phi(t_{n}),$
(32)

under the assumption (B), we get from the standard Fourier projection
properties[31]

$\left\|I_{M}\Phi^{n}-\Phi(t_{n})\right\|_{L^{2}}\leq\left\|I_{M}\Phi^{n}-P_{M}\Phi(t_{n})\right\|_{L^{2}}+C_{2}h^{m},\
0\leq n\leq\frac{T/\varepsilon}{\tau}.$ (33)

Thus, it suffices to consider the error function ${\bf e}^{n}\in Y_{M}$ at
$t_{n}$ as

${\bf e}^{n}:={\bf e}^{n}(x)=I_{M}\Phi^{n}-P_{M}\Phi(t_{n}),\quad 0\leq
n\leq\frac{T/\varepsilon}{\tau}.$ (34)

Since $\Phi^{0}_{j}=\Phi_{0}(x_{j})$, it is obvious from (B) that $\|{\bf
e}^{0}\|_{L^{2}}\leq C_{3}h^{m}$. From the local truncation error (25), we
have the error equation for ${\bf e}^{n}$ ($0\leq
n\leq\frac{T/\varepsilon}{\tau}-1$),

$\displaystyle{\bf
e}^{n+1}=I_{M}\Phi^{n+1}-P_{M}\mathcal{S}_{\tau}(\Phi(t_{n}))+P_{M}(\mathcal{E}^{n}).$
(35)

Recall the semi-discretization (10) and the full-discretization (15), we have

$\begin{split}&I_{M}\Phi^{n+1}=e^{\frac{\tau{\mathbf{T}}}{2}}(I_{M}\Phi^{(2)}),\
I_{M}\Phi^{(2)}=I_{M}\left(e^{\varepsilon\tau\mathbf{V}}\Phi^{(1)}\right),\
I_{M}\Phi^{(1)}=e^{\frac{\tau{\mathbf{T}}}{2}}(I_{M}\Phi^{n}),\\\
&P_{M}(\mathcal{S}_{\tau}(\Phi(t_{n})))=e^{\frac{\tau{\mathbf{T}}}{2}}(P_{M}\Phi^{\langle
2\rangle}),\ \Phi^{\langle
2\rangle}:=e^{\varepsilon\tau\mathbf{V}}\Phi^{\langle 1\rangle},\
\Phi^{\langle
1\rangle}:=e^{\frac{\tau{\mathbf{T}}}{2}}\Phi(t_{n}).\end{split}$

In view of the facts that $I_{M}$ and $P_{M}$ are identical on $Y_{M}$ and
$e^{\tau{\mathbf{T}}/2}$ preserves the $H^{k}$-norm ($k\geq 0$), using the
Taylor expansion
$e^{\varepsilon\tau{\mathbf{V}}}=1+\varepsilon\tau{\mathbf{V}}\int^{1}_{0}e^{\varepsilon\tau\theta{\mathbf{V}}}d\theta$
and assumptions (A) and (B), we have

$\displaystyle\left\|I_{M}\Phi^{n+1}-P_{M}(\mathcal{S}_{\tau}(\Phi(t_{n})))\right\|_{L^{2}}=\left\|I_{M}\Phi^{(2)}-P_{M}\Phi^{\langle
2\rangle}\right\|_{L^{2}},$ (36) $\displaystyle\left\|P_{M}\Phi^{\langle
2\rangle}-I_{M}\Phi^{\langle 2\rangle}\right\|_{L^{2}}\leq
C_{4}\varepsilon\tau h^{m}.$ (37)

In addition, noticing $\mathbf{V}=-i(V(x)I_{2}-A_{1}(x)\sigma_{1})$, by direct
computation and Parseval’s identity, we have

$\displaystyle\left\|I_{M}\Phi^{(2)}-I_{M}\Phi^{\langle
2\rangle}\right\|_{L^{2}}$
$\displaystyle=\sqrt{h\sum_{j=0}^{M-1}\left|\Phi^{(2)}_{j}-\Phi^{\langle
2\rangle}(x_{j})\right|^{2}}=\sqrt{h\sum_{j=0}^{M-1}\left|\Phi^{(1)}_{j}-\Phi^{\langle
1\rangle}(x_{j})\right|^{2}}$
$\displaystyle=\left\|I_{M}\Phi^{(1)}-I_{M}\Phi^{\langle
1\rangle}\right\|_{L^{2}}=\left\|I_{M}\Phi^{n}-P_{M}\Phi(t_{n})\right\|_{L^{2}}$
$\displaystyle=\left\|{\bf e}^{n}\right\|_{L^{2}}.$ (38)

From (26), (27) and the assumption (B), it is clear that there exists
$C_{5}>0$ such that $\|\mathcal{E}^{n}\|_{L^{2}}\leq C_{5}\varepsilon\tau^{3}$
for $0\leq n\leq\frac{T/\varepsilon}{\tau}-1$. Taking the $L^{2}$-norm on both
sides of (35) and combining the above estimates, we obtain for $0\leq
n\leq\frac{T/\varepsilon}{\tau}-1$,

$\displaystyle\left\|{\bf e}^{n+1}\right\|_{L^{2}}$
$\displaystyle\leq\left\|I_{M}\Phi^{(2)}-P_{M}\Phi^{\langle
2\rangle}\right\|_{L^{2}}+\left\|\mathcal{E}^{n}\right\|_{L^{2}}$
$\displaystyle\leq\left\|I_{M}\Phi^{(2)}-I_{M}\Phi^{\langle
2\rangle}\right\|_{L^{2}}+\left\|P_{M}\Phi^{\langle
2\rangle}-I_{M}\Phi^{\langle
2\rangle}\right\|_{L^{2}}+\left\|\mathcal{E}^{n}\right\|_{L^{2}}$
$\displaystyle\leq\left\|{\bf
e}^{n}\right\|_{L^{2}}+C_{6}\left(\varepsilon\tau
h^{m}+\varepsilon\tau^{3}\right),$ (39)

where $C_{6}=\max\\{C_{4},C_{5}\\}$. Thus, we arrive at

$\left\|{\bf e}^{n+1}\right\|_{L^{2}}\leq C_{6}\varepsilon
t_{n+1}\left(h^{m}+\tau^{2}\right)+C_{3}h^{m},\quad 0\leq
n\leq\frac{T/\varepsilon}{\tau}-1,$ (40)

which completes the proof of Theorem 3.4 by taking $C_{0}=C_{2}+C_{3}$ and
$C_{1}=C_{6}$. $\square$

##  4 Improved uniform error bounds

In this section, we establish the improved uniform error bounds for the time-
splitting methods applied to the Dirac equation (4) up to the long-time
$T_{\varepsilon}$ under the assumptions (A) and (B).

###  4.1 For semi-discretization

######  Theorem 4.1.

Let $\Phi^{[n]}(x)$ be the numerical approximation obtained from the semi-
discretized second-order time-splitting (10) for the Dirac equation (7). Under
the assumptions (A) and (B), for $0<\tau_{0}<1$ sufficiently small and
independent of $\varepsilon$, when
$0<\tau<\alpha\frac{\pi(b-a)\tau_{0}}{\sqrt{\tau_{0}^{2}(b-a)^{2}+4\pi^{2}(1+\tau_{0})^{2}}}$
for a fixed constant $\alpha\in(0,1)$, we have the following improved uniform
error bound for any $\varepsilon\in(0,1]$

$\|\Phi(t_{n},x)-\Phi^{[n]}\|_{L^{2}}\lesssim\varepsilon\tau^{2}+\tau^{m}_{0},\quad
0\leq n\leq\frac{T/\varepsilon}{\tau}.$ (41)

In particular, if the exact solution is smooth, e.g., $\Phi(t,x)\in
L^{\infty}([0,T_{\varepsilon}];(H^{\infty}_{\rm per}(\Omega))^{2})$, the last
term $\tau_{0}^{m}$ decays exponentially fast and could be ignored practically
for small enough $\tau_{0}$, where the improved uniform error bound for
sufficiently small $\tau$ will be

$\|\Phi(t_{n},x)-\Phi^{[n]}\|_{L^{2}}\lesssim\varepsilon\tau^{2},\quad 0\leq
n\leq\frac{T/\varepsilon}{\tau}.$ (42)

######  Proof 4.2.

From the local truncation error (25) and the error equation (28), we have

${\bf
e}^{[n+1]}=\mathcal{S}_{\tau}(\Phi^{[n]})-\mathcal{S}_{\tau}(\Phi(t_{n}))+\mathcal{E}^{n}=e^{\tau\mathbf{T}}{\bf
e}^{[n]}+W^{n}+\mathcal{E}^{n},\quad n\geq 0,$ (43)

where $W^{n}:=W^{n}(x)$ ($n=0,1,\ldots$) is given by

$W^{n}(x)=\varepsilon\tau
e^{\frac{\tau\mathbf{T}}{2}}{\mathbf{V}}\int^{1}_{0}e^{\varepsilon\tau\theta\mathbf{V}}d\theta\;e^{\frac{\tau\mathbf{T}}{2}}{\bf
e}^{[n]}.$ (44)

Under the assumption (A), we have the following estimate

$\left\|W^{n}(x)\right\|_{L^{2}}\lesssim\varepsilon\tau\left\|{\bf
e}^{[n]}\right\|_{L^{2}}.$ (45)

Based on (43), we arrive at

${\bf e}^{[n+1]}=e^{(n+1)\tau\mathbf{T}}{\bf
e}^{[0]}+\sum_{k=0}^{n}e^{(n-k)\tau\mathbf{T}}\left(W^{k}+\mathcal{E}^{k}\right),\quad
0\leq n\leq\frac{T/\varepsilon}{\tau}-1.$ (46)

Combining (26), (27) and (45), noticing ${\bf e}^{[0]}={\bm 0}$, we have the
estimates for $0\leq n\leq\frac{T/\varepsilon}{\tau}-1$,

$\left\|{\bf
e}^{[n+1]}\right\|_{L^{2}}\lesssim\varepsilon\tau^{2}+\varepsilon\tau\sum_{k=0}^{n}\left\|{\bf
e}^{[k]}\right\|_{L^{2}}+\left\|\sum_{k=0}^{n}e^{(n-k)\tau\mathbf{T}}\mathcal{F}(\Phi(t_{k}))\right\|_{L^{2}}.$
(47)

In order to obtain the improved uniform error bounds, we shall employ the
regularity compensation oscillation (RCO) technique[2, 3] to deal with the
last term in the RHS of the inequality (47).

The key idea is a summation-by-parts procedure combined with spectrum cut-off
and phase cancellation. First, we employ the spectral projection on
$\Phi(t_{k})$ such that only finite Fourier modes of $\Phi(t_{k})$ need to be
considered and the projection error could be controlled by the regularity of
$\Phi(t_{k})$. Then, we apply the summation-by parts formula to the low
Fourier modes such that the phase could be cancelled for small $\tau$ (the
terms of the type $\sum_{k=0}^{n}e^{(n-k)\tau\mathbf{T}}$) and an extra order
of $\varepsilon$ could be gained from the terms like
$\mathcal{F}(\Phi(t_{k}))-\mathcal{F}(\Phi(t_{k+1}))$.

From the Dirac equation (7) and the assumption (A), we notice that
$\partial_{t}\Phi(t,x)-\mathbf{T}\Phi(t,x)=O(\varepsilon)$. In order to gain
an extra order of $\varepsilon$, it is natural to introduce the ‘twisted
variable’ as

$\Psi(t,x)=e^{-t\mathbf{T}}\Phi(t,x),\quad t\geq 0,$ (48)

which satisfies the equation

$\partial_{t}\Psi(t,x)=\varepsilon
e^{-t\mathbf{T}}\left(\mathbf{V}e^{t\mathbf{T}}\Psi(t,x)\right),\quad t\geq
0.$ (49)

Noticing $\mathbf{T}=-\sigma_{1}\partial_{x}-i\sigma_{3}$, under the
assumptions (A) and (B), we can prove that

$\|\Psi\|_{L^{\infty}([0,T_{\varepsilon}];(H_{\rm
per}^{m}(\Omega))^{2})}\lesssim
1,\quad\|\partial_{t}\Psi\|_{L^{\infty}([0,T_{\varepsilon}];(H_{\rm
per}^{m}(\Omega))^{2})}\lesssim\varepsilon$ (50)

and

$\left\|\Psi(t_{n+1})-\Psi(t_{n})\right\|_{H^{m}}\lesssim\varepsilon\tau,\quad
0\leq n\leq\frac{T/\varepsilon}{\tau}-1.$ (51)

The RCO technique will be used to force $\partial_{t}\Psi$ to appear with the
gain of order $\varepsilon$ for the summation-by-parts procedure in the last
term $\sum_{k=0}^{n}e^{(n-k)\tau\mathbf{T}}\mathcal{F}(\Phi(t_{k}))$, where
the small $\tau$ is required to control the accumulation of the frequency of
the type $e^{(n-k)\tau\mathbf{T}}$.

Step 1. As introduced in Ref. BCF1, we choose the cut-off parameter
$\tau_{0}\in(0,1)$ and $M_{0}=2\lceil 1/\tau_{0}\rceil\in\mathbb{Z}^{+}$
($\lceil\cdot\rceil$ is the ceiling function) with $1/\tau_{0}\leq
M_{0}/2<1+1/\tau_{0}$. Under the assumptions (A) and (B), we have the
following estimate

$\left\|P_{M_{0}}\mathcal{F}(P_{M_{0}}\Phi(t_{n}))-\mathcal{F}(\Phi(t_{n}))\right\|_{L^{2}}\lesssim\varepsilon\tau\tau_{0}^{m}.$
(52)

Based on the above estimates, (47) would imply for $0\leq n\leq
T_{\varepsilon}/\tau-1$,

$\left\|{\bf
e}^{[n+1]}\right\|_{L^{2}}\lesssim\tau_{0}^{m}+\varepsilon\tau^{2}+\varepsilon\tau\sum_{k=0}^{n}\left\|{\bf
e}^{[k]}\right\|_{L^{2}}+\left\|\mathcal{L}^{n}\right\|_{L^{2}},$ (53)

where

$\mathcal{L}^{n}=\sum_{k=0}^{n}e^{-(k+1)\tau\mathbf{T}}P_{M_{0}}\mathcal{F}\left(e^{t_{k}\mathbf{T}}(P_{M_{0}}\Psi(t_{k}))\right).$
(54)

Step 2. Now, we concentrate on $\mathcal{L}^{n}$, which represents the
summation of low Fourier modes. For $l\in\mathcal{T}_{M_{0}}$, define the
index set $\mathcal{I}_{l}^{M_{0}}$ associated to $l$ as

$\mathcal{I}_{l}^{M_{0}}=\left\\{(l_{1},l_{2})\ |\ l_{1}+l_{2}=l,\
l_{1}\in\mathbb{Z},\,l_{2}\in\mathcal{T}_{M_{0}}\right\\}.$ (55)

Following the notations in (16), we denote

$\Pi_{l}^{+}=Q_{l}\,\text{diag}(1,0)(Q_{l})^{T},\quad\Pi_{l}^{-}=Q_{l}\,\text{diag}(0,1)(Q_{l})^{T},$
(56)

where $\Pi_{l}^{\pm}$ are the projectors onto the eigenspaces of $\Gamma_{l}$
corresponding to the eigenvalues $\pm\delta_{l}$, respectively. Moreover, we
have $(\Pi_{l}^{\pm})^{T}=\Pi_{l}^{\pm}$, $\Pi_{l}^{+}+\Pi_{l}^{-}=I_{2}$,
$(\Pi_{l}^{\pm})^{2}=\Pi_{l}^{\pm}$, $\Pi_{l}^{\pm}\Pi_{l}^{\mp}={\bf 0}$. By
direct computation, we have

$e^{t{\bm
T}}P_{M_{0}}\Psi(t_{k})=\sum\limits_{l\in\mathcal{T}_{M_{0}}}\left(e^{-it\delta_{l}}\Pi_{l}^{+}+e^{it\delta_{l}}\Pi_{l}^{-}\right)\widehat{\Psi}_{l}(t_{k})e^{i\mu_{l}(x-a)}.$
(57)

According to the definition of $\mathcal{F}$ in (22), the expansion below
follows

$\displaystyle e^{-(k+1)\tau{\bm
T}}P_{M_{0}}\left(e^{(\tau-s)\mathbf{T}}\mathbf{V}e^{(t_{k}+s)\mathbf{T}}P_{M_{0}}\Psi(t_{k}))\right)$
$\displaystyle=\sum\limits_{l\in\mathcal{T}_{M_{0}}}\sum\limits_{(l_{1},l_{2})\in\mathcal{I}_{l}^{M_{0}}}\sum\limits_{\nu_{j}=\pm,j=1,2}\mathcal{G}_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}(s)e^{i\mu_{l}(x-a)},$
(58)

where $\mathcal{G}_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}(s)$ is a function of
$s$ defined as

$\mathcal{G}_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}(s)=e^{i(t_{k}+s)\delta^{\nu_{1},\nu_{2}}_{l,l_{2}}}\Pi_{l}^{\nu_{1}}\widehat{\mathbf{V}}_{l_{1}}\Pi_{l_{2}}^{\nu_{2}}\widehat{\Psi}_{l_{2}}(t_{k}),$
(59)

with
$\delta_{l,l_{2}}^{\nu_{1},\nu_{2}}=\nu_{1}\delta_{l}-\nu_{2}\delta_{l_{2}}$.
Then, the remainder term $\mathcal{L}^{n}$ in (53) reads

$\mathcal{L}^{n}=\varepsilon\sum\limits_{k=0}^{n}\sum\limits_{l\in\mathcal{T}_{M_{0}}}\sum\limits_{(l_{1},l_{2})\in\mathcal{I}_{l}^{M_{0}}}\sum\limits_{\nu_{j}=\pm,j=1,2}\lambda_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}e^{i\mu_{l}(x-a)},$
(60)

where the coefficients $\lambda_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}$ given by

$\displaystyle\lambda_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}$
$\displaystyle=\tau\mathcal{G}_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}(\tau/2)-\int_{0}^{\tau}\mathcal{G}_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}(s)\,ds$
$\displaystyle=r_{l,l_{2}}^{\nu_{1},\nu_{2}}e^{it_{k}\delta^{\nu_{1},\nu_{2}}_{l,l_{2}}}c_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}},$
(61)

and

$\displaystyle
c_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}=\Pi_{l}^{\nu_{1}}\widehat{\mathbf{V}}_{l_{1}}\Pi_{l_{2}}^{\nu_{2}}\widehat{\Psi}_{l_{2}}(t_{k}),$
(62) $\displaystyle r_{l,l_{2}}^{\nu_{1},\nu_{2}}=\tau
e^{i\tau\delta^{\nu_{1},\nu_{2}}_{l,l_{2}}/2}-\int^{\tau}_{0}e^{is\delta^{\nu_{1},\nu_{2}}_{l,l_{2}}}\;ds=O(\tau^{3}(\delta^{\nu_{1},\nu_{2}}_{l,l_{2}})^{2}).$
(63)

We only need to consider the case $\delta^{\nu_{1},\nu_{2}}_{l,l_{2}}\neq 0$
as $r_{l,l_{2}}^{\nu_{1},\nu_{2}}=0$ if
$\delta^{\nu_{1},\nu_{2}}_{l,l_{2}}=0$. For $l\in\mathcal{T}_{M_{0}}$ and
$(l_{1},l_{2})\in\mathcal{I}_{l}^{M_{0}}$, we have

$|\delta^{\nu_{1},\nu_{2}}_{l,l_{2}}|\leq
2\delta_{M_{0}/2}=2\sqrt{1+\mu_{{M_{0}}/2}^{2}}<2\sqrt{1+\frac{4\pi^{2}(1+\tau_{0})^{2}}{\tau_{0}^{2}(b-a)^{2}}},$
(64)

which implies for
$0<\tau\leq\alpha\frac{\pi(b-a)\tau_{0}}{\sqrt{\tau_{0}^{2}(b-a)^{2}+4\pi^{2}(1+\tau_{0})^{2}}}$
with $0<\alpha,\tau_{0}<1$, there holds
$\frac{\tau}{2}|\delta_{l,l_{2}}^{\nu_{1},\nu_{2}}|\leq\alpha\pi$. Denoting
$S^{\nu_{1},\nu_{2}}_{n,l,l_{2}}=\sum\limits_{k=0}^{n}e^{it_{k}\delta_{l,l_{2}}^{\nu_{1},\nu_{2}}}$
($n\geq 0$), for
$0<\tau\leq\alpha\frac{\pi(b-a)\tau_{0}}{\sqrt{\tau_{0}^{2}(b-a)^{2}+4\pi^{2}(1+\tau_{0})^{2}}}$,
we obtain

$|S^{\nu_{1},\nu_{2}}_{n,l,l_{2}}|\leq\frac{1}{|\sin(\tau\delta_{l,l_{2}}^{\nu_{1},\nu_{2}}/2)|}\leq\frac{C}{\tau|\delta_{l,l_{2}}^{\nu_{1},\nu_{2}}|},\quad
C=\frac{2\alpha\pi}{\sin(\alpha\pi)},\quad\forall n\geq 0.$ (65)

Using the summation-by-parts, we find that

$\displaystyle\sum_{k=0}^{n}\lambda_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}=$
$\displaystyle\
r_{l,l_{2}}^{\nu_{1},\nu_{2}}\left[\sum_{k=0}^{n-1}S^{\nu_{1},\nu_{2}}_{k,l,l_{2}}(c_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}-c_{k+1,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}})+S_{n,l,l_{2}}^{\nu_{1},\nu_{2}}c_{n,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}\right],$
(66)

with

$c_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}-c_{k+1,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}=\Pi_{l}^{\nu_{1}}\widehat{\mathbf{V}}_{l_{1}}\Pi_{l_{2}}^{\nu_{2}}\left(\widehat{\Psi}_{l_{2}}(t_{k})-\widehat{\Psi}_{l_{2}}(t_{k+1})\right).$
(67)

Combining (63), (65), (66), and (67), we have

$\left|\sum_{k=0}^{n}\lambda_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}\right|\lesssim\tau^{2}|\delta_{l,l_{2}}^{\nu_{1},\nu_{2}}|\left|\widehat{\mathbf{V}}_{l_{1}}\right|\left[\sum\limits_{k=0}^{n-1}\left|\widehat{\Psi}_{l_{2}}(t_{k})-\widehat{\Psi}_{l_{2}}(t_{k+1})\right|+\left|\widehat{\Psi}_{l_{2}}(t_{n})\right|\right].$
(68)

Step 3. Now, we are ready to give the improved estimates. For
$l\in\mathcal{T}_{M_{0}}$ and $(l_{1},l_{2})\in\mathcal{I}_{l}^{M_{0}}$,
simple calculations show ($l=l_{1}+l_{2}$)

$|\delta_{l,l_{2}}^{\nu_{1},\nu_{2}}|\lesssim\prod_{j=1}^{2}(1+\mu_{l_{j}}^{2})^{1/2}.$
(69)

Based on (60), (68) and (69), using Cauchy inequality, we have

$\displaystyle\left\|\mathcal{L}^{n}\right\|^{2}_{L^{2}}$ (70)
$\displaystyle=\varepsilon^{2}\sum\limits_{l\in\mathcal{T}_{N_{0}}}\big{|}\sum\limits_{(l_{1},l_{2})\in\mathcal{I}_{l}^{N_{0}}}\sum\limits_{\nu_{j}=\pm,j=1,2}\sum\limits_{k=0}^{n}\lambda_{k,l,l_{1},l_{2}}^{\nu_{1},\nu_{2}}\big{|}^{2}$
$\displaystyle\lesssim\varepsilon^{2}\tau^{4}\bigg{\\{}\sum_{l\in\mathcal{T}_{M_{0}}}\bigg{(}\sum\limits_{(l_{1},l_{2})\in\mathcal{I}_{l}^{N_{0}}}\left|\widehat{\mathbf{V}}_{l_{1}}\right|\left|\widehat{\Psi}_{l_{2}}(t_{n})\right|\prod_{j=1}^{2}(1+\mu_{l_{j}}^{2})^{1/2}\bigg{)}^{2}$
$\displaystyle\quad+n\sum\limits_{k=0}^{n-1}\bigg{[}\sum_{l\in\mathcal{T}_{M_{0}}}\bigg{(}\sum\limits_{(l_{1},l_{2})\in\mathcal{I}_{l}^{N_{0}}}\left|\widehat{\mathbf{V}}_{l_{1}}\right|\left|\widehat{\Psi}_{l_{2}}(t_{k})-\widehat{\Psi}_{l_{2}}(t_{k+1})\right|\prod_{j=1}^{2}(1+\mu_{l_{j}}^{2})^{1/2}\bigg{)}^{2}\bigg{]}\bigg{\\}}.$

In order to estimate each term in the above inequality, we use the auxiliary
function
$\Theta(x)=\sum_{l\in\mathbb{Z}}(1+\mu_{l}^{2})^{1/2}\left|\widehat{\Psi}_{l}(t_{n})\right|e^{i\mu_{l}(x-a)}$,
where $\Theta(x)\in\left(H_{\rm per}^{m-1}(\Omega)\right)^{2}$ which can be
verified from the assumption (B) and we can prove
$\|\Theta(x)\|_{H^{s}}\lesssim\|\Psi(t_{n})\|_{H^{s+1}}$ ($s\leq m-1$).
Similarly, introduce the function
$\mathbf{U}(x)=\sum_{l\in\mathbb{Z}}(1+\mu_{l}^{2})^{1/2}\left|\widehat{\mathbf{V}}_{l}\right|e^{i\mu_{l}(x-a)}$,
where $\mathbf{U}(x)\in{W^{m-1,\infty}_{\rm per}}(\Omega)$ as can be derived
from the assumption (A). We can prove directly that
$\|\mathbf{U}(x)\|_{W^{0,\infty}}\lesssim\|\mathbf{V}(x)\|_{W^{2,\infty}}$.
Expanding

$\mathbf{U}(x)\Theta(x)=\sum\limits_{l\in\mathbb{Z}}\sum\limits_{l_{1}+l_{2}=l}\prod_{j=1}^{2}(1+\mu_{l_{j}}^{2})^{1/2}\left|\widehat{\mathbf{V}}_{l_{1}}\right|\left|\widehat{\Psi}_{l_{2}}(t_{n})\right|e^{i\mu_{l}(x-a)},$

we can obtain

$\displaystyle\sum_{l\in\mathcal{T}_{M_{0}}}\bigg{(}\sum\limits_{(l_{1},l_{2})\in\mathcal{I}_{l}^{M_{0}}}\left|\widehat{\mathbf{V}}_{l_{1}}\right|\left|\widehat{\Psi}_{l_{2}}(t_{n})\right|\prod_{j=1}^{2}(1+\mu_{l_{j}}^{2})^{1/2}\bigg{)}^{2}$
$\displaystyle\leq\
\|\mathbf{U}(x)\Theta(x)\|^{2}_{L^{2}}\lesssim\|\mathbf{V}(x)\|_{W^{2,\infty}}^{2}\|\Psi(t_{n})\|_{H^{1}}^{2},$
(71)

which together with the assumption (A), (50) and (51) gives

$\displaystyle\left\|\mathcal{L}^{n}\right\|^{2}_{L^{2}}$
$\displaystyle\lesssim\varepsilon^{2}\tau^{4}\|\mathbf{V}(x)\|_{W^{2,\infty}}^{2}\bigg{(}\|\Psi(t_{n})\|_{H^{1}}^{2}+n\sum\limits_{k=0}^{n-1}\|\Psi(t_{k})-\Psi(t_{k+1})\|_{H^{1}}^{2}\bigg{)}$
$\displaystyle\lesssim\varepsilon^{2}\tau^{4}+n^{2}\varepsilon^{4}\tau^{6}\lesssim\varepsilon^{2}\tau^{4}.$
(72)

for $0\leq n\leq\frac{T/\varepsilon}{\tau}-1$, where the same trick is applied
to the rest terms. Combining (53) and (72), we have

$\|{\bf
e}^{n+1}\|_{L^{2}}\lesssim\tau_{0}^{m}+\varepsilon\tau^{2}+\varepsilon\tau\sum\limits_{k=0}^{n}\|{\bf
e}^{k}\|_{L^{2}},\quad 0\leq n\leq\frac{T/\varepsilon}{\tau}-1.$ (73)

Discrete Gronwall’s inequality would yield $\|{\bf
e}^{n+1}\|_{L^{2}}\lesssim\varepsilon\tau^{2}+\tau_{0}^{m}$ $(0\leq
n\leq\frac{T/\varepsilon}{\tau}-1)$, and the proof for the improved uniform
error bound (41) in Theorem 42 is completed.

######  Remark 4.3.

$\tau_{0}\in(0,1)$ is a parameter introduced in analysis and the requirement
on $\tau$ (essentially $\tau\lesssim\tau_{0}$) enables the improved estimates
on the low Fourier modes $|l|\leq 1/\tau_{0}$, where the error constant
depends on $\alpha$. $\tau_{0}$ can be arbitrary as long as the assumed
relation between $\tau$ and $\tau_{0}$ holds, i.e. $\tau_{0}$ could be either
fixed, or dependent on $\tau$, e.g.
$\tau_{0}=\frac{\sqrt{16\pi^{2}+(b-a)^{2}}}{\alpha(b-a)\pi}\tau$.

######  Remark 4.4.

The improved uniform error bounds are established for the second-order Strang
splitting method. Under appropriate assumptions of the exact solution, the
improved uniform error bounds can be extended to the first-order Lie-Trotter
splitting and the fourth-order PRK splitting method with improved uniform
error bounds at $\varepsilon\tau$ and $\varepsilon\tau^{4}$, respectively.

###  4.2 For full-discretization

For the TSFP method (15), we could establish the following improved uniform
error bounds.

######  Theorem 4.5.

Let $\Phi^{n}$ be the approximation obtained from the TSFP (15) for the Dirac
equation (7). Under the assumptions (A) and (B), for $0<\tau_{0}<1$
sufficiently small and independent of $\varepsilon$, when
$0<\tau\leq\alpha\frac{\pi(b-a)\tau_{0}}{\sqrt{\tau_{0}^{2}(b-a)^{2}+4\pi^{2}(1+\tau_{0})^{2}}}$
for a fixed constant $\alpha\in(0,1)$, we have

$\left\|\Phi(t_{n},x)-I_{M}\Phi^{n}\right\|_{L^{2}}\lesssim
h^{m}+\varepsilon\tau^{2}+\tau^{m}_{0},\quad 0\leq
n\leq\frac{T/\varepsilon}{\tau}$ (74)

for any $0<\varepsilon\leq 1$. In particular, if the exact solution is smooth,
e.g., $\Phi(t,x)\in L^{\infty}([0,T_{\varepsilon}];(H^{\infty}_{\rm
per}(\Omega))^{2})$, the improved uniform error bounds for sufficiently small
$\tau$ will be

$\left\|\Phi(t_{n},x)-I_{M}\Phi^{n}\right\|_{H^{1}}\lesssim
h^{m}+\varepsilon\tau^{2},\quad 0\leq n\leq\frac{T/\varepsilon}{\tau}.$ (75)

######  Proof 4.6.

From the error estimates in previous sections, we have for $0\leq
n\leq\frac{T/\varepsilon}{\tau}$,

$\left\|\Phi(t_{n},x)-\Phi^{[n]}\right\|_{L^{2}}\lesssim\varepsilon\tau^{2}+\tau^{m}_{0},\quad\left\|\Phi^{[n]}-P_{M}\Phi^{[n]}\right\|_{L^{2}}\lesssim
h^{m}.$ (76)

Since
$\Phi(t_{n},x)-I_{M}\Phi^{n}=\Phi(t_{n},x)-\Phi^{[n]}+\Phi^{[n]}-P_{M}\Phi^{[n]}+P_{M}\Phi^{[n]}-I_{M}\Phi^{n}$,
we derive

$\left\|\Phi(t_{n},x)-I_{M}\Phi^{n}\right\|_{L^{2}}\leq\left\|P_{M}\Phi^{[n]}-I_{M}\Phi^{n}\right\|_{L^{2}}+C_{1}\left(h^{m}+\varepsilon\tau^{2}+\tau^{m}_{0}\right),$
(77)

where $C_{1}$ is a constant independent of $h$, $\tau$, $n$, $\varepsilon$ and
$\tau_{0}$. As a result, it remains to establish the estimates on the error
function ${\widetilde{\bf e}}^{n}:={\widetilde{\bf e}}^{n}(x)\in Y_{M}$ given
as

${\widetilde{\bf e}}^{n}(x):=P_{M}\Phi^{[n]}-I_{M}\Phi^{n},\quad 0\leq
n\leq\frac{T/\varepsilon}{\tau}.$

From (15) and (10), we get

$\displaystyle
I_{M}\Phi^{n+1}=e^{\frac{\tau{\mathbf{T}}}{2}}\left(I_{M}\left(e^{\varepsilon\tau\mathbf{V}}e^{\frac{\tau{\mathbf{T}}}{2}}(I_{M}\Phi^{n})\right)\right),$
$\displaystyle
P_{M}\Phi^{[n+1]}=e^{\frac{\tau{\mathbf{T}}}{2}}\left(P_{M}\left(e^{\varepsilon\tau\mathbf{V}}e^{\frac{\tau{\mathbf{T}}}{2}}(P_{M}\Phi^{[n]})\right)\right),$

which lead to

$\displaystyle\widetilde{\bf e}^{n+1}=e^{\tau\mathbf{T}}\widetilde{\bf
e}^{n}+W^{n}(x),$ (78)

where

$W^{n}(x)=e^{\frac{\tau\mathbf{T}}{2}}\left[P_{M}\left((e^{\varepsilon\tau\mathbf{V}}-1)e^{\frac{\tau{\mathbf{T}}}{2}}(P_{M}\Phi^{[n]})\right)-I_{M}\left((e^{\varepsilon\tau\mathbf{V}}-1)e^{\frac{\tau{\mathbf{T}}}{2}}(I_{M}\Phi^{n})\right)\right].$

Similar to the error estimates in Ref. BCJY, we have the following error
bounds

$\left\|W^{n}(x)\right\|_{L^{2}}\lesssim\varepsilon\tau\left(h^{m}+\left\|\widetilde{\bf
e}^{n}\right\|_{L^{2}}\right),$ (79)

Thus, we could obtain

$\left\|\widetilde{\bf e}^{n+1}\right\|_{L^{2}}\leq\left\|\widetilde{\bf
e}^{n}\right\|_{L^{2}}+C_{2}\varepsilon\tau\left(h^{m}+\left\|\widetilde{\bf
e}^{n}\right\|_{L^{2}}\right),\quad 0\leq n\leq\frac{T/\varepsilon}{\tau}-1,$
(80)

where $C_{2}$ is a constant independent of $h$, $\tau$, $n$ and $\varepsilon$.
Since $\widetilde{\bf e}^{0}=P_{M}\Phi_{0}-I_{M}\Phi_{0}$, we have
$\left\|\widetilde{\bf e}^{0}\right\|_{L^{2}}\lesssim h^{m}$ and discrete
Gronwall’s inequality implies $\left\|\widetilde{\bf
e}^{n+1}\right\|_{L^{2}}\lesssim h^{m}$ for $0\leq
n\leq\frac{T/\varepsilon}{\tau}-1$. Combining the above estimates with (77),
we derive

$\left\|\Phi(t_{n},x)-I_{M}\Phi^{n}\right\|_{L^{2}}\lesssim
h^{m}+\varepsilon\tau^{2}+\tau^{m}_{0},\quad 0\leq
n\leq\frac{T/\varepsilon}{\tau},$

which shows the improved uniform error bound (74) and the proof for Theorem 75
is completed.

##  5 Numerical results

In this section, we present numerical results of the TSFP method for the long-
time dynamics of the Dirac equation with $O(\varepsilon)$-potentials up to the
long-time $T_{\varepsilon}=T/\varepsilon$.

###  5.1 For $\varepsilon=1$ with $T\gg 1$ regime

First, we show an example to confirm that the uniform error bound grows
linearly with respect to the time $t$. We take $\Omega=(0,1)$, the
electromagnetic potentials

$V(x)=x^{2}(x-1)^{2}+1,\quad A_{1}(x)=x^{2}(x-1)^{2}+1,\quad x\in\Omega.$ (81)

and the $(H^{2}_{\rm per}(\Omega))^{2}$ initial data

$\phi_{1}(x)=\phi_{2}(x)=\frac{1}{2}x^{2}(1-x)^{2}+3,\quad x\in\Omega.$ (82)

The regularity is enough to ensure the uniform and improved error bounds. The
‘exact’ solution $\Phi(t,x)$ is obtained numerically by the TSFP (15) with a
very fine mesh size $h_{e}=1/128$ and time step size $\tau_{e}=10^{-4}$. To
quantify the error, we introduce the following error functions:

$e(t_{n})=\left\|\Phi(t_{n},x)-I_{N}\Phi^{n}\right\|_{L^{2}},\quad
e_{\max}(t_{n})=\max_{0\leq q\leq n}e(t_{q}).$

In the rest of the paper, the spatial mesh size is always chosen sufficiently
small and thus spatial errors can be ignored when considering the long-time
error growth and/or the temporal errors.

Figure 1: Long-time temporal errors of the TSFP (15) for the Dirac equation
(4) with $\varepsilon=1$ and different time step $\tau$.

Figure 2: Long-time temporal errors of the PRK4 method for the Dirac equation
(4) with $\varepsilon=1$ and different time step $\tau$.

Fig. 1 depicts the long-time temporal errors of the TSFP method for the Dirac
equation (4) with $\varepsilon=1$ and different time step $\tau$, which shows
that the uniform errors grows linearly with respect to the time. In addition,
for a given accuracy bound, the time to exceed the error bar is quadruple when
the time step is half, which also confirms the linear growth. For comparisons,
Fig. 2 plots the long-time errors of the PRK4 method, which indicates that
higher order time-splitting methods could get better accuracy with the same
time step size as well as longer time simulations within a given accuracy
bound.

###  5.2 For $\varepsilon\to 0$ with fixed $T$ regime

Next, we report the convergence test for the TSFP method (15) for the Dirac
equation (4) with the electromagnetic potentials (81) and the initial data
(82).

Figure 3: Long-time temporal errors of the TSFP (15) for the Dirac equation
(4) with $\tau=0.1$ and different $\varepsilon$.

Figure 4: Temporal convergence rates of the TSFP (15) for the Dirac equation
(4) over long-time dynamics at $t=1/\varepsilon$: convergence rates in $\tau$
(a), and convergence rates in $\varepsilon$ (b).

Figure 5: Spatial convergence rates of the TSFP (15) for the Dirac equation
(4) over long-time dynamics at $t=1/\varepsilon$: convergence rates in $\tau$
(a), and convergence rates in $\varepsilon$ (b).

Fig. 3 plots the long-time errors of the TSFP method for the Dirac equation
(4) with the fixed time step $\tau$ and different $\varepsilon$, which
confirms the improved uniform error bound at $O(\varepsilon\tau^{2}$) up to
the long-time at $O(1/\varepsilon)$. Figs. 4 & 5 exhibit the temporal and
spatial errors of the TSFP (15) for the Dirac equation (4) at
$t=1/\varepsilon$. Fig. 4(a) shows the second-order convergence of the TSFP
method in time. Each line in Fig. 4(b) gives the global errors with a fixed
time step $\tau$ and verifies that the global error performs like
$O(\varepsilon\tau^{2})$ up to the long-time at $O(1/\varepsilon)$. Each line
in Fig. 5(a) shows the spectral accuracy of the TSFP method in space and Fig.
5(b) verifies the spatial errors are independent of the small parameter
$\varepsilon$ in the long-time regime.

###  5.3 Comparisons of different temporal discretizaitons

In this subsection, we compare the long-time temporal errors of the time-
splitting methods with the finite difference method (FDM) and the exponential
wave integrator (EWI) method[18, 19]. In order to compare the temporal errors,
we adopt the Fourier pseudospectral method in space combined with each
temporal discretization and choose a fine mesh size such that the spatial
errors are neglected.

Figure 6: Temporal convergence rates of different temporal discretizations for
the Dirac equation (4) over long-time dynamics at $t=1/\varepsilon$:
convergence rates in $\tau$ (a), and convergence rates in $\varepsilon$ (b).

Fig. 6(a) displays the temporal errors for the fixed $\varepsilon=1$ with
different time step $\tau$. For the three second-order schemes, the second-
order (Strang) time-splitting method obtains smaller temporal errors than the
other two methods with the same time step. The fourth-order time-splitting
method (PRK4) not only has higher order convergent rate but also gives much
smaller errors than the other three methods with the same time step. Fig. 6(b)
shows the long-time temporal errors of these methods for different
$\varepsilon$ with with a fixed time step $\tau$. The splitting methods have
improved uniform error bounds like $O(\varepsilon\tau^{2})$ up to the long-
time at $O(1/\varepsilon)$. The EWI method has uniform error bounds, while the
long-time temporal errors of the finite difference method depend on the
parameter $\varepsilon$ and behave like $O(1/\varepsilon)$. As a result, time-
splitting methods perform much better than FDM and EWI in the long-time
simulations.

##  6 Conclusions

Improved uniform error bounds for the time-splitting methods for the long-time
dynamics of the Dirac equation with small electromagnetic potentials were
rigorously established. With the help of the unitary property of the solution
flow in $L^{2}(\Omega)$, the linear growth of the uniform error bound for the
time-splitting methods was strictly proven. By employing the regularity
compensation oscillation (RCO) technique, the improved uniform error bounds
were proved to be $O(\varepsilon\tau^{2})$ and $O(h^{m}+\varepsilon\tau^{2})$
up to the long-time at $O(1/\varepsilon)$ for the semi-discretization and
full-discretization, respectively. Numerical results were shown to validate
our error bounds and to demonstrate that they are sharp. Finally, comparisons
of different time discretizations were presented to illustrate the superior
property of the time-splitting methods for the numerical simulation of the
long-time dynamics of the Dirac equation.

## Acknowledgment

This work was partially supported by the Ministry of Education of Singapore
grant MOE2019-T2-1-063 (R-146-000-296-112, W. Bao & Y. Feng).

## References

  * [1] W. Bao and Y. Cai, Mathematical theory and numerical methods for Bose-Einstein condensation, Kinet. Relat. Mod. 6 (2013) 1–135.
  * [2] W. Bao, Y. Cai and Y. Feng, Improved uniform error bounds for the time-splitting methods for the long-time dynamics of the Schrödinger/nonlinear Schrödinger equation, arXiv: 2109.08940.
  * [3] W. Bao, Y. Cai and Y. Feng, Improved uniform error bounds on time-splitting methods for long-time dynamics of the nonlinear Klein–Gordon equation with weak nonlinearity, arXiv: 2109.14902.
  * [4] W. Bao, Y. Cai, X. Jia and Q. Tang, Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime, J. Sci. Comput. 71 (2017) 1094–1134.
  * [5] W. Bao, Y. Cai, X. Jia and J. Yin, Error estimates of numerical methods for the nonlinear Dirac equation in the nonrelativistic limit regime, Sci. China Math. 59 (2016) 1461–1494.
  * [6] W. Bao, Y. Cai and J. Yin, Super-resolution of time-splitting methods for the Dirac equation in the nonrelativistic regime, Math. Comp. 89 (2020) 2141–2173.
  * [7] W. Bao, Y. Cai and J. Yin, Uniform error bounds of time-splitting methods for the nonlinear Dirac equation in the nonrelativistic limit regime, SIAM J. Numer. Anal. 59 (2021) 1040–1066.
  * [8] W. Bao and J. Yin, A fourth-order compact time-splitting Fourier pseudospectral method for the Dirac equation, Res. Math. Sci. 6 (2019) article 11.
  * [9] S. Blanes and P. C. Moan, Pratical symplectic partitioned Runge–Kutta and Runge–Kutta–Nyström methods, J. Comput. Appl. Math. 142 (2002) 313–330.
  * [10] J. W. Braun, Q. Su and R. Grobe, Numerical approach to solve the time-dependent Dirac equation, Phys. Rev. A 59 (1999) 604–612.
  * [11] W. Chen, X. Li and D. Liang, Energy-conserved splitting FDTD methods for Maxwell’s equations, Numer. Math. 108 (2008) 445–485.
  * [12] W. Chen, X. Li and D. Liang, Energy-conserved splitting finite-difference time-domain methods for Maxwell’s equations in three dimensions, SIAM J. Numer. Anal. 48 (2010) 1530–1554.
  * [13] P. A. M. Dirac, The quantum theory of the electron, Proc. R. Soc. Lond. A 117 (1928) 610–624.
  * [14] P. A. M. Dirac, Principles of Quantum Mechanics (Oxford University Press, 1958).
  * [15] G. Dujardin and E. Faou, Long time behavior of splitting methods applied to the linear Schrödinger equation, C. R. Acad. Sci. Paris 344 (2007) 89–92.
  * [16] G. Dujardin and E. Faou, Normal form and long time analysis of splitting schemes for the linear Schrödinger equation with small potential, Numer. Math. 108 (2007) 223–262.
  * [17] M. Esteban and E. Séré, Existence and multiplicity of solutions for linear and nonlinear Dirac problems, Partial Differ. Equ. Appl. 12 (1997) 107–112.
  * [18] Y. Feng and J. Yin, Spatial resolution of different discretizations over long-time for the Dirac equation with small potentials, arXiv: 2105.10468.
  * [19] Y. Feng, Z. Xu, and J. Yin, Uniform error bounds of exponential wave integrator methods for the long-time dynamics of the Dirac equation with small potentials, Appl. Numer. Math. 172 (2022) 50–66.
  * [20] F. Fillion-Gourdeau, E. Lorin and A. D. Bandrauk, Numerical solution of the time-dependent Dirac equation in coordinate space without fermion-doubling, Comput. Phys. Commun. 183 (2012) 1403–1415.
  * [21] L. Gauckler and C. Lubich, Splitting integrators for nonlinear Schrödinger equations over long times, Found. Comput. Math. 10 (2010) 275–302.
  * [22] F. Gesztesy, H. Grosse and B. Thaller, A rigorous approach to relativistic corrections of bound state energies for spin-1/2 particles, Ann. Inst. Henri Poincaré Phys. Theor. 40 (1984) 159–174.
  * [23] L. Gosse, A well-balanced and asymptotic-preserving scheme for the one-dimensional linear Dirac equation, BIT 55 (2015) 433–458.
  * [24] L. Gross, The Cauchy problem for the coupled Maxwell and Dirac equations, Commun. Pure Appl. Math. 19 (1966) 1–15.
  * [25] B.-Y. Guo, J. Shen and C.-L. Xu, Spectral and pseudospectral approximations using Hermite functions: Application to the Dirac equation, Adv. Comput. Math. 19 (2003) 35–55.
  * [26] E. Hairer, C. Lubich and G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations (Springer-Verlag, 2002).
  * [27] C. Lubich, On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations, Math. Comp. 77 (2008) 2141–2153.
  * [28] Y. Ma and J. Yin, Error bounds of the finite difference time domain methods for the Dirac equation in the semiclassical regime, J. Sci. Comput. 81 (2019) 1801–1822.
  * [29] R. I. McLachlan and G. R. W. Quispel, Splitting methods, Acta Numer. 11 (2002) 341–434.
  * [30] G. R. Mocken and C. H. Keitel, Quantum dynamics of relativistic electrons, J. Comput. Phys. 199 (2004) 558-588.
  * [31] J. Shen, T. Tang and L.-L. Wang, Spectral Methods: Algorithms, Analysis and Applications (Springer-Verlag, 2011).
  * [32] G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506–517.
  * [33] B. Thaller, The Dirac Equation (Springer-Verlag, 1992).
  * [34] H. F. Trotter, On the product of semi-groups of operators, Proc. Amer. Math. Soc. 10 (1959) 545–551.
  * [35] J. Yin, A fourth-order compact time-splitting method for the Dirac equation with time-dependent potentials, J. Comput. Phys. 430 (2021) 110109.