# Improving Security in McAdams Coefficient-Based Speaker Anonymization by
Watermarking Method

Candy Olivia Mawalim and Masashi Unoki Japan Advanced Institute of Science and
Technology,
1-1 Asahidai, Nomi, Ishikawa 923–1292 Japan
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Speaker anonymization aims to suppress speaker individuality to protect
privacy in speech while preserving the other aspects, such as speech content.
One effective solution for anonymization is to modify the McAdams coefficient.
In this work, we propose a method to improve the security for speaker
anonymization based on the McAdams coefficient by using a speech watermarking
approach. The proposed method consists of two main processes: one for
embedding and one for detection. In embedding process, two different McAdams
coefficients represent binary bits “0” and “1”. The watermarked speech is then
obtained by frame-by-frame bit inverse switching. Subsequently, the detection
process is carried out by a power spectrum comparison. We conducted objective
evaluations with reference to the VoicePrivacy 2020 Challenge (VP2020) and of
the speech watermarking with reference to the Information Hiding Challenge
(IHC) and found that our method could satisfy the blind detection,
inaudibility, and robustness requirements in watermarking. It also
significantly improved the anonymization performance in comparison to the
secondary baseline system in VP2020.

##  1 Introduction

Speech technology has been developed rapidly in recent years, with
applications such as a virtual assistant that receives speech commands as its
input. On the one hand, this technology helps us to live more conveniently. On
the other hand, it can cause privacy issues because speech encapsulates
personal identifiable information (PII) such as age, gender, health, and
emotional state [1]. This issue has been escalating due to the increasing
availability of advanced speech synthesis technology. At this stage, with only
a few speech samples spoken by a target speaker, we could easily clone a new
utterance spoken by the speaker [2].

To preserve speech privacy, speaker anonymization or de-identification is used
as a preprocessing method before the speech is publicly distributed. Research
on speaker anonymization constitutes a relatively new area that has emerged
following the efforts to suppress speaker identity to fool speaker
identification systems. In prior studies, the Kaldi phone has been used for
voice transformation [3, 4] and could successfully fool the earlier speaker
identification system based on a Gaussian mixture model (GMM). In addition,
another approach based on cepstral frequency warping has been used for speech
transformation for the purpose of de-identification [5].

Although a few solutions have been proposed for anonymization, there is no
formal definition or task in place yet, which causes ambiguity in determining
the level of anonymization performance. Hence, the VoicePrivacy 2020 Challenge
(VP2020) was initiated with the goal of fostering the development of speaker
anonymization techniques by defining the task, metrics, scenarios, and a
benchmark for the initial solutions (baseline). In this challenge, an
anonymization system suppressed PII while maintaining the intelligibility and
naturalness of a speech signal [1].Two baseline systems were provided: (1) a
primary baseline based on a state-of-the-art x-vector embedding and neural
waveform modelling [6], and (2) a secondary baseline based on the manipulation
of pole locations derived by linear predictive coding (LPC) using the McAdams
coefficient [7].

As a result of the VP2020, several systems have been proposed for speaker
anonymization [8]. Most of these systems focused on improving the primary
baseline [9, 10, 11], since it achieved superior results in comparison to the
secondary baseline. In this paper, we focus on improving the security of the
secondary baseline, which is much less complex than the primary baseline and
requires no training data. Our proposed method is based on a speech
watermarking approach. With speech watermarking, it is possible to embed a
secret message (also called a watermark) [12] that can be used to verify the
authentication of the speech. We speculate that by embedding a watermark in
anonymized speech, we should be able to identify the originality of the speech
or use it to prevent speech spoofing.

In Section 2 of this paper, we provide a brief overview of related studies.
Section 3 introduces our proposed method. Our experimental setting,
evaluations, and results are reported in Section 4. We conclude in Section 5
with a brief summary and mention of future work.

##  2 McAdams coefficient-based Speaker Anonymization

As one solution for preserving the privacy of speech data, VP2020 [1] was
initiated with the goal of fostering the development of speaker anonymization
techniques. Such techniques are ideally expected to suppress the leakage of
PII while maintaining the linguistic information of the speech signal. In this
challenge, four requirements were determined for the speaker anonymization
technique: (1) the output had to be a speech waveform, (2) it must maximize
the suppression of speaker individuality information, (3) it must preserve
speech naturalness and intelligibility, and (4) it must ensure the distinction
of voices of different speakers. This challenge also provided two different
baseline systems111https://github.com/Voice-Privacy-Challenge/Voice-Privacy-
Challenge-2020, datasets, and protocols for evaluating the speaker
anonymization performance.

The primary baseline system was developed on the basis of x-vector embeddings
and a neural source filter-based waveform model [6]. The idea behind this
system is to modify the speech identity information (x-vector) after
disentangling it from the speech content (phoneme posteriorgram and the
fundamental frequency). The secondary baseline system (shown in Fig. 1) was
developed on the basis of modifications to the McAdams coefficient [7]. The
McAdams coefficient is related to the adjustment of harmonic frequency
distributions, which affects the perception of timbre [13]. Although the
results of the secondary baseline were not as good as the first baseline, it
requires no training data and is much less complex.

Figure 1: Block diagram of speaker anonymization based on McAdams coefficient
[7]. “LP coeff.” is referred to as linear prediction coefficients. “LPC” is
referred to as linear predictive coding. “$\phi$” is the angle of poles with a
non-zero imaginary part. “$\alpha$” is the McAdams coefficient.

The McAdams coefficient proposed in [13] is a parameter derived on the basis
of the additive synthesis method in music signal processing [14]. This method
is applied to timbre generation by resynthesizing multiple harmonic
cosinusoidal oscillations, as

$y_{syn}(t)=\sum_{h=1}^{H}r_{h}(t)\cos(2\pi(hf_{0})^{\alpha}t+\Phi_{h}),$ (1)

where $y_{syn}(t)$ is the synthesized signal by combining harmonic
consinusoidal oscillations with inverse Fourier series, $h$ is the harmonic
index, $r_{h}(t)$ is the amplitude, $\Phi_{h}$ is the phase, and $\alpha$ is
the McAdams coefficient [13]. Prior work on speaker anonymization [7], has
shown that the McAdams coefficient can transform the spectral envelope of
speech signals and affect timbre perception. To examine this in more detail,
we show the spectral envelopes of the frequency response by various McAdams
coefficients in Fig. 2. We can see that the farther away the McAdams
coefficient is from the original speech ($\alpha=1$), the greater the shift of
spectral envelope we can obtain. The degree of this shifting affects our
perception in perceiving the speech formants [13].

Figure 2: Spectral envelopes of speech frame with
$\alpha\in\\{1.0,0.9,0.7,0.5\\}$ [7].

##  3 Proposed Method

In this study, we propose a method based on the secondary baseline to enhance
the security of anonymized speech by using a watermarking approach. Speech
watermarking aims to protect the security in a speech signal by imperceptibly
embedding within it a particular message, such as a signature that indicates
the speech’s ownership. Speech watermarking should fulfill at least three
requirements: inaudibility (not perceivable by the human auditory system),
blindness (detection without the availability of original signal), and
robustness against common signal processing operations. The trade-off between
inaudibility and robustness has been the most pressing issue in existing
speech watermarking techniques [12]. In the present work, we investigate the
effectiveness of using a watermarking approach for a McAdams coefficient-based
speech anonymization method with regard to these requirements.

Generally, speech watermarking consists of two main processes: embedding and
detection. Figure 3 shows the block diagram of our embedding process. As the
first step, we generated the anonymized signals from the original signal
($x(n)$) using two different McAdams coefficients ($\alpha_{0}$ and
$\alpha_{1}$). The anonymization procedure follows the steps in Fig. 1.

Figure 3: Block diagram of embedding process. $\alpha_{0}$ and $\alpha_{1}$
are the McAdams coefficients for representing binary bit “0” and “1”.
$a_{0}(n)$ and $a_{1}(n)$ are the output anonymized speech in time domain.

Firstly, the original speech was segmented into speech frames with frame
length depending on the watermarking payload. Subsequently, the speech frame
was analyzed using linear predictive coding (LPC) with an order of 20
($M=20$). The mathematical form of the LPC is characterized by the following
differential equation:

$s(n)=\sum_{i=1}^{M}c(i)s(n-i)+e(n)$ (2)

where $s(n)$ is a speech signal, $c(i)$ corresponds to the filter coefficient
in $i$-th order, $M$ is the maximum order of the prediction (in this study
$M=20$), and $e(n)$ is the prediction error. The corresponding transfer
function ($H(z)$) for Eq. (2) is represented by twentieth-order all-pole
autoregressive filters, which is given by:

$H(z)=\frac{1}{1-\sum_{i=1}^{M}c(i)z^{-i}}$ (3)

From LPC analysis, we obtained the linear prediction coefficients (LP
coefficients) and residuals. These LP coefficients ($c(i)$) were then used to
derive the pole positions. The derived poles were comprised of complex poles
(poles with non-zero-valued imaginary terms) and real poles (poles with zero-
valued imaginary terms). The shift of complex poles position ($\phi^{\alpha}$)
was resulting in the angle shifting to either clockwise or counter-clockwise
of the complex positions [7]. The effect of this angle shifting in the
spectral envelope is shown in Fig. 2. After the modification of McAdams
coefficients $\alpha$, the modified complex poles and the real poles were
converted back to LP coefficients. The anonymized speech frame was
resynthesized from these LP coefficients and the original residuals.

After obtaining the anonymized signals, we normalized them to be in similar
relative loudness (fixed a target peak level in decibel relative to full scale
(dBFS)) and range of frequency components (using a bandpass filter (BPF)). The
cut-off frequencies for the BPF were 125 Hz and 4 kHz. Finally, we constructed
the watermarked signal ($y(n)$) by frame-by-frame concatenation of the
anonymized signals obtained by bit inverse according to watermarked bit-
stream.

We found that the anonymized signals from different McAdams coefficients
carried different amounts of power spectrum, specifically in the lower
frequency components. Using a higher McAdams coefficient results in a higher
power spectrum in the lower frequency. On the basis of this characteristic, we
determined a power threshold for the blind detection process (shown in Fig.
4). The detection process was conducted by comparing the power spectrum
obtained by fast Fourier transform (FFT) of the watermarked signal
($|Y(\omega)|$) in a specific frequency range with the designated threshold
$\theta$.

Figure 4: Block diagram of blind detection process. “BPF” stands for the band-
pass filtering. “FFT” stands for fast Fourier transform. $|Y(\omega)|$ is the
power spectrum of the watermarked signal $y(n)$ obtained by FFT. $\theta$ is
the power spectrum threshold for blind detection process. $w^{\prime}(k)$ is
the detected watermark bit of the $k$-th frame.

