# Linguistic Collapse:
Neural Collapse in (Large) Language Models

Robert Wu
University of Toronto
<EMAIL_ADDRESS>
&Vardan Papyan
University of Toronto
<EMAIL_ADDRESS>

###### Abstract

Neural collapse ($\mathcal{NC}$) is a phenomenon observed in classification
tasks where top-layer representations collapse into their class means, which
become equinorm, equiangular and aligned with the classifiers. These behaviors
– associated with generalization and robustness – would manifest under
specific conditions: models are trained towards zero loss, with noise-free
labels belonging to balanced classes, which do not outnumber the model’s
hidden dimension. Recent studies have explored $\mathcal{NC}$ in the absence
of one or more of these conditions to extend and capitalize on the associated
benefits of ideal geometries. Language modeling presents a curious frontier,
as training by token prediction constitutes a classification task where none
of the conditions exist: the vocabulary is imbalanced and exceeds the
embedding dimension; different tokens might correspond to similar contextual
embeddings; and large language models (LLMs) in particular are typically only
trained for a few epochs. This paper empirically investigates the impact of
scaling the architectures and training of causal language models (CLMs) on
their progression towards $\mathcal{NC}$. We find that $\mathcal{NC}$
properties that develop with scaling are linked to generalization. Moreover,
there is evidence of some relationship between $\mathcal{NC}$ and
generalization independent of scale. Our work therefore underscores the
generality of $\mathcal{NC}$ as it extends to the novel and more challenging
setting of language modeling. Downstream, we seek to inspire further research
on the phenomenon to deepen our understanding of LLMs – and neural networks at
large – and improve existing architectures based on $\mathcal{NC}$-related
properties.

Figure 1: Simultaneous development of the four neural collapse
($\mathcal{NC}$) [1] properties in 230 causal language models trained on
TinyStories [2], alongside improvement in generalization (i.e. validation
performance). Left to right: $\mathcal{NC}_{1}$) within-class (representation)
variability collapse; $\mathcal{GNC}_{2}$) hyperspherical uniformity of class
means; $\mathcal{UNC}_{3}$) uniform duality between class means and
corresponding classifiers; and $\mathcal{NC}_{4}$) agreement between token
(maximum a prior) classifiers and implicit nearest-class center classifiers.
Colored by model size and annotated with coefficient of determination
($R^{2}$).

##  1 Introduction

###  1.1 Neural Collapse

A learning phenomenon known as neural collapse ($\mathcal{NC}$) emerges during
the terminal phase of training neural networks with cross-entropy (CE) loss
for classification [1]. It was originally characterized as the co-occurrence
of the following properties in a model’s top-layer features (also known as
last-layer representations or embeddings) and classifiers:

  1. ($\mathcal{NC}_{1}$)

Within-class variability collapse: Top-layer features collapse to their class
means.

  2. ($\mathcal{NC}_{2}$)

Convergence to a simplex ETF: Class means tend towards equinorm and
equiangular vectors when centered about the global average. The resulting
geometry – known as a simplex equiangular tight frame (ETF) – maximizes
pairwise angles and distances.

  3. ($\mathcal{NC}_{3}$)

Convergence to self-duality: Top-layer classifier vectors converge to their
corresponding class means, and thus also form a simplex ETF.

  4. ($\mathcal{NC}_{4}$)

Nearest decision rule: Linear classifiers approximate a nearest-class center
(NCC) classifier: top-layer embeddings predict the class with the closest mean
(implied by $\mathcal{NC}$1-3).

These behaviors, often associated with improved generalization and robustness
(among other benefits, such as those discussed in §1.4), traditionally
manifest under the following conditions:

  1. C1)

Few classes: The number of classes is upper-bounded by the embedding dimension
plus one: $C\leq d+1$; this is required to construct a perfect simplex ETF.

  2. C2)

Balanced classes: The number of samples is equal across classes:
$N_{c}=N_{c^{\prime}},\forall c\neq c^{\prime}$.

  3. C3)

Noise-free labels: Identical (or very similar) embeddings should belong to the
same class.

  4. C4)

Terminal phase of training (TPT): The model is trained past zero error,
towards zero loss.

Absent these conditions, one does not typically anticipate $\mathcal{NC}$.
Since then, follow-up studies have extended beyond and proposed techniques of
quantifying or achieving $\mathcal{NC}$ (discussed in Section 2).

###  1.2 (Large) Language Models

$\mathcal{NC}$ is a phenomenon observed specifically in classification tasks.
While not traditionally thought of as classifiers, language models – including
large language models (LLMs) – learn to model aleatoric uncertainty, which can
be viewed as stochastic token prediction [3]. For instance, masked language
models (MLMs) such as BERT [4] predict one or several masked tokens within an
input sequence based on the surrounding context. Likewise, autoregressive or
causal language models (CLMs) such as GPT [5] predict the next token in a
sequence given the context of previous tokens. Most of these models are
essentially (pre-)trained by token classification on their vocabularies. This
parallel – also drawn by [6] – raises a few natural questions:

  1. 1.

Does the pre-training stage of a language model exhibit $\mathcal{NC}$?

  2. 2.

How do scaling and other training configurations influence $\mathcal{NC}$ in
(L)LMs?

  3. 3.

To what extent is $\mathcal{NC}$ in (L)LMs correlated with their
generalization abilities?

  4. 4.

Do such correlations, between $\mathcal{NC}$ and improved generalization,
persist independent of the (potential confounders of) model size and training?

To address these, we first examine the specific settings of training CLMs as
they are opposed ($\neg$) to the prerequisites (C1-4, §1.1) for $\mathcal{NC}$
to manifest.

  1. $\neg$C1)

Many classes: The unique tokens found in language modeling vocabularies are
vast, usually numbering in the tens of thousands and far exceeding the hidden
dimension [7].

  2. $\neg$C2)

Imbalanced classes: The distribution of tokens in natural language is
typically very imbalanced [8, 9], as is the case in TinyStories [2], the
dataset we use (Appendix Figure 4). It has an average of 16K samples per class
but a standard deviation of over 32K.

  3. $\neg$C3)

Ambiguous contexts: There may exist very similar or even identical contexts
followed by different tokens in the natural language data [10]. For instance,
over half of the sequences in TinyStories [2] lead with “Once upon a time”,
but only three-quarters of these follow with a comma (“Once upon a time,”). In
other words, there is almost certainly some ambiguity to contend with in our
context embeddings.

  4. $\neg$C4)

Undertraining: Most practical language models (including LLMs) are not trained
for more than a few epochs [11, 12]. Further optimization typically renders
diminishing returns in improving evaluation performance [13] long before any
possible TPT.

###  1.3 Contributions

We train a suite of Transformer-based [14] CLMs111Our models range from 3.4M
(small) to 205M (large), so we inclusively use “CLM” instead of “LLM”. across
a grid of model widths, depths, and training epochs on the TinyStories dataset
[2] to assess the degrees to which $\mathcal{NC}$ properties develop and how
they relate to generalization performance. Our findings (summarized in Figure
1) reveal:

  * •

Emergence of $\mathcal{NC}$ with scale: As we scale model size and training,
the properties of $\mathcal{NC}$ emerge; within-class variability
($\mathcal{NC}_{1}$) and interference are reduced while hyperpsherical
uniformity ($\mathcal{GNC}_{2}$), uniform duality ($\mathcal{UNC}_{3}$), and
classifier agreement ($\mathcal{NC}_{4}$) improve.

  * •

Progression towards hyperspherical uniformity: Class means, while unable to
form a simplex ETF ($\mathcal{NC}_{2}$), nonetheless tend towards uniform
dispersion on a hypersphere, a geometry theorized by [15] and formalized by
[16] as hyperspherical uniformity ($\mathcal{GNC}$2).

  * •

Tendency towards uniform duality: Direct alignment (self-duality) between
class means and classifiers ($\mathcal{NC}_{3}$) does not appear to develop
with scale. Instead, the variability of (mis)alignment across classes
decreases with width and training, suggesting its minimization – which we term
uniform duality ($\mathcal{UNC}_{3}$) – may be more cohesive with
$\mathcal{NC}$.

  * •

Correlation between $\mathcal{NC}$ and generalization: The developments of
$\mathcal{NC}$ properties are correlated with improved validation performance.
We show these correlations to persist even when fixing the (potential
confounders of) model architecture and training by simply varying the random
seeds for initialization and data shuffling. This suggests that $\mathcal{NC}$
is not simply a side-effect of training, but possibly a factor of
generalization in its own right.

###  1.4 Significance

Recently, methods building on $\mathcal{NC}$ have found use in deep learning
at large. We highlight areas such as federated learning [17], graph neural
networks [18], incremental/continual learning [19, 20, 21], meta-learning [21,
22], out-of-distribution detection [23, 24, 25], privacy [26, 27], learning
theory [28, 29, 30, 31, 32, 33] transfer learning [34, 35, 36, 37, 38, 39],
and even LLM prompting [40]. With our results, we aim to extend insights from
such contributions and other related works to the language modeling domain and
ultimately assist researchers in improving and interpreting their (L)LMs.

##  2 Related Works

$\mathcal{NC}$ was initially observed in image classification tasks such as on
CIFAR-10/100 [41] and ImageNet [42]. Since then, the phenomenon has been
further studied both theoretically and empirically [15, 32, 37, 43, 44, 45,
46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
65], with several works venturing into settings without some of the
traditional prerequisites (C1-4, §1.1):

#### A large number of classes

$\mathcal{NC}$ established the simplex ETF as an optimal configuration.
However, a perfect simplex ETF ($\mathcal{NC}_{2}$) requires that the number
of classes $C$ not exceed $d+1$ where $d$ is the embedding dimension. This
requirement that $d$ be sufficiently large is not practical when the classes
number beyond the thousands. GPT-2 [66], for instance, has a vocabulary of
over 50K tokens. 222Following [67], one might describe such a setting
($C>d+1$) as a model “in superposition”.

In such a scenario, one might still expect class means to be uniformly
distributed on a $d$-dimensional hypersphere [15]. [16] formalize this as
hyperspherical uniformity within a generalized neural collapse
($\mathcal{GNC}$) framework, which [6] then empirically demonstrate. These
latter two works mention language modeling as applicable for $\mathcal{GNC}$;
[6] even frame token prediction as a classification task, just as we do. We
remark however that both earlier works simulate a large number of classes by
drastically shrinking the embedding dimension. In contrast, we study next-
token prediction, using the full class space (vocabulary) with an imbalanced
token distribution and ambiguous samples.

#### Class imbalance

$\mathcal{NC}$ traditionally assumed that classes were sample-balanced. Since
then, follow-up works have investigated the effect of class imbalance on the
prevalence of $\mathcal{NC}$. [46] studied the layer-peeled model (LPM) and
discovered that minority collapse occurs when classes are imbalanced across
two equally-sized groups of classes; a threshold for minority collapse was
later characterized by [32]. [51] showed that $\mathcal{NC}$ still occurs in
such an LPM when the classifier is initialized as a fixed ETF. [68] introduced
simplex-encoded-label interpolation (SELI) but noted that more severe
imbalance worsens even this geometry. Recently, feature regularization has
been employed to induce $\mathcal{NC}$ and improve generalization in class-
imbalanced settings [57, 60, 61].

#### Multi-label classification

In some problems, one might encounter mixed or multi-label samples, be they
natural or due to noise or augmentation [69, 70]. $\mathcal{NC}$ was also
recently studied for such data by [63], who observed that multi-label class
means arrive at the average of their individual labels’ means. They also
devise an augmented CE loss function to accommodate such samples.

However, most of our ambiguous samples are not soft- nor multi-label, but
rather identical (or very similar) context samples with different hard token
labels ($\neg$C3, §1.2). This effectively rules out the prospect of achieving
zero classification error and potentially introduces irreducible noise.

A recent study showed that memorization of noisy labels likely leads to
degradation of $\mathcal{NC}$ [59].

#### Early stages of training

[52] studied $\mathcal{NC}$ in small ResNet [71] models that had not yet
converged (similar to most LLMs). They show that within-class variability
drops and “plateaus” ($\mathcal{NC}$1) earlier than other $\mathcal{NC}$
metrics, a result that we also observe (§4.1, Figures 5, 14).

#### Natural language processing (NLP)

An earlier study reported that the ratio of within-class to between-class
covariances of word embeddings increases with model depth [72, 73], seemingly
at odds with $\mathcal{NC}_{1}$. It can however be distinguished from
literature studying layer-wise $\mathcal{NC}$ in that it does not center the
mean representation vectors (i.e. subtract the global mean).

[16] fine-tuned BERT [4] using their hyperspherical uniformity gap (HUG)
objective on binary classification tasks from the GLUE benchmark [74]. [75]
conducted a tangentially-related investigation of convolutional neural
networks for few-class semi-supervised clustering in which they identify
$\mathcal{NC}$. Our work is distinct from these in several ways: a) our class
space far exceeds our embedding dimension ($C\gg d+1$) because we classify on
the full token vocabulary; b) we analyze embeddings at a token-level
granularity rather than at the sequence level; and c) our next token
prediction is causal (context-dependent) as opposed to the per-sample
independence of their few-category classification.

