# DeBERTinha: A Multistep Approach to Adapt DebertaV3 XSmall for Brazilian
Portuguese Natural Language Processing Tasks

Israel Campiotti1<EMAIL_ADDRESS>Matheus Rodrigues2<EMAIL_ADDRESS>Yuri Albuquerque3<EMAIL_ADDRESS>
Rafael Azevedo4<EMAIL_ADDRESS>Alyson Andrade
<EMAIL_ADDRESS>
(1University of Campinas, 2Federal University of Pernambuco, 3University of
São Paulo,
4Pontifical Catholic University of Paraná)

###### Abstract

This paper presents an approach for adapting the DebertaV3 XSmall model pre-
trained in English for Brazilian Portuguese natural language processing (NLP)
tasks. A key aspect of the methodology involves a multi-step training process
to ensure the model is effectively tuned for the Portuguese language. Initial
datasets from Carolina and BrWac are preprocessed to address issues like
emojis, HTML tags, and encodings. A Portuguese-specific vocabulary of 50,000
tokens is created using SentencePiece. Rather than training from scratch, the
weights of the pre-trained English model are used to initialize most of the
network, with random embeddings, recognizing the expensive cost of training
from scratch. The model is fine-tuned using the replaced token detection task
in the same format of DebertaV3 training. The adapted model, called
DeBERTinha, demonstrates effectiveness on downstream tasks like named entity
recognition, sentiment analysis, and determining sentence relatedness,
outperforming BERTimbau-Large in two tasks despite having only 40M parameters.

##  1 Introduction

In this study, we leverage the advancements introduced by the DebertaV3 model,
which currently offers only a large Portuguese variant with 900M parameters.
From the perspective of DebertaV3, we aim to harness the merits of its
architectural excellence and craft a Portuguese version. When compared to
other Portuguese encoder Transformers, like BERTimbau, DeBERTinha holds a
distinctive advantage with its considerably reduced parameter count. Despite
its reduced size, DeBERTinha manages to outperform BERTimbau in two out of
four selected NLP tasks, all while being approximately five times more
compact. The substantial advantage of a diminished parameter count becomes
especially evident in corporate applications, as it demands fewer
computational resources for predictive tasks and offers the benefit of reduced
response times.

Our construction steps begin with the acquisition of essential datasets,
including the Carolina and BrWac datasets, which serve as the foundation for
pre-training. These datasets, however, present challenges such as the presence
of emojis, HTML tags, and non-standard characters. To address these issues, we
employ emoji removal, html tag cleansing using ftfy, and the normalization of
utf8 encodings. To create a tailored vocabulary for Portuguese, we utilize the
Carolina dataset and employ the SentencePiece algorithm, resulting in 50
thousand tokens. These tokens form the basis for creating our tokenizer, a
crucial component for text processing. Subsequently, we merge the Carolina and
BrWac datasets to create a unified dataset for pre-training. Tokenization is a
pivotal step in NLP, and our next step involves tokenizing the merged dataset.
To ensure compatibility with model input limitations, we divide texts into
smaller pieces, each containing a maximum of 510 tokens. This is further
complemented by the addition of special tokens, [CLS] at the beginning and
[SEP] at the end, resulting in a 512-token structure.

One unique aspect of our approach involves the initialization of the model.
While we start with random embeddings, the rest of the model benefits from the
weights of the DebertaV3 XSmall model pre-trained in English. This approach is
based on the recognition that training from scratch in Portuguese can be
prohibitively expensive.

For model pre-training, we leverage the combinations of Masked Language
Modeling (MLM) and Replaced Token Detection (RTD) losses with a hyperparameter
$\lambda$ within the same mathematical framework explained in [1]. To ensure
the preservation of model progress, we save checkpoints for both the generator
and discriminator models. These checkpoints play a crucial role in subsequent
fine-tuning tasks.

In the latter part of this article, we demonstrate the versatility of our
adapted model by using the saved weights for Named Entity Recognition (NER)
and classification tasks. Our experimentation extends to the ASSIN2 dataset,
where we evaluate the model’s ability to determine the relatedness of two
sentences. Our entire multi-step approach can therefore offer a simple
methodology for adapting pre-trained models to new languages, with a specific
focus on Brazilian Portuguese. Through this work, we contribute to the growing
body of research in the field of NLP, offering insights into effective cross-
lingual adaptation and showcasing the versatility of DebertaV3 XSmall in
addressing various NLP tasks in the Portuguese language context.

##  2 Related Work

When tackling tasks like replace token detection, where the model needs to
accurately predict missing words or phrases within a given context, language-
specific models tend to reach better results than multi language models.
Similarly, [2] exhibited superior performance across tasks and settings by
training on Catalan, a moderately under-resourced language. In [3], a model
with 10 billion parameters trained primarily on a Chinese corpus showcased
significant outperformance on 54 Chinese NLP tasks. Exploring the feasibility
of training monolingual Transformer-based language models in French, [4]
demonstrated comparable results using a 4GB web-crawled dataset, rivaling
outcomes from larger datasets that typically exceed 130GB.

While multilingual models yield impressive results, they tend to have larger
sizes. In specific instances, such as Portuguese, their performance can fall
short compared to monolingual counterparts [5], especially for high-resource
languages. Unfortunately, this limitation has constrained the availability of
these cutting-edge models primarily to English in the monolingual context,
causing inconvenience by impeding their practical application in NLP systems.
Furthermore, it hampers exploration into their language modeling capabilities,
particularly concerning morphologically rich languages like Portuguese, [6]
address those details in the context of large language models. This
specialization allows those monolingual models to capture the nuances, idioms,
and domain-specific terms unique to a language more effectively. In replace
token detection, having a precise vocabulary is crucial for accurately
predicting missing tokens.

In our current study, we adopt a specific methodology to leverage the
capabilities of the pre-trained language model known as DeBERTaV3, as outlined
in the research [1]. This particular model represents an advancement over the
original DeBERTa model introduced by [7]. Notably, this improvement is
achieved by departing from the conventional masked language modeling (MLM)
approach and instead embracing replaced token detection (RTD). This RTD
methodology is rooted in the concept elucidated in [8], where the core
technique revolves around the utilization of two transformer encoders – a
generator and a discriminator – within the framework of a replaced token
detection task.

Within the basic conceptional structure of DeBERTaV3, a noteworthy
architectural feature involves the sharing of embeddings between the generator
and the discriminator. This departure from the conventional model setup, where
embeddings are handled independently, is achieved through a method referred to
as gradient-disentangled embedding sharing (GDES). In this approach, the
generator and the discriminator share embeddings, but a critical distinction
is that the flow of gradients from the discriminator to the generator
embeddings is intentionally halted.

##  3 Data sets

To train language models that can compete with other state-of-the-art models,
a substantial volume of data is indispensable. In the context of Brazilian
Portuguese, we can draw attention to BERTimbau [5], which harnessed the brWAC
dataset [9], yielding a hefty 17.5 GB of raw text after rigorous processing.
Recent works, such as Albertina PT-BR [10], have also made use of this same
dataset.

In this particular endeavor, we exploit two distinct datasets during the pre-
training phase, brWAC and Carolina. In line the well established research
practices, we made use of the brWAC dataset, comprising a web crawl of
Brazilian web pages, culminating in an approximately 17.5 GB corpus of
processed text. Additionally, we incorporated the Carolina dataset [11], with
823 million tokens and encompassing two million text entries. To provide a
contextual background, Carolina dataset encompasses a diverse array of
domains, including wikis, university documents, social media interactions, the
legislative branch, and various other public domain sources. In particular, we
selectively utilized only the document body, disregarding titles, metadata,
and other information associated with the text. We employed the emoji and ftfy
[12] libraries to respectively handle emojis and HTML-related content, while
also organizing the data into a single file format with each line containing
510 tokens. In total, the raw pre-processed texts amounts to 33 GB of data.

##  4 DeBERTinha a Brazilian Portuguese model

Prior research efforts [5, 6, 13] have highlighted the benefits of training
Transformers models for specific languages. While multilingual models have
shown improvements, monolingual training and tokenization methods continue to
hold significance in achieving state-of-the-art performance on monolingual
datasets. Consequently, our focus is on training the Deberta-V3 model for the
Portuguese language.

Deberta-V3 has demonstrated superior performance compared to other encoder
architectures, primarily attributed to its pre-training methodology and
extensive vocabulary. Although there has been prior work in training Deberta
for Portuguese [10], we identified opportunities for enhancement:

  * •

We trained a Portuguese tokenizer using text from the Carolina dataset. Given
that our dataset was smaller than the one used for training the English
tokenizer, we opted for a vocabulary size of 50 thousand tokens.

  * •

We maintained fixed examples with a size of 512 tokens, fully utilizing the
maximum context window of the Deberta model.

  * •

Our goal was on training the Deberta-V3 XSmall version, which consists of only
40M parameters. Additionally, we continued training on the RTD task using the
weights available for the English model on the Hugging Face model hub.

Due to resource constraints, the training process was divided into two phases.
The first phase involved training for one epoch using a batch size of 1664 on
8x 80GB A100 GPUs, while the second phase encompassed an additional two epochs
with a batch size of 288 on 8x 32GB V100 GPUs. The cumulative training time
amounted to 12.5 hours, incurring a cost of cloud computing of near four
hundred dollars only. This calculation includes the 2.5 hours required for
loading the dataset on each machine.

##  5 Experimental Setup

In this section we report the experimental results obtained by DeBERTinha in 4
different tasks: ASSIN2-RTE and ASSIN2-STS [14], LeNERBR [15], HateBR [16]. We
used a Google Colab instance with one 16GB T4 GPU to finetune on each task for
an average time of 30min. Each training uses AdamW as optimizer with a
learning rate of 0.00005, trained for a maximum of 20 epochs with early
stopping. We compare the results of our DeBERTinha, that contains 40M
paramenters, against the baseline models BERTimbau-Large (335M) and Albertina-
Large (900M).

ASSIN2 contains sentence pairs of premises and hypothesis and for each pair
there are two annotated targets: a binary classification, indicating whether
premise and hypthosis are related or not, this task is denoted as ASSIN2-RTE;
a similarity score between 0 and 5, indicating how similar the premise and
hypothesis are, this task is named ASSIN2-STS. Accuracy and Pearson
correlation are used as metrics for ASSIN2-RTE and ASSIN2-STS, respectively.