##  4 Experiments

This section describes the experimental setting in our study and reports the
objective evaluation results.

###  4.1 Experimental Setting

We used LibriSpeech [15] and VCTK [16] datasets (development and testing sets)
that were provided in VP2020. LibriSpeech is an English speech corpus designed
for automatic speech recognition (ASR) research sampled at 16 kHz [15]. VCTK
is an English speech corpus that contains 109 native speakers with various
accents and was designed for text-to-speech (TTS) research sampled at 48 kHz.
The development and training data of both datasets in VP2020 consists of more
than 20,000 utterances from almost 200 speakers. The sampling rate of the
speech data is set to 16 kHz. Similar to the secondary baseline [7], we do not
need any training data for our proposed method. The McAdams coefficient used
to represent bit “0” was 0.6 ($\alpha_{0}=0.6$) and bit “1” ($\alpha_{1}=0.8$)
was 0.8. We evaluated the speaker verifiability of our proposed method by
using a pretrained automatic speaker verification system (ASVeval) and the
intelligibility by using a pretrained automatic speech recognition system
(ASReval), similar to the protocol in VP2020. The suggested requirement for a
reasonable payload of an audio watermarking is 72 bits per 30 seconds (around
2–3 bps) [17]. In this study, we used a 4-bps watermarked signal to evaluate
the anonymization performance.

(a) Female utterances
(b) Male utterances

Figure 5: Speaker verifiability evaluation in equal-error rate (EER). “f”
stands for female utterances, whereas “m” stands for male utterances. There
are three scenarios considered, i.e., (1) o-o (both original data for
enrollment and trial), (2) o-a (original data for enrollment and anonymized
data for trial), and (3) a-a (both anonymized data for enrollment and trial).

For speech watermarking, we also evaluated the speech quality and robustness
of our method with a total 100 randomly selected utterances from the
LibriSpeech and VCTK datasets. Since the original signal was not available, we
used MOSNet, the pretrained mean opinion score (MOS) predictor proposed in
[18]. MOSNet is an objective evaluation tool based on deep learning approach
for predicting human MOS ratings in a voice conversion system. Subsequently,
we evaluated the robustness of our propsed method as suggested in [19] by
calculating the bit error rate (BER) and balanced F1-score during normal (no
attack) operations along with several signal processing operations, including
noise addition, resampling, requantization, compression, and speech codecs. We
also examined the security level by calculating the false acceptance rate
(FAR) and false rejection rate (FRR). The maximum acceptable BER threshold as
the robustness indication is $10\%$ [19]. We embedded random binary streams
with payloads of 2, 4, 8, 16, and 32 bps and varied the detection threshold in
the order of lower to higher payloads (0.15, 0.09, 0.05, 0.02, and 0.01,
respectively). Due to the space limitation, the results reported in this paper
are mainly in mean value.

###  4.2 Results

Our intent with this evaluation is mainly to investigate the effectiveness and
reliability of the proposed method in anonymizing the PII of the speaker and
watermarking the speech. For the speaker anonymization, we conducted our
evaluation using ASVeval to check the speaker verifiability performance in
three cases: original enrolls and original trials case (o-o), original enrolls
and anonymized trials case (o-a), and anonymized enrolls and anonymized trials
case (a-a). The results are shown in Fig. 5, where the top and bottom graphs
shows the results of female and male utterances, respectively. As we can see,
the proposed method (P) improved the speaker similarity in comparison to the
secondary baseline in VP2020 (B2): the EER improved significantly to more than
35% in the o-a case and to more than 20% in the a-a case for both female and
male utterances.

Table 1: MOSNet evaluation results. | payload (bps) | MOS
---|---|---
original | - | 3.15 ± 0.49
anonymized | - | 2.70 ± 0.18
| 2 | 2.73 ± 0.20
| 4 | 2.73 ± 0.21
| 8 | 2.70 ± 0.19
| 16 | 2.67 ± 0.18
proposed method | 32 | 2.60 ± 0.18
Figure 6: Speech intelligibility evaluation in word error rate (WER). “ASReval
ori” and “ASReval anon” denote ASR trained on original and anonymized data,
respectively. “Libri” is referred to as LibriSpeech dataset. “VCTK” is
referred to as VCTK dataset. “ori”, “B2”, “P” are original speech, the output
of the secondary baseline in VP2020, the output of our proposed method,
respectively. Figure 7: Robustness test results in terms of BER (bit error
rate) in nine cases: (a) normal, (b) AWGN, (c) resample-8, (d) resample-24,
(e) requant-8, (f) requant-24, (g) mp3, (h) flv, and (i) G723.1. Figure 8:
Robustness test results in nine cases: (a) normal, (b) AWGN, (c) resample-8,
(d) resample-24, (e) requant-8, (f) requant-24, (g) mp3, (h) flv, and (i)
G723.1. For metrics were used for the robusness evaluation, including F1
(F1-score), FAR (false acceptance rate), and FRR (false rejection rate).

For the speech intelligibility, we conducted an objective evaluation by using
ASReval trained on original data and anonymized data. The results are shown in
Fig. 6, where the blue and orange bars indicate the results obtained from
ASReval trained on original and anonymized data, respectively. In contrast to
the speech verifiability, our proposed method caused a noticeable degradation
to intelligibility, especially for the VCTK dataset (the degradation is
increased from 11.79% to 54.78%). This degradation was mainly occurred due to
the anonymization method as we can obtain almost similar average WER with the
stochastic approach proposed in [7]. The further shift to McAdams
coefficients, the more WER occurred. To improve the intelligibility, we could
retrain the ASReval using the anonymized data as in [7]. For example, the WER
is shown by the blue bar for “VCTK (P)” could be reduced to 11.91%.

Inaudibility and robustness are essential requirements in speech watermarking.
To evaluate the inaudibility of the proposed method, we used MOSNet [18] to
derive the MOS of the watermarked signal. Table 1 shows the MOSNet results of
original signal, anonymized signal with McAdams coefficients $\alpha=0.8$, and
the output signal of our proposed method with various payloads. We can see
that there was a speech quality degradation (MOS degraded from 3.15 to 2.70)
caused by the McAdams coefficient-based anonymization method, while in
contrast, the proposed method could maintain a similar MOS even with a
relatively high payload.

Figure 9: Detected watermarks from image embedding in 4-bps watermarked signal
after several operations: (a) original, (b) normal, (c) resample-8, (d)
resample-24, (e) requant-8, (f) requant-24, (g) AWGN, (h) MP3, and (i) G723.1.

We carried out a robustness test by calculating the detection accuracy from
the output speech after several common signal processing operations. Figures 7
and 8 show the robustness test results. We examined nine cases: no attack
(normal), addition of white Gaussian noise (AWGN), downsampling to 8 kHz
(resample-8), upsampling to 24 kHz (resample-24), bit compression to 8 bits
(requantize-8), bit extension to 24 bits (requantize-24), MP3 compression
(MP3), flash video format (flv), and G723.1 codec. For AWGN processing, the
signal to noise ratio used is 40 decibel (dB). Meanwhile, the bitrate range
for MP3 compression was from 220-260 kbps (kilo bits per second). The bitrate
of G723.1 codec was 5.3 kbps with algebraic code-excited linear prediction
(ACELP) algorithm. As we can see, our proposed method was robust against other
operations (the BER was similar to the normal case), except for the conversion
to video codec (flv). The results here demonstrate that our method is suitable
for watermarking purposes, since the BER for 4 bps satisfied the robustness
criteria ($\mathrm{BER}<10\%$). The results also suggest that the security
level indicated by the FAR, FRR, and F1-score is adequate for payloads up to 4
bps.

In summary, our proposed method could improve significantly the speaker
verifiability by ASVeval. Although intelligibility was degraded, the WER could
reduced by simply retraining the ASReval. The evaluation for speech
watermarking also confirmed that our proposed method could resulted in
inaudible and robust watermark. As an example, Figure 9 shows the example of
original image watermarked in the speech with the detected image after several
operations. We could easily perceive the 2021 logo in the detected watermarks
in almost all processing (BER is approximately 7%).

##  5 Conclusion and Future Work

In this paper, we proposed a technique to improve the security of McAdams
coefficient-based speaker anonymization by using a watermarking approach. By
adding a watermark into the anonymized signal, we could ensure the originality
of the speech (one of the anti-spoofing countermeasures). The performance of
the proposed method was evaluated objectively based on the standard in speaker
anonymization (VP2020 protocol) and watermarking. The results showed that our
method is suitable for watermarking, as it could satisfy the blind, inaudible,
and robustness requirements. In addition, it significantly improved the
speaker dissimilarity in ASVeval. The limitation of our method in comparison
to the baseline system is that it caused degradation in speech
intelligibility, although this could be resolved by retraining the ASReval
with the anonymized speech.

In spite of these promising results, we acknowledge that the evaluation,
particularly in consideration to human perception such as the speech quality
test, might not be sufficient. As future work, we will explore more suitable
ways to evaluate the watermarking technique for speaker anonymization. We also
intend to optimize the proposed method, as its detection rate and payload are
relatively limited.

##  6 Acknowledgements

This work was supported by a Grant-in-Aid for Scientific Research (B) (No.
17H01761), JSPS KAKENHI Grant (No.20J20580), Fund for the Promotion of Joint
International Research (Fostering Joint International Research (B))(20KK0233),
and KDDI foundation (Research Grant Program).

## References

  * [1] N. Tomashenko, B. M. L. Srivastava, X. Wang, E. Vincent, A. Nautsch, J. Yamagishi, N. Evans, J. Patino, J.-F. Bonastre, P.-G. Noé, and M. Todisco, “The VoicePrivacy 2020 Challenge evaluation plan,” 2020.
  * [2] S. Ö. Arik, J. Chen, K. Peng, W. Ping, and Y. Zhou, “Neural voice cloning with a few samples,” _CoRR_ , vol. abs/1802.06006, 2018. [Online]. Available: http://arxiv.org/abs/1802.06006
  * [3] Q. Jin, A. R. Toth, A. W. Black, and T. Schultz, “Is voice transformation a threat to speaker identification?” in _2008 IEEE International Conference on Acoustics, Speech and Signal Processing_. IEEE, 2008, pp. 4845–4848.
  * [4] Q. Jin, A. R. Toth, T. Schultz, and A. W. Black, “Speaker de-identification via voice transformation,” in _2009 IEEE Workshop on Automatic Speech Recognition & Understanding_. IEEE, 2009, pp. 529–533.
  * [5] C. Magariños, P. Lopez-Otero, L. Docio-Fernandez, E. Rodriguez-Banga, D. Erro, and C. Garcia-Mateo, “Reversible speaker de-identification using pre-trained transformation functions,” _Computer Speech & Language_, vol. 46, pp. 36–52, 2017.
  * [6] F. Fang, X. Z. Wang, J. Yamagishi, I. Echizen, M. Todisco, N. Evans, and J.-F. Bonastre, “Speaker anonymization using x-vector and neural waveform models,” 2019.
  * [7] J. Patino, N. Tomashenko, M. Todisco, A. Nautsch, and N. Evans, “Speaker anonymisation using the McAdams coefficient,” _arXiv preprint arXiv:2011.01130_ , 2020.
  * [8] “The VoicePrivacy 2020 Challenge, challenge setup and results,” 2020. [Online]. Available: https://www.voiceprivacychallenge.org/docs/1___VoicePrivacy_