A more related work studied feature collapse in individual word
representations [76], but the authors note that their analysis is limited to
shallow NLP modeling on more rigid and symmetrical text data.

##  3 Methodology

In this section, we describe the training setup for our CLMs (§3.1),
procedures for collecting top-layer embeddings (§3.2, 3.7), and measurements
of $\mathcal{NC}$ and generalization (§3.4, 3.5, 3.6, 3.7, 3.8).

###  3.1 Dataset and Language Models

TinyStories is a synthetic333 TinyStories was developed and evaluated as a
faithful emulation of large language modeling [2], so we chose it for
experimentation to train hundreds of CLMs and analyze embeddings at low cost.
See Appendix B. dataset generated by GPT-3.5 and GPT-4 [2] with token
vocabulary $\mathbb{V}=\llbracket 1,29233\rrbracket$. Token prediction can
therefore be framed as classification among $C=$29,233$$ token
classes.444Although the GPT-2 [77] tokenizer has over 50K tokens, only a
subset vocabulary appears in TinyStories. Following the GPT-style regime[5],
raw sequences are packed into $S$ chunks of size $T$, providing $N=S(T-1)$
token samples for training.555We cannot use the first ground-truth nor last
predicted token in any chunk. Details are listed in Appendix B.

We use 30 CLM architectures based on GPT Neo [78], configured similarly to
[2]. They vary in width (embedding dimension)
$d\in\\{64,128,256,512,768,1024\\}$ and depth (number of self-attention
layers) $L\in\\{1,2,4,8,12\\}$. Our models were trained by teacher-
forcing666Parallelized training using sequences’ ground-truth labels for
context as opposed to predicted tokens. using CE loss. For each architecture,
we trained multiple models for 1, 3, and 10 epochs. We albated over weight
decay factors $\beta=0.0005$ [50] and $\beta=0.1$ [77]. Further details are
listed in Appendices A, C.

###  3.2 Context Embeddings

Base CLMs perform next-token prediction: given a sequence of tokens
$\boldsymbol{x}_{1:t}\in\mathbb{V}^{t}$, a top-layer context embedding
$\boldsymbol{h}(\boldsymbol{x}_{1:t})\in\mathbb{R}^{d}$ is used to predict the
next token $x_{t+1}^{\prime}\in\mathbb{V}$ where $1\leq t\leq T$. A classifier
for class $c$ with weights $\boldsymbol{w}_{c}$ and bias777Similar to many
CLMs [78, 79, 80, 81], our classifiers do not include additive biases, so
$b_{c}=0$. $b_{c}$ would make maximum a prior (MAP) predictions as

$x_{t+1}^{\prime}:=\operatorname*{argmax}_{c\in\mathbb{V}}\left\langle\boldsymbol{w}_{c},\boldsymbol{h}(\boldsymbol{x}_{1:t})\right\rangle+b_{c}.$
(1)

#### Class embedding means

For each token class $c$, we are interested in the mean embedding
$\boldsymbol{\mu}_{c}\in\mathbb{R}^{d}$ across sequences $s$ and their
contexts $\boldsymbol{x}_{1:t}^{(s)}$, where the next token $x_{t+1}^{(s)}$ is
ground-truth ($t<T$) and equal to $c$. To center these means, we compute their
unweighted 888Different from most literature where classes were balanced and
weighting is already equal. average $\boldsymbol{\bar{\mu}}\in\mathbb{R}^{d}$.
Precisely, we take

$\boldsymbol{\mu}_{c}:=\frac{1}{N_{c}}\sum_{s=1}^{S}\sum_{t=1}^{T-1}\boldsymbol{h}\left(\boldsymbol{x}_{1:t}^{(s)}\right)\mathbb{I}\left(x_{t+1}^{(s)}=c\right),\qquad\boldsymbol{\bar{\mu}}:=\frac{1}{C}\sum_{c=1}^{C}\boldsymbol{\mu}_{c},$
(2)

where $N_{c}$ is the number of samples of class $c$ and $\mathbb{I}$ is the
(binary) indicator function.

#### Class embedding variances

In a second pass, we accumulate the unbiased sample variances999This sample
variance is computed across all unnormalized dimensional entries.:

$\sigma_{c}^{2}:=\frac{1}{N_{c}-1}\sum_{s=1}^{S}\sum_{t=1}^{T-1}\left\|\boldsymbol{h}\left(\boldsymbol{x}_{1:t}^{(s)}\right)-\boldsymbol{\mu}_{c}\right\|_{2}^{2}\mathbb{I}\left(x_{t+1}^{(s)}=c\right).$
(3)

###  3.3 Homogeneity and Variability

For some collapse measures (such as $\mathcal{(G)NC}_{2}$ and
$\mathcal{NC}_{3}$), we are primarily interested in the variation rather than
the average of pair-wise relationships. To that end, we also include in our
analysis the coefficient of variation (CoV) of several measurements, which is
the ratio of their standard deviations to their means:
$\sigma(\cdot)/\mu(\cdot)$. We can interpret this as a normalized measure of
variability.

###  3.4 Signal-to-Noise Ratio – $\mathcal{NC}_{1}$

The ability to disambiguate between classes depends on the ratio of within-
class to between-class variabilities. Building upon foundational works [82,
83], $\mathcal{NC}$ originally measured variability through an inverse signal-
to-noise ratio (SNR), whose minimization constitutes within-class variability
collapse ($\mathcal{NC}_{1}$). We employ the class-distance normalized
variance (CDNV) proposed by [34]:

$\text{CDNV}_{c,c^{\prime}}:=\frac{\sigma_{c}^{2}+\sigma_{c^{\prime}}^{2}}{2\|\boldsymbol{\mu}_{c}-\boldsymbol{\mu}_{c^{\prime}}\|_{2}^{2}},\quad\forall
c\neq c^{\prime}.$ (4)

These pair-wise measures constitute the off-diagonal101010The main diagonal of
CDNVs is undefined (due to the zero denominator) and ignored. entries of a
symmetric matrix in $\mathbb{R}^{C\times C}$, whose average we use as an
inverse SNR. Within-class variability collapse is then re-characterized by the
minimization of this quantity: $\text{CDNV}_{c,c^{\prime}}\rightarrow
0,\forall c\neq c^{\prime}$. This alternative convergence is faithful to
$\mathcal{NC}_{1}$ but more robust and numerically stable [34].

###  3.5 Geometric Structures – $\mathcal{(G)NC}_{2}$

The separability of our representations also depends on the geometric
structures found in our embeddings. [1] characterize $\mathcal{NC}_{2}$ as
convergence to a simplex equiangular tight frame (ETF) [84, 85].

#### Equinormness

Such a near-orthonormal configuration firstly implies that class means are
equinorm:

$\log\|\boldsymbol{\mu}_{c}-\boldsymbol{\bar{\mu}}\|_{2}-\log\|\boldsymbol{\mu}_{c^{\prime}}-\boldsymbol{\bar{\mu}}\|_{2}\rightarrow
0,\quad\forall c\neq c^{\prime}.$ (5)

We measure CoV in the logarithms of class mean norms to assess “equinormness”.

#### Equiangularity

$\mathcal{NC}_{2}$ also entails that class means are equiangular about their
center $\boldsymbol{\bar{\mu}}$: pair-wise distances and angles between their
class means should be maximized and similar. Following [1], we measure
interference (sometimes known as similarity or coherence [86, 87]). Its
minimization,

$\left\langle\frac{\boldsymbol{\mu}_{c}-\boldsymbol{\bar{\mu}}}{\|\boldsymbol{\mu}_{c}-\boldsymbol{\bar{\mu}}\|_{2}},\frac{\boldsymbol{\mu}_{c^{\prime}}-\boldsymbol{\bar{\mu}}}{\|\boldsymbol{\mu}_{c^{\prime}}-\boldsymbol{\bar{\mu}}\|_{2}}\right\rangle\rightarrow\frac{-1}{C-1},\quad\forall
c\neq c^{\prime},$ (6)

together with equinormness (Equation 5) constitute convergence to a simplex
ETF. Although this geometry is not ultimately attainable since there are too
many classes ($C>d+1$), it can still be meaningful to measure a model’s
tendency towards one. As with CDNV noise (Equation 4), pair-wise interferences
form off-diagonal111111The main diagonal of interferences is equal to
$\boldsymbol{1}$ (maximal coherence or self-similarity). entries in a
symmetric matrix in $\mathbb{R}^{C\times C}$. The minimization of CoV in
interferences therefore expresses the degree of “equiangularity”.

#### Hypershperical uniformity

A relaxation from the ETF is hyperspherical uniformity ($\mathcal{GNC}_{2}$),
with mean vectors $\boldsymbol{\mu}_{c}$ uniformly distributed on the
$d$-dimensional hypersphere [15, 16]. We gauge the degree of this uniformity
with pair-wise interactions under a logarithmic inverse distance kernel:
121212Following [16], we employ the logarithmic inverse distance kernel
$K_{\log}(\boldsymbol{a},\boldsymbol{b})=\log\|\boldsymbol{a}-\boldsymbol{b}\|^{-1}$
for its ability to emphasize gaps between small distances while scaling down
the effect of larger distances.

$\log\left\|\frac{\boldsymbol{\mu}_{c}-\boldsymbol{\bar{\mu}}}{\|\boldsymbol{\mu}_{c}-\boldsymbol{\bar{\mu}}\|_{2}}-\frac{\boldsymbol{\mu}_{c^{\prime}}-\boldsymbol{\bar{\mu}}}{\|\boldsymbol{\mu}_{c^{\prime}}-\boldsymbol{\bar{\mu}}\|_{2}}\right\|^{-1},\;\forall
c\neq c^{\prime}.$ (7)

###  3.6 Alignment Between Classifiers and Class Mean Embeddings –
$\mathcal{(U)NC}_{3}$

The linear classifiers $\\{\boldsymbol{w}_{c}\\}_{c=1}^{C}$ lie in a dual
vector space to that of the class means
$\\{\boldsymbol{\mu}_{c}\\}_{c=1}^{C}$. While convergence to self-duality
($\mathcal{NC}_{3}$) was previously measured as distances [1] between class
means classifiers (Equation 11), we instead simplify ($\mathcal{NC}_{3}$) as
class-wise similarities:

$\left\langle\frac{\boldsymbol{w}_{c}}{\|\boldsymbol{w}_{c}\|_{2}},\frac{\boldsymbol{\mu}_{c^{\prime}}-\boldsymbol{\bar{\mu}}}{\|\boldsymbol{\mu}_{c^{\prime}}-\boldsymbol{\bar{\mu}}\|_{2}}\right\rangle\rightarrow
1,\quad\forall c\in\mathbb{V}.$ (8)

We define uniform duality ($\mathcal{UNC}_{3}$) to be the minimization of the
CoV of these similarities.

###  3.7 Agreement of the Classifiers – $\mathcal{NC}_{4}$

Finally, $\mathcal{NC}_{4}$ is described as the simplification of the linear
classifier’s MAP predictions (Equation 1) for test131313We used the validation
set since it was not used for training or hyperparameter optimization. points
to those of the implicit nearest-class center (NCC) classifier:

$\operatorname*{argmax}_{c\in\mathbb{V}}\left\langle\boldsymbol{w}_{c},\boldsymbol{h}\right\rangle+b_{c}\rightarrow\operatorname*{argmin}_{c\in\mathbb{V}}\|\boldsymbol{h}-\boldsymbol{\mu}_{c}\|_{2}.$
(9)

We calculate141414In practice, we use an equivalent decomposition (Eq. 12).
agreement as the proportion of samples on which the classifiers agree:
151515Note that agreement is not equivalent to accuracy.

$\frac{1}{N_{\text{val}}}\sum_{s=1}^{S_{\text{val}}}\sum_{t=1}^{T_{\text{val}}-1}\mathbb{I}\left({x_{t+1}^{\prime}}^{(s)}=\operatorname*{argmin}_{c\in\mathcal{V}}\left\|\boldsymbol{h}\left(x_{1:t}^{(s)}\right)-\boldsymbol{\mu_{c}}\right\|_{2}\right).$
(10)

###  3.8 Probing a Relationship Between $\mathcal{NC}$ and Generalization

To isolate the effect of $\mathcal{NC}$ on generalization independent of model
scaling and training (if it exists), we selected a two-layer $768$-wide
architecture of which to train twenty more instances with weight decay
$\beta=0.0005$, each with a different data shuffling seed. We then followed
the remainder of the pipeline described above to collect and analyze
embeddings with respect to $\mathcal{NC}$. Finally, we performed a permutation
test with $10^{4}$ trials to determine the statistical significance of any
correlation between $\mathcal{NC}$ and generalization that remains when we
hold constant all factors but shuffling seeds.

##  4 Results

In this section, we present the results from our empirical study on scaling
and generalization:

  1. ($\mathcal{NC}_{1}$)