LeNERBR is a dataset for the Named Entity Recognition (NER) task in Portuguese
legal text. It contains seven different classes: Organization, Person, Time,
Local, Legislation, Jurisprudence and Other. The Other label is given to every
token that does not fall into any of the previious categories. We use the
standard BIO representation and only predict scores for the first token of
each word. We use F1 score to asses performance on this task.

HateBR is a dataset for classifying whether a tweet is inappropriate/malicious
or not. We use Accuracy as our metric for this task.

From the results shown in Table 1 we see that our DeBERTinha model surpasses
BERTimbau-Large in two out of the four datasets: for the ASSIN2-RTE DeBERTinha
achieves 89.99% Accuracy against the 89.13% of the BERTimbau-Large; in LeNERBR
our model achieves a F1 score of 90.19%, while BERTimbau-Large achieves
90.15%. In the other two datasets DeBERTinha achieves over 97% of the
BERTimbau-Large performance.

Both models still fall behind Albertina’s 91.30% Accuracy on ASSIN2-RTE and
86.76% Pearson Correlation on ASSIN2-STS. However, Albertina is 22.5 times
bigger than DeBERTinha and 2.6 times larger than BERTimbau-Large. Because
training Albertina demands heavy computational resources we do not train it on
the other datasets and report the results taken directly from the original
paper.

Model/Dataset | LeNERBR | ASSIN2-RTE | ASSIN2-STS | HateBR
---|---|---|---|---
BERTimbau-Large | 90.15 | 89.13 | 85.31 | 93.60
Albertina | - | 91.30 | 86.76 | -
DeBERTinha | 90.19 | 89.99 | 84.75 | 91.28
Table 1: Results from BERTimbau-Large, Albertina and DeBERTinha on 4 commonly
used Portuguese datasets.

##  6 Conclusion

This work has shown promising results for adapting a pre-trained language
model to Portuguese through a carefully designed multi-step methodology. The
resulting DeBERTinha model demonstrated competitive performance across several
NLP tasks compared to larger baseline models, highlighting the effectiveness
of the proposed approach.

On the ASSIN2-RTE task for sentence relatedness classification, DeBERTinha
achieved an accuracy of 89.99%, outperforming BERTimbau-Large which scored
89.13%. For the named entity recognition task on the LeNERBR dataset,
DeBERTinha achieved an F1-score of 90.19%, slightly higher than BERTimbau-
Large’s score of 90.15%. While Albertina still achieved better results on some
tasks, DeBERTinha was able to attain over 97% of the performance of BERTimbau-
Large, despite having only 40M parameters compared to BERTimbau-Large’s 335M
parameters.

This shows that the multi-step adaptation process leveraging an existing pre-
trained model can produce a more specialized Portuguese model that rivals or
exceeds the performance of much larger baseline models, demonstrating an
effective approach for resource-constrained scenarios. The methodology
introduced in this work contributes to research on cross-lingual transfer
learning and model adaptation techniques. Overall, the results validate
DeBERTinha as a promising lightweight model for Brazilian Portuguese NLP
applications.

For future work we aim at assessing DeBERTinha’s performance in a bigger range
of Portuguese NLP tasks, such as information retrieval, and extend its context
to 1024 and 2048 tokens.

##  7 Acknowledgements

We would like to express our heartfelt gratitude to Letrus for their generous
funding of the cloud services related to the processing and training of our
language model. Their support was of a central importance in the execution of
this work.

## References

  * [1] Pengcheng He, Jianfeng Gao, and Weizhu Chen.  Debertav3: Improving deberta using electra-style pre-training with gradient-disentangled embedding sharing.  arXiv preprint arXiv:2111.09543, 2021.
  * [2] Jordi Armengol-Estapé, Casimiro Pio Carrino, Carlos Rodriguez-Penagos, Ona de Gibert Bonet, Carme Armentano-Oller, Aitor Gonzalez-Agirre, Maite Melero, and Marta Villegas.  Are multilingual models the best choice for moderately under-resourced languages? a comprehensive assessment for catalan.  arXiv preprint arXiv:2107.07903, 2021.
  * [3] Yu Sun, Shuohuan Wang, Shikun Feng, Siyu Ding, Chao Pang, Junyuan Shang, Jiaxiang Liu, Xuyi Chen, Yanbin Zhao, Yuxiang Lu, et al.  Ernie 3.0: Large-scale knowledge enhanced pre-training for language understanding and generation.  arXiv preprint arXiv:2107.02137, 2021.
  * [4] Louis Martin, Benjamin Muller, Pedro Javier Ortiz Suárez, Yoann Dupont, Laurent Romary, Éric Villemonte de La Clergerie, Djamé Seddah, and Benoît Sagot.  Camembert: a tasty french language model.  arXiv preprint arXiv:1911.03894, 2019.
  * [5] Fábio Souza, Rodrigo Nogueira, and Roberto Lotufo.  Bertimbau: Pretrained bert models for brazilian portuguese.  In Ricardo Cerri and Ronaldo C. Prati, editors, Intelligent Systems, pages 403–417, Cham, 2020. Springer International Publishing.
  * [6] Ramon Pires, Hugo Abonizio, Thales Rogério, and Rodrigo Nogueira.  Sabiá: Portuguese large language models.  arXiv preprint arXiv:2304.07880, 2023.
  * [7] Pengcheng He, Xiaodong Liu, Jianfeng Gao, and Weizhu Chen.  Deberta: Decoding-enhanced bert with disentangled attention.  arXiv preprint arXiv:2006.03654, 2020.
  * [8] Kevin Clark, Minh-Thang Luong, Quoc V Le, and Christopher D Manning.  Electra: Pre-training text encoders as discriminators rather than generators.  arXiv preprint arXiv:2003.10555, 2020.
  * [9] Jorge A Wagner Filho, Rodrigo Wilkens, Marco Idiart, and Aline Villavicencio.  The brwac corpus: A new open resource for brazilian portuguese.  In Proceedings of the Eleventh International Conference on Language Resources and Evaluation (LREC 2018), 2018.
  * [10] João Rodrigues, Luís Gomes, João Silva, António Branco, Rodrigo Santos, Henrique Lopes Cardoso, and Tomás Osório.  Advancing neural encoding of portuguese with transformer albertina pt.  arXiv preprint arXiv:2305.06721, 2023.
  * [11] Marcelo Finger, Maria Clara Paixão de Sousa, Cristiane Namiuti, Vanessa Martins do Monte, Aline Silva Costa, Felipe Ribas Serras, Mariana Lourenço Sturzeneker, Raquel de Paula Guets, Renata Morais Mesquita, Guilherme Lamartine de Mello, Maria Clara Ramos Morales Crespo, Maria Lina de Souza Jeannine Rocha, Patrícia Brasil, Mariana Marques da Silva, and Mayara Feliciano Palma.  Carolina: The open corpus for linguistics and artificial intelligence.  https://sites.usp.br/corpuscarolina/corpus, 2022.  Version 1.1 (Ada).
  * [12] Robyn Speer.  ftfy.  Zenodo, 2019.  Version 5.5.
  * [13] Diedre Carmo, Marcos Piau, Israel Campiotti, Rodrigo Nogueira, and Roberto Lotufo.  Ptt5: Pretraining and validating the t5 model on brazilian portuguese data, 2020.
  * [14] Livy Real, Erick Fonseca, and Hugo Goncalo Oliveira.  The assin 2 shared task: a quick overview.  In International Conference on Computational Processing of the Portuguese Language, pages 406–412. Springer, 2020.
  * [15] Pedro H. Luz de Araujo, Teófilo E. de Campos, Renato R. R. de Oliveira, Matheus Stauffer, Samuel Couto, and Paulo Bermejo.  LeNER-Br: a dataset for named entity recognition in Brazilian legal text.  In International Conference on the Computational Processing of Portuguese (PROPOR), Lecture Notes on Computer Science (LNCS), pages 313–323, Canela, RS, Brazil, September 24-26 2018. Springer.
  * [16] Francielle Alves Vargas, Isabelle Carvalho, Fabiana Rodrigues de Góes, Fabrício Benevenuto, and Thiago Alexandre Salgueiro Pardo.  Hatebr: A large expert annotated corpus of brazilian instagram comments for offensive language and hate speech detection, 2022.


# Numerical study of $\delta$-function current sheets arising from resonant
magnetic perturbations

Yi-Min Huang<EMAIL_ADDRESS>Department of Astrophysical Sciences,
Princeton University, Princeton, New Jersey 08544, USA  Princeton Plasma
Physics Laboratory, Princeton, New Jersey 08543, USA  Stuart R. Hudson
Princeton Plasma Physics Laboratory, Princeton, New Jersey 08543, USA  Joaquim
Loizu  École Polytechnique Fédérale de Lausanne, Swiss Plasma Center, CH-1015
Lausanne, Switzerland  Yao Zhou  Princeton Plasma Physics Laboratory,
Princeton, New Jersey 08543, USA  Institute of Natural Sciences, School of
Physics and Astronomy, and MOE-LSC, Shanghai Jiao Tong University, Shanghai
200240, China  Amitava Bhattacharjee  Department of Astrophysical Sciences,
Princeton University, Princeton, New Jersey 08544, USA  Princeton Plasma
Physics Laboratory, Princeton, New Jersey 08543, USA

###### Abstract

General three-dimensional toroidal ideal magnetohydrodynamic equilibria with a
continuum of nested flux surfaces are susceptible to forming singular current
sheets when resonant perturbations are applied. The presence of singular
current sheets indicates that, in the presence of non-zero resistivity,
magnetic reconnection will ensue, leading to the formation of magnetic islands
and potentially regions of stochastic field lines when islands overlap.
Numerically resolving singular current sheets in the ideal MHD limit has been
a significant challenge. This work presents numerical solutions of the Hahm-
Kulsrud-Taylor (HKT) problem, which is a prototype for resonant singular
current sheet formation. The HKT problem is solved by two codes: a Grad-
Shafranov (GS) solver and the SPEC code. The GS solver has built-in nested
flux surfaces with prescribed magnetic fluxes. The SPEC code implements multi-
region relaxed magnetohydrodynamics (MRxMHD), whereby the solution relaxes to
a Taylor state in each region while maintaining force balance across the
interfaces between regions. As the number of regions increases, the MRxMHD
solution appears to approach the ideal MHD solution assuming a continuum of
nested flux surfaces. We demonstrate agreement between the numerical solutions
obtained from the two codes through a convergence study.