challenge_setup_and_results_N_Tomashenko.pdf

  * [9] C. O. Mawalim, K. Galajit, J. Karnjana, and M. Unoki, “X-vector singular value modification and statistical-based decomposition with ensemble regression modeling for speaker anonymization system,” in _Interspeech 2020_. ISCA, 2020, pp. 1703–1707. [Online]. Available: https://doi.org/10.21437/Interspeech.2020-1887
  * [10] H. Turner, G. Lovisotto, and I. Martinovic, “Speaker anonymization with distribution-preserving x-vector generation for the VoicePrivacy challenge 2020.” [Online]. Available: https://www.voiceprivacychallenge.org/docs/Oxford.pdf
  * [11] F. M. Ezpinoza-Cuadros, J. M. Perero-Codosero, J. Anton-Martin, and L. A. Hernandez-Gomez, “Speaker de-identification system using autoencoders and adversarial training.” [Online]. Available: https://www.voiceprivacychallenge.org/docs/Sigma.pdf
  * [12] G. Hua, J. Huang, Y. Q. Shi, J. Goh, and V. L. Thing, “Twenty years of digital audio watermarking—a comprehensive review,” _Signal Processing_ , vol. 128, pp. 222–242, 2016. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0165168416300263
  * [13] S. Mcadams, “Spectral fusion, spectral parsing and the formation of auditory images,” _Ph. D. Thesis, Stanford_ , 1984.
  * [14] C. Dodge and T. A. Jerse, “Computer music: Synthesis, composition, and performance,” 1997.
  * [15] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, “Librispeech: An ASR corpus based on public domain audio books,” in _2015 IEEE ICASSP_ , 2015, pp. 5206–5210.
  * [16] C. Veaux, J. Yamagishi, and K. Macdonald, “CSTR VCTK corpus: English multi-speaker corpus for CSTR voice cloning toolkit,” 2017.
  * [17] STEP2001, “News release, final selection of technology toward the global spread of digital audio watermarks,” 2001. [Online]. Available: http://www.jasrac.or.jp/ejhp/release/2001/0629.html
  * [18] C.-C. Lo, S.-W. Fu, W.-C. Huang, X. Wang, J. Yamagishi, Y. Tsao, and H.-M. Wang, “MOSNet: Deep learning based objective assessment for voice conversion,” 2019.
  * [19] IHC Committee, “IHC evaluation criteria and competition.” [Online]. Available: https://www.ieice.org/iss/emm/ihc/IHC_criteriaVer6.pdf


# Exactly energy-conserving electromagnetic Particle-in-Cell method in
curvilinear coordinates

J. Croonen Centre for mathematical Plasma Astrophysics, Department of
Mathematics, KU Leuven, Celestijnenlaan 200B, B-3001 Leuven, Belgium  L.
Pezzini Centre for mathematical Plasma Astrophysics, Department of
Mathematics, KU Leuven, Celestijnenlaan 200B, B-3001 Leuven, Belgium  Royal
Observatory of Belgium, Solar-Terrestrial Centre of Excellence, Ringlaan 3,
1180 Uccle, Belgium  F. Bacchini Centre for mathematical Plasma Astrophysics,
Department of Mathematics, KU Leuven, Celestijnenlaan 200B, B-3001 Leuven,
Belgium  Royal Belgian Institute for Space Aeronomy, Solar-Terrestrial Centre
of Excellence, Ringlaan 3, 1180 Uccle, Belgium  G.Lapenta Centre for
mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven,
Celestijnenlaan 200B, B-3001 Leuven, Belgium

###### Abstract

In this paper, we introduce and discuss an exactly energy-conserving Particle-
in-Cell method for arbitrary curvilinear coordinates. The flexibility provided
by curvilinear coordinates enables the study of plasmas in complex-shaped
domains by aligning the grid to the given geometry, or by focusing grid
resolution on regions of interest without overresolving the surrounding,
potentially uninteresting domain. We have achieved this through the
introduction of the metric tensor, the Jacobian matrix, and contravariant
operators combined with an energy-conserving fully implicit solver. We
demonstrate the method’s capabilities using a Python implementation to study
several one- and two-dimensional test cases: the electrostatic two-stream
instability, the electromagnetic Weibel instability, and the geomagnetic
environment modeling (GEM) reconnection challenge. The test results confirm
the capability of our new method to reproduce theoretical expectations (e.g.
instability growth rates) and the corresponding results obtained with a
Cartesian uniform grid when using curvilinear grids. Simultaneously, we show
that the method conserves energy to machine precision in all cases.

Computational Methods (1965) — Plasma Physics (2089) — Plasma Astrophysics
(1261)

††journal: ApJS

##  1 Introduction

Advanced plasma simulation methods have become a key tool for studying a wide
variety of systems and phenomena throughout the fields of plasma physics and
astrophysics. The difficulty of recreating the relevant plasma conditions in a
laboratory environment make simulations irreplaceable for obtaining otherwise
unachievable key insight. Among plasma simulations techniques, Particle-in-
Cell (PIC) methods stand out as particularly useful when studying plasmas
where kinetic effects are important. The state of the art of PIC methods is
constantly evolving with new capabilities and improvements being developed.
Recently, the astro- and plasma communities have steered their interest toward
more advanced PIC methods, especially involving flexible (i.e. adaptive) grids
and high-stability algorithms, to allow for long-term simulations of systems
where Cartesian (uniform) grids might be ill-suited (for example Tokamak
simulations using magnetic coordinates as shown in e.g. Jolliet et al. 2007,
or accretion disks around astrophysical compact objects using spherical Kerr-
Schild coordinates as shown in e.g. Crinquand et al. 2022). A nonuniform,
adaptive grid could be particularly beneficial in several cases, e.g. i) when
the physical shape of the domain of interest is nontrivial, such that a
significant part of an encompassing Cartesian box would be occupied by empty
space or otherwise uninteresting regions where plasma behavior is trivial or
quasi-static; and ii) when a small region of interest is embedded in a much
larger domain, and the latter might be significantly overresolved (in terms of
spatial and temporal resolution) due to physical requirements only imposed by
the small region of interests. In both cases, employing a standard uniform
grid can substantially increase the computational cost of the simulation, in
the worst case making certain studies entirely unfeasible. The introduction of
a nonuniform grid is one way to relax these limitations, by adapting the grid
to the physical setup in question and thus assigning the computational
resources more efficiently (see e.g. Fichtl et al. 2012; Delzanno et al. 2013;
Chacón & Chen 2016; Stanier & Chacón 2022).

Several methods have been developed to use curvilinear grids, but most are
implemented with one specific grid in mind, and thus lack flexibility (see
e.g. Ringle 2011; Gonzalez-Herrero et al. 2019). Others are fully general,
i.e. they work with an arbitrary grid, but come with other trade-offs such as:
i) using a reduced set of equations, rather than a full electromagnetic
implementation, thus limiting their applicability to the constraints of these
approximations (e.g. an electrostatic model in Fichtl et al. 2012; Delzanno et
al. 2013, or the Vlasov-Darwin approximation in Chacón & Chen 2016); or ii)
introducing large deviations from energy conservation in long simulations, due
to discretization errors that are intrinsic of the PIC implementation with
explicit methods. Energy conservation is a universal physical property, of
particular interest in setups where energy is converted from one form to
another. For example, when modeling instabilities or reconnection events,
preserving energy is extremely important to avoid numerical heating or cooling
which could significantly influence the plasma dynamics and potentially
invalidate the results (e.g. Markidis & Lapenta 2011a).

In this work, we provide a method combining general curvilinear grids, through
the introduction of coordinate transformations, with a fully electromagnetic-
PIC implementation using an exactly energy-conserving, Jacobian-free Newton-
Krylov solver. The method was designed with the architecture of the ECsim code
(Lapenta 2017) in mind, to allow for easy implementation in a production-ready
infrastructure. To the best of our knowledge, this may represent the first
production-scale code which combines the aforementioned properties
(curvilinear grids and exact energy conservation).

This paper is organized as follows: the mathematical foundation of the method
is described in detail in Section 2; in Section 3, a series of one- and two-
dimensional tests are discussed to assess the validity of the new method; the
main results, conclusions, potential applications, and next steps in this line
of research are discussed in the final Section 4.

##  2 The Curvilinear PIC Method

The method described in this paper combines a fully implicit PIC
implementation (Chen et al. 2011; Markidis & Lapenta 2011b) using finite
differences and a Jacobian-free Newton-Krylov solver, with a curvilinear grid
through the introduction of the metric tensor, Jacobian matrix, and covariant
differentiation operators. As we will show in Section 3, the method was
successfully implemented and tested in Python.

###  2.1 Governing equations

The PIC method is a first-principles approach for fully kinetic plasma
simulations. It describes the plasma as a coupled system between freely moving
particles existing in a position-velocity phase space, and a discrete
representation of electric and magnetic fields on a computational grid.

To evolve the fields in time we use Maxwell’s equations, here in CGS units:

$\mbox{\boldmath{$\nabla$}}\cdot\boldsymbol{E}=4\pi\rho,$ (1)
$\mbox{\boldmath{$\nabla$}}\cdot\boldsymbol{B}=0,$ (2)
$\frac{1}{c}\frac{\partial\boldsymbol{B}}{\partial
t}=-\mbox{\boldmath{$\nabla$}}\times\boldsymbol{E},$ (3)
$\frac{1}{c}\frac{\partial\boldsymbol{E}}{\partial
t}=-\frac{4\pi}{c}\boldsymbol{I}+\mbox{\boldmath{$\nabla$}}\times\boldsymbol{B},$
(4)