Within-class variability is reduced across model scale (more so by width than
depth) and training (up to 6 epochs), and is tightly correlated with
validation performance.

  2. ($\mathcal{NC}_{2}$)

Equinormness/equiangularity shows some improvement with scale, training, and
performance. Hyperspherical uniformity ($\mathcal{GNC}_{2}$) also improves but
more clearly and consistently.

  3. ($\mathcal{NC}_{3}$)

Our models fail to achieve self-duality: class means do not align with
classifiers. However, uniform duality ($\mathcal{UNC}3$) is correlated with
model width, training, and performance.

  4. ($\mathcal{NC}_{4}$)

Larger or more trained models exhibit closer agreement between their linear
and implicit NCC classifiers. Agreement is also associated with validation
performance.

We find that $\mathcal{NC}$ is generally promoted by model size and training
and correlated with generalization (validation performance). We also discern
some of this correlation independent of scale (§4.5).

###  4.1 Within-Class Variability Collapse – $\mathcal{NC}_{1}$

Scaling our models dramatically reduces normalized variance, which is further
aided by more training epochs and stronger weight decay (Appendix Figs. 5, 6).
These noise reductions tightly associate with generalization (Fig. 1, left,
“$\mathcal{NC}_{1}$”). The relationship is strongest at larger model scale.

###  4.2 Geometric Structures – $\mathcal{(G)NC}_{2}$

Figure 2: Validation loss shows a negligible correlation with equiangularity
($\mathcal{NC}_{2}$, left) but a strong one with hyperspherical uniformity
($\mathcal{GNC}_{2}$, right). In other words, $\mathcal{GNC}_{2}$ develops
with scale and correlates well with generalization, much better than
$\mathcal{NC}_{2}$.

#### Equinormness

Mean norms grow with model width and training (Appendix Fig. 7), and subtly
with depth (Appendix Fig. 8). Despite this growth, the variation of norms
consistently decreases with scale (Appendix Figs. 9, 10). Both trends
correlate with improved generalization (Appendix Fig. 25).

#### Equiangularity

Scaling model dimensions reduces average interference (Appendix Figs. 11, 12)
down to an empirical plateau of approximately $10^{-3}$, which is more
apparent in less-trained models. However, the variation of interference rises
and persists when scaling (Appendix Figs. 13, 14), suggesting that
interference becomes more concentrated between some pairs of classes. These
results could be due to various factors, including but not limited to
unfriendly conditions of language modeling (§1.2) or the impossibility of a
simplex ETF (§3.5).

We only see a rough correlation between average interference and
generalization and almost none with its variation, equiangularity (Appendix
Fig. 26). In other words, the limited trends we observed toward a simplex ETF
do not appear to be associated with generalization when $C>d+1$.

#### Hyperspherical uniformity

Logarithmic distances drop more gradually and consistently with scale
(Appendix Figs. 15, 16), implying this quantity is more robust or may not face
the same geometric barriers seen in conventional interference (Appendix Figs.
13, 14). Variation of these logarithmic distances is also consistently reduced
with scale (Appendix Fig. 17, 18).

And finally, generalization has much stronger correlations with logarithmic
distances than it has with regular interference (Fig. 2), validating the
applicability of $\mathcal{GNC}$ [16] when $C>d+1$.

###  4.3 Classifier (Mis)alignment and Duality – $\mathcal{(U)NC}_{3}$

Model width ($d$) is correlated faintly with average similarity between class
means and their respective classifiers (Appendix Fig. 19), but strongly with
variation in similarity (Appendix Fig. 21). The relationships to
generalization follow the same pattern (Fig. 3), suggesting that uniform
duality ($\mathcal{UNC}_{3}$) might serve as a better $\mathcal{NC}$ property
than self-duality ($\mathcal{NC}_{3}$) overall; we discuss this in §4.5.

Figure 3: Validation loss shows a negligible relationship with self-duality
($\mathcal{NC}_{3}$, left) and some correlation with uniform duality
($\mathcal{UNC}_{3}$, right). In other words, $\mathcal{UNC}_{3}$ develops
with scale and correlates with generalization well and much better than
$\mathcal{NC}_{3}$.

###  4.4 Classifier Agreement – $\mathcal{NC}_{4}$

The linear and NCC classifiers agree on far more examples than random chance.
Scaling encourages agreement (Appendix Figs. 23, 24). Increasing width for
certain depths happens to plateau or even regress the agreement rate, but this
limitation is overcome with further training. And finally, agreement is a
strong indicator of generalization (Fig. 1, right, $\mathcal{NC}_{4}$).

###  4.5 Neural Collapse’s Relationship with Generalization

Table 1 presents the correlation scores of $\mathcal{NC}$ metrics with
generalization alongside their associated p-values from the permutation tests
described in §3.8. The majority of the correlations are statistically
significant ($p<0.05$) independent of scale, affirming that $\mathcal{NC}$ is
correlated with generalization.

###  4.6 Duality of Duality

The dichotomy of self-duality ($\mathcal{NC}_{3}$) and uniform duality
($\mathcal{UNC}_{3}$) is rather conspicuous. Our main experiments find that
$\mathcal{NC}_{3}$ does not consistently develop with scale while
$\mathcal{UNC}_{3}$ does (Fig. 3). However, within fixed scale, the opposite
is true, implying that $\mathcal{UNC}_{3}$ may be confounded by model capacity
while $\mathcal{NC}_{3}$ is a subtle and fine-grained indicator of
generalization.

###  4.7 The Effect of Weight Regularization

Our models trained with either weight decay factor exhibited very similar
patterns in the emergence of $\mathcal{NC}$ (or lack thereof), but the more
aggressive factor $\beta=0.1$ typically used to train LLMs [77] resulted in
stronger development of $\mathcal{NC}$ properties (Appendices D, E, F, G, H,
I, J, K, L, M).

Table 1: Permutation test of $\mathcal{NC}$ measurements with respect to validation loss. Twenty-one (21) identical two-layer $768$-wide models were trained with different data shuffling seeds and permuted with $10^{4}$ trials. The $p$-values below $0.05$ (bolded) show those properties to be statistically significant. Property | Measurement |  $R^{2}$ Correlation ($\uparrow$) |  $p$-value ($\downarrow$)
---|---|---|---
$\mathcal{NC}$1 | Within-Class Variability Collapse | $0.192011$ | $\mathbf{0.0485}$
$\mathcal{NC}$2 | Equinormness | $0.026174$ | $0.4870$
$\mathcal{NC}$2 | Equiangularity | $0.218574$ | $\mathbf{0.0317}$
$\mathcal{GNC}$2 | Hyperspherical Uniformity | $0.487935$ | $\mathbf{0.0002}$
$\mathcal{NC}$3 | Self-Duality | $0.322210$ | $\mathbf{0.0063}$
$\mathcal{UNC}$3 | Uniform Duality | $0.000036$ | $0.9784$
$\mathcal{NC}$4 | Classifier Agreement | $0.490278$ | $\mathbf{0.0001}$

##  5 Limitations

While to the best of our knowledge, no previous work has studied realistic
stochastic token prediction, it is possible that the quantities that we
measure are not perfectly suited for $\mathcal{NC}$ in language modeling.

Our work is on causal language modeling in its most basic form. We did not
conduct experiments into encoding, multi-modal, or instruction-tuned models.

The models that we used in our permutation test (§4.5, Table 1) are only of
single small architecture. Therefore, our results on scale-independent
$\mathcal{NC}$ do not necessarily translate to larger models.

##  6 Discussion

#### Deep neural collapse ($\mathcal{DNC}$)

A large body of work has established that properties resembling $\mathcal{NC}$
evolve as a function of model depth [33, 52, 58, 59, 88, 89, 90, 91, 92, 93,
94, 95, 96, 97, 98, 99, 100, 101]. Our results mostly concur, but encourage
further exploration of $\mathcal{DNC}$ in language modeling.

#### Learning to collapse

Given the evidence for the development of $\mathcal{NC}$ and associated
generalization under various loss functions [16, 48, 59, 97, 102, 103, 104] in
other domains, NLP researchers may also benefit from analyzing, adapting or
training towards $\mathcal{NC}$ with respect to alternative objectives.

#### Compression

CLMs appear to reduce noise, demonstrating their ability to encapsulate
different textual contexts into single embeddings. This subscribes to a theory
that LLMs are essentially a compression of internet content by aleatoric
modeling [105]. However, owing to contextual ambiguity, we hypothesize that CE
loss and noise (Equation 4) cannot be fully eliminated by current practices.

#### Interpretability

Our formulae (Section 3) and results (Section 4) expose the pair-wise token
class interactions in noise, interference, classifier duality, and classifier
agreement in the top-level features. Similarly to works in other domains [18,
54, 70, 106], these $\mathcal{NC}$ metrics can serve as a form of low-level
interpretability to aid understanding certain behaviors of (L)LMs.

#### Fairness

Foundation models are ubiquitous for their comprehensive capabilities and
adaptability. As previous work discussed class imbalance [32, 57, 60, 61], our
work may extend these strategies to measure and perhaps promote fairness in
foundation LLMs, some of which are designed to be multilingual or
multicultural. While $\mathcal{NC}$ itself would not lead to unfairness, its
potential interpretability may, in theory, enable an adversarial agent to
measure (un)fairness as they manipulate a LLM.

##  7 Conclusion

In this paper, we adapt the neural collapse ($\mathcal{NC}$) framework to the
language modeling domain, where models are undertrained and next-tokens are
variably drawn from numerous and imbalanced token classes. We leverage the
canonical and more recent metrics to demonstrate that $\mathcal{NC}$ emerges
as we scale the size and training of hundreds of causal language models. Our
results show a correlation between $\mathcal{NC}$ and generalization, a
relationship that persists even when the model scale is fixed.

With these novel insights, we hope to inspire new efforts to study
$\mathcal{NC}$ (and related behaviors) in language modeling and capitalize on
their previously established benefits to improve and better understand the
increasingly complex (large) language models and their training processes.