##  I Introduction

Ideal magnetohydrodynamics (MHD) permits solutions with singular current
sheets.[1] General three-dimensional (3D) ideal MHD equilibria with a
continuum of nested flux surfaces, as often assumed by stellarator equilibrium
solvers such as VMEC[2] and NSTAB,[3] are susceptible to the formation of
singular current sheets at rational surfaces.[4, 5, 6] Nominally two-
dimensional (2D) systems such as tokamaks can also develop singular current
sheets when subjected to resonant magnetic perturbations (RMPs). The formation
of ideal MHD singular current sheets has significant practical implications.
With a finite resistivity or other non-ideal effects that enable magnetic
reconnection, magnetic field lines surrounding the ideal singular current
sheets will break and reconnect, thereby releasing magnetic energy;
consequently, the magnetic field will evolve into a field with magnetic
islands and possibly regions of stochastic field lines if islands overlap.[7,
8] The sites of ideal MHD singular current sheets, therefore, serve as an
indicator of where magnetic reconnection will occur. The intensities of the
current sheets also measure the amount of energy available for reconnection.

Figure 1: A sketch of the Hahm-Kulsrud-Taylor problem. The in-plane components
of the magnetic field reverse directions at the mid-plane (the dashed line).
The upper and lower boundaries are shaped by mirror-symmetric sinusoidal
perturbations. In response to the perturbation, a singular current sheet
develops at the mid-plane.

A prototype for singular current sheet formation driven by RMPs is the Hahm-
Kulsrud-Taylor (HKT) problem, [9, 10, 11] shown in Figure 1. This 2D problem
has a magnetized plasma enclosed by two conducting walls in slab geometry.
Before the conducting walls are perturbed, the initial magnetic field is a
smooth function of space. The in-plane component points along the $y$
direction and $B_{y}$ reverses direction at the mid-plane (the dashed line in
Figure 1). A non-uniform $B_{z}$ component renders the magnetic field force-
free. We then impose a sinusoidal perturbation with an up-down symmetry to the
conducting walls and look for a new ideal equilibrium that is consistent with
the boundary perturbation while conserving magnetic fluxes between flux
surfaces. In this new equilibrium, a singular current sheet will develop at
the mid-plane, which is a flux surface that resonates with the boundary
perturbation.

The primary objectives of this work are (1) investigating the nature of ideal
singular current sheets by constructing numerical solutions of the HKT problem
as an example and (2) assessing the accuracy of the numerical solutions via
convergence tests. We limit ourselves to the case of vanishing plasma pressure
in this paper, which results in a Dirac $\delta$-function current singularity.
For more general cases with a non-vanishing pressure gradient, a
Pfirsch–Schlüter current density that diverges algebraically towards the
resonant surfaces could arise in addition to the $\delta$-function current
singularities.[6] We leave the Pfirsch–Schlüter current singularity to a
future study.

This study employs two numerical codes: (1) a flux conserving Grad-Shafranov
(GS) solver [12] and (2) the Stepped Pressure Equilibrium Code (SPEC).[13]

The GS solver assumes a continuum of nested flux surfaces, which precludes
magnetic island formation. By prescribing the toroidal (i.e., out-of-plane)
and poloidal (i.e., in-plane) fluxes, the geometry of flux surfaces determines
the magnetic field. The geometry of the flux surfaces is described by a
mapping from coordinate space to physical space. The numerical implementation
discretizes the mapping with a Chebyshev-Fourier pseudospectral method,[14,
15] where the residual MHD force $\boldsymbol{J}\times\boldsymbol{B}-\nabla p$
is calculated on a set of collocation points. Here, we use standard notations
for the magnetic field ($\boldsymbol{B}$), the electric current density
($\boldsymbol{J}$), and the plasma pressure ($p$). The collocation points are
uniformly spaced along the Fourier ($y$) direction and correspond to the
interior Chebyshev–Lobato points along the $x$ direction. The mapping is
iteratively updated by an energy descent algorithm, similar to that of VMEC,
until the residual MHD force is below a threshold.

Previously, numerical solutions of the ideal HKT problem from the GS solver
have been tested, showing agreement with the solutions of a fully Lagrangian
solver[16] and analytic solutions obtained with an asymptotic boundary-layer
analysis.[17] However, the accuracy of the GS solution has not been fully
assessed and quantified. The Lagrangian solver does not yield a converged
magnetic field at the resonant surface and therefore cannot facilitate
quantification of errors. The boundary-layer analytic solution also cannot be
used to assess the accuracy of the GS solution, because it is approximate and
not exact.

To further assess the accuracy of the GS solution, it is not sufficient to
rely on self-convergence. Even if the GS solution converges as the resolution
increases, there is no guarantee that the converged solution is correct. To
address this issue, we employ SPEC as an independent solver to benchmark the
GS solver. Another motivation for employing SPEC in this study is that SPEC
can handle a much broader class of 3D configurations. If SPEC can obtain
approximate solutions to the ideal HKT problem, it can potentially be applied
to more complicated 3D problems involving multiple resonant surfaces.

The SPEC code solves for multi-region relaxed magnetohydrodynamic (MRxMHD)
equilibria.[13] MRxMHD does not assume a continuum of nested flux surfaces.
Instead, the physical domain is divided into nested regions. In each region,
the magnetic field relaxes to a Taylor state,[18] i.e., a Beltrami field
satisfying the condition $\nabla\times\boldsymbol{B}=\mu\boldsymbol{B}$, where
$\mu$ is a constant, while conserving magnetic helicity as well as the
poloidal and the toroidal magnetic fluxes. Force-balance conditions are
enforced across the interfaces between adjacent regions. Within each MRxMHD
region, formation of magnetic islands and stochastic field line regions is
allowed;111However, note that stochastic field line regions are not possible
for the 2D HKT problem even when magnetic reconnection is allowed; only
magnetic islands are possible. and the interfaces between MRxMHD regions serve
as ideal flux surfaces that prevent the magnetic field from relaxing to a
global Taylor state. MRxMHD can be viewed as a bridge between Taylor’s
relaxation theory and ideal MHD. When there is only one region in the entire
domain, MRxMHD is equivalent to Taylor’s relaxation. On the other hand, in the
limit of an infinite number of regions such that the ideal interfaces become a
continuum, it has been shown that MRxMHD approaches ideal MHD under some
conditions.[20, 21]

If these conditions hold, we expect SPEC solutions to approach the ideal MHD
solution as the number of regions increases. Hence, we should be able to use
SPEC solutions with a large number of regions to benchmark the GS solutions.
However, Loizu _et al._[22] previously studied a similar problem of imposing
an $m=2$, $n=1$ perturbation on a cylindrical screw pinch with SPEC and
concluded that a minimal finite jump, approximately proportional to the
perturbation amplitude, in the rotational transform across the resonant
surface is a _sine qua non_ condition for the existence of a solution. Because
the HKT problem has a continuous rotational transform, the _sine qua non_
condition raises the question of whether the solutions previously obtained
with the fully Lagrangian solver and the GS solver can also be obtained by
SPEC. As it turns out, in this study we find that SPEC can actually obtain
solutions to the HKT problem without requiring a discontinuous rotational
transform; therefore, the previous interpretation of the _sine qua non_
condition as a necessary condition for the existence of a solution is not
valid.

This paper is organized as follows. In Sec. II, we briefly describe the flux
preserving formulation of the GS equation and review the linear and nonlinear
solutions of the HKT problem. In Sec. III, we present numerical solutions and
convergence tests from the GS solver. In Sec. IV, the numerical solutions and
convergence tests from SPEC are presented. In Sec. V, we further examine the
nature of the singular solution and discuss possible reasons for why SPEC
failed to find solutions in Ref. [22] when the _sine qua non_ condition was
not satisfied, as well as the correct interpretation of the _sine qua non_
condition. We present a case when the _sine qua non_ condition is marginally
satisfied to demonstrate how that affects the nature of the solution. Finally,
we conclude and discuss future perspectives in Sec. VI.

##  II Grad–Shafranov Formulation of the Hahm-Kulsrud-Taylor Problem

Two-dimensional MHD equilibria in Cartesian geometry satisfy the Grad-
Shafranov equation

$\nabla^{2}\psi=-\frac{dP}{d\psi},$ (1)

where

$P=p+\frac{B_{z}^{2}}{2}$ (2)

is a function of $\psi$. Here, the Cartesian coordinate $z$ is the direction
of translational symmetry. The flux function $\psi$ determines the
perpendicular components of the magnetic field through the relation

$\boldsymbol{B}_{\perp}=\boldsymbol{\hat{z}}\times\nabla\psi.$ (3)

Both the out-of-plane component $B_{z}$ and the plasma pressure $p$ are
functions of $\psi$. The component $B_{z}$ is determined by the conservation
of magnetic flux. In this study, we set $p$ equal to zero.

With the magnetic fluxes prescribed, the magnetic field is determined by the
geometry of the flux surfaces. We can label the flux surfaces with an
arbitrary variable, and a convenient choice is to use the initial positions
$x_{0}$ of flux surfaces before the boundary perturbation is imposed. The flux
surfaces are described by a mapping from $\left(x_{0},y\right)$ to
$\left(x,y\right)$ via a function $x\left(x_{0},y\right)$. Using the chain
rule, we can express the partial derivatives with respect to the Cartesian
coordinates in terms of the partial derivatives with respect to the
coordinates $\left(x_{0},y\right)$:

$\left(\frac{\partial}{\partial x}\right)_{y}=\dfrac{1}{\partial x/\partial
x_{0}}\frac{\partial}{\partial x_{0}},$ (4) $\left(\frac{\partial}{\partial
y}\right)_{x}=\frac{\partial}{\partial y}-\dfrac{\partial x/\partial
y}{\partial x/\partial x_{0}}\frac{\partial}{\partial x_{0}}.$ (5)