where $\boldsymbol{E}$ and $\boldsymbol{B}$ are the electric and magnetic
field respectively, $\rho$ is the charge density, $c$ is the speed of light,
and $\boldsymbol{I}$ is the current density, which is computed from the
particle velocities (see Section 2.3 below). Note that the PIC method only
requires Faraday’s law (3) and Ampère’s law (4) to form a closed system of
equations. It can easily be shown that Gauss’s law for magnetism (2) will
always be satisfied if the spatial discretization on the computational grid
mimics the continuous vector identity
$\mbox{\boldmath{$\nabla$}}\cdot\mbox{\boldmath{$\nabla$}}\times\mbox{\boldmath{$V$}}=0$
for a generic vector $V$ (Lapenta, 2017). Gauss’s law for $E$ (1), however, is
not automatically satisfied in fully implicit PIC. While it is possible to
construct a method that also satisfies this equation (e.g. Chen et al. 2011;
Chen & Tóth 2019), we are currently only concerned with presenting a first
implementation of our method, which conserves energy exactly in curvilinear
coordinates, and that can than be further improved upon in future work.
Moreover, in our tests we did not detect numerical artifacts linked to charge
conservation, implying that numerical errors are limited.

Computational particles in the PIC method are updated using the standard
equations of motion,

$\frac{\mathrm{d}\boldsymbol{x}}{\mathrm{d}t}=\boldsymbol{v},$ (5)
$\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t}=\frac{q}{m}\left(\boldsymbol{E}+\frac{1}{c}\boldsymbol{v}\times\boldsymbol{B}\right),$
(6)

where $\boldsymbol{x}$ and $\boldsymbol{v}$ are the particle position and
velocity, $q$ the particle charge, and $m$ the particle mass.

In the PIC paradigm, the sources for Maxwell’s equations ($\boldsymbol{I}$ and
$\rho$) are collected from the particles onto the computational grid via
interpolation; electromagnetic fields are interpolated at the particle
positions to obtain the Lorentz force needed to update the particles. Because
of this, particle and field equations are in general nonlinearly coupled. In
the method described in this paper a fully implicit nonlinear iterative solver
is used to solve this set of equations, which are discretized as described in
Section 2.3.

###  2.2 Curvilinear PIC: Coordinate transformation

Figure 1:  Diagram of the coordinate transformation between physical and
logical space. The Jacobian and its inverse are used to transform
contravariant vector components. In physical space the Euclidean metric is
used, while a nontrivial metric is used in logical space.

To work with a curved grid we consider an arbitrary mapping

$\boldsymbol{\xi}=f(\boldsymbol{x}).$ (7)

The space described by the new coordinates $\boldsymbol{\xi}$ will be called
the “logical” space and the space described by the original set of coordinates
$\boldsymbol{x}$ will be called the “physical” space. The mapping $f$ is
chosen such that the curvilinear grid in physical space becomes a regular
orthogonal grid in logical space. This is depicted in figure 1. Next, we
introduce the Jacobian and its inverse,

$j_{\mu}^{i}=\frac{\partial x_{i}}{\partial\xi_{\mu}},$ (8)
${j^{-1}}_{i}^{\mu}=\frac{\partial\xi_{\mu}}{\partial x_{i}},$ (9)

and the associated metric defined as

$g_{\mu\nu}=j^{i}_{\mu}j^{i}_{\nu}.$ (10)

Here, Greek indices are used to denote $\boldsymbol{\xi}$ coordinates in the
logical space and Latin indices to denote $\boldsymbol{x}$ coordinates in
physical space. In this way, the difficulty of working with a nonuniform
curvilinear grid has been circumvented by transforming to the logical grid, at
the cost of a non-trivial metric tensor and Jacobian matrix. Calculating the
metric and Jacobian only has to be done once at initialization for each grid
point and can subsequently be stored in memory. Note that, in our
implementation, the grid is static and therefore the Jacobian and metric are
time-invariant.

Due to the introduction of a nontrivial metric the vector operators have to be
changed from their usual Cartesian implementation to a covariant
implementation. The only differential operator appearing in our PIC equations
is the curl operator, e.g. for a generic (covariant) vector field
$V_{\kappa}$,

$(\mbox{\boldmath{$\nabla$}}\times\mbox{\boldmath{$V$}})^{\mu}\rightarrow\frac{1}{J}\epsilon^{\mu\nu\kappa}\partial_{\nu}V_{\kappa},$
(11)

where $J=|j|$ is the Jacobian determinant, and $\epsilon^{\mu\nu\kappa}$ is
the Levi-Civita symbol. Furthermore, the metric tensor can be used to raise
and lower indices, e.g. $V_{\kappa}=g_{\kappa\lambda}V^{\lambda}$.

For the particle equations of motion we keep the Cartesian description from
eqs. (5)–(6). We avoid solving these equations in logical space, since this
would require the Jacobian and metric tensor to be known at each particle
position at every timestep. Given the computational complexity of deriving
these values from the initial mapping function, this was not considered a
practical operation to be carried out at every iteration for every particle.
This operation is in principle feasible if all mapping quantities are
analytically known, but this imposes additional constraints on the mapping
function (i.e. invertibility and differentiability) which we want to avoid for
generality. Calculating the equations of motion in Cartesian coordinates
instead introduces the need to convert particle positions and vector-field
components between physical and logical space. This strategy is still less
computationally demanding than evolving particles in the logical space.
Particle positions are transformed from Cartesian physical space to general
logical space simply using the original mapping function $f$ in eq. (7).
Contravariant vector components can be transformed using the Jacobian and
inverse Jacobian matrices as

$V^{i}=j_{\mu}^{i}V^{\mu},$ (12) $V^{\mu}={j^{-1}}_{i}^{\mu}V^{i}.$ (13)

This is necessary in the interpolation of the current density and the fields,
as we show in the following Section. For further details or interest in the
topic we recommend the book Liseikin (1999).

###  2.3 Discretization

The notation in this paper uses superscript parentheses to denote discrete
time steps, and subscript parentheses to denote discrete spatial coordinates.
All quantities are known at integer time steps; when the method requires half-
integer times (see below), the relevant quantities will be linearly
interpolated between subsequent time levels. Field quantities are defined at
discrete locations on the computational grid, e.g. in two spatial dimensions,
$V_{(\xi,\eta)}^{(n)}$ denotes a quantity $V$ at the grid location
$(\xi,\eta)$ in the logical space and discrete time level $n$. This allows us
to avoid confusion with indices of vector components. When the exact spatial
indices are not relevant they are replaced by subscript $g$ to indicate a
quantity known at a generic grid point, or subscript $p$ for a quantity known
at a particle’s position.

To discretize the governing equations, we use a finite-difference scheme.
Magnetic-field components are placed on grid nodes, with integer spatial
indices, while electric field components are placed on the grid centers, with
half-integer spatial indices111This choice conforms to the “colocated” grid
discretization that has been used for the iPic3D and ECSIM family of codes,
see e.g. Markidis et al. (2010); Lapenta (2017).. Maxwell’s equations (3)–(4)
require spatial derivatives of the magnetic field to be colocated with the
electric field, and conversely the spatial derivatives of the electric field
to be colocated with the magnetic field. This requires the quantities to be
averaged along the direction perpendicular to that of the derivative before
taking finite differences. This allows for a centered finite-difference
spatial derivative which is second-order accurate; in two dimensions, the
discrete derivatives for a generic quantity $V$ along the $\xi$ or $\eta$
directions are

$D_{\xi}V=\frac{(V_{(\xi+1,\eta)}+V_{(\xi+1,\eta+1)})-(V_{(\xi,\eta)}+V_{(\xi,\eta+1)})}{2\Delta\xi},$
(14)
$D_{\eta}V=\frac{(V_{(\xi,\eta+1)}+V_{(\xi+1,\eta+1)})-(V_{(\xi,\eta)}+V_{(\xi+1,\eta)})}{2\Delta\eta},$
(15)

where the spatial indices can be integers or half-integers depending on
whether the specific quantity exists on nodes or centers respectively. Time
derivatives are discretized using a second-order accurate central-difference
scheme. Using these discrete derivatives, we can now write the discrete
Maxwell’s equations (3)–(4):

$\frac{B^{\mu}{}_{g}^{(n+1)}-B^{\mu}{}_{g}^{(n)}}{\Delta
t}=-\frac{c}{J_{g}}\epsilon^{\mu\nu\kappa}(D_{\nu}E_{\kappa})_{g}^{(n+1/2)},$
(16) $\begin{split}\frac{E^{\mu}{}_{g}^{(n+1)}-E^{\mu}{}_{g}^{(n)}}{\Delta
t}=&-4\pi I^{\mu}{}_{g}^{(n+1/2)}\\\
&+\frac{c}{J_{g}}\epsilon^{\mu\nu\kappa}(D_{\nu}B_{\kappa})_{g}^{(n+1/2)}.\end{split}$
(17)

The current density on each grid point must be gathered from the particles as

$I^{\mu}{}_{g}=\frac{({j^{-1}}_{i}^{\mu})_{g}}{J_{g}\Delta\xi\Delta\eta}\sum_{p}w_{pg}q_{p}v^{i}{}_{p},$
(18)

where $w_{pg}$ is the interpolation function linking particle $p$ with the
grid point $g$. $w_{pg}$ is a function of the particle position in logical
coordinates; this can be found using the mapping function
$\boldsymbol{\xi}=f(\boldsymbol{x})$. The inverse-Jacobian factor in the
equation above transforms the current from physical to logical components as
described by eq. (13).

Since particles are evolved in the physical space the discretization is
straightforward,

$\frac{x^{i}{}_{p}^{(n+1)}-x^{i}{}_{p}^{(n)}}{\Delta
t}=v^{i}{}_{p}^{(n+1/2)},$ (19)
$\frac{v^{i}{}_{p}^{(n+1)}-v^{i}{}_{p}^{(n)}}{\Delta
t}=\frac{q_{p}}{m_{p}}(E^{i}+\epsilon^{ijk}v_{j}B_{k})^{(n+1/2)}_{p}.$ (20)

Note that the fields in these equations are located at the particle positions,
which requires an interpolation: for a generic field $V$,

$V^{i}{}_{p}=\sum_{g}w_{gp}(j_{\mu}^{i}V^{\mu})_{g},$ (21)

where $w_{gp}$ is the interpolation function from grid points to particle
positions. By construction we have chosen $w_{gp}=w_{pg}$, which is a
requirement for energy conservation (see Appendix A). As described by eq.
(12), we added the Jacobian $j_{\mu}^{i}$ to convert $V^{\mu}$ from
curvilinear to Cartesian coordinates, as required by the Cartesian equations
of motion.