## References

  * [1] Vardan Papyan, X. Y. Han, and David L. Donoho.  Prevalence of neural collapse during the terminal phase of deep learning training.  Proceedings of the National Academy of Sciences, 117(40):24652–24663, sep 2020.
  * [2] Ronen Eldan and Yuanzhi Li.  Tinystories: How small can language models be and still speak coherent english?, 2023.
  * [3] Emily M. Bender, Timnit Gebru, Angelina McMillan-Major, and Shmargaret Shmitchell.  On the dangers of stochastic parrots: Can language models be too big?  In Proceedings of the 2021 ACM Conference on Fairness, Accountability, and Transparency, FAccT ’21, page 610–623, New York, NY, USA, 2021. Association for Computing Machinery.
  * [4] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding, 2019.
  * [5] Alec Radford, Karthik Narasimhan, Tim Salimans, Ilya Sutskever, et al.  Improving language understanding by generative pre-training.  2018\.
  * [6] Jiachen Jiang, Jinxin Zhou, Peng Wang, Qing Qu, Dustin Mixon, Chong You, and Zhihui Zhu.  Generalized neural collapse for a large number of classes, 2023.
  * [7] Zhilin Yang, Zihang Dai, Ruslan Salakhutdinov, and William W. Cohen.  Breaking the softmax bottleneck: A high-rank rnn language model, 2018\.
  * [8] Claude Elwood Shannon.  A mathematical theory of communication.  The Bell system technical journal, 27(3):379–423, 1948.
  * [9] P Sargant Florence.  Human behaviour and the principle of least effort.  The Economic Journal, 60(240):808–810, 1950.
  * [10] Daniel Jurafsky and James H Martin.  Speech and language processing: An introduction to natural language processing, computational linguistics, and speech recognition, 2009.
  * [11] Jared Kaplan, Sam McCandlish, Tom Henighan, Tom B. Brown, Benjamin Chess, Rewon Child, Scott Gray, Alec Radford, Jeffrey Wu, and Dario Amodei.  Scaling laws for neural language models, 2020.
  * [12] Jordan Hoffmann, Sebastian Borgeaud, Arthur Mensch, Elena Buchatskaya, Trevor Cai, Eliza Rutherford, Diego de Las Casas, Lisa Anne Hendricks, Johannes Welbl, Aidan Clark, Tom Hennigan, Eric Noland, Katie Millican, George van den Driessche, Bogdan Damoc, Aurelia Guy, Simon Osindero, Karen Simonyan, Erich Elsen, Jack W. Rae, Oriol Vinyals, and Laurent Sifre.  Training compute-optimal large language models, 2022.
  * [13] Niklas Muennighoff, Alexander M. Rush, Boaz Barak, Teven Le Scao, Aleksandra Piktus, Nouamane Tazi, Sampo Pyysalo, Thomas Wolf, and Colin Raffel.  Scaling data-constrained language models, 2023.
  * [14] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  Advances in neural information processing systems, 30, 2017.
  * [15] Jianfeng Lu and Stefan Steinerberger.  Neural collapse with cross-entropy loss, 2021.
  * [16] Weiyang Liu, Longhui Yu, Adrian Weller, and Bernhard Schölkopf.  Generalizing and decoupling neural collapse via hyperspherical uniformity gap, 2023.
  * [17] Zexi Li, Xinyi Shang, Rui He, Tao Lin, and Chao Wu.  No fear of classifier biases: Neural collapse inspired federated learning with synthetic and fixed classifier.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5319–5329, 2023.
  * [18] Vignesh Kothapalli, Tom Tirer, and Joan Bruna.  A neural collapse perspective on feature evolution in graph neural networks.  Advances in Neural Information Processing Systems, 36, 2024.
  * [19] Yibo Yang, Haobo Yuan, Xiangtai Li, Jianlong Wu, Lefei Zhang, Zhouchen Lin, Philip Torr, Dacheng Tao, and Bernard Ghanem.  Neural collapse terminus: A unified solution for class incremental learning and its variants, 2023.
  * [20] Qinhao Zhou, Xiang Xiang, and Jing Ma.  Hierarchical task-incremental learning with feature-space initialization inspired by neural collapse.  Neural Processing Letters, 55(8):10811–10827, 2023.
  * [21] Hang Ran, Weijun Li, Lusi Li, Songsong Tian, Xin Ning, and Prayag Tiwari.  Learning optimal inter-class margin adaptively for few-shot class-incremental learning via neural collapse-based meta-learning.  Information Processing & Management, 61(3):103664, 2024.
  * [22] Saaketh Medepalli and Naren Doraiswamy.  On the role of neural collapse in meta learning models for few-shot learning, 2023.
  * [23] Jarrod Haas, William Yolland, and Bernhard Rabus.  Linking neural collapse and l2 normalization with improved out-of-distribution detection in deep neural networks, 2023.
  * [24] Mouïn Ben Ammar, Nacim Belkhir, Sebastian Popescu, Antoine Manzanera, and Gianni Franchi.  Neco: Neural collapse based out-of-distribution detection, 2024.
  * [25] Jiawei Zhang, Yufan Chen, Cheng Jin, Lei Zhu, and Yuantao Gu.  Epa: Neural collapse inspired robust out-of-distribution detector, 2024\.
  * [26] Donghao Li, Yang Cao, and Yuan Yao.  Neuromixgdp: A neural collapse-inspired random mixup for private data release.  In Conference on Parsimony and Learning, pages 480–514. PMLR, 2024\.
  * [27] Chendi Wang, Yuqing Zhu, Weijie J Su, and Yu-Xiang Wang.  Neural collapse meets differential privacy: Curious behaviors of noisySGD with near-perfect representation learning, 2024.
  * [28] Tolga Ergen, Arda Sahiner, Batu Ozturkler, John Pauly, Morteza Mardani, and Mert Pilanci.  Demystifying batch normalization in relu networks: Equivalent convex optimization models and implicit regularization, 2022.
  * [29] Tolga Ergen and Mert Pilanci.  Revealing the structure of deep neural networks via convex duality.  In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 3004–3014. PMLR, 18–24 Jul 2021.
  * [30] Mariia Seleznova, Dana Weitzner, Raja Giryes, Gitta Kutyniok, and Hung-Hsu Chou.  Neural (tangent kernel) collapse.  In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 16240–16270. Curran Associates, Inc., 2023.
  * [31] Matus Telgarsky.  Feature selection with gradient descent on two-layer networks in low-rotation regimes, 2022.
  * [32] Wanli Hong and Shuyang Ling.  Neural collapse for unconstrained feature model under cross-entropy loss with imbalanced data, 2023.
  * [33] Gerard Ben Arous, Reza Gheissari, Jiaoyang Huang, and Aukosh Jagannath.  High-dimensional sgd aligns with emerging outlier eigenspaces, 2023.
  * [34] Tomer Galanti, András György, and Marcus Hutter.  On the role of neural collapse in transfer learning, 2022.
  * [35] Like Hui, Mikhail Belkin, and Preetum Nakkiran.  Limitations of neural collapse for understanding generalization in deep learning, 2022.
  * [36] Tomer Galanti, András György, and Marcus Hutter.  Improved generalization bounds for transfer learning via neural collapse.  In First Workshop on Pre-training: Perspectives, Pitfalls, and Paths Forward at ICML 2022, 2022.
  * [37] Vignesh Kothapalli.  Neural collapse: A review on modelling principles and generalization.  Transactions on Machine Learning Research, 2023.
  * [38] Zijian Wang, Yadan Luo, Liang Zheng, Zi Huang, and Mahsa Baktashmotlagh.  How far pre-trained models are from neural collapse on the target dataset informs their transferability.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5549–5558, 2023.
  * [39] Anonymous.  Understanding and improving transfer learning of deep models via neural collapse.  Submitted to Transactions on Machine Learning Research, 2024.  Under review.
  * [40] Didi Zhu, Zexi Li, Min Zhang, Junkun Yuan, Yunfeng Shao, Jiashuo Liu, Kun Kuang, Yinchuan Li, and Chao Wu.  Understanding prompt tuning for v-l models through the lens of neural collapse, 2023.
  * [41] Alex Krizhevsky and Geoffrey Hinton.  Learning multiple layers of features from tiny images.  Technical Report 0, University of Toronto, Toronto, Ontario, 2009.
  * [42] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee, 2009.
  * [43] Dustin G. Mixon, Hans Parshall, and Jianzong Pi.  Neural collapse with unconstrained features, 2020.
  * [44] Tomaso Poggio and Qianli Liao.  Explicit regularization and implicit bias in deep network classifiers trained with the square loss, 2020.
  * [45] Weinan E and Stephan Wojtowytsch.  On the emergence of simplex symmetry in the final and penultimate layers of neural network classifiers, 2021.
  * [46] Cong Fang, Hangfeng He, Qi Long, and Weijie J Su.  Exploring deep neural networks via layer-peeled model: Minority collapse in imbalanced training.  Proceedings of the National Academy of Sciences, 118(43):e2103091118, 2021.
  * [47] Zhihui Zhu, Tianyu Ding, Jinxin Zhou, Xiao Li, Chong You, Jeremias Sulam, and Qing Qu.  A geometric analysis of neural collapse with unconstrained features, 2021\.
  * [48] X.Y. Han, Vardan Papyan, and David L. Donoho.  Neural collapse under MSE loss: Proximity to and dynamics on the central path.  In International Conference on Learning Representations, 2022.
  * [49] Can Yaras, Peng Wang, Zhihui Zhu, Laura Balzano, and Qing Qu.  Neural collapse with normalized features: A geometric analysis over the riemannian manifold.  Advances in neural information processing systems, 35:11547–11560, 2022.
  * [50] Akshay Rangamani and Andrzej Banburski-Fahey.  Neural collapse in deep homogeneous classifiers and the role of weight decay.  In ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4243–4247, 2022.
  * [51] Yibo Yang, Shixiang Chen, Xiangtai Li, Liang Xie, Zhouchen Lin, and Dacheng Tao.  Inducing neural collapse in imbalanced learning: Do we really need a learnable classifier at the end of deep neural network?, 2022.
  * [52] Tom Tirer and Joan Bruna.  Extended unconstrained features model for exploring deep neural collapse.  In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 21478–21505. PMLR, 17–23 Jul 2022.
  * [53] Peng Wang, Huikang Liu, Can Yaras, Laura Balzano, and Qing Qu.  Linear convergence analysis of neural collapse with unconstrained features.  In OPT 2022: Optimization for Machine Learning (NeurIPS 2022 Workshop), 2022.
  * [54] Jinxin Zhou, Chong You, Xiao Li, Kangning Liu, Sheng Liu, Qing Qu, and Zhihui Zhu.  Are all losses created equal: A neural collapse perspective, 2022.
  * [55] Akshay Rangamani, Marius Lindegaard, Tomer Galanti, and Tomaso A Poggio.  Feature learning in deep classifiers through intermediate neural collapse.  In International Conference on Machine Learning, pages 28729–28745. PMLR, 2023.
  * [56] Tom Tirer, Haoxiang Huang, and Jonathan Niles-Weed.  Perturbation analysis of neural collapse.  In International Conference on Machine Learning, pages 34301–34329. PMLR, 2023.
  * [57] Hien Dang, Tho Tran, Stanley Osher, Hung Tran-The, Nhat Ho, and Tan Nguyen.  Neural collapse in deep linear networks: From balanced to imbalanced data, 2023.
  * [58] Peter Súkeník, Marco Mondelli, and Christoph H Lampert.  Deep neural collapse is provably optimal for the deep unconstrained features model.  Advances in Neural Information Processing Systems, 36, 2024.
  * [59] Duc Anh Nguyen, Ron Levie, Julian Lienen, Eyke Hüllermeier, and Gitta Kutyniok.  Memorization-dilation: Modeling neural collapse under noise.  In The Eleventh International Conference on Learning Representations, 2023.
  * [60] Zhisheng Zhong, Jiequan Cui, Yibo Yang, Xiaoyang Wu, Xiaojuan Qi, Xiangyu Zhang, and Jiaya Jia.  Understanding imbalanced semantic segmentation through neural collapse.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 19550–19560, 2023.
  * [61] Xuantong Liu, Jianfeng Zhang, Tianyang Hu, He Cao, Yuan Yao, and Lujia Pan.  Inducing neural collapse in deep long-tailed learning.  In International Conference on Artificial Intelligence and Statistics, pages 11534–11544. PMLR, 2023.
  * [62] Peifeng Gao, Qianqian Xu, Peisong Wen, Huiyang Shao, Zhiyong Yang, and Qingming Huang.  A study of neural collapse phenomenon: Grassmannian frame, symmetry and generalization, 2023.
  * [63] Mufan Bill Li, Mihai Nica, and Daniel M. Roy.  The neural covariance sde: Shaped infinite depth-and-width networks at initialization, 2023.
  * [64] Zhanxuan Hu, Yichen Wang, Hailong Ning, Yonghang Tai, and Feiping Nie.  Neural collapse inspired semi-supervised learning with fixed classifier.  Information Sciences, 667:120469, 2024.
  * [65] Gao Peifeng, Qianqian Xu, Yibo Yang, Peisong Wen, Huiyang Shao, Zhiyong Yang, Bernard Ghanem, and Qingming Huang.  Towards demystifying the generalization behaviors when neural collapse emerges, 2024.
  * [66] Alec Radford, Jeff Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever.  Language models are unsupervised multitask learners.  2019\.
  * [67] Nelson Elhage, Tristan Hume, Catherine Olsson, Nicholas Schiefer, Tom Henighan, Shauna Kravec, Zac Hatfield-Dodds, Robert Lasenby, Dawn Drain, Carol Chen, Roger Grosse, Sam McCandlish, Jared Kaplan, Dario Amodei, Martin Wattenberg, and Christopher Olah.  Toy models of superposition.  Transformer Circuits Thread, 2022.  https://transformer-circuits.pub/2022/toy_model/index.html.
  * [68] Christos Thrampoulidis, Ganesh Ramachandra Kini, Vala Vakilian, and Tina Behnia.  Imbalance trouble: Revisiting neural-collapse geometry.  Advances in Neural Information Processing Systems, 35:27225–27238, 2022.
  * [69] Nagarajan Natarajan, Inderjit S Dhillon, Pradeep K Ravikumar, and Ambuj Tewari.  Learning with noisy labels.  Advances in neural information processing systems, 26, 2013.
  * [70] Quinn Fisher, Haoming Meng, and Vardan Papyan.  Pushing boundaries: Mixup’s influence on neural collapse, 2024.
  * [71] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [72] David Mimno and Laure Thompson.  The strange geometry of skip-gram with negative sampling.  In Martha Palmer, Rebecca Hwa, and Sebastian Riedel, editors, Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing, pages 2873–2878, Copenhagen, Denmark, September 2017. Association for Computational Linguistics.
  * [73] Kawin Ethayarajh.  How contextual are contextualized word representations? comparing the geometry of bert, elmo, and gpt-2 embeddings, 2019.
  * [74] Alex Wang, Amanpreet Singh, Julian Michael, Felix Hill, Omer Levy, and Samuel R. Bowman.  Glue: A multi-task benchmark and analysis platform for natural language understanding, 2019.
  * [75] Jia Hui Feng, Edmund M-K Lai, and Weihua Li.  A study of neural collapse for text classification.  In International Conference on Deep Learning Theory and Applications, pages 126–142. Springer, 2023.
  * [76] Thomas Laurent, James H. von Brecht, and Xavier Bresson.  Feature collapse, 2023.
  * [77] Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, et al.  Language models are few-shot learners.  Advances in neural information processing systems, 33:1877–1901, 2020.
  * [78] Sid Black, Leo Gao, Phil Wang, Connor Leahy, and Stella Biderman.  GPT-Neo: Large Scale Autoregressive Language Modeling with Mesh-Tensorflow, March 2021.  If you use this software, please cite it using these metadata.
  * [79] Susan Zhang, Stephen Roller, Naman Goyal, Mikel Artetxe, Moya Chen, Shuohui Chen, Christopher Dewan, Mona Diab, Xian Li, Xi Victoria Lin, Todor Mihaylov, Myle Ott, Sam Shleifer, Kurt Shuster, Daniel Simig, Punit Singh Koura, Anjali Sridhar, Tianlu Wang, and Luke Zettlemoyer.  Opt: Open pre-trained transformer language models, 2022.
  * [80] Hugo Touvron, Thibaut Lavril, Gautier Izacard, Xavier Martinet, Marie-Anne Lachaux, Timothée Lacroix, Baptiste Rozière, Naman Goyal, Eric Hambro, Faisal Azhar, Aurelien Rodriguez, Armand Joulin, Edouard Grave, and Guillaume Lample.  Llama: Open and efficient foundation language models, 2023.
  * [81] Albert Q. Jiang, Alexandre Sablayrolles, Arthur Mensch, Chris Bamford, Devendra Singh Chaplot, Diego de las Casas, Florian Bressand, Gianna Lengyel, Guillaume Lample, Lucile Saulnier, Lélio Renard Lavaud, Marie-Anne Lachaux, Pierre Stock, Teven Le Scao, Thibaut Lavril, Thomas Wang, Timothée Lacroix, and William El Sayed.  Mistral 7b, 2023.
  * [82] Ronald A Fisher.  The use of multiple measurements in taxonomic problems.  Annals of eugenics, 7(2):179–188, 1936.
  * [83] C Radhakrishna Rao.  The utilization of multiple measurements in problems of biological classification.  Journal of the Royal Statistical Society. Series B (Methodological), 10(2):159–203, 1948.
  * [84] Thomas Strohmer and Robert W Heath Jr.  Grassmannian frames with applications to coding and communication.  Applied and computational harmonic analysis, 14(3):257–275, 2003\.
  * [85] Shayne FD Waldron.  An introduction to finite tight frames.  Springer, 2018.
  * [86] David L Donoho, Michael Elad, and Vladimir N Temlyakov.  Stable recovery of sparse overcomplete representations in the presence of noise.  IEEE Transactions on information theory, 52(1):6–18, 2005.
  * [87] Joel A Tropp.  Just relax: Convex programming methods for identifying sparse signals in noise.  IEEE transactions on information theory, 52(3):1030–1051, 2006\.
  * [88] Vardan Papyan.  Traces of class/cross-class structure pervade deep learning spectra, 2020\.
  * [89] Christopher R. Hoyt and Art B. Owen.  Probing neural networks with t-sne, class-specific projections and a guided tour, 2021.
  * [90] John Zarka, Florentin Guth, and Stéphane Mallat.  Separation and concentration in deep networks, 2021.
  * [91] Tomer Galanti, Liane Galanti, and Ido Ben-Shaul.  On the implicit bias towards minimal depth of deep neural networks, 2022\.
  * [92] Ido Ben-Shaul and Shai Dekel.  Nearest class-center simplification through intermediate layers.  In Alexander Cloninger, Timothy Doster, Tegan Emerson, Manohar Kaul, Ira Ktena, Henry Kvinge, Nina Miolane, Bastian Rieck, Sarah Tymochko, and Guy Wolf, editors, Proceedings of Topological, Algebraic, and Geometric Learning Workshops 2022, volume 196 of Proceedings of Machine Learning Research, pages 37–47. PMLR, 25 Feb–22 Jul 2022.
  * [93] Hangfeng He and Weijie J. Su.  A law of data separation in deep learning, 2022.
  * [94] Liam Parker, Emre Onal, Anton Stengel, and Jake Intrater.  Neural collapse in the intermediate hidden layers of classification neural networks, 2023.
  * [95] Akshay Rangamani, Marius Lindegaard, Tomer Galanti, and Tomaso A Poggio.  Feature learning in deep classifiers through intermediate neural collapse.  In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, Proceedings of the 40th International Conference on Machine Learning, volume 202 of Proceedings of Machine Learning Research, pages 28729–28745. PMLR, 23–29 Jul 2023.
  * [96] Wojciech Masarczyk, Mateusz Ostaszewski, Ehsan Imani, Razvan Pascanu, Piotr Mił oś, and Tomasz Trzcinski.  The tunnel effect: Building data representations in deep neural networks.  In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine, editors, Advances in Neural Information Processing Systems, volume 36, pages 76772–76805. Curran Associates, Inc., 2023.
  * [97] Daniel Beaglehole, Peter Súkeník, Marco Mondelli, and Mikhail Belkin.  Average gradient outer product as a mechanism for deep neural collapse, 2024.
  * [98] Peng Wang, Xiao Li, Can Yaras, Zhihui Zhu, Laura Balzano, Wei Hu, and Qing Qu.  Understanding deep representation learning via layerwise feature compression and discrimination, 2024.
  * [99] Sicong Wang, Kuo Gai, and Shihua Zhang.  Progressive feedforward collapse of resnet training, 2024.
  * [100] Connall Garrod and Jonathan P. Keating.  Unifying low dimensional observations in deep learning through the deep linear unconstrained feature model, 2024.
  * [101] Emanuele Zangrando, Piero Deidda, Simone Brugiapaglia, Nicola Guglielmi, and Francesco Tudisco.  Neural rank collapse: Weight decay and small within-class variability yield low-rank bias, 2024.
  * [102] Mengjia Xu, Akshay Rangamani, Qianli Liao, Tomer Galanti, and Tomaso Poggio.  Dynamics in deep classifiers trained with the square loss: Normalization, low rank, neural collapse, and generalization bounds.  Research, 6:0024, 2023.
  * [103] Tong Liang and Jim Davis.  Inducing neural collapse to a fixed hierarchy-aware frame for reducing mistake severity.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 1443–1452, 2023.
  * [104] Guglielmo Bonifazi, Iason Chalas, Gian Hess, and Jakub Łucki.  Can we understand plasticity through neural collapse?, 2024.
  * [105] Ilya Sutskever, Aug 2023.
  * [106] Li Guo, Keith Ross, Zifan Zhao, George Andriopoulos, Shuyang Ling, Yufeng Xu, and Zixuan Dong.  Cross entropy versus label smoothing: A neural collapse perspective, 2024\.
  * [107] Neil Burgess, Jelena Milanovic, Nigel Stephens, Konstantinos Monachopoulos, and David Mansell.  Bfloat16 processing for neural networks.  In 2019 IEEE 26th Symposium on Computer Arithmetic (ARITH), pages 88–91, 2019.