Here, the subscripts of the partial derivatives on the left-hand side indicate
the coordinates that are held fixed; the partial derivatives on the right-hand
side are with respect to the $\left(x_{0},y\right)$ coordinates. Hereafter,
partial derivatives are taken to be with respect to the $\left(x_{0},y\right)$
coordinates by default, unless otherwise indicated by the subscripts.

Using these relations, the Cartesian components of the in-plane magnetic field
are given by

$B_{x}=-\left(\frac{\partial\psi}{\partial y}\right)_{x}=\dfrac{\partial
x/\partial y}{\partial x/\partial x_{0}}\frac{d\psi}{dx_{0}}$ (6)

and

$B_{y}=\left(\frac{\partial\psi}{\partial x}\right)_{y}=\dfrac{1}{\partial
x/\partial x_{0}}\frac{d\psi}{dx_{0}}.$ (7)

The out-of-plane component $B_{z}$ is determined by conservation of magnetic
flux as

$B_{z}\left(x_{0}\right)=\frac{B_{z0}\left(x_{0}\right)}{\left\langle\dfrac{\partial
x}{\partial x_{0}}\right\rangle}.$ (8)

Here, $B_{z0}$ is the initial $z$-component of the magnetic field; the flux
surface average $\left\langle f\right\rangle$ is defined as

$\left\langle
f\right\rangle\equiv\frac{1}{L}\int_{0}^{L}f\left(x_{0},y\right)dy$ (9)

for an arbitrary function $f\left(x_{0},y\right)$, with
$y\in$$\left[0,L\right]$ being the domain of the system along the $y$
direction. The out-of-plane component of the current density is given by

$J_{z}=\nabla^{2}\psi=\left(\frac{d\psi}{dx_{0}}\right)^{-1}\frac{\partial}{\partial
x_{0}}\left(\frac{B_{x}^{2}+B_{y}^{2}}{2}\right)-\frac{\partial
B_{x}}{\partial y},$ (10)

and the GS equation can be written as

$H=-\left(\frac{d\psi}{dx_{0}}\right)^{-1}\frac{\partial}{\partial
x_{0}}\left(\frac{B_{x}^{2}+B_{y}^{2}}{2}+P\right)+\frac{\partial
B_{x}}{\partial y}=0.$ (11)

The residual MHD force is given by

$\boldsymbol{F}\equiv H\nabla\psi.$ (12)

To obtain the solution, we can use $F_{x}=HB_{y}$ to push the flux surfaces
along the $x$ direction, subjected to a friction force to damp the energy
until the system settles down to an equilibrium.

For the HKT problem, we consider an initial force-free equilibrium

$\boldsymbol{B}_{0}=x_{0}\boldsymbol{\hat{y}}+\sqrt{B_{0}^{2}-x_{0}^{2}}\boldsymbol{\hat{z}}$
(13)

in the domain $x_{0}\in\left[-a,a\right]$ and $y\in\left[0,L\right]$, where
the $y$ direction is assumed to be periodic. The corresponding in-plane flux
function is $\psi=x_{0}^{2}/2$. We impose a sinusoidal perturbation on the
boundary that deforms $x=\pm a$ to
$x=\pm\left(a+\delta\cos\left(ky\right)\right)$ and let the system evolve
under the constraints of ideal MHD to a new equilibrium.

For a small boundary perturbation, we may linearize the GS equation in terms
of the displacements of the flux surfaces $\xi\left(x_{0},y\right)\equiv
x\left(x_{0},y\right)-x_{0}$. To the leading order in $\xi$, the magnetic
field components are

$B_{x}\simeq\frac{\partial\xi}{\partial y}\frac{d\psi}{dx_{0}},$ (14)
$B_{y}\simeq\left(1-\partial\xi/\partial x_{0}\right)\frac{d\psi}{dx_{0}},$
(15)

and

$B_{z}\simeq B_{z0}\left(1-\left\langle\frac{\partial\xi}{\partial
x_{0}}\right\rangle\right).$ (16)

The linearized GS equation now reads

$\frac{\partial}{\partial
x_{0}}\left(\left(\frac{d\psi}{dx_{0}}\right)^{2}\frac{\partial\xi}{\partial
x_{0}}+B_{z0}^{2}\left\langle\frac{\partial\xi}{\partial
x_{0}}\right\rangle\right)+\frac{\partial^{2}\xi}{\partial
y^{2}}\left(\frac{d\psi}{dx_{0}}\right)^{2}=0.$ (17)

For the HKT problem with $\psi=x_{0}^{2}/2$ and the boundary condition
$\xi\left(\pm a,y\right)=\pm\delta\cos\left(ky\right)$, if we adopt the ansatz
$\xi=\bar{\xi}(x_{0})\cos(ky)$, then $\left\langle\partial\xi/\partial
x_{0}\right\rangle=0$ and the linearized GS equation reduces to

$\frac{d^{2}}{dx_{0}^{2}}\left(x_{0}\bar{\xi}\right)-k^{2}x_{0}\bar{\xi}=0.$
(18)

The general solution of Eq. (18) is a linear superposition of two independent
solutions

$\bar{\xi}=c_{1}\frac{\sinh\left(k\left|x_{0}\right|\right)}{x_{0}}+c_{2}\frac{\cosh\left(kx_{0}\right)}{x_{0}},$
(19)

and the boundary condition $\bar{\xi}(\pm a)=\pm\delta$ requires

$\delta=\frac{c_{1}\sinh\left(ka\right)+c_{2}\cosh\left(ka\right)}{a}.$ (20)

We can immediately see that the linear solution is problematic near the
resonant surface at $x_{0}=0$. The divergence of
$\cosh\left(kx_{0}\right)/x_{0}$ at $x_{0}=0$ suggests that the coefficient
$c_{2}$ must be set to zero, and the boundary condition (20) then determines
the coefficient $c_{1}=a\delta/\sinh\left(ka\right)$. However, the limit that
$\lim_{x_{0}\to 0}\sinh\left(k\left|x_{0}\right|\right)/x_{0}=k$ yields
$x\simeq x_{0}+\left(ka\delta/\sinh\left(ka\right)\right)\cos\left(ky\right)$
in the vicinity of $x_{0}=0$, leading to overlap of flux surfaces when
$\left|x_{0}\right|\leq ka\delta/\sinh\left(ka\right)$, which amounts to a
physical inconsistency and is unpermitted. Therefore, within an inner region
$\left|x_{0}\right|\lesssim\mathcal{O}\left(ka\delta/\sinh\left(ka\right)\right)$,
the linear solution is not valid and we must consider the nonlinear solution.

The nonlinear solution of the inner region was first derived by Rosenbluth,
Dagazian, and Rutherford (hereafter RDR) for the ideal internal kink
instability [23] and was later adapted to the HKT problem.[24, 17] Because
$d\psi/dx_{0}\to 0$ in the inner region, the dominant balance of the GS
equation (11) is approximately given by

$\frac{\partial}{\partial
x_{0}}\left(\frac{B_{y}^{2}}{2}+P\left(x_{0}\right)\right)=0;$ (21)

here, we have neglected $B_{x}^{2}$ compared to $B_{y}^{2}$ in Eq. (11) by
assuming $\left|\partial x/\partial y\right|\ll 1$. Integrating Eq. (21)
yields

$B_{y}=\dfrac{1}{\partial x/\partial
x_{0}}\frac{d\psi}{dx_{0}}=\operatorname{sgn}\left(\frac{d\psi}{dx_{0}}\right)\sqrt{f(x_{0})+g(y)},$
(22)

where

$f\left(x_{0}\right)=-2P\left(x_{0}\right)+\text{const}$ (23)

and $g(y)$ is an arbitrary function that will be determined later by
asymptotic matching to the outer solution; the
$\operatorname{sgn}\left(d\psi/dx_{0}\right)$ factor comes from the
requirement that $\partial x/\partial x_{0}>0$ must be satisfied to avoid
overlapping flux surfaces. Without loss of generality, we are free to set
$f(0)=0$, and the tangential discontinuity of $B_{y}$ at $x_{0}=0$ is

$\left.B_{y}\right|_{0^{\pm}}=\pm\sqrt{g(y)}.$ (24)

Using the flux function $\psi=x_{0}^{2}/2$ for the HKT problem and integrating
Eq. (22) one more time yields the inner solution of RDR

$x_{RDR}\left(x_{0},y\right)=\int_{0}^{x_{0}}\frac{\left|x^{\prime}\right|}{\sqrt{f(x^{\prime})+g(y)}}dx^{\prime}.$
(25)

Note that the functions $f\left(x_{0}\right)$ and $g\left(y\right)$ are not
independent, but are related through Eq. (23). Here, $P=B_{z}^{2}/2$ and
$B_{z}$ is determined by Eq. (8) with $x\left(x_{0},y\right)$ replaced by
$x_{RDR}\left(x_{0},y\right)$. The resulting relation is cumbersome. To
simplify the problem, we further assume that
$\left|B_{z}\right|\gg\left|B_{\perp}\right|$ (i.e., in the so-called reduced
MHD regime) and replace the constraint (23) by the incompressible constraint
$\left\langle\partial x/\partial x_{0}\right\rangle=1$, yielding

$\left\langle\frac{1}{\sqrt{f(x_{0})+g(y)}}\right\rangle=\frac{1}{|x_{0}|}.$
(26)

Once $g(y)$ is obtained, the constraint (26) then determines
$f\left(x_{0}\right)$.

The function $g(y)$ can be obtained via asymptotic matching to the linear
solution in the outer region. The readers are referred to Ref. [17] for
further detail of adapting the matching method of RDR for the internal kink
mode to the HKT problem. Here, we simply quote the relevant results. The
function $g(y)$ can be obtained by numerically solving an integral
equation,[25] but a good analytic approximation for $g(y)$ is

$g(y)\simeq\frac{4c_{1}^{2}k^{2}}{3}\sin^{8}(ky/2),$ (27)

where $c_{1}$ is the coefficient of the outer solution given by

$c_{1}=\frac{12\sinh(ka)-\sqrt{144\sinh^{2}(ka)-84\delta
k^{2}\cosh(ka)}}{7k^{2}\cosh(ka)}.$ (28)

##  III Numerical Solutions of the Grad-Shafranov Equation