To solve for equations (16)–(17) and (19)–(20), a nonlinear Newton-Krylov
iterative solver is used. Convergence of this solver is typically case-based,
and is not guaranteed when spatiotemporal resolution is insufficient to
capture relevant physical phenomena. As a rule of thumb for the tests in the
next Sections, the temporal and spatial resolutions were constrained to
satisfy $c\Delta t/\Delta x<1$, with a typical value for this ratio usually
set to $\sim 0.25$. This was sufficient in all test cases to achieve absolute
and relative errors below $10^{-14}$, typically within $\sim$10 iterations.
Note that $\Delta x$ may vary throughout our nonuniform grids, in which case
the smallest value of $\Delta x$ in the grid must be considered for this
constraint.

###  2.4 Discrete energy

To verify the exactly energy-conserving nature of this method, a proper
definition for the relevant energies must be established. The electric- and
magnetic-field energies are respectively defined as

$U_{\mathrm{E}}=\sum_{g}\frac{\Delta\xi\Delta\eta}{8\pi}(JE_{\mu}E^{\mu})_{g},$
(22)
$U_{\mathrm{B}}=\sum_{g}\frac{\Delta\xi\Delta\eta}{8\pi}(JB_{\mu}B^{\mu})_{g},$
(23)

which combined give the total electromagnetic-field energy
$U_{\mathrm{F}}=U_{\mathrm{E}}+U_{\mathrm{B}}$. The particle kinetic energy is
defined as

$U_{\mathrm{K}}=\sum_{p}\frac{m(v_{i}v^{i})_{p}}{2}.$ (24)

The sum of the kinetic and electromagnetic field energy is the total energy in
the system: $U=U_{\mathrm{K}}+U_{\mathrm{F}}$. In Appendix A we show
analytically that, using these definitions, the total energy will be conserved
exactly (i.e. to machine precision).

##  3 Test cases

To validate the methods developed in this paper, we considered three different
test cases: the one-dimensional, electrostatic two-stream instability; the
one-dimensional, electromagnetic Weibel (or filamentation) instability
(Weibel, 1959; Bret et al., 2010); and the two-dimensional Geospace
Environmental Modeling (GEM) challenge, i.e. a paradigmatic magnetic-
reconnection setup (Birn et al., 2001).

###  3.1 Two-stream Instability

The two-stream instability represents a classical one-dimensional test for PIC
codes. The initial setup consists of two counterstreaming electron beams with
uniform number density $n_{0}/2$. We initialize the particles in each beam
according to a drifting Maxwellian with thermal speed
$v_{\mathrm{th}}/c=0.001$. The drift-velocity component is along the only
spatial direction of the system $x$, and is added to both beams with opposite
sign, i.e. $v_{\mathrm{d}}/c=\pm 0.2$. The analytic solution of the dispersion
relation predicts that the system is unstable if
$kv_{\mathrm{d}}/\omega_{\mathrm{p}}<1$ (Goldston & H., 1995), with $k$ the
wave number and $\omega_{\mathrm{p}}$ the plasma frequency. This results in an
unstable, purely electrostatic evolution where charge bunching will generate a
self-reinforcing electric-field perturbation. During the linear instability
phase, the electric-field energy (eq. (22)) in the fastest-growing unstable
mode will grow as

$U_{\mathrm{E}}\sim\exp\left(\omega_{\mathrm{p}}t\right).$ (25)

To test the capability of our code to handle curvilinear coordinates, we
employed two nonuniform, one-dimensional grids alongside a standard Cartesian
grid used as a reference. The first curvilinear grid has a small sinusoidal
perturbation of one period along the domain, with mapping function

$\xi(x)=x+\varepsilon L_{x}\sin\left(\frac{2\pi x}{L_{x}}\right),$ (26)

where $L_{x}$ the length of the domain and $\varepsilon$ parametrizes the
strength of the perturbation, which we set to $\varepsilon=0.02$. The second
grid has a perturbation generated by a hyperbolic tangent, producing a grid
with a $\operatorname{sech}^{2}$ profile of cell density. A high density peak
is present near the center of the domain; the mapping function is

$\xi(x)=\frac{L_{x}-2\varepsilon L_{x}}{L_{x}}x-\varepsilon
L_{x}\tanh\left[\frac{1}{w}\left(\frac{L_{x}}{2}-x\right)\right]+\varepsilon
L_{x},$ (27)

where $w$ determines the width of the high-density region, which we set to $4$
with $\varepsilon=0.04$. These two perturbed grids are constructed to ensure
the total system length $L_{x}$ is maintained. This is necessary to ensure a
fair comparison to the Cartesian reference grid. The parameters used in the
sinusoidal and hyperbolic grids are chosen such that the ratio of the cell-
width between the smallest and largest cells is $\sim 0.75$. This provides a
meaningful deviation from the Cartesian case, while still allowing the same
time steps without breaking the constraint $c\Delta t/\Delta x<1$ mentioned in
Section 2.3. Figure 2 depicts the three grids used in this test for a
resolution of 64 cells. The gray-scale is an indication for local cell
density, with darker colors indicating higher density. The grids used in the
simulations had $N_{x}=2048$ cells with a domain length of
$32c/\omega_{\mathrm{p}}$, time step $\Delta
t=0.00390625\omega_{\mathrm{p}}^{-1}$ and $ppc=144$ particles per cell.

Figure 2: The three one-dimensional grids used in the study of the two-stream
and Weibel instabilities: from top to bottom, Cartesian (i.e. the usual grid
employed in standard PIC), sinusoidal, and hyperbolic-tangent mapping
functions. They are shown here for a 64-cell grid with the grays-cale
indicating the local grid density. Figure 3: Top left: Evolution of the
electric-field energy for the different grids as well as the theoretical
growth rate associated with the linear-growth phase of the two-stream
instability. Bottom left: Evolution of the relative error in the total energy
for the different grids. Right: Snapshots of the $(x,v_{x})$ phase space
during the nonlinear stage for the different grids.

Our results for this test are shown in Figure 3. The top left panel shows the
evolution in time of the electric-field energy for the three distinct grids in
one plot. It shows good agreement between the three cases as well as close
adherence to the theoretical value of the growth rate outlined in eq. (25)
during the linear-growth phase of the instability. The bottom left panel plots
the time history of the relative error in the total energy, which confirms
energy conservation to machine precision for all three cases. Lastly, the
three panels on the right show a snapshot of the position-velocity phase space
$(x,v_{x})$, where we observe the formation of holes that are typical for the
nonlinear phase of this instability. These results give a good indication that
the method is capable of faithfully reproducing electrostatic effects in a
nonuniform, one-dimensional setup while maintaining exact energy conservation.

###  3.2 Weibel Instability

To test our method in a more general electromagnetic case, we study the Weibel
(i.e. filamentation) instability. The setup is almost equivalent to the two-
stream instability: $N_{x}=2048$, $L_{x}=32c/\omega_{\mathrm{p}}$, $\Delta
t=0.00390625\omega_{\mathrm{p}}^{-1}$ and $ppc=144$ uniformly distributed and
divided between the two electron populations. The difference here is the
direction of the beam drift velocity, which is now perpendicular to the
domain’s direction i.e. $v_{y}=v_{d}=\pm 0.2c$, and similarly, the thermal
velocity only has a component along the perpendicular direction as well, with
magnitude $v_{\mathrm{th}}=0.001c$. The Weibel instability is excited by
current separation generated by the Lorentz force in a perturbed magnetic
field when acting on the two counterstreaming beams. The resulting current
filaments will reinforce the initial magnetic-field perturbations. This leads
to an exchange of energy between particles and magnetic fields. The magnetic-
field energy in fastest-growing unstable mode will evolve as

$U_{\mathrm{B}}\sim\exp(2\,v_{\mathrm{d}}\,t)$ (28)

during the linear phase of the instability.

Figure 4 shows the results of the Weibel instability tests. The upper-left
panel shows the evolution in time of the magnetic-field energy for all grids.
For all cases, we observe a slight discrepancy between the simulated and
theoretical value of the growth rate during the linear instability phase.
However, the results of the different grids are in strong agreement with each
other. The slight discrepancy with the theoretical growth rate could be
attributed to different factors: for example, the finite spatial resolution
prohibits resolving large wave numbers which would otherwise increase the
observed growth rate, potentially agreeing better with the theoretical maximum
growth rate that corresponds to $k\to\infty$; in addition, noise induced by
the motion of the limited number of particles drowns out small waves that
might otherwise contribute to the instability. At any rate, the discrepancy is
very small and the results appear robust in terms of agreement between
different cases we run. The lower-left panel shows again the evolution of the
relative energy error, indicating conservation of energy to machine precision
for all three grid setups. The panels on the right-hand side of the Figure
show the current filaments forming and growing throughout the domain as time
increases. Considering these results, the method appears capable of correctly
reproducing fully electromagnetic effects while conserving energy to machine
precision, at least in this simple one-dimensional test case.

Figure 4: Top left: Evolution of the magnetic-field energy for the different
grids as well as the theoretical growth rate associated with the linear-growth
phase of the Weibel instability. Bottom left: Evolution of the relative error
in the total energy for the different grids. Right: Formation and growth of
current filaments in the domain with evolving time for the different grids.

###  3.3 GEM reconnection challenge

The last test case presented in this paper is the Geospace Environmental
Modeling (GEM) reconnection challenge. This is a well-known two- (or
three-)dimensional setup to study magnetic-reconnection events (Birn et al.,
2001). We adopt a two-dimensional $xy$ geometry where the $B_{x}$ field
undergoes a sign reversal along the $y$-axis, creating a current sheet in the
inversion region. For simplicity, we adopt periodic boundary conditions,
initializing two successive inversions to respect periodicity. The initial
setup of the magnetic field is

$B_{x}=B_{0}\left[1-\tanh\left(\frac{y-L_{y}/4}{\delta}\right)+\tanh\left(\frac{y-3L_{y}/4}{\delta}\right)\right].$
(29)

To respect the
$c\mbox{\boldmath{$\nabla$}}\times\mbox{\boldmath{$B$}}=4\pi\mbox{\boldmath{$I$}}$
condition, we initialize a drifting electron and positron population in the
current sheets with a number density of

$n=n_{0}\left[\operatorname{sech}^{2}\left(\frac{y-L_{y}/4}{\delta}\right)+\operatorname{sech}^{2}\left(\frac{y-3L_{y}/4}{\delta}\right)\right].$
(30)

Here $B_{0}$ and $n_{0}$ are the peak magnetic field strength and density in
the current sheet respectively, and $\delta$ is the width of the current
sheet, which we choose equal to $0.5c/\omega_{\mathrm{p}}$.