##  Appendix A Architectural Details

Table 2: Sample architectural configuration for a $12$-layer $1024$-dimensional causal language model (CLM) based on [2] and GPT-Neo [78]. Shallower models have configurations with attention_layers and attention_types truncated. Narrower models adjust hidden_size to their width ($d$). All other configuration values are the same across models. Setting | Value
---|---
activation_function | gelu_new
architectures | GPTNeoForCausalLM
attention_dropout | 0
attention_layers | global, local, global, local, …
attention_types | [[global, local], 6]
bos_token_id | 50256
embed_dropout | 0
eos_token_id | 50256
gradient_checkpointing | false
hidden_size | 1024
initializer_range | 0.02
intermediate_size | null
layer_norm_epsilon | 1e-05
max_position_embeddings | 2048
model_type | gpt_neo
num_heads | 16
num_layers | 12
resid_dropout | 0
summary_activation | null
summary_first_dropout | 0.1
summary_proj_to_labels | true
summary_type | cls_index
summary_use_proj | true
torch_dtype | float32
transformers_version | 4.28.1
use_cache | true
vocab_size | 50257
window_size | 256

##  Appendix B Dataset

TinyStories is a synthetic dataset of short children’s stories generated by
GPT-3.5 and GPT-4 [2].

  * •

The dataset is released with the CDLA-Sharing-1.0 licence.

  * •

Our models were trained and evaluated using the first version of TinyStories
[2].

  * •

The $2,141,709$ stories are split into $2,119,719$ train and $21,990$
validation stories.

  * •

Their experimental setup [2] called for the GPT-2 [66] tokenizer, of which
only a subset vocabulary $\mathcal{V}$ of size $29,233$ appears in
TinyStories.

  * •

Following the GPT-style regime for training/evaluation [5], raw sequences
(stories) from the train set are packed into $229,367$ ($S$) chunks of $2048$
($T$) tokens each. This setup provides $469,514,249$ ($N$) ground-
truth161616$N=S(T-1)$ as we cannot use the first ground-truth nor the last
predicted token in any chunk. token samples for training.

  * •

Dataloading and preprocessing (including chunking) each used two workers.

Figure 4: The $500$ most frequent classes from TinyStories [2] exhibit
significant sample imbalance.

##  Appendix C Optimization

Training was performed using an adaptation of an open-source causal language
modeling script from Huggingface:
https://github.com/huggingface/transformers/blob/main/examples/pytorch/language-
modeling/run_clm.py

  * •

Each model was trained on a single NVIDIA A100 (40GB) GPU for up to 8 hours
per epoch.

  * •

Learning rates were set by a linear schedule based on the number of steps with
no warm-up.

  * •

Training was performed in bfloat16 [107] mixed precision.

  * •

The results presented in this work are from two sets of models trained with
weight decay $\beta=0.0005$ [50] and $\beta=0.1$ [77]. A previous set of
models was trained without weight decay and the results are very similar to
$\beta=0.0005$.

Table 3: Batch sizes used to train models on a single NVIDIA A100 (40GB) GPU. Width ($d$) and depth ($L$) correspond to hidden_size and length of attention_layers, respectively, in Table 2. Depth $(L)\downarrow$ Width $(d)\rightarrow$ | $64$ | $128$ | $256$ | $512$ | $768$ | $1024$
---|---|---|---|---|---|---
$1$-layer | $16$ | $16$ | $16$ | $16$ | $16$ | $16$
$2$-layer | $16$ | $16$ | $16$ | $16$ | $16$ | $8$
$4$-layer | $8$ | $8$ | $8$ | $8$ | $8$ | $8$
$8$-layer | $8$ | $8$ | $8$ | $4$ | $4$ | $4$
$12$-layer | $4$ | $4$ | $4$ | $4$ | $4$ | $4$

##  Appendix D Within-Class Variability Collapse with Scale – $\mathcal{NC}$1

Figure 5: Average (logarithmic) class-distance normalized variance (CDNV) is
reduced ($\mathcal{NC}$1) when scaling width ($d$) and across training for 1
(left) through 10 (right) epochs with weight decays $\beta=0.0005$ (top) and
$\beta=0.1$ (bottom).

Figure 6: Average (logarithmic) class-distance normalized variance (CDNV) is
reduced ($\mathcal{NC}$1) when scaling depth ($L$) and across training for 1
(left) through 10 (right) epochs with weight decays $\beta=0.0005$ (top) and
$\beta=0.1$ (bottom).

##  Appendix E Mean Norms Growth with Scale – (Related to $\mathcal{NC}$2)

Figure 7: Class embedding mean norms grow when scaling width ($d$) and across
training for 1 (left) through 10 (right) epochs with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom).

Figure 8: Class embedding mean norms grow when scaling depth ($L$) and across
training for 1 (left) through 10 (right) epochs with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom).

##  Appendix F Equinormness with Scale – $\mathcal{NC}$2

Figure 9: Variation in embedding norms decreases ($\mathcal{NC}$2) when
scaling width ($d$) in models trained for 1 (left) through 10 (right) with
weight decays $\beta=0.0005$ (top) and $\beta=0.1$ (bottom). Note that the
degree of equinormness eventually plateaus for sufficiently deep and trained
models.

Figure 10: Variation in embedding norms decreases ($\mathcal{NC}$2,) when
scaling depth ($L$) in models trained for 1 (left) through 10 (right) with
weight decays $\beta=0.0005$ (top) and $\beta=0.1$ (bottom).

##  Appendix G Interference with Scale – (Related to $\mathcal{NC}$2)

Figure 11: Average interference decreases (to some extent) when scaling width
($d$) in models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom).

Figure 12: Average interference decreases (to some extent) when scaling depth
($L$) in models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom).

##  Appendix H Equiangularity with Scale – (Against $\mathcal{NC}$2)

Figure 13: Variation in interference roughly increases when scaling width
($d$) in models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom). Note this trend is against
equiangularity, affirming the traditional $\mathcal{NC}$2 to be less useful
than $\mathcal{GNC}$2 [16].

Figure 14: Variation in interference increases when scaling depth ($L$) in
models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom). Note this trend is against
equiangularity, affirming the traditional $\mathcal{NC}$2 to be less useful
than $\mathcal{GNC}$2 [16].

##  Appendix I Logarithmic Distances with Scale – (Related to
$\mathcal{GNC}$2)

Figure 15: Average logarithmic distance decreases when scaling width ($d$) in
models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom).

Figure 16: Average logarithmic distance decreases when scaling depth ($L$) in
models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom).

##  Appendix J Hyperspherical Uniformity with Scale – $\mathcal{GNC}$2

Figure 17: Variation in logarithmic distances decreases when scaling width
($d$) in models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom). This consistent trend towards
hyperspherical uniformity affirms that $\mathcal{GNC}$2 [16] is more useful
than $\mathcal{NC}$2.