Figure 2: A selection of flux surfaces obtained from the highest resolution
calculation that we use as the reference for convergence tests in this study.

Now we present numerical solutions to the HKT problem obtained by the GS
solver. We set the free parameters of this problem to $B_{0}=10$, $a=1/2$,
$\delta=0.1$, and $k=2\pi/L=1$. We assume a mirror symmetry of the solution
and solve in only half of the domain $x_{0}\in[0,1/2]$.

We perform two sets of numerical calculations. The first set employs a direct
Chebyshev-Fourier pseudospectral discretization of the GS equation. For the
second set, we take advantage of the knowledge of RDR’s analytic solution and
express the geometry of flux surfaces as
$x(x_{0},y)=x_{\text{RDR}}(x_{0},y)+\tilde{x}(x_{0},y)$. Here, to calculate
the analytic solution $x_{RDR}$, we adopt the analytic approximation (27) for
$g(y)$ and numerically solve the incompressible constraint (26) to obtain
$f\left(x_{0}\right)$. We then numerically integrate Eq. (25) to obtain
$x_{RDR}$. We rewrite the GS equation in terms of the deviation $\tilde{x}$
from the RDR solution and implement a special version of the GS solver for
this formulation. Because the analytic solution accounts for most of the
singular behavior near the resonant surface, the accuracy of the second set of
solutions is substantially improved.

We use $N_{y}=512$ collocation points along the $y$ direction to ensure that
most of the numerical errors are due to the discretization along the $x_{0}$
direction. We then test the convergence of the numerical solution by
increasing the number of Chebyshev collocation points $N_{x}$. We perform
calculations with $N_{x}=8$, 16, 32, 64, and 128. The Chebyshev collocation
points cluster near the edges of the domain, with the shortest distance
between the collocation points scales as $1/N_{x}^{2}$. For $N_{x}=128$, the
closest collocation point is at $x_{0}=7.5\times 10^{-5}$. Due to the lack of
a perfectly precise solution for the convergence test, we take the most
accurate numerical solution available as a substitute. For that purpose, the
$N_{x}=128$ solution from the second set (with the subtraction of the RDR
solution) serves as the reference.

Our primary diagnostics for the convergence test are: (a) the discontinuity of
magnetic field at the resonant surface $\left.B_{y}\right|_{0^{+}}$
($\left.B_{y}\right|_{0^{-}}=-\left.B_{y}\right|_{0^{+}}$ from symmetry); and
(b) the geometry of a selection of flux surfaces. For the latter, we use the
flux surfaces labeled by $x_{0}=1/4$, $1/8$, $1/16$,$\ldots$, $1/1024$. This
set of flux surfaces is shown in Fig. 2. We quantify the errors of a solution
by the $L_{2}$ norms of the differences of relevant quantities relative to the
reference solution. Specifically, we use

$\left\|\left.\Delta
B_{y}\right|_{0^{+}}\right\|_{2}\equiv\left\langle\left(\left.B_{y}\right|_{0^{+}}-\left.B_{y}^{ref}\right|_{0^{+}}\right){}^{2}\right\rangle^{1/2}$
(29)

and

$\left\|\Delta
x\right\|_{2}\equiv\left\langle\left(x-x^{ref}\right){}^{2}\right\rangle^{1/2},$
(30)

where the flux surface average is defined in Eq. (9).

The calculation of $B_{y}$ using Eq. (7) fails at the resonant surface,
because both the denominator and the numerator approach zero. To obtain
$\left.B_{y}\right|_{0^{+}}$, we perform a polynomial extrapolation using the
barycentric formula[26] with values of $B_{y}$ on all the collocation points
other than $x_{0}=0$. Additionally, because the flux surfaces of choice for
the convergence test do not coincide with the Chebyshev collocation points, we
have to perform a polynomial interpolation to determine their geometry.

Figure 3: The magnetic field $B_{y}|_{0^{+}}$ at the lower boundary of the
computational domain obtained from the GS solver without RDR subtraction. For
reference, the dotted line shows the RDR solution.

Figure 3 shows $\left.B_{y}\right|_{0^{+}}$ from the GS solver with increasing
grid resolutions. For $N_{x}=8$, we can see that $\left.B_{y}\right|_{0^{+}}$
becomes negative near $y=0$ and $y=2\pi$. This is a numerical error due to
discretization and extrapolation, as the true solution should remain positive
and only becomes zero at $y=0$ and $y=2\pi$. As $N_{x}$ increases, the
solution quickly converges and the curves are virtually on top of one another
when $N_{x}\geq 16$. The values near $y=0$ and $y=2\pi$ remain slightly
negative, but the magnitude rapidly decreases as $N_{x}$ increases. For the
second set of solutions with RDR subtraction, the curves virtually overlap
with each other for all the cases we have done (not shown). The dotted line in
Fig. 3 shows the RDR solution. We can see that although the RDR solution is a
good approximation, there is a visible difference between the RDR solution and
the converged GS solution.

Figure 4 shows the convergence of $\left.B_{y}\right|_{0^{+}}$ errors for both
sets of solutions. We can see that applying RDR subtraction reduces the errors
by approximately one order of magnitude, but the overall convergence rates are
similar for both sets of solutions. Likewise, the convergence of flux surface
errors is shown in Figure 5 for both sets of solutions. Evidently, flux
surfaces closer to the resonant surface are more difficult to solve
accurately. Again, the RDR subtraction reduces the errors by approximately an
order of magnitude, but the overall convergence rate remains similar.

Note that the data points for $N_{x}=128$ with RDR subtraction are missing in
Figures 4 and 5, because that solution serves as the reference. The
convergence tests provide a base for estimating the errors of the reference
solution. Because the same reference solution will also be used for the
convergence test of SPEC solutions, it is important to ensure that the
reference solution is sufficiently accurate. By extrapolating the trends in
Fig. 4 and Fig. 5(a), we estimate the reference solution’s error of
$B_{y}|_{0^{+}}$ to be smaller than $10^{-6}$, error of the flux surface
labeled by $x_{0}=1/4$ smaller than $10^{-9}$, and error of the flux surface
labeled by $x_{0}=1/1024$ smaller than $10^{-6}$.

Figure 4: Convergence of $\left.B_{y}\right|_{0^{+}}$ errors from the GS
solvers with and without subtracting the RDR solution.  Figure 5: Convergence
of flux surface errors from the GS solvers with and without subtracting the
RDR solution.

##  IV SPEC solutions

Now we continue with the SPEC solutions to the HKT problem. Here we also
present the results from two sets of numerical calculations. For the first
set, the initial positions of interfaces between volumes are uniformly spaced
before the boundary perturbation is imposed. We start from the number of
volumes $N_{vol}=2$, then increase to $N_{vol}=4$, $8$, up to $N_{vol}=128$.
For the second set of calculations, we explore the possible advantages of
packing more volumes near the resonant surface. Because the best strategy for
packing volumes is not a priori clear, we adopt a procedure of refining only
the nearest volume to the resonant surface to see how SPEC performs under this
extreme scenario of local refinement. The procedure goes as follows: We start
from $N_{vol}=2$. At each level of refinement, the volume adjacent to the
resonant surface is divided into two equal volumes. In this way, we go up to
an “effective” $N_{vol}=512$, meaning that the smallest volume is $1/512$ of
the domain, while the actual number of volumes is $N_{vol}=10$. The interfaces
between the volumes for the highest resolution case of the second set exactly
correspond to the flux surfaces we use for convergence tests shown in Fig. 2.

We test the convergence of the two sets of SPEC solutions as the number of
volumes increases, using the highest resolution GS solution as the reference.
The number of Fourier harmonics along the $y$ direction is 48 for all the SPEC
calculations presented here.

When SPEC finds a solution, it is not guaranteed that the ideal interfaces
between volumes will not overlap with each other. Overlapping ideal interfaces
are not permitted on physical grounds, but they do occasionally occur in SPEC
solutions, especially for those interfaces close to the resonant surface, and
this will cause the SPEC algorithm to crash. Because SPEC uses Newton’s method
to find the solution, having a good initial guess is crucial. A useful
approach to overcome the problem of overlapping ideal interfaces is to start
from a small boundary perturbation, find the solution, then use the solution
as the initial guess for a slightly increased boundary perturbation. This
process is repeated until the full amplitude of boundary perturbation is
reached.

Figure 6: Convergence of $B_{y}|_{0^{+}}$ of SPEC solutions as the number of
uniformly spaced volumes increases. Figure 7: Ideal interfaces between MRxMHD
volumes (red) and samples of Poincaré plot in each volume (black). The left-
hand side shows the case with eight volumes, and the right-hand side shows the
case with sixteen volumes.

Figure 6 shows $\left.B_{y}\right|_{0^{+}}$ from SPEC using uniformly-spaced
volumes. Similar to the GS solutions shown in Figure 3, the values of
$\left.B_{y}\right|_{0^{+}}$ in SPEC solutions also become negative near $y=0$
and $y=2\pi$, but the magnitude rapidly decreases as the number of volumes
increases. The reason for negative $\left.B_{y}\right|_{0^{+}}$ is the
presence of residual magnetic islands near the resonant surface,[11] as we can
see in Figure 7. Here, the red lines are the ideal interfaces and the black
dots represent samples of the Poincaré plot from field line tracing. The left-
hand-side of the figure shows the $N_{vol}=8$ case, while the right-hand-side
shows the $N_{vol}=16$ case. The Poincaré plot reveals the residual islands in
the lower left and the lower right corners. The size of the island decreases
as $N_{vol}$ increases from $8$ to $16$. This trend continues as $N_{vol}$
further increases, resulting in the decrease of the magnitude of negative
$\left.B_{y}\right|_{0^{+}}$.

Figure 8: Convergence of the $\left.B_{y}\right|_{0^{+}}$ errors of SPEC
solutions for cases of packed volume and uniformly-spaced volumes. Figure 9:
Convergence of the flux surface errors of SPEC solutions. Panel (a) shows the
cases of uniformly-spaced volumes, and panel (b) shows the cases of packed
volumes near the resonant surface. Note that $N_{vol}$ for the packed cases
corresponds to the “effective” number of volumes as discussed in the text, not
the actual number of volumes. Because the outer region never gets refined for
the packed cases, we do not expect the solutions to approach the ideal MHD
solution even in the limit of $N_{vol}\to\infty$.