At initialization, the thermal pressure of the particles in the current sheet
must be able to withstand the magnetic pressure to ensure the structure
neither collapses nor expands, i.e. for pressure equilibrium

$\frac{B^{2}_{0}}{8\pi}=2nk_{\mathrm{b}}T,$ (31)

with $k_{\mathrm{b}}$ the Boltzmann constant and $T$ the plasma temperature.
The factor $2$ is introduced to include both electron and positron
contribution to the density. By equating $k_{\mathrm{b}}T$ with
$mv^{2}_{\mathrm{th}}$, a constraint for the thermal speed in the current
sheet can be found,

$\frac{v_{\mathrm{th}}}{c}=\frac{B}{c\sqrt{16\pi nm}}.$ (32)

Note that we do not manually perturb the initial equilibrium, so that tearing
modes eventually causing reconnection are only excited by random particle
noise.

As in the 1D cases, a Cartesian grid is used to establish a baseline to which
we compare nonuniform-grid results. The Cartesian grid has dimensions
$L_{x}=16c/\omega_{\mathrm{p}}$ by $L_{y}=32c/\omega_{\mathrm{p}}$ with
$N_{x}=160$ by $N_{y}=320$ cells and $ppc=16$ particles per cell. The time
step was set to $\Delta t=0.025\omega_{\mathrm{p}}^{-1}$ for $nt=30000$ steps,
for a total time of $750\omega_{\mathrm{p}}^{-1}$. This satisfies the
convergence constraint $c\Delta t/\Delta x<1$.

To test the behavior of our code with curvilinear coordinates, we construct a
nonuniform grid such that we achieve the same density of computational cells
(i.e. the same numerical resolution) near the current sheets, while
significantly reducing the grid resolution in the upstream regions. This
naturally results in fewer total grid cells than in the uniform Cartesian
case. This grid setup enables us to resolve the reconnection region with the
same precision of the Cartesian setup while greatly reducing the total number
of cells in regions of no physical interest. The mapping function for this
grid was created using piecewise linear functions that were connected using
hyperbolic tangents as an approximation for the Heaviside step function,

$\xi(x)=x,$ (33)
$\eta(y)=\frac{L_{y}}{p_{\mathrm{max}}}\left(y-\sum_{n=1}^{4}g_{n}(y)\right),$
(34)

where $g_{n}(y)=((r_{n}y+p_{n})-(r_{n-1}y+p_{n-1}))\tanh(s(y-b_{n}))$. Here,
$r_{n}$ is the cell-density ratio with respect to the low-density upstream
regions. This is set to $1$ in the upstream regions and to a desired larger
ratio for the current-sheet regions. $s$ determines the sharpness of the
transition, with larger values more closely resembling the Heaviside function
and lower values creating a smoother transition region. This value was set to
$s=5$ for all grids. The $b_{n}$ terms determine the break points i.e. the
$y$-coordinate where the piecewise function switches from one linear function
to the next. These breakpoints therefore determine the points of transition
between high and low density regions. These are set at $L_{y}/8$, $3L_{y}/8$,
$5L_{y}/8$ and $7L_{y}/8$ respectively. The $p_{n}$ factors are the offsets of
the individual linear functions, which are fully constrained by the previous
parameters. Finally, $p_{\mathrm{max}}$ is the maximum value reached by the
piecewise function, and is introduced such that $\eta\in[0,L_{y}]$.

We employed two nonuniform grids, the first with a density ratio $r=5$ and
number of cells in the $y$-direction $N_{y}=192$, and the second with ratio
$r=10$ and $N_{y}=176$. The number of cells in the $x$-direction remains
unchanged from the Cartesian case ($N_{x}=160$). By design, the curvilinear
grid has the same physical size $L_{x}=16c/\omega_{\mathrm{p}}$ and
$L_{y}=32c/\omega_{\mathrm{p}}$ of the Cartesian grid. Furthermore, the time
step and number of iterations were deliberately set to be identical in both
cases to again satisfy the convergence constrained, since the smallest $\Delta
x$ is also unchanged.

Our results are shown in Figure 5. The left panel shows the Cartesian grid in
the top half of the domain and the curvilinear grid with density ratio $r=10$
in the bottom half. In the middle panel we show a comparison of the out-of-
plane current $\log_{10}|J_{z}|$ and the in-plane magnetic-field lines between
the Cartesian grid (top) and the curvilinear grid with density ratio $r=10$
(bottom) at time $t=750\omega_{\mathrm{p}}^{-1}$. Although it is clear that
the upstream regions have considerably lower resolution in the curvilinear
case, the reconnection regions of interest look very similar. Since the exact
position of the formed X-points and magnetic islands is due to the random
initialization of the particle properties, they appear in different places in
the top and bottom current sheet in the domain. In the top right panel, we
plot total magnetic and kinetic energies for all cases, to show the energy
exchange between magnetic fields and particles over time. The initial magnetic
field contains the majority of the energy, and as the magnetic fields
reconfigure into a more relaxed state via reconnection, the particles gain
energy and get accelerated. The bottom right panel shows the evolution in time
of the relative error in the total energy, which is again conserved to machine
precision, supporting our findings from the 1D tests and the theory in
Appendix A.

Figure 5: Left: Comparison of the uniform Cartesian grid (top) and curvilinear
grid with $r=10$ (bottom) which has a reduced resolution in the upstream
regions, while maintaining the same resolution in the current sheets that was
used in the Cartesian grid. The magnetic islands formed in the top and bottom
halves are not collocated since their positions are randomized by the starting
values of the particle properties. Middle: Comparison of a snapshot of the
out-of-plane current and the in-plane magnetic-field lines for the
aforementioned grids. Top right: Evolution of the magnetic-field energy and
kinetic energy for the different grids. Bottom right: Evolution of the
relative error in the total energy for the different grids.

##  4 Conclusions

In this paper, a method is presented to adapt a fully implicit PIC method to
work with arbitrary coordinate systems, including curvilinear coordinates,
while maintaining exact (i.e. to machine precision) energy conservation. While
previous methods have been published for reduced equations, like the Darwin
approximation (Chacón & Chen 2016), this is the first time to our knowledge
when this is achieved for a PIC with the full Maxwell equations. This is
achieved by changing the coordinates to a system where the grid becomes
logically rectangular. To change coordinates we introduce the metric tensor,
Jacobian matrix, and contravariant curl operator into Maxwell’s equations for
the PIC algorithm. The particle equations of motion are solved in Cartesian
coordinates, to avoid the complexity of calculating the metric and Jacobian at
each particle position and at each iteration.

To validate the method, we have presented three different test cases. The two-
stream instability and Weibel instability were tested in one spatial dimension
and the GEM reconnection challenge was tested in two dimensions. The tests
runs with nonuniform (curvilinear) grids were compared to the results obtained
with a uniform Cartesian grid, taken as a baseline. Two nonuniform grids were
tested in the one-dimensional cases, both generated by perturbing the baseline
Cartesian grid: the first with a sinusoidal perturbation and the second with a
hyperbolic-tangent perturbation. The nonuniform grid in the GEM challenge was
designed to greatly reduce the total number of grid cells, while maintaining
the same spatial resolution of the uniform grid around the current sheets, to
properly resolve reconnection dynamics while underresolving the unimportant
upstream regions. Comparisons were performed on global quantities such as
individual energy components and the total energy of the system.

We demonstrated that our new method achieves energy conservation to machine
precision in all three test cases with both Cartesian and curvilinear grids.
The two-stream instability modeled with curvilinear grids shows a growth rate
of the electric-field energy which closely matches both the (Cartesian
uniform) baseline and the theoretical growth rate. The Weibel instability also
shows excellent alignment of the magnetic-field energy to the Cartesian
baseline. In our nonuniform-grid runs, we observed that the 2D GEM challenge
produced reconnection events inside the current-sheet region as expected of
this setup, despite the nonuniform grid having only $55-60\%$ the total number
of cells with respect to the uniform Cartesian run. Reconnection can be
observed in the in-plane magnetic-field line rearrangement, which shows the
formation of magnetic islands, and in the exchange between magnetic and
particle energy, showing a conversion of magnetic-field energy into kinetic
energy as particles get accelerated during the reconnection process.

Although our method possesses desirable properties (exact energy conservation
and grid adaptability to complex geometries), ample ground is left for further
improvements, which we will pursue in future work. Unlike the energy, the
total charge is not conserved, since
$\mbox{\boldmath{$\nabla$}}\cdot\mbox{\boldmath{$E$}}=4\pi\rho$ is not
enforced in our implementation. However, as shown in other works (e.g. Chen &
Tóth 2019), a charge-conserving adaptation of fully implicit PIC methods is
possible (Chen et al., 2011), and could be implemented in our code in the
future. Another shortcoming is represented by the convergence of the nonlinear
iteration, on which our method is based, which is not guaranteed; from our
experiments, we observed that a sufficiently high spatiotemporal resolution
must be enforced to avoid nonconvergence issues, which could break the
temporal iteration. In many practical cases however, this will not pose a
problem since a fine spatiotemporal resolution is required anyway, to resolve
the physics of interest. In our tests, convergence (i.e. an absolute error
below a tolerance $\sim 10^{-14}$) was always attained with $\sim 10$
nonlinear iterations; hence, the computational cost in our example simulations
was approximately a factor 10 larger than a standard explicit PIC method,
since each nonlinear iteration roughly corresponds to one explicit time cycle.
The greater computational cost of our method is however counterbalanced by
numerous advantages, as discussed below. In addition, speeding up convergence
is in principle possible with more advanced strategies: for example, Chen et
al. (2014) have shown that by preconditioning the nonlinear solver the average
number of iterations can be drastically reduced, and the spatiotemporal
resolution constraints can be relaxed.