Figure 18: Variation in logarithmic distances decreases when scaling depth
($L$) in models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom). This consistent trend towards
hyperspherical uniformity affirms that $\mathcal{GNC}$2 [16] is more useful
than $\mathcal{NC}$2.

##  Appendix K Self-Duality with Scale – (Against $\mathcal{NC}$3)

Self-duality ($\mathcal{NC}3$) was originally the convergence of classifiers
to class means after normalization:

$\left\|\frac{\boldsymbol{w}_{c}}{\|\boldsymbol{w}_{c}\|_{2}}-\frac{\boldsymbol{\mu}_{c}-\boldsymbol{\bar{\mu}}}{\|\boldsymbol{\mu}_{c}-\boldsymbol{\bar{\mu}}\|_{2}}\right\|_{2}\rightarrow
0,\quad\forall c.$ (11)

Instead, we use class-wise similarity (Equation 8) and its variation
($\mathcal{UNC}$3).

Figure 19: Average classifier alignment increases $\mathcal{NC}$3 when
training for 1 (left) through 10 (right) with weight decays $\beta=0.0005$
(top) and $\beta=0.1$ (bottom). However, we see no meaningful trend when
scaling width $d$, suggesting that $\mathcal{NC}$3 does not coalesce with
language modeling training.

Figure 20: Average classifier alignment increases $\mathcal{NC}$3 when
training for 1 (left) through 10 (right) with weight decays $\beta=0.0005$
(top) and $\beta=0.1$ (bottom). However, we see no meaningful trend when
scaling depth $L$, suggesting that $\mathcal{NC}$3 does not coalesce with
language modeling training.

##  Appendix L Uniformity Duality with Scale – $\mathcal{UNC}$3

Figure 21: Variation in classifier alignment decreases when scaling width
($d$) in models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom).

Figure 22: Variation in classifier alignment increases when scaling depth
($L$) in models trained for 1 (left) through 10 (right) with weight decays
$\beta=0.0005$ (top) and $\beta=0.1$ (bottom). This negative trend of
$\mathcal{UNC}$3 in more learnt models (right) suggests that the link of
$\mathcal{(U)NC}$3 with scale and performance is still weak.

##  Appendix M Classifier Agreement – $\mathcal{NC}$4

For computational reasons, we compute Equations 9, 10 using a simple
decomposition:

$\operatorname*{argmin}_{c\in\mathbb{V}}\|\boldsymbol{h}_{b}-\boldsymbol{\mu}_{c}\|_{2}=\operatorname*{argmin}_{c\in\mathbb{V}}\left(\|\boldsymbol{h}_{b}\|^{2}+\|\boldsymbol{\mu}_{c}\|^{2}-2\boldsymbol{h}_{b}^{\top}\boldsymbol{\mu}_{c}\right)$
(12)

where $b\in[1,B]$ and $c\in\mathbb{V}$ with batch size $B$ and vocabulary
$\mathbb{V}$.

Figure 23: Classifier agreement improves when scaling width ($d$) in models
trained for 1 (left) through 10 (right) with weight decays $\beta=0.0005$
(top) and $\beta=0.1$ (bottom).

Figure 24: Classifier agreement improves when scaling depth ($L$) in models
trained for 1 (left) through 10 (right) with weight decays $\beta=0.0005$
(top) and $\beta=0.1$ (bottom).

##  Appendix N Correlations with Generalization Performance

Figure 25: Validation loss shows some correlation with average mean norms
(left) and their variations ($\mathcal{NC}_{2}$) (right). Figure 26:
Validation loss shows weak correlation with average interference and almost
none with its variation (i.e. equiangularity, $\mathcal{NC}_{2}$) (right).
Figure 27: Validation loss shows some correlation with average kernel
distances and with its variation (i.e. hyperspherical uniformity,
$\mathcal{GNC}_{2}$) (right).


* (113) X. Yin, Y. Zhu, and J. Hu, “A comprehensive survey of privacy-preserving federated learning: A taxonomy, review, and future directions,” _ACM Computing Surveys (CSUR)_ , vol. 54, no. 6, pp. 1–36, 2021.
  * (114) T. Orekondy, S. J. Oh, Y. Zhang, B. Schiele, and M. Fritz, “Gradient-leaks: Understanding and controlling deanonymization in federated learning,” _arXiv preprint arXiv:1805.05838_ , 2018.
  * (115) Y. Qu, L. Gao, T. H. Luan, Y. Xiang, S. Yu, B. Li, and G. Zheng, “Decentralized privacy using blockchain-enabled federated learning in fog computing,” _IEEE Internet of Things Journal_ , vol. 7, no. 6, pp. 5171–5183, 2020.
  * (116) L. Su and J. Xu, “Securing distributed gradient descent in high dimensional statistical learning,” _Proceedings of the ACM on Measurement and Analysis of Computing Systems_ , vol. 3, no. 1, pp. 1–41, 2019.
  * (117) D. Enthoven and Z. Al-Ars, “An overview of federated deep learning privacy attacks and defensive strategies,” _Federated Learning Systems: Towards Next-Generation AI_ , pp. 173–196, 2021.


# YUNIC: A Multi-Dimensional Particle-In-Cell Code for Laser-Plasma
Interaction

Huai-Hang Song  Beijing National Laboratory for Condensed Matter Physics,
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China  Wei-
Min Wang  Department of Physics and Beijing Key Laboratory of Opto-electronic
Functional Materials and Micro–nano Devices, Renmin University of China,
Beijing 100872, China  Yu-Tong Li  Beijing National Laboratory for Condensed
Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing
100190, China

###### Abstract

For simulating laser-plasma interactions, we developed a parallel, multi-
dimensional, fully relativistically particle-in-cell (PIC) code, named yunic.
The core algorithm is introduced, including field solver, particle pusher,
field interpolation, and current interpolation. In addition to the classical
electromagnetic interaction in plasmas, nonlinear Compton scattering and
nonlinear Breit-Wheeler pair production are also implemented based on Monte-
Carlo methods to study quantum electrodynamics (QED) processes. We benchmark
yunic against theories and other PIC codes through several typical cases.
Yunic can be applied in varieties of physical scenes, from relativistic laser-
plasma interactions to astrophysical plasmas and strong-field QED physics.

###### pacs:

###### Contents

  1. I Introduction
  2. II Standard particle-in-cell algorithm
    1. II.1 Normalization
    2. II.2 Field solver
    3. II.3 Particle pusher
    4. II.4 Field interpolation
    5. II.5 Current interpolation
    6. II.6 Algorithm structure
  3. III Benchmarks in several typical cases
  4. IV QED modules
  5. V Conclusion

##  I Introduction

The particle-in-cell (PIC) method Birdsall and Langdon (1991), belonging to a
kinetic description, can accurately simulate the collective plasma behavior,
from linear to relativistically nonlinear processes. Compared to
magnetohydrodynamics simulation, PIC simulation utilizing the quasi-particle
concept can resolve plasma dynamics on smaller spatial and temporal scales.
Meanwhile, it requires much less computational expense than that by directly
solving Vlasov-Boltzmann equations. Since it was initially developed in 1970s
Dawson (1983), PIC method has become one of the most powerful and
indispensable tools in various plasma areas, particularly in laser-plasma
interactions Gibbon (2005); Macchi (2013). In this manuscript, we introduce a
recently developed PIC code yunic, and demonstrate a few typical benchmarks to
validate this code. Besides the classical plasma dynamics, yunic is also
capable of simulating extremely laser-plasma interactions in the quantum
electrodynamics (QED) regime, including spin and polarization effects in the
processes of nonlinear Compton scattering and nonlinear Breit-Wheeler pair
production Baier  _et al._ (1998).

##  II Standard particle-in-cell algorithm

The core idea of PIC method is to solve Maxwell’s equations on the discrete
spatial grid [Sec. II.2], while pushing quasi-particles in the free space
[Sec. II.3]. The currents generated by moving charged particles should be
interpolated to the spatial grid as sources to solve field equations [Sec.
II.5], and fields should also be interpolated back to an arbitrary particle
position to push them [Sec. II.4]. Hence, the discrete fields and non-discrete
particles are self-consistently connected. The main equations for solving
collisionless plasma problems are as follows (in Gaussian units):

$\displaystyle\nabla\cdot{\bf E}=4\pi\rho,$ (1) $\displaystyle\nabla\cdot{\bf
B}=0,$ (2) $\displaystyle\nabla\times{\bf E}=-\frac{1}{c}\frac{\partial\bf
B}{\partial t},$ (3) $\displaystyle\nabla\times{\bf
B}=\frac{1}{c}\left(4\pi{\bf J}+\frac{\partial\bf E}{\partial t}\right),$ (4)
$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot{\bf J}=0,$ (5)
$\displaystyle\frac{d{\bf p}}{dt}=q({\bf E}_{p}+\frac{\bf p}{\gamma
mc}\times{\bf B}_{p}),$ (6) $\displaystyle\frac{d{\bf r}_{p}}{dt}=\frac{\bf
p}{\gamma m},$ (7)

where Eqs. (1)-(4) are Maxwell’s equations, Eq. (5) is the charge conservation
equation, Eqs. (6) and (7) are Newton-Lorentz equations.

###  II.1 Normalization

In the PIC code, after defining a reference frequency $\omega_{r}$, it is
convenient to normalize time $t$, length $L$, electric field $\bf E$, magnetic
field $\bf B$, particle velocity $\bf v$, momentum $\bf p$, number density
$n$, current density $\bf J$, charge $q$ and mass $m$ to following quantities,
respectively:

$\displaystyle t_{r}$ $\displaystyle=$
$\displaystyle\frac{1}{\omega_{r}},\quad L_{r}=\frac{c}{\omega_{r}},\quad
E_{r}=\frac{m_{e}\,c\,\omega_{r}}{e},\quad
B_{r}=\frac{m_{e}\,c\,\omega_{r}}{e},\quad v_{r}=c,$ (8) $\displaystyle p_{r}$
$\displaystyle=$ $\displaystyle m_{e}\,c,\quad
n_{r}=\frac{m_{e}\,\omega_{r}^{2}}{4\pi e^{2}},\quad J_{r}=c\,e\,n_{r},\quad
q_{r}=e,\quad m_{r}=m_{e},$

where $m_{e}$ is the electron rest mass, and $e$ is the electron charge. In
laser-plasma interactions, $\omega_{r}$ is usually chosen as the laser
frequency $\omega_{L}$, and it could also be chosen as the plasma oscillating
frequency $\omega_{p}$ in other interactions. For convenience, we use
normalized quantities below based on Eq. (8).

###  II.2 Field solver

Figure 1:  The staggered Yee grid in (a) one dimensional (1D), (b) two
dimensional (2D), and (c) three dimensional (3D) spaces for solving Maxwell’s
equations, respectively.

The electromagnetic fields ${\bf E}$ and ${\bf B}$ are self-consistently
evolved by solving Maxwell’s equations. One only needs to solve two curl
equations [Eqs. (3) and (4)] and a charge conservation equation [Eq. (5)],
because the other two divergence equations [Eqs. (1) and (2)] are
automatically satisfied with time if they hold initially Eastwood (1991);
Villasenor and Buneman (1992). The details of how to realize the charge
conservation are presented in Sec. II.5.

The finite-difference time-domain (FDTD) method Taflove and Hagness (2000) is
adopted for numerically solving Maxwell’s equations. Equations. (3) and (4)
can be written as the following discrete forms with the 2nd-order accuracy:

$\displaystyle{\bf E}^{n+1/2}-{\bf E}^{n-1/2}$ $\displaystyle=$
$\displaystyle\Delta t(\nabla\times{\bf B}^{n}-{\bf J}^{n}),$ (9)
$\displaystyle{\bf B}^{n+1}-{\bf B}^{n}$ $\displaystyle=$
$\displaystyle-\Delta t\nabla\times{\bf E}^{n+1/2}.$ (10)

Here, the leapfrog scheme is adopted in time and the famous staggered Yee grid
Kane Yee (1966) is employed in space, as illustrated in Fig. 1.

###  II.3 Particle pusher

We proceed to solve the motion of charged particles in the electromagnetic
field by discretizing Eqs. (6) and (7) as following:

$\displaystyle\frac{{\bf u}^{n+1}-{\bf u}^{n}}{\Delta t}$ $\displaystyle=$
$\displaystyle q_{m}({\bf E}_{p}^{n+1/2}+\frac{1}{{\bf\gamma}^{n+1/2}}{\bf
u}^{n+1/2}\times{\bf B}_{p}^{n+1/2}),$ (11) $\displaystyle\frac{{\bf
r}_{p}^{n+3/2}-{\bf r}_{p}^{n+1/2}}{\Delta t}$ $\displaystyle=$
$\displaystyle\frac{{\bf u}^{n+1}}{\gamma^{n+1}},$ (12)

where ${\bf u}={\bf p}/m$ and $q_{m}=q/m$.

Figure 2:  Diagram of Boris algorithm.