Figure 8 shows the convergence of the errors of $\left.B_{y}\right|_{0^{+}}$
as $N_{vol}$ increases, for both sets of SPEC solutions. Here, for the cases
of packed volumes, $N_{vol}$ corresponds to the “effective” number of volumes
as discussed above. We can see that the errors from both sets are nearly
identical for the same $N_{vol}$, even though the volumes far away from the
resonant surface are much coarser for the packed cases. The errors
approximately scale as $N_{vol}^{-1.5}$ for both uniform and packed cases.
This finding suggests that $\left.B_{y}\right|_{0^{+}}$ may not strongly
depend on the accuracy in the outer region.

On the other hand, the effects of inadequate resolution in the outer region
are evident in the convergence of flux surface errors, shown in Figure 9.
Here, the errors consistently scale as $N_{vol}^{-3}$ for cases of uniform
volumes. For cases of packed volumes, although the errors initially decrease
as $N_{vol}^{-3}$, the trend eventually flattens as $N_{vol}$ further
increases. It is possible that the stalling of convergence in the outer region
may eventually affect the convergence of $\left.B_{y}\right|_{0^{+}}$ for the
packed cases. We can see that the last point of packed cases in Figure 8
exhibits some deviation from the $N_{vol}^{-1.5}$ scaling. Another possible
reason for the deviation is that SPEC solutions only use 48 Fourier modes,
which may not be sufficient to accurately represent the flux surfaces near the
resonant surface.

The results of packed-volume solutions show the effectiveness of local
refinement, even when using the extreme refinement scenario adopted here. A
better strategy in practical applications would be to refine over the entire
domain while placing more volumes near the resonant surface. Nonetheless, our
results indicate that a good approximation of the resonant singular current
density may be obtained even with relatively coarse volumes away from the
resonant surface.

##  V Discussion

###  V.1 Nature of the singular solution

The agreement between the solutions of the GS solver and SPEC suggests that
both codes are approaching the true solution of the HKT problem as the
resolution (or number of volumes) increases. Now we further examine the nature
of the singular solution.

The finite tangential discontinuity $\left.B_{y}\right|_{0^{\pm}}$ arises from
a continuous initial magnetic field through the compression of the space
between flux surfaces, which is evident from the flux surfaces shown in Figure
2. As we can infer from the RDR solution (25), for flux surfaces sufficiently
close to the resonant surface such that the condition

$f\left(x_{0}\right)\ll g(y)$ (31)

is satisfied, we have

$x\left(x_{0},y\right)\simeq\frac{x_{0}^{2}}{\sqrt{g(y)}}.$ (32)

Because $f(0)=0$ and $g(y)\simeq\left(4c_{1}^{2}k^{2}/3\right)\sin^{8}(\pi
y/L)$, the condition (31) will eventually be satisfied for sufficiently small
$x_{0}$ for all $y$ except at $y=0$ and $y=L$, but the transition to the
quadratic mapping $x\sim x_{0}^{2}$ occurs at different $x_{0}$ for different
$y$. To compensate for the strong compression of the quadratic mapping, the
“downstream” regions of flux surfaces near $y=0$ and $y=L$ have to bulge
outward to maintain approximate incompressibility.

Now we show that the flux surfaces sufficiently close to the resonant surface
satisfy a similarity relation near the downstream region after a proper
rescaling. To reveal the rescaling rules, we first need to establish the
behavior of $f\left(x_{0}\right)$ near $x_{0}=0$. When the function $g(y)$ is
known, the function $f(x_{0})$ can be obtained by solving Eq. (26). Because
$f(x_{0})\to 0$ in the limit $|x_{0}|\to 0$, the function
$1/\sqrt{f(x_{0})+g(y)}$ is localized near $y=0$, $L.$ Hence, in this limit we
can approximate $g\left(y\right)$ by its leading order Taylor expansion,
yielding

$\displaystyle\left\langle\frac{1}{\sqrt{f(x_{0})+g(y)}}\right\rangle$
$\displaystyle\simeq$
$\displaystyle\frac{2}{L}\int_{0}^{\infty}\frac{dy}{\sqrt{f(x_{0})+\left(4c_{1}^{2}k^{2}/3\right)(ky/2)^{8}}}$
$\displaystyle=$
$\displaystyle\frac{2}{\pi^{3/2}}\dfrac{\Gamma(3/8)\Gamma(9/8)}{f(x_{0})^{3/8}\left(4c_{1}^{2}k^{2}/3\right)^{1/8}},$
(33)

where $\Gamma$ is the gamma function.[27] Plugging Eq. (33) into Eq. (26)
yields the leading order behavior of $f(x_{0})$ in the limit
$\left|x_{0}\right|\to 0$:

$f(x_{0})\simeq\left[c_{f}\left|x_{0}\right|\right]^{8/3},$ (34)

where

$\displaystyle c_{f}$
$\displaystyle\equiv\frac{2(3/4)^{1/8}}{\pi^{3/2}}\Gamma(3/8)\Gamma(9/8)\left(c_{1}k\right)^{-1/4}$
$\displaystyle\simeq 0.7735\left(c_{1}k\right)^{-1/4}.$ (35)

Without loss of generality, here we consider $x_{0}\geq 0$. Applying the
leading order approximations of $f\left(x_{0}\right)$ and $g\left(y\right)$
near $x_{0}=0$ and $y=0$ to the RDR solution (25) yields

$x\left(x_{0},y\right)\simeq\int_{0}^{x_{0}}\frac{x^{\prime}}{\sqrt{d_{1}x^{\prime
8/3}+d_{2}y^{8}}}dx^{\prime},$ (36)

where $d_{1}$ and $d_{2}$ are some constants. With a change of variables
$\zeta=x^{\prime}/x_{0}$, equation (36) can be rewritten as

$\frac{x\left(x_{0},y\right)}{x_{0}^{2/3}}\simeq\int_{0}^{1}\frac{\zeta}{\sqrt{d_{1}\zeta^{8/3}+d_{2}\left(y/x_{0}^{1/3}\right)^{8}}}d\zeta.$
(37)

Equation (37) suggests that if we rescale $x$ and $y$ to $x/x_{0}^{2/3}$ and
$y/x_{0}^{1/3}$, the flux surfaces near $(x,y)=(0,0)$ will approximately
coincide. This similarity relation is borne out by our numerical solutions,
shown in Figure 10 for a selection of flux surfaces before and after
rescaling.

The similarity relation implies that the heights of the bulged flux surfaces
in the downstream region scale as $x_{0}^{2/3}$ and the widths scale as
$x_{0}^{1/3}$; the enclosed volumes scale as $x_{0}$, to be consistent with
the incompressible constraint. Therefore, in the limit of $x_{0}\to 0$, the
width of the bulged region becomes narrower and narrower. Because the enclosed
volumes scale as $\sim x_{0}$, the geometry of flux surfaces $x(x_{0},y)$ may
be viewed as approaching a Dirac $\delta$-function $\sim x_{0}\delta(y)$.

Figure 10: The nested flux surfaces near the resonant surface exhibit a
similarity relation. Panel (a) shows a selection of flux surfaces near the
lower-left corner of Figure 2. After rescaling by $x\to x/x_{0}^{2/3}$ and
$y\to y/x_{0}^{1/3}$, the flux surfaces become nearly identical, as shown in
panel (b). This figure uses the reference solution obtained by the GS solver
with RDR subtraction. Using the SPEC solution with packed volumes or the GS
solver without RDR subtraction yields nearly identical curves. Figure 11: The
distorted mesh shows the Lagrangian mapping to the final equilibrium from a
rectangular uniform mesh before imposing the boundary perturbation.

Examining the solution from a Lagrangian perspective provides further insight
to its singular nature. The Lagrangian formulation of ideal MHD describes a
state in terms of the mapping from the initial positions $\boldsymbol{x}_{0}$
of fluid elements to their final positions $\boldsymbol{x}$. The magnetic
field at $\boldsymbol{x}$ is determined by the initial magnetic field
$\boldsymbol{B}_{0}$ at $\boldsymbol{x}_{0}$ and the mapping
$\boldsymbol{x}\left(\boldsymbol{x}_{0}\right)$ via the relation [28, 29]

$\boldsymbol{B}=\dfrac{\boldsymbol{B}_{0}\cdot\dfrac{\partial\boldsymbol{x}}{\partial\boldsymbol{x}_{0}}}{\mathcal{J}},$
(38)

where
$\mathcal{J}=\det\left(\partial\boldsymbol{x}/\partial\boldsymbol{x}_{0}\right)$
is the Jacobian of the mapping.

Although our GS solver is not fully Lagrangian because the mesh can move along
the $x$ direction but not along the $y$ direction, we can reconstruct the full
Lagrangian mapping of fluid elements from the initial to the final state once
the solution is obtained. For each fluid element labeled by
$\left(x_{0},y\right)$ in the final solution, we need to find its initial
position $\left(x_{0},y_{0}\right)$. This “inverse” Lagrangian mapping can be
expressed as a function $y_{0}\left(x_{0},y\right)$. From the conservation of
magnetic flux through an infinitesimal fluid element

$B_{z0}\left(x_{0}\right)dx_{0}\left[\frac{\partial y_{0}}{\partial
y}dy\right]=B_{z}\left(x_{0}\right)\left[\frac{\partial x}{\partial
x_{0}}dx_{0}\right]dy$ (39)

and using Eq. (8) to relate $B_{z0}$ and $B_{z}$, we can calculate

$\frac{\partial y_{0}}{\partial y}=\dfrac{\partial x/\partial
x_{0}}{\left\langle\partial x/\partial x_{0}\right\rangle}$ (40)

and integrate it along each constant-$x_{0}$ contour to obtain
$y_{0}(x_{0},y)$.

Figure 11 visualizes how a rectangular uniform mesh in the initial state is
deformed by the Lagrangian mapping in the final state. We can see that the
mapping is highly distorted near the resonant surface. All the vertical mesh
lines in the initial state now converge towards the lower corners in the final
state, and the single point $\left(x_{0},y_{0}\right)=\left(0,L/2\right)$ is
stretched to an entire line of the lower boundary. This result strongly
suggests that the solution we find here for the HKT problem can only be
approached, but cannot be reached by ideal MHD evolution described via smooth,
diffeomorphic Lagrangian mapping.[30]