By combining a curvilinear grid with an exactly energy-conserving PIC method,
we have achieved a unique set of features with numerous potential
applications. In general, nonuniform grids are especially well suited for
setups with a small region of interest within a much larger domain, or setups
where the bulk of the plasma is contained within a spatial region with a
nontrivial shape. Meanwhile, energy conservation is important in many problems
where energy conversion from one type to another occurs and has significant
impact on the plasma dynamics, e.g. reconnection events or instabilities (see
Markidis & Lapenta 2011a and references therein). In several cases, the
combination of energy conservation and custom grids can be of high interest:
in fusion plasma physics, our method could be used in simulations of the
scrape-off layer (i.e. the plasma layer close to the separatrix) of a Tokamak
plasma (e.g. Ohtsuka et al. 1978; Xu & Cohen 1998; Krasheninnikov 2001;
Fundamenski et al. 2007). The nonuniform grid can be tailored to match the
geometry of the boundary layer, greatly reducing the total computational cost.
Another example is the expanding solar wind (e.g. Innocenti et al. 2019, 2020;
Bott et al. 2021; Micera et al. 2021): here, a large domain, which expands in
time as the solar wind travels, makes it difficult to employ a standard
Cartesian box while maintaining sufficient resolution throughout the domain.
By using a grid that mimics the shape of the expanding solar wind, with a high
resolution closer to the Sun and decreasing grid density further out (where
less resolution is required), the computational cost of this problem can be
significantly reduced. Conversely, when considering compression-driven
dynamics in plasmas, a compressing box has been used in the past in the
context of simulations of the intergalactic medium (e.g. Sironi & Narayan
2015), which could similarly be replaced with our nonuniform grids. Other
solar wind-related phenomena whose modeling could benefit from a custom grid
are switchbacks. Switchbacks are magnetic structures in the solar wind
characterized by a typical S-shape of the plasma stream, and may be related to
interchange reconnection events and particle scattering and acceleration (e.g.
Drake et al. 2021; Wyper et al. 2022). With our method, a curved grid can be
fitted to this S-shaped stream, reducing the total grid resolution and cutting
out the surrounding plasma that might be of little interest. This might open
up the possibility to a fully kinetic study of the driving phenomena which was
thus far unfeasible. In higher-energy plasma scenarios, our method could also
find applicability e.g. for general-relativistic plasma simulations. Recent
works have simulated plasma accretion around compact objects such as black
holes (e.g. Parfrey et al. 2019; Bransgrove et al. 2021; Crinquand et al.
2022; El Mellah et al. 2022; Galishnikova et al. 2023), where plasma flows in
structures such as accretion disks and jets. These works present large-scale
simulations utilizing explicit codes, usually implemented with one specific
(four-dimensional) metric corresponding to the particular physical case in
question. With the method presented in this work (modified by adding a time
component of the metric and relativistic effects; see e.g. Bacchini et al.
2019; Bacchini 2023), we can simultaneously avoid the aforementioned downsides
of explicit codes (i.e. numerical instability and lack of energy conservation)
and generalize such simulations to any arbitrary metric tensor that might be
of interest, opening the possibility for more physical cases to be studied.

We foresee several potential developments for the future. The logical next
step in this line of work is to implement the method in a production-ready,
parallel infrastructure such as those currently employed for implicit-PIC
simulations (e.g. iPic3D, ECsim). This will open the door to larger-scale
simulations outside the capabilities of our current test implementation in
Python. The most obvious candidate would be the ECsim code (Lapenta 2017),
since it shares (by design) a very similar discretization of our governing PIC
equations. This step will be carried out in the future. It would also be of
great interest to consider novel architectures such as GPUs. The traditional
bottlenecks in PIC codes are the particle gathering and deposition steps. With
the increasing popularity of GPU-accelerated codes, these bottlenecks could be
mitigated, since the gathering and deposition steps are well suited to
parallelization (e.g. Joseph et al. 2011; Decyk & Singh 2014; Vasileska et al.
2021). The field solver, on the other hand, whether implemented on CPUs or
GPUs, requires global communication of grid quantities, which hinders scaling
behavior to a large number of cores, so it becomes imperative to further
reduce the cost of the field-solution step through other means. With the
method described in this paper, this can be achieved by optimizing the grid
and thus reducing the total number of cells, which is the main driving factor
of the field solver’s computational cost.

##  Appendix A Formal proof of energy conservation

Throughout this Appendix the time-level notation has been shortened for
improved readability by dropping the ubiquitous $n$, thus only writing the
offset from this arbitrary time level. For example $(1/2)$ is equivalent to
$(n+1/2)$.

Starting from Maxwell’s equations (3)–(4), we contract both sides with
$B_{\mu}{}_{g}^{(1/2)}$ and $E_{\mu}{}_{g}^{(1/2)}$ respectively:

$\frac{c}{J_{g}}(B_{\mu}\epsilon^{\mu\nu\kappa}\partial_{\nu}E_{\kappa})_{g}^{(1/2)}=-(B_{\mu}\partial_{t}B^{\mu})_{g}^{(1/2)},$
(A1)
$\frac{c}{J_{g}}(E_{\mu}\epsilon^{\mu\nu\kappa}\partial_{\nu}B_{\kappa})_{g}^{(1/2)}=(E_{\mu}\partial_{t}E^{\mu})_{g}^{(1/2)}+4\pi(E_{\mu}I^{\mu})_{g}^{(1/2)}.$
(A2)

Subtracting eq. (A1) from eq.(A2) and rearranging terms gives

$\begin{split}&(E_{\mu}\partial_{t}E^{\mu})_{g}^{(1/2)}+(B_{\mu}\partial_{t}B^{\mu})_{g}^{(1/2)}\\\
=\,\,&\frac{c}{J_{g}}(E_{\mu}\epsilon^{\mu\nu\kappa}\partial_{\nu}B_{\kappa})_{g}^{(1/2)}-\frac{c}{J_{g}}(B_{\mu}\epsilon^{\mu\nu\kappa}\partial_{\nu}E_{\kappa})_{g}^{(1/2)}-4\pi(E_{\mu}I^{\mu})_{g}^{(1/2)}.\end{split}$
(A3)

On the left-hand side, the time derivatives are discretized as described in
Section 2.3 and half-integer timestep values are substituted with a temporal
linear interpolation, producing

$\begin{split}&\frac{E_{\mu}{}_{g}^{(1)}+E_{\mu}{}_{g}^{(0)}}{2}\frac{E^{\mu}{}_{g}^{(1)}-E^{\mu}{}_{g}^{(0)}}{\Delta
t}+\frac{B_{\mu}{}_{g}^{(1)}+B_{\mu}{}_{g}^{(0)}}{2}\frac{B^{\mu}{}_{g}^{(1)}-B^{\mu}{}_{g}^{(1)}}{\Delta
t}\\\
=\,\,&\frac{c}{J_{g}}(E_{\mu}\epsilon^{\mu\nu\kappa}\partial_{\nu}B_{\kappa})_{g}^{(1/2)}-\frac{c}{J_{g}}(B_{\mu}\epsilon^{\mu\nu\kappa}\partial_{\nu}E_{\kappa})_{g}^{(1/2)}-4\pi(E_{\mu}I^{\mu})_{g}^{(1/2)}.\end{split}$
(A4)

Summing over all grid points and multiplying with
$J_{g}\Delta\xi\Delta\eta/4\pi$ while working out the products on the left-
hand side results in

$\begin{split}&\sum_{g}\frac{J_{g}\Delta\xi\Delta\eta}{\Delta
t}\left(\frac{(E_{\mu}E^{\mu})_{g}^{(1)}+(E_{\mu}E^{\mu})_{g}^{(0)}}{8\pi}+\frac{(B_{\mu}B^{\mu})_{g}^{(1)}+(B_{\mu}B^{\mu})_{g}^{(0)}}{8\pi}\right)\\\
=\,\,&\sum_{g}\frac{c\Delta\xi\Delta\eta}{4\pi}\left(E_{\mu}\epsilon^{\mu\nu\kappa}\partial_{\nu}B_{\kappa}-B_{\mu}\epsilon^{\mu\nu\kappa}\partial_{\nu}E_{\kappa}\right)_{g}^{(1/2)}-\sum_{g}J_{g}\Delta\xi\Delta\eta(E_{\mu}I^{\mu})_{g}^{(1/2)}.\end{split}$
(A5)

On the left-hand side we can now recognize the definition of the field
energies as described by eqs. (22)–(23). The first term on the right-hand side
is the divergence of the Poynting-flux vector; this discrete integral
(expressed by the sum over all grid points $g$) vanishes under periodic
boundary conditions. This can be easily shown for the one-dimensional case by
expanding the sum,

$\begin{split}\sum_{i=0}^{N}&-E_{\eta(i+1/2)}(B_{\zeta(i+1)}-B_{\zeta(i)})+E_{\zeta(i+1/2)}(B_{\eta(i+1)}-B_{\eta(i)})\\\
&+B_{\eta(i)}(E_{\zeta(i+1/2)}-E_{\zeta(i-1/2)})-B_{\zeta(i)}(E_{\eta(i+1/2)}-E_{\eta(i-1/2)})=0,\end{split}$
(A6)

where every term can indeed be canceled out if the 0-th and $N$-th grid points
are identical. This argument can easily be extended to two and three
dimensions as well. The right-hand side in eq. (LABEL:eq:apx_long) thus
reduces to the second term and can be used to denote the change in field
energy between two subsequent time levels,

$\frac{\Delta U_{\mathrm{F}}}{\Delta
t}=-\sum_{g}J_{g}\Delta\xi\Delta\eta(E_{\mu}I^{\mu})_{g}^{(1/2)},$ (A7)

where the remaining term represents energy exchange between fields and
particles.

We now have to verify that this term is identical to the energy-exchange term
obtained from the particle equations. To do so we start from eq. (20) and
contract both sides with $v_{i}{}_{p}^{(1/2)}$,

$m_{p}v_{i}{}_{p}^{(1/2)}\frac{v^{i}{}_{p}^{(1)}-v^{i}{}_{p}^{(0)}}{\Delta
t}=q_{p}\left(v_{i}E^{i}+v_{i}\epsilon^{ijk}v_{j}B_{k}\right)_{p}^{(1/2)}.$
(A8)

In the last term on the right-hand side, we can recognize
$v_{i}\epsilon^{ijk}v_{j}B_{k}=0$ since this is the dot product of $v$ with
the cross product between itself and $B$. Like before, we rewrite the half-
integer timesteps as temporal averages in the left-hand side, and we take the
sum over all particles $p$ to get

$\sum_{p}\frac{m}{\Delta
t}\frac{(v_{i}v^{i})_{p}^{(1)}-(v_{i}v^{i})_{p}^{(0)}}{2}=\sum_{p}q_{p}(v_{i}E^{i})_{p}^{(1/2)}.$
(A9)

We can recognize the definition of the kinetic energy as shown in eq. (24) and
use the expression for interpolated quantities (eq. (21)) in the right hand
side:

$\frac{\Delta U_{\mathrm{K}}}{\Delta
t}=\sum_{p}q_{p}v_{i}{}_{p}^{(1/2)}\sum_{g}w_{gp}(j_{\nu}^{i}E^{\nu})_{g}^{(1/2)}.$
(A10)

Since the interpolation weight from grid to particles $w_{gp}$ is the same as
the weight from particles to grid $w_{pg}$, it is possible to switch the order
of the summations,