Boris algorithm Boris (1970) is employed to push charged particles due to its
advantage of long term accuracy Qin  _et al._ (2013), which splits the
electric and magnetic forces by defining ${\bf u}^{+}$ and ${\bf u}^{-}$,

$\displaystyle{\bf u}^{n}$ $\displaystyle=$ $\displaystyle{\bf
u}^{-}-\frac{q_{m}\Delta t}{2}{\bf E}_{p}^{n+1/2},$ (13) $\displaystyle{\bf
u}^{n+1}$ $\displaystyle=$ $\displaystyle{\bf u}^{+}+\frac{q_{m}\Delta
t}{2}{\bf E}_{p}^{n+1/2}.$ (14)

Substituting Eqs. (13) and (14) into Eq. (11), one can obtain a rotation
equation of $\bf u^{+}$ and $\bf u^{-}$ about the magnetic field ${\bf
B}_{p}$, i.e., $2\gamma^{n+1/2}({\bf u}^{+}-{\bf u}^{-})=q_{m}{\Delta t}({\bf
u}^{+}+{\bf u}^{-})\times{\bf B}_{p}^{n+1/2}$, and then solve it by the
following implementation:

$\displaystyle{\bf u}^{\prime}$ $\displaystyle=$ $\displaystyle{\bf
u}^{-}+{\bf u}^{-}\times{\bf t},$ (15) $\displaystyle{\bf u}^{+}$
$\displaystyle=$ $\displaystyle{\bf u}^{-}+{\bf u}^{\prime}\times\frac{2{\bf
t}}{1+t^{2}}$ (16)

where ${\bf t}=q_{m}\Delta t{\bf B}_{p}^{n+1/2}/{(2\gamma^{n+1/2})}$. The
diagram of Boris algorithm is illustrated in Fig. 2.

If contributions of electric field and magnetic field to the Lorentz force
nearly cancel out, e.g., in the rest frame of an ultrarelativistic beam Vay
(2007), Boris algorithm might lead to large errors, and Vay algorithm Vay
(2008) is more suitable.

###  II.4 Field interpolation

Figure 3:  (a) Shape functions $S_{n}(\xi)$ of 2nd-order ($n=2$) and 4th-order
($n=4$) used in PIC for field and current interpolations. (b) (c) The area of
the shadow $W_{n}(i-i_{c})$ represents the weight of the grid field acting on
the particle for $n=2$ and $n=4$, respectively. Normally, weight function
$W_{2}(i-i_{c})$ crosses three grid cells, and $W_{4}(i-i_{c})$ crosses five
grid cells.

In Sec. II.3, we have discussed how to push a charged particle provided we
have known the electric field ${\bf E}_{p}$ and magnetic field ${\bf B}_{p}$
at the particle position ${\bf r}_{p}(x_{p},y_{p},z_{p})$. In this section, we
discuss how to obtain ${\bf E}_{p}$ and ${\bf B}_{p}$ through the
interpolation. Actually, it depends on the interpolation shape function
$S_{n}(\xi)$ we choose, where $n$ is the interpolation order. The function
$S_{n}(\xi)$ defines the shape and smoothness of quasi-particles and also
determines the simulation accuracy Birdsall and Langdon (1991). Taking 1D as
an example, its 2nd-order and 4th-order forms shown in Fig. 3(a) are given by
Abe  _et al._ (1986)

$\displaystyle S_{2}(\xi)$ $\displaystyle=$
$\displaystyle\begin{dcases}1-|\xi|&\text{if}\,\,|\xi|\leq 1,\\\
0&\text{otherwise},\end{dcases}$ (17) $\displaystyle S_{4}(\xi)$
$\displaystyle=$
$\displaystyle\begin{dcases}\frac{2}{3}-|\xi|^{2}+\frac{1}{2}|\xi|^{3}&\text{if}\,\,|\xi|\leq
1,\\\ \frac{4}{3}(1-\frac{1}{2}|\xi|)^{3}&\text{if}\,\,1<|\xi|\leq 2,\\\
0&\text{otherwise}.\end{dcases}$ (18)

For a particle located at $x_{p}$, its acting field contributed by grid point
$i$ can be expressed as $F(i)*W_{n}(i-i_{c})$, where its weight function is
defined by $W_{n}(i-i_{c})=\int_{i-1/2}^{i+1/2}S_{n}(\xi-x_{p}/\Delta x)d\xi$
and $i_{c}$ is the grid point nearest to the particle, as shown in Figs. 3(b)
and 3(c). The expression of $W_{n}(i-i_{c})$ after the integration can be
found from APPENDIX A of Abe  _et al._ (1986). In 3D, the electromagnetic
field acting on the particle can be calculated through following
interpolations under the condition that the fields are constant over each
cell:

$\displaystyle{\bf E}({\bf
B})_{p}(x_{p},y_{p},z_{p})=\sum_{i,j,k}W_{n}(i-i_{c})*W_{n}(j-j_{c})*W_{n}(k-k_{c})*{\bf
E}({\bf B})(i,j,k).$ (19)

Notice that different field components generally have different weights since
the grid is staggered [see Fig. 1].

###  II.5 Current interpolation

As we mentioned in Sec. II.2, one needs to ensure the charge conservation in
order to avoid solving Poisson’s equation [Eq. (1)] Villasenor and Buneman
(1992), since the local computation of the former is much simpler and
computationally cheaper than the global computation of the latter. In yunic,
Esirkepov algorithm Esirkepov (2001) is adopted to ensure the charge
conservation in the current calculation. By assuming the particle trajectory
over one time step is a straight line, the current flux in 3D is decomposed
into twelve segments along $x$, $y$, and $z$ axes, respectively, i.e.,

$\displaystyle{\rm J}_{x}^{n+1}(i+1,j,k)-{\rm J}_{x}^{n+1}(i,j,k)$
$\displaystyle=$ $\displaystyle-q\frac{\Delta x}{\Delta t}W_{x}(i,j,k),$ (20)
$\displaystyle{\rm J}_{y}^{n+1}(i,j+1,k)-{\rm J}_{y}^{n+1}(i,j,k)$
$\displaystyle=$ $\displaystyle-q\frac{\Delta y}{\Delta t}W_{y}(i,j,k),$ (21)
$\displaystyle{\rm J}_{z}^{n+1}(i,j,k+1)-{\rm J}_{z}^{n+1}(i,j,k)$
$\displaystyle=$ $\displaystyle-q\frac{\Delta z}{\Delta t}W_{z}(i,j,k),$ (22)

where $W_{x}(i,j,k)$, $W_{y}(i,j,k)$ and $W_{z}(i,j,k)$ can be found from Eq.
(31) in Esirkepov (2001).

This method is easily extended to an arbitrary high-order shape function
$S_{n}(\xi)$. Note that $S_{n}(\xi)$ employed in the current interpolation
should be the same as that in the field interpolation [Sec. II.4] to eliminate
the self-force.

Figure 4:  Flow chart of the core algorithm in yunic.

###  II.6 Algorithm structure

Yunic is written in C++ language and massively parallelized by MPI. The
simulation results are output in parallel by MPI-IO and then
analyzed/visualized by Python/Matplotlib. Its three versions aimed at
different spatial dimensions (1D, 2D, and 3D) are constructed separately for
efficiency. Yunic employs a modified algorithm that originally adopted in psc
Ruhl  by Hartmut Ruhl, and latter also adopted in epoch Arber  _et al._
(2015), which is slightly different from the common algorithm as described in
Secs. II.2-II.5. The modified algorithm updates the fields and particle
positions at both full-time steps and half-time steps, as sketched in Fig. 4.
Hence, the drawback of the leapfrog method is overcome and one can obtain the
information of fields and particles at the same time, which is important in
some cases, e.g., for simulating QED processes.

##  III Benchmarks in several typical cases

Here, we benchmark our PIC code yunic against the open-source PIC code smilei
Derouillat  _et al._ (2018), including 1D, 2D and 3D versions. At these
presented cases, the simulation results of two PIC codes are in good
agreement.

(a) _1D PIC simulation: hole boring_. A linearly (circularly) polarized laser
with a wavelength of $\lambda_{L}=1~{}\mu\rm m$ incidents from the left
boundary at $t=0$ ps. The laser has a 5-laser-period rise time before a
constant normalized intensity of
$a_{0}=eE_{L}/m_{e}c\omega_{L}=10{\sqrt{2}}~{}(10)$. An uniform overdense
plasma with an electron density of $n_{p}=50n_{c}$ is located at $5~{}\mu{\rm
m}<x<15~{}\mu{\rm m}$, where $n_{c}=m_{e}\omega_{L}^{2}/4\pi e^{2}$. The
computational domain has a size of $20\lambda_{L}$ with 5120 cells in the $x$
direction. Each cell contains 100 electrons and 100 protons. Absorbing
boundaries are used for both particles and fields. The 4th-order interpolation
is employed. The comparison of ion charge density $n_{i}$ at different times
are shown in Fig. 5.

Figure 5:  Spatial distributions of ion density $n_{i}$ at different times,
(a)(d) $t=0.17$ ps, (b)(e) $t=0.33$ ps, and (c)(f) $t=0.5$ ps. (a)-(c) and
(d)-(f) correspond to circularly and linearly polarized lasers, respectively.
Figure 6:  Time evolution of the relative increase of total energy $\Delta
E/E_{0}$ with different interpolation orders and particles per cell (ppc). The
colored thick lines are the simulation results by yunic, and black thin lines
correspond to those by smilei. Figure 7:  (a)(d) Electron density $n_{e}$ of
the background plasma. (b)(e) Longitudinal electric field $E_{x}$. (c)(f) On-
axis distribution of $E_{x}$. (a)-(c) and (d)-(f) are driven by electron and
positron beams, respectively.

(b) _2D PIC simulation: self-heating_. A plasma with an initial temperature of
1 keV has an uniform electron density of $n_{p}=100n_{c}$. Each plasma
wavelength ($\lambda_{p}=2\pi/k_{p}$) contains 320 grids. Periodic boundaries
are applied for both particles and fields. No current smoothing or other
additional algorithms are used to control self-heating. The comparison of
relative energy increasing $\Delta E/E_{0}$ due to self-heating are shown in
Fig. 6.

(c) _3D PIC simulation: wakefield driven by electron or positron beam_. A
1-GeV drive electron (positron) beam has a bi-Gaussian density profile with a
transverse size $k_{p}\sigma_{r}=0.5$, bunch length $k_{p}\sigma_{x}=0.5$, and
peak density $n_{b}/n_{p}=4$. The computational domain has a size of
$4\lambda_{p}\times 2\lambda_{p}\times 2\lambda_{p}$ in $x\times y\times z$
directions, sampled by $512\times 192\times 192$ cells. Each cell contains 8
electrons and 8 protons for the uniform background plasma and 4 electrons or
positrons for the drive beam. For the field initialization of an
ultrarelativistic charged beam, yunic first solves the Poisson’s equation in
the beam’s rest frame, and then applies Lorentz transformation to obtain its
self-generated fields in the laboratory frame Massimo  _et al._ (2016). The
comparison of background plasma density $n_{e}$ and excited longitudinal
electric field $E_{x}$ are shown in Fig. 7.

##  IV QED modules

The available petawatt and next-generation 10-petawatt and 100-petawatt laser
systems can provide an extreme field density of $I_{0}=10^{22-25}\rm
W/cm^{2}$. To explore interactions of ultraintense lasers with plasmas, yunic
have implemented QED modules with Monte-Carlo methods Elkina  _et al._ (2011);
Ridgers  _et al._ (2014); Gonoskov  _et al._ (2015) to calculate the photon
emission via nonlinear Compton scattering and the electron-positron pair
production via nonlinear Breit-Wheeler process.

The spin- and polarized-averaged differential probability of the photon
emission is written as Baier  _et al._ (1998)

$\displaystyle\frac{d^{2}W_{\rm rad}}{dudt}=\frac{\alpha
m^{2}c^{4}}{\sqrt{3}\pi\hbar\varepsilon_{e}}\left[\frac{u^{2}-2u+2}{1-u}K_{2/3}(y)-\int_{y}^{\infty}K_{1/3}(x)dx\right],$
(23)

where $K_{\nu}(y)$ is the second-kind $\nu$-order modified Bessel function,
$y=2u/[3(1-u)\chi_{e}]$, $u=\varepsilon_{\gamma}/\varepsilon_{e}$,
$\varepsilon_{e}$ is the electron energy before the photon emission,
$\varepsilon_{\gamma}$ is the emitted photon energy, and $\alpha\approx 1/137$
is the fine structure constant. Quantum parameter
$\chi_{e}=(e\hbar/m_{e}^{3}c^{4})|F_{\mu\nu}p^{\nu}|$ presents the field
experienced by the electron in its rest frame normalized to the Schwinger
critical field $E_{\rm cr}=1.3\times 10^{16}$ V/cm or $B_{\rm cr}=4.4\times
10^{13}$ G. The laser-plasma interaction enters the QED-dominated region if
$\chi_{e}\gtrsim 1$.