Now we discuss some limitations of our present methods in tackling the
$\delta$-function singularities. At first sight, the pseudospectral method
employed by the GS solver may seem ill-suited for problems with
discontinuities. However, note that the primary variable to describe the
solution is the geometry of flux surfaces represented by the mapping
$x\left(x_{0},y\right)$, which is not discontinuous. Although the magnetic
field does become discontinuous, we only take derivatives on the total
pressure $B^{2}/2+p$, which is continuous, when evaluating the residual force.
For that reason, the pseudospectral method does not perform poorly because of
the discontinuous magnetic field. In this study, because we assume a mirror
symmetry and solve for half of the domain, the magnetic field discontinuity is
not present within the computational domain and therefore does not pose a
problem. However, the GS solver works fine even when we do not assume the
symmetry, provided that the collocation points do not fall on (or very close
to) the resonant surface.

Although the mapping $x\left(x_{0},y\right)$ is continuous, it appears to
become non-differentiable when the rotational transform is continuous. Even
though the exact form of $x\left(x_{0},y\right)$ is not known, we may use the
RDR solution, Eq. (25), as a proxy. The function $x_{RDR}\left(x_{0},y\right)$
is infinitely differentiable along the $y$ direction over the entire domain,
and is infinitely differentiable along $x$ everywhere except at the point
$(x_{0},y)=(0,0)$ (and also $(0,L)$ because of the periodicity). Because
$x\sim x_{0}^{2/3}$ when $y=0$, the partial derivative of $x_{RDR}$ along
$x_{0}$ diverges at $(x_{0},y)=(0,0)$ as $\partial x_{RDR}/\partial
x_{0}|_{y=0}\sim x_{0}^{-1/3}$. This singular behavior leads to the non-
smoothness of the flux surfaces near the resonant surface. Consequently, the
convergence rate of the GS solver is algebraic with respect to the number of
collocation points (see Figures 4 and 5), as opposed to an exponential
convergence we usually expect from a pseudospectral method. In contrast, when
applying to a problem that satisfies the _sine qua non_ condition in Sec. V.2,
the GS solver can achieve much more rapid convergence (see Figures 13 and 15).

We can appreciate the non-smoothness of HKT flux surfaces near the resonant
surface through the similarity relation we discussed earlier. Because the
width of the bulged region scales as $x_{0}^{1/3}$, when we increase the
resolution along the $x$ direction, we need to increase the resolution along
$y$ direction as well to resolve the localized structure. For the GS solver,
since the closest Chebyshev collocation point to the resonant surface has
$x_{0}\propto 1/N_{x}^{2}$, roughly speaking, the resolution in $y$ needs to
scale as $N_{y}\sim N_{x}^{2/3}$ to resolve the localized structures. A
similar requirement also applies to SPEC when the number of volumes increases.
Therefore, the Fourier representation employed by both solvers is inefficient
for ideal flux surfaces near the resonant surface. A possible remedy is to
employ alternative basis functions for the flux surfaces. The version of GS
solver with RDR subtraction effectively uses the RDR solution as one of the
basis functions. However, although subtracting the RDR solution significantly
improves the accuracy of the GS solver, it does not completely remove the
effect of the singularity and the convergence rate remains similarly
algebraic. The convergence rate could potentially be further improved by
adopting a more accurate $g\left(y\right)$ in the RDR solution (25), either by
numerically solving the RDR integral equation [25] or by dynamically solving
$g\left(y\right)$ as a part of the solver.

Note that this singular behavior of flux surfaces near the resonant surface
only arises when we try to obtain the ideal MHD solution. SPEC, which
implements MRxMHD, is not an ideal MHD equilibrium solver with nested flux
surfaces by design. By changing the number of volumes, SPEC allows a
transition from Taylor relaxation to ideal MHD. When modeling non-ideal
plasmas that allow magnetic islands and regions of stochastic field lines with
MRxMHD, an active area of research is to understand where the ideal interfaces
should be placed and when an ideal interface should be removed.[21] The
presence of a strong current sheet on an ideal interface is an indication that
the interface should be removed. If we remove the ideal interface at $x_{0}=0$
and allow reconnection, the singular behavior of flux surfaces may no longer
be a problem.

###  V.2 Reinterpreting the _sine qua non_ condition in Loizu _et al._
(2015)[22]

We mention in the Introduction that Loizu _et al._[22] previously studied an
$m=2$, $n=1$ perturbation on a cylindrical screw pinch with SPEC and concluded
that a minimal finite jump in the rotational transform is necessary for the
existence of a solution. This finding motivated Loizu _et al._ to call the
minimal finite jump a _sine qua non_ condition. However, in the present study,
we show that SPEC actually can find solutions for the HKT problem, which has a
continuous rotational transform, provided that Newton’s method is initialized
with care. Therefore, the previous interpretation of the _sine qua non_
condition by Loizu _et al._ is incorrect. To further clarify the issue, it is
instructive to discuss the _sine qua non_ condition in the context of the HKT
problem. A similar discussion can also be found in Sec. 3.3 of Ref. [31].

Instead of a continuous initial magnetic field, let us now suppose that the
initial field has a finite discontinuity at $x_{0}=0$:

$B_{y0}=x_{0}\pm b.$ (41)

Here, we take the plus sign for $x_{0}>0$ and the minus sign for $x_{0}<0$.
The discontinuity parameter $b$ provides a finite jump in the rotational
transform. In the limit $b\to 0$, the original HKT problem is recovered.

For this modified HKT problem, the linearized GS equation (17) becomes

$\frac{d^{2}}{dx_{0}^{2}}\left(\left(x_{0}\pm
b\right)\bar{\xi}\right)-k^{2}\left(x_{0}\pm b\right)\bar{\xi}=0.$ (42)

With boundary conditions $\bar{\xi}(0)=0$ and $\bar{\xi}(\pm a)=\pm\delta$,
the solution is

$\bar{\xi}=\frac{(a+b)\delta}{\sinh(ka)}\frac{\sinh\left(kx_{0}\right)}{x_{0}\pm
b}.$ (43)

The geometry of perturbed flux surfaces up to the linear order is given by
$x=x_{0}+\bar{\xi}\cos(ky)$. To prevent overlapping of flux surfaces requires
$\partial x/\partial x_{0}>0$, which amounts to

$\left|\frac{d\bar{\xi}}{dx_{0}}\right|=\frac{(a+b)\delta}{\sinh(ka)}\left|k\frac{\cosh\left(kx_{0}\right)}{x_{0}\pm
b}-\frac{\sinh\left(kx_{0}\right)}{(x_{0}\pm b)^{2}}\right|<1$ (44)

for the linear perturbation. Since

$\left|\frac{d\bar{\xi}}{dx_{0}}\right|=\frac{(a+b)k\delta}{b\sinh(ka)}\left|1\mp\frac{2x_{0}}{b}+O\left(x_{0}^{2}\right)\right|$
(45)

in the vicinity of $x_{0}=0$, it is sufficient to ensure that
$\left|d\bar{\xi}/dx_{0}\right|<1$ at $x_{0}=0$. That leads to the _sine qua
non_ condition for the HKT problem:

$b>b_{\min}=\frac{ka\delta}{\sinh\left(ka\right)-k\delta}.$ (46)

The _sine qua non_ condition ensures that the flux surfaces of the linear
solution do not overlap. However, not satisfying the _sine qua non_ condition
does not imply the nonexistence of a solution; it simply means that a
nonlinear solution must be sought non-perturbatively. The RDR solution
demonstrates how an approximate nonlinear solution can be obtained through a
boundary layer analysis and asymptotic matching. The previous
misinterpretation of the _sine qua non_ condition as the necessary condition
for the existence of a solution further led to an erroneous claim that the RDR
solution has a discontinuous rotational transform.[25] This latter mistake has
been corrected by Zhou _et al._[17]

Let us now examine how the _sine qua non_ condition affects the solution. For
the same boundary perturbation with $k=1$, $a=0.5$, and $\delta=0.1$ as
before, the _sine qua non_ condition (46) gives $b_{\min}\simeq 0.119$. In
what follows, we consider the case $b=0.12$ such that the _sine qua non_
condition is marginally satisfied. We numerically calculate the solution with
the GS solver and perform exactly the same convergence tests as before.

Figure 12: The magnetic field discontinuity $B_{y}|_{0^{+}}$ obtained by the
GS solver with different number of collocation points $N_{x}$, for the case
$b=0.12$.  Figure 13: Self convergence test of the $B_{y}|_{0^{+}}$ error for
the case $b=0.12$. The $N_{x}=128$ solution serves as the reference. The error
is dominated by round-off errors when $N_{x}\geq 32$. Figure 14: A zoom-in
view of a selection of flux surfaces for the case $b=0.12$. This figure can be
compared with Fig. 10(a) for the case $b=0$.

Figure 12 shows the magnetic field discontinuity $B_{y}|_{0^{+}}$ for
different number of collocation points $N_{x}$. We can see that all the curves
are visually indistinguishable, even with a resolution as low as $N_{x}=8$.
Note that $B_{y}|_{0^{+}}>0$ everywhere because of the discontinuous
rotational transform.

Figure 13 shows a self convergence test for the $B_{y}|_{0^{+}}$ error. Here,
we use the $N_{x}=128$ solution as the reference. We can see that the error
reaches a level below $10^{-10}$ at $N_{x}=32$. Further increasing the
resolution does not lower the error, suggesting that the error is dominated by
round-off errors when $N_{x}\geq 32$.

Figure 15: Self convergence test of the flux surface errors for the case
$b=0.12$. The $N_{x}=128$ solution serves as the reference. The errors are
dominated by round-off errors when $N_{x}\geq 32$.

Figure 14 shows a zoom-in view of a selection of flux surfaces. These flux
surfaces correspond to the same flux surfaces shown in Fig. 10(a) for the case
with a continuous rotational transform. Comparing these two figures, we can
see that with the discontinuous rotational transform, the space between flux
surfaces is no longer strongly squeezed; consequently, the flux surfaces also
do not strongly bulge out in the “downstream” region near the lower-left
corner. Figure 15 shows a self convergence test of flux surface errors, using
the $N_{x}=128$ case as the reference. Again, the errors appear to be limited
by round-off errors when $N_{x}\geq 32$.