$\frac{\Delta U_{\mathrm{K}}}{\Delta
t}=\sum_{g}\sum_{p}q_{p}v_{i}{}_{p}^{(1/2)}w_{pg}(j_{\nu}^{i}E^{\nu})_{g}^{(1/2)},$
(A11)

and by adding the missing factors we can recover the right-hand side of eq.
(18) for the current density,

$\frac{\Delta U_{\mathrm{K}}}{\Delta
t}=\sum_{g}J_{g}\Delta\xi\Delta\eta(j_{\nu}^{i}E^{\nu}j_{\mu}^{i})_{g}^{(1/2)}\frac{({j^{-1}}_{i}^{\mu})_{g}}{J_{g}\Delta\xi\Delta\eta}\sum_{p}w_{pg}q_{p}v^{i}{}_{p}^{(1/2)}.$
(A12)

Substituting and rearranging the Jacobian factors results in

$\frac{\Delta U_{\mathrm{K}}}{\Delta
t}=\sum_{g}J_{g}\Delta\xi\Delta\eta(j_{\nu}^{i}j_{\mu}^{i}E^{\nu}I^{\mu})_{g}^{(1/2)},$
(A13)

where we now recognize the metric $j_{\nu}^{i}j_{\mu}^{i}=g_{\nu\mu}$ and use
it to lower the index of the electric field, such that

$\frac{\Delta U_{\mathrm{K}}}{\Delta
t}=\sum_{g}J_{g}\Delta\xi\Delta\eta(E_{\mu}I^{\mu})_{g}^{(1/2)},$ (A14)

where the energy-exchange term with the field energy is recovered. Combining
the change in kinetic energy (A14) and field energy (A7) gives the change in
total energy of the system between time steps,

$\frac{\Delta U}{\Delta t}=\frac{\Delta U_{\mathrm{K}}}{\Delta t}+\frac{\Delta
U_{\mathrm{F}}}{\Delta
t}=\sum_{g}J_{g}\Delta\xi\Delta\eta(E_{\mu}I^{\mu})_{g}^{(1/2)}-\sum_{g}J_{g}\Delta\xi\Delta\eta(E_{\mu}I^{\mu})_{g}^{(1/2)}=0.$
(A15)

Since this vanishes exactly between any two timesteps, the method conserves
the total energy $U$ exactly (to machine precision) independently of the
timestep, grid spacing, and any other numerical parameter.

## Acknowledgements

This work was supported by the KU Leuven Bijzonder Onderzoeksfonds (BOF) under
the C1 project TRACESpace. J.C. acknowledges support from the European Union
Horizon 2020 project DEEP-SEA (Grant agreement 955606). L.P. acknowledges
support from a PhD grant awarded by the Royal Observatory of Belgium. F.B.
acknowledges support from the FED-tWIN programme (profile Prf-2020-004,
project “ENERGY”) issued by BELSPO. G.L. acknowledges support from the
European Research Council (ERC) Advanced Grant “TerraVirtualE”. The
computational resources and services used in this work were provided by the
VSC (Flemish Supercomputer Center), funded by the Research Foundation –
Flanders (FWO) and the Flemish Government – department EWI.

## References

  * Bacchini (2023) Bacchini, F. 2023, ApJS, 268, 60
  * Bacchini et al. (2019) Bacchini, F., Amaya, J., & Lapenta, G. 2019, J. Phys: Conf. Ser., 1225, 012011
  * Birn et al. (2001) Birn, J., Drake, J. F., Shay, M. A., et al. 2001, Journal of Geophysical Research: Space Physics, 106, 3715, doi: https://doi.org/10.1029/1999JA900449
  * Bott et al. (2021) Bott, A. F. A., Arzamasskiy, L., Kunz, M. W., Quataert, E., & Squire, J. 2021, The Astrophysical Journal Letters, 922, L35, doi: 10.3847/2041-8213/ac37c2
  * Bransgrove et al. (2021) Bransgrove, A., Ripperda, B., & Philippov, A. 2021, Phys. Rev. Lett., 127, 055101, doi: 10.1103/PhysRevLett.127.055101
  * Bret et al. (2010) Bret, A., Gremillet, L., & Dieckmann, M. E. 2010, Physics of Plasmas, 17, 120501, doi: 10.1063/1.3514586
  * Chacón & Chen (2016) Chacón, L., & Chen, G. 2016, Journal of Computational Physics, 316, 578, doi: https://doi.org/10.1016/j.jcp.2016.03.070
  * Chen et al. (2011) Chen, G., Chacón, L., & Barnes, D. C. 2011, Journal of Computational Physics, 230, 7018
  * Chen et al. (2014) Chen, G., Chacón, L., Leibs, C., Knoll, D., & Taitano, W. 2014, Journal of Computational Physics, 258, 555, doi: 10.1016/j.jcp.2013.10.052
  * Chen & Tóth (2019) Chen, Y., & Tóth, G. 2019, Journal of Computational Physics, 386, 632, doi: 10.1016/j.jcp.2019.02.032
  * Crinquand et al. (2022) Crinquand, B., Cerutti, B., Dubus, G., Parfrey, K., & Philippov, A. 2022, Phys. Rev. Lett., 129, 205101, doi: 10.1103/PhysRevLett.129.205101
  * Decyk & Singh (2014) Decyk, V. K., & Singh, T. V. 2014, Computer Physics Communications, 185, 708, doi: https://doi.org/10.1016/j.cpc.2013.10.013
  * Delzanno et al. (2013) Delzanno, G. L., Camporeale, E., Moulton, J. D., et al. 2013, IEEE Transactions on Plasma Science, 41, 3577, doi: 10.1109/TPS.2013.2290060
  * Drake et al. (2021) Drake, J., Agapitov, O., Swisdak, M., et al. 2021, Astronomy and astrophysics (Berlin), 650, A2
  * El Mellah et al. (2022) El Mellah, I., Cerutti, B., Crinquand, B., & Parfrey, K. 2022, Astronomy and Astrophysics, 663, A169, doi: 10.1051/0004-6361/202142847
  * Fichtl et al. (2012) Fichtl, C. A., Finn, J. M., & Cartwright, K. L. 2012, Computational Science and Discovery, 5, 014011, doi: 10.1088/1749-4699/5/1/014011
  * Fundamenski et al. (2007) Fundamenski, W., Garcia, O., Naulin, V., et al. 2007, Nuclear Fusion, 47, 417, doi: 10.1088/0029-5515/47/5/006
  * Galishnikova et al. (2023) Galishnikova, A., Philippov, A., Quataert, E., et al. 2023, Phys. Rev. Lett., 130, 115201, doi: 10.1103/PhysRevLett.130.115201
  * Goldston & H. (1995) Goldston, R. J., & H., R. P. 1995, Introduction to plasma physics, 1st edn. (IOP Publishing Ltd)
  * Gonzalez-Herrero et al. (2019) Gonzalez-Herrero, D., Micera, A., Boella, E., Park, J., & Lapenta, G. 2019, Computer Physics Communications, 236, 153, doi: https://doi.org/10.1016/j.cpc.2018.10.026
  * Innocenti et al. (2020) Innocenti, M. E., Boella, E., Tenerani, A., & Velli, M. 2020, The Astrophysical Journal Letters, 898, L41, doi: 10.3847/2041-8213/aba591
  * Innocenti et al. (2019) Innocenti, M. E., Tenerani, A., & Velli, M. 2019, The Astrophysical Journal, 870, 66, doi: 10.3847/1538-4357/aaf1be
  * Jolliet et al. (2007) Jolliet, S., Bottino, A., Angelino, P., et al. 2007, Computer Physics Communications, 177, 409, doi: https://doi.org/10.1016/j.cpc.2007.04.006
  * Joseph et al. (2011) Joseph, R. G., Ravunnikutty, G., Ranka, S., D’Azevedo, E., & Klasky, S. 2011, in 2011 IEEE International Parallel and Distributed Processing Symposium, 395–406, doi: 10.1109/IPDPS.2011.46
  * Krasheninnikov (2001) Krasheninnikov, S. 2001, Physics Letters A, 283, 368, doi: https://doi.org/10.1016/S0375-9601(01)00252-3
  * Lapenta (2017) Lapenta, G. 2017, Journal of Computational Physics, 334, 349, doi: https://doi.org/10.1016/j.jcp.2017.01.002
  * Liseikin (1999) Liseikin, V. D. 1999, Grid generation methods, Vol. 1 (Springer)
  * Markidis & Lapenta (2011a) Markidis, S., & Lapenta, G. 2011a, Journal of Computational Physics, 230, 7037, doi: 10.1016/j.jcp.2011.05.033
  * Markidis & Lapenta (2011b) —. 2011b, Journal of Computational Physics, 230, 7037
  * Markidis et al. (2010) Markidis, S., Lapenta, G., & Rizwan-Uddin. 2010, Mathematics and Computers in Simulation, 80, 7
  * Micera et al. (2021) Micera, A., Zhukov, A. N., López, R. A., et al. 2021, The Astrophysical Journal, 919, 42, doi: 10.3847/1538-4357/ac1067
  * Ohtsuka et al. (1978) Ohtsuka, H., Kimura, H., Shimomura, S., et al. 1978, Plasma Physics, 20, 749, doi: 10.1088/0032-1028/20/8/003
  * Parfrey et al. (2019) Parfrey, K., Philippov, A., & Cerutti, B. 2019, Phys. Rev. Lett., 122, 035101, doi: 10.1103/PhysRevLett.122.035101
  * Ringle (2011) Ringle, R. 2011, International journal of mass spectrometry, 303, 42
  * Sironi & Narayan (2015) Sironi, L., & Narayan, R. 2015, The Astrophysical Journal, 800, 88, doi: 10.1088/0004-637X/800/2/88
  * Stanier & Chacón (2022) Stanier, A., & Chacón, L. 2022, Journal of Computational Physics, 459, 111144, doi: https://doi.org/10.1016/j.jcp.2022.111144
  * Vasileska et al. (2021) Vasileska, I., Bogdanović, L., & Kos, L. 2021, in 2021 44th International Convention on Information, Communication and Electronic Technology (MIPRO), 291–294, doi: 10.23919/MIPRO52101.2021.9596959
  * Weibel (1959) Weibel, E. S. 1959, Physical review letters, 2, 83
  * Wyper et al. (2022) Wyper, P. F., DeVore, C. R., Antiochos, S. K., et al. 2022, The Astrophysical Journal Letters, 941, L29, doi: 10.3847/2041-8213/aca8ae
  * Xu & Cohen (1998) Xu, X. Q., & Cohen, R. H. 1998, Contributions to Plasma Physics, 38, 158, doi: https://doi.org/10.1002/ctpp.2150380124