Figure 8:  (Upper row) The spectra of emitted photons with
$\varepsilon_{e}=510$ MeV and (a) $\chi_{e}=0.49$, (b) $\chi_{e}=4.9$, and (c)
$\chi_{e}=49$, respectively. (Lower row) The spectra of generated positrons
with $\varepsilon_{\gamma}=510$ MeV and (d) $\chi_{\gamma}=0.49$, (e)
$\chi_{\gamma}=4.9$, and (f) $\chi_{\gamma}=49$, respectively.

To calculate Eq. (23), two different Monte-Carlo methods have been implemented
in yunic. Both methods require two uniformly distributed random numbers
$r_{1}$ and $r_{2}$ to simulate the stochastic photon emission, where
$r_{1},r_{2}\in(0,1)$. Method I Elkina  _et al._ (2011): First, a random
number $r_{1}$ is generated to compare with the total radiation probability
$W_{\rm rad}$; if $r_{1}>W_{\rm rad}$, a photon is emitted, and its energy
ratio $u_{0}$ is determined by the other random number $r_{2}$ according to
$\int_{u_{\rm min}}^{u_{0}}dW_{\rm rad}/du=r_{2}W_{\rm rad}$ with a low-energy
cutoff $u_{\rm min}$. Method II Elkina  _et al._ (2011); Gonoskov  _et al._
(2015): a photon with an energy ratio of $u_{0}=r_{1}^{3}$ is emitted if
$3r_{1}^{2}dW_{\rm rad}(u_{0})/du>r_{2}$. The photon spectrum calculated by
two methods are both in good agreement with the theoretical spectrum of Eq.
(23), as shown in Figs. 8(a)-(c).

Similarly, the spin- and polarized-averaged differential probability of the
pair production is written as Baier  _et al._ (1998)

$\displaystyle\frac{d^{2}W_{\rm pairs}}{d\varepsilon_{+}dt}$ $\displaystyle=$
$\displaystyle\frac{\alpha
m^{2}c^{4}}{\sqrt{3}\pi\hbar\varepsilon_{\gamma}^{2}}\left[\frac{\varepsilon_{+}^{2}+\varepsilon_{-}^{2}}{\varepsilon_{+}\varepsilon_{-}}K_{2/3}(y)+\int_{y}^{\infty}K_{1/3}(x)dx\right],$
(24)

where
$y=2\varepsilon_{\gamma}^{2}/(3\chi_{\gamma}\varepsilon_{+}\varepsilon_{-})$,
and $\varepsilon_{\gamma}$, $\varepsilon_{-}$ and $\varepsilon_{+}$ are the
energies of the parent $\gamma$ photon, newly created electron and positron,
respectively. Another parameter
$\chi_{\gamma}=(e\hbar^{2}/m_{e}^{3}c^{4})|F_{\mu\nu}k^{\nu}|$ characterizes
the pair production. Two Monte-Carlo methods for the photon emission discussed
above are also employed to calculate the pair production in the similar way,
which are in good agreement with the theory, as shown in Figs. 8(d)-(f).

Figure 9:  Energy spectra of (a)(d) electrons $\Phi_{-}$, (b)(e) photons
$\Phi_{\gamma}$, and (c)(f) positrons $\Phi_{+}$. (a)-(c) and (d)-(f)
correspond to strengths of the external magnetic field $B_{0}=0.001B_{\rm cr}$
and $B_{0}=0.009B_{\rm cr}$, respectively.

Now, we benchmark our QED module in the process of electron-positron cascades
with just Method II, where the photon emission, quantum radiation reaction,
and pair production are self-consistently included. In Fig. 9, we consider the
same simulation setups as Fig. 2 and Fig. 4 of Ridgers  _et al._ (2014), where
an electron bunch of an initial Lorentz factor $\gamma_{0}=1000$ are moving
under a perpendicularly external magnetic field of a strength of
$B_{0}=0.001B_{\rm cr}$ or $B_{0}=0.009B_{\rm cr}$. The energy spectra of
electrons, photons and positrons both agree well with the results of Ridgers
_et al._ (2014).

In Fig. 10, we present a typical cascade case to further test our QED module.
An electron with an initial Lorentz factor $\gamma_{0}=2\times 10^{5}$ is
moving under a perpendicularly external magnetic field of a strength of
$B_{0}=0.2B_{\rm cr}$. The total number of electrons and positrons of energies
above 100 MeV are countered. The simulation result is shown in Fig. 10, which
is averaged over 2000 simulation runs with different random seeds. Our
simulation result is in good agreement with those in Anguelov and Vankov
(1999); Elkina  _et al._ (2011); Gonoskov  _et al._ (2015).

The electron/positron spin and $\gamma$-photon polarization are also
implemented into the QED module of yunic Song  _et al._ (2019, 2021a, 2021b)
based on spin- and polarization resolved photon emission and pair production
probabilities Baier  _et al._ (1998); Li  _et al._ (2019, 2020a, 2020b), which
have been benchmarked in Song  _et al._ (2021a). Hence, yunic can be used to
investigate spin and polarization related effects in laser-plasma
interactions, which are ignored in the previously employed QED-PIC.

Figure 10:  The number of pairs of energies above 200 MeV, where $t_{\rm
rad}=1.16\times 10^{-16}$ s.

##  V Conclusion

In summary, we have introduced a multi-dimensional PIC code, named yunic. Its
core algorithm is described and several benchmarks are preformed. Yunic was
employed to investigate the low-frequency whistler waves excited by
relativistic laser pulses Song  _et al._ (2020) and the spin and polarization
effects on the nonlinear Breit-Wheeler pair production in laser-plasma
interactions Song  _et al._ (2021b). In the future, yunic will continue to be
used to explore interesting phenomena or physical processes of application
prospects and new physical modules may be added for specific problems.

###### Acknowledgements.

This work was supported by the National Key R&D Program of China (Grant No.
2018YFA0404801), National Natural Science Foundation of China (Grant Nos.
11775302 and 11721091), the Strategic Priority Research Program of Chinese
Academy of Sciences (Grant Nos. XDA25050300, XDA25010300).

## References

  * Birdsall and Langdon (1991) Charles K Birdsall and A Bruce Langdon,  _Plasma Physics via Computer Simulation_ (Institute of Physics, Bristol, 1991).
  * Dawson (1983) John M. Dawson, “Particle simulation of plasmas,” Rev. Mod. Phys. 55, 403–447 (1983).
  * Gibbon (2005) Paul Gibbon,  _Short Pulse Laser Interactions with Matter_ (Imperial College Press, London, 2005).
  * Macchi (2013) A. Macchi,  _A Superintense Laser-Plasma Interaction Theory Primer_ (Springer, 2013).
  * Baier  _et al._ (1998) V. N. Baier, V.M. Katkov, and V. M. Strakhovenko,  _Electromagnetic Processes at High Energies in Oriented Single Crystals_ (World Scientific, Singapore, 1998).
  * Eastwood (1991) James W. Eastwood, “The virtual particle electromagnetic particle-mesh method,” Comput. Phys. Commun. 64, 252–266 (1991).
  * Villasenor and Buneman (1992) John Villasenor and Oscar Buneman, “Rigorous charge conservation for local electromagnetic field solvers,” Comput. Phys. Commun. 69, 306–316 (1992).
  * Taflove and Hagness (2000) A. Taflove and S. C. Hagness,  _Computational Electromagnetics: The Finite-Difference Time-Domain Method_ (Artech House, Boston, 2000).
  * Kane Yee (1966) Kane Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966).
  * Boris (1970) J. P. Boris, “Proceedings of 4th conference on numerical simulation of plasmas,” (Naval Research Laboratory, Washington D. C., 1970) p. 3–67.
  * Qin  _et al._ (2013) Hong Qin, Shuangxi Zhang, Jianyuan Xiao, Jian Liu, Yajuan Sun, and William M. Tang, “Why is Boris algorithm so good?” Phys. Plasmas 20, 084503 (2013).
  * Vay (2007) J.-L. Vay, “Noninvariance of space- and time-scale ranges under a Lorentz transformation and the implications for the study of relativistic interactions,” Phys. Rev. Lett. 98, 130405 (2007).
  * Vay (2008) J.-L. Vay, “Simulation of beams or plasmas crossing at relativistic velocity,” Phys. Plasmas 15, 056701 (2008).
  * Abe  _et al._ (1986) Hirotade Abe, Natsuhiko Sakairi, Ryohei Itatani, and Hideo Okuda, “High-order spline interpolations in the particle simulation,” J. Comput. Phys. 63, 247–267 (1986).
  * Esirkepov (2001) T. Zh. Esirkepov, “Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor,” Comput. Phys. Commun. 135, 144 – 153 (2001).
  * (16) Hartmut Ruhl,  _Classical Particle Simulations with the PSC code_.
  * Arber  _et al._ (2015) T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, “Contemporary particle-in-cell approach to laser-plasma modelling,” Plasma Phys. Control. Fusion 57, 113001 (2015).
  * Derouillat  _et al._ (2018) J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai, J. Dargent, C. Riconda, and M. Grech, “Smilei : A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation,” Comput. Phys. Commun. 222, 351–373 (2018).
  * Massimo  _et al._ (2016) F. Massimo, A. Marocchino, and A.R. Rossi, “Electromagnetic self-consistent field initialization and fluid advance techniques for hybrid-kinetic PWFA code Architect,” Nucl. Instrum. Methods A 829, 378–382 (2016), 2nd European Advanced Accelerator Concepts Workshop - EAAC 2015.
  * Elkina  _et al._ (2011) N. V. Elkina, A. M. Fedotov, I. Yu. Kostyukov, M. V. Legkov, N. B. Narozhny, E. N. Nerush, and H. Ruhl, “QED cascades induced by circularly polarized laser fields,” Phys. Rev. ST Accel. Beams 14, 054401 (2011).
  * Ridgers  _et al._ (2014) C. P. Ridgers, J. G. Kirk, R. Duclous, T. G. Blackburn, C. S. Brady, K. Bennett, T. D. Arber, and A. R. Bell, “Modelling gamma-ray photon emission and pair production in high-intensity laser–matter interactions,” J. Comput. Phys. 260, 273–285 (2014).
  * Gonoskov  _et al._ (2015) A. Gonoskov, S. Bastrakov, E. Efimenko, A. Ilderton, M. Marklund, I. Meyerov, A. Muraviev, A. Sergeev, I. Surmin, and E. Wallin, “Extended particle-in-cell schemes for physics in ultrastrong laser fields: Review and developments,” Phys. Rev. E 92, 023305 (2015).
  * Anguelov and Vankov (1999) V. Anguelov and H. Vankov, “Electromagnetic showers in a strong magnetic field,” J. Phys. G: Nucl. Part. Phys. 25, 1755–1764 (1999).
  * Song  _et al._ (2019) Huai-Hang Song, Wei-Min Wang, Jian-Xing Li, Yan-Fei Li, and Yu-Tong Li, “Spin-polarization effects of an ultrarelativistic electron beam in an ultraintense two-color laser pulse,” Phys. Rev. A 100, 033407 (2019).
  * Song  _et al._ (2021a) Huai-Hang Song, Wei-Min Wang, Yan-Fei Li, Bing-Jun Li, Yu-Tong Li, Zheng-Ming Sheng, Li-Ming Chen, and Jie Zhang, “Spin and polarization effects on the nonlinear Breit-Wheeler pair production in laser-plasma interaction,” (2021a), arXiv:2102.05882 [physics.plasm-ph] .
  * Song  _et al._ (2021b) Huai-Hang Song, Wei-Min Wang, and Yu-Tong Li, “Generation of polarized positron beams via collisions of ultrarelativistic electron beams,” (2021b), arXiv:2103.10417 [physics.acc-ph] .
  * Li  _et al._ (2019) Yan-Fei Li, Rashid Shaisultanov, Karen Z. Hatsagortsyan, Feng Wan, Christoph H. Keitel, and Jian-Xing Li, “Ultrarelativistic electron-beam polarization in single-shot interaction with an ultraintense laser pulse,” Phys. Rev. Lett. 122, 154801 (2019).
  * Li  _et al._ (2020a) Yan-Fei Li, Rashid Shaisultanov, Yue-Yue Chen, Feng Wan, Karen Z. Hatsagortsyan, Christoph H. Keitel, and Jian-Xing Li, “Polarized ultrashort brilliant multi-GeV $\gamma$ rays via single-shot laser-electron interaction,” Phys. Rev. Lett. 124, 014801 (2020a).
  * Li  _et al._ (2020b) Yan-Fei Li, Yue-Yue Chen, Wei-Min Wang, and Hua-Si Hu, “Production of highly polarized positron beams via helicity transfer from polarized electrons in a strong laser field,” Phys. Rev. Lett. 125, 044802 (2020b).
  * Song  _et al._ (2020) Huai-Hang Song, Wei-Min Wang, Jia-Qi Wang, Yu-Tong Li, and Jie Zhang, “Low-frequency whistler waves excited by relativistic laser pulses,” Phys. Rev. E 102, 053204 (2020).