The numerical calculations shown here all use $N_{y}=512$. However, because
the flux surfaces no longer have strongly localized geometric structures, the
same accuracy can be achieved with a lot fewer grid points along $y$. We find
that a similar accuracy can be achieved with $N_{x}=32$ and $N_{y}=64$.

This example demonstrates that the _sine qua non_ condition significantly
alters the nature of the solution, which may have contributed to why SPEC had
no problem finding solutions when the condition was satisfied in Ref. [22].
Because the space between flux surfaces is no longer strongly squeezed, the
Newton’s solver of SPEC is less likely to have overlapped ideal interfaces.
Moreover, because the SPEC calculations in Ref. [22] only use a small number
of Fourier harmonics (the toroidal mode number $n\leq 6$ and the poloidal mode
number $m\leq 3$), the fact that the flux surfaces do not develop localized
structures also helps.

##  VI Conclusions and Future Perspectives

In conclusion, we have demonstrated that with the increase of resolution or
the number of volumes, the GS solver and SPEC both appear to approach the
solution of the ideal HKT problem with a $\delta$-function singularity. Our
result is also the first to show that SPEC can obtain approximate solutions of
the ideal HKT problem without requiring a discontinuous rotational transform
across the resonant surface in the initial condition.

In the previous calculation by Loizu et al.,[22] the sine qua non condition
originated from a breakdown of the linear solution near the resonant surface,
which was misinterpreted as a lack of a solution. This misunderstanding,
compounded with the fact that SPEC uses a Newton method that may fail to find
the solution without a carefully chosen initial guess, led to the erroneous
conclusion that a finite threshold of discontinuous rotational transform is
necessary for the existence of a solution. However, as we have discussed in
Sec. II, the breakdown of the linear solution does not imply a lack of a
solution, but rather that a nonlinear solution must be sought. Furthermore, by
carefully initiating Newton’s method, we have demonstrated that SPEC can
obtain the solution. The present study also calls for reconsideration of the
previous study by Loizu et al.,[22] as well as the benchmark study between
SPEC and VMEC on the same problem. [32]

In future work, it would be prudent to implement in SPEC a steepest descent
algorithm for the energy functional, which should be beneficial for tackling
this and similar problems. For instance, one could first use the more robust
steepest descent algorithm to obtain an approximate solution, then switch to
Newton’s method for more rapid convergence to the final solution.

Following up this work, several further investigations will be pursued in the
future. Some of the present approaches could be adapted to singular current
sheets arising from the ideal internal kink instability [23, 33] and more
general 3D magnetic resonant perturbations. In addition to the
$\delta$-function singularities, the algebraically divergent Pfirsch–Schlüter
current in the presence of a pressure gradient should also be investigated.
Recent studies have shown that ideal current singularities on resonant
surfaces may be eliminated by modifying the plasma boundary.[34, 35] This new
approach could also be investigated with SPEC. Finally, the tendency to form
current sheets is thwarted in real plasmas by non-ideal effects, which will
tend to drive magnetic reconnection, forming magnetic islands or regions of
stochastic field lines when island overlap occurs. An important question of
practical significance is whether the sizes of saturated islands or regions of
stochastic field lines can be predicted from the intensity of current
singularities. [36, 5, 37, 38, 39] If such a relationship can be established,
we may use singularity intensities as a proxy for the sizes of magnetic
islands (or regions of stochastic field lines) in stellarator optimization.

###### Acknowledgements.

This research was supported by the U.S. Department of Energy under Contract
No. DE-AC02-09CH11466 and by a grant from the Simons Foundation/SFARI (560651,
AB). Part of this work has been carried out within the framework of the
EUROfusion Consortium and has received funding from the Euratom research and
training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053.
The views and opinions expressed herein do not necessarily reflect those of
the European Commission. YZ was partly sponsored by Shanghai Pujiang Program
under Grant No. 21PJ1408600. Part of the numerical calculations were performed
with computers at the National Energy Research Scientific Computing Center.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Parker [1994] E. N. Parker,  _Spontaneous Current Sheets in Magnetic Fields_ (Oxford University Press, Inc., 1994).
  * Hirshman and Whitson [1983] S. P. Hirshman and J. C. Whitson, Physics of Fluids 26, 3553 (1983).
  * Garabedian [2002] P. R. Garabedian, Proceedings of the National Academy of Sciences 99, 10257 (2002).
  * Grad [1967] H. Grad, Phys. Fluids 10, 137 (1967).
  * Bhattacharjee  _et al._ [1995] A. Bhattacharjee, T. Hayashi, C. C. Hegna, N. Nakajima, and T. Sato, Phys. Plasmas 2, 883 (1995).
  * Helander [2014] P. Helander, Reports on Progress in Physics 77, 087001 (2014).
  * Biskamp [1993] D. Biskamp,  _Nonlinear Magnetohydrodynamics_ (Cambridge University Press, 1993).
  * Biskamp [2000] D. Biskamp,  _Magnetic Reconnection in Plasmas_ (Cambridge University Press, 2000).
  * Hahm and Kulsrud [1985] T. S. Hahm and R. M. Kulsrud, Phys. Fluids 28, 2412 (1985).
  * Wang and Bhattacharjee [1992] X. Wang and A. Bhattacharjee, Phys. Fluids B 4, 1795 (1992).
  * Dewar  _et al._ [2017] R. L. Dewar, S. R. Hudson, A. Bhattacharjee, and Z. Yoshida, Phys. Plasmas 24, 042507 (2017).
  * Huang  _et al._ [2009] Y.-M. Huang, A. Bhattacharjee, and E. G. Zweibel, Astrophys. J. Lett. 699, L144 (2009).
  * Hudson  _et al._ [2012] S. R. Hudson, R. L. Dewar, G. Dennis, M. J. Hole, M. McGann, G. von Nessi, and S. Lazerson, Physics of Plasmas 19, 112502 (2012).
  * Fornberg [1995] B. Fornberg,  _A Practical Guide to Pseudospectral Methods_ (Cambridge University Press, 1995).
  * Trefethen [2000] L. N. Trefethen,  _Spectral Methods in Matlab_ (SIAM Philadelphia, 2000).
  * Zhou  _et al._ [2016] Y. Zhou, Y.-M. Huang, H. Qin, and A. Bhattacharjee, Phys. Rev. E 93, 023205 (2016).
  * Zhou  _et al._ [2019] Y. Zhou, Y.-M. Huang, A. H. Reiman, H. Qin, and A. Bhattacharjee, Physics of Plasmas 26, 022103 (2019).
  * Taylor [1974] J. B. Taylor, Phys. Rev. Lett. 33, 1139 (1974).
  * Note [1] However, note that stochastic field line regions are not possible for the 2D HKT problem even when magnetic reconnection is allowed; only magnetic islands are possible.
  * Dennis  _et al._ [2013] G. R. Dennis, S. R. Hudson, R. L. Dewar, and M. J. Hole, Physics of Plasmas 20, 032509 (2013).
  * Qu  _et al._ [2021] Z. S. Qu, S. R. Hudson, R. L. Dewar, J. Loizu, and M. J. Hole, Plasma Physics and Controlled Fusion 63, 125007 (2021).
  * Loizu  _et al._ [2015] J. Loizu, S. R. Hudson, A. Bhattacharjee, S. Lazerson, and P. Helander, Physics of Plasmas 22, 090704 (2015).
  * Rosenbluth  _et al._ [1973] M. N. Rosenbluth, R. Y. Dagazian, and P. H. Rutherford, Phys. Fluids 16, 1894 (1973).
  * Boozer and Pomphrey [2010] A. H. Boozer and N. Pomphrey, Physics of Plasmas 17, 110707 (2010).
  * Loizu and Helander [2017] J. Loizu and P. Helander, Phys. Plasmas 24, 040701 (2017).
  * Berrut and Trefethen [2004] J.-P. Berrut and L. N. Trefethen, SIAM Review 46, 501 (2004).
  * Abramowitz and Stegun [1972] M. Abramowitz and I. A. Stegun, eds.,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ , 10th ed. (National Bureau of Standards, 1972).
  * Newcomb [1962] W. A. Newcomb, Nuclear Fusion Supplement, Part 2 , 451 (1962).
  * Zhou  _et al._ [2014] Y. Zhou, H. Qin, J. W. Burby, and A. Bhattacharjee, Phys. Plasmas 21, 102109 (2014).
  * Pfefferlé  _et al._ [2020] D. Pfefferlé, L. Noakes, and Y. Zhou, Plasma Physics and Controlled Fusion 62, 074004 (2020).
  * Zhou [2017] Y. Zhou,  _Variational Integration for Ideal Magnetohydrodynamics and Formation of Current Singularities_ , Ph.D. thesis, Princeton University (2017), arXiv:1708.08523 .
  * Lazerson  _et al._ [2016] S. A. Lazerson, J. Loizu, S. Hirshman, and S. R. Hudson, Physics of Plasmas 23, 012507 (2016).
  * Park  _et al._ [1980] W. Park, D. A. Monticello, R. B. White, and S. C. Jardin, Nucl. Fusion 20, 1181 (1980).
  * Mikhailov  _et al._ [2019] M. Mikhailov, J. Nührenberg, and R. Zille, Nuclear Fusion 59, 066002 (2019).
  * Kim  _et al._ [2020] E. Kim, G. B. McFadden, and A. J. Cerfon, Plasma Physics and Controlled Fusion 62, 044002 (2020).
  * Cary and Hanson [1991] J. R. Cary and J. D. Hanson, Physics of Fluids B: Plasma Physics 3, 1006 (1991).
  * Loizu  _et al._ [2020] J. Loizu, Y.-M. Huang, S. R. Hudson, A. Baillod, A. Kumar, and Z. S. Qu, Physics of Plasmas 27, 070701 (2020).
  * Geraldini  _et al._ [2021] A. Geraldini, M. Landreman, and E. Paul, Journal of Plasma Physics 87 (2021), 10.1017/s0022377821000428.
  * Rodríguez and Bhattacharjee [2021] E. Rodríguez and A. Bhattacharjee, Physics of Plasmas 28, 092506 (2021).