# Unsupervised topological learning
for identification of atomic structures

Sébastien Becker1,2, Emilie Devijver2, Rémi Molinier3, and Noël Jakse1 1Univ.
Grenoble Alpes, CNRS, Grenoble INP, SIMaP, F-38000 Grenoble, France
2Univ. Grenoble Alpes, CNRS, Grenoble INP, LIG, F-38000 Grenoble, France
3Univ. Grenoble Alpes, CNRS, IF, F-38000 Grenoble,France

###### Abstract

We propose an unsupervised learning methodology with descriptors based on
Topological Data Analysis (TDA) concepts to describe the local structural
properties of materials at the atomic scale. Based only on atomic positions
and without a priori knowledge, our method allows for an autonomous
identification of clusters of atomic structures through a Gaussian mixture
model. We apply successfully this approach to the analysis of elemental Zr in
the crystalline and liquid states as well as homogeneous nucleation events
under deep undercooling conditions This opens the way to deeper and autonomous
study of complex phenomena in materials at the atomic scale.

atomic structures, molecular dynamics, homogeneous nucleation, unsupervised
learning, topological data analysis

##  I Introduction

Knowledge of microscopic events driving properties of matter is of fundamental
and technological interest. More precisely, understanding the structure at the
atomic scale is crucial for the design of new materials, which can play a role
in solving contemporary economical and social issues such as the reduction of
energy consumption or new drug design [1]. However, detailed information on
such a small spatial scale is still out of reach experimentally.

To access physical and chemical properties at the atomic scale, simulation
tools such as molecular dynamics (MD) represent dedicated in-depth means [2].
Thanks to the increase of computational power, such simulations can nowadays
easily reach millions of atoms evolving through several nanoseconds. Although
production of large amounts of valuable data has been streamlined, the
processing of the so-called “big data” is far from being obvious and has led
to the fourth paradigm [3] of science in many scientific fields in the last
two decades. Likewise, in condensed matter physics and materials science, new
protocols and tools based on supervised learning and neural networks [4, 5, 6,
7] as well as unsupervised learning [8, 9, 10, 11, 12], or even combinations
of these approaches [13] are currently built at an accelerated rate.

In this paper, we propose an unsupervised method to analyze structural
information, where descriptors from atomic positions are constructed using
persistent homology (PH), a classical TDA tool [14]. PH has been successfully
applied to MD simulations to analyze the medium-range structural environments
[15] in amorphous solids [16, 17, 18], ice [19], and complex molecular liquids
[20]. However, to the best of our knowledge, PH has never been used as a
descriptor to encode local atomic structures, whereas it provides topological
features at different resolutions and for different levels of homology.
Similar local atomic structures, with respect to those topological features,
are clustered using a Gaussian mixture model (GMM). To illustrate the
potential of this approach, we analyze here previous simulations of
homogeneous crystal nucleation of elemental Zirconium (Zr) [21].

Section II is dedicated to a presentation of our method illustrated on pure Zr
simulations, while Section III presents and discusses the physical results
highlighted on these simulations. Finally, Section IV draw the conclusions.

##  II Unsupervised learning based on topological descriptors

We present in this section the complete methodology behind our unsupervised
protocol, giving all the steps necessary for its implementation.

###  II.1 Data production

Our unsupervised learning approach is illustrated on pure Zr simulations. All
the simulations mentioned in this paper have been performed with the lammps
code [22], in the isobaric-isothermal ensemble with the Nosé-Hoover thermostat
and barostat [23] and interatomic interaction described by the MEAM potential
[24] developed in [21]. For a detailed description of the procedure of the
simulation, we refer the reader to Ref. [21], and only the essential
characteristics are recalled here. The phase space trajectory is obtained by
integrating Newton’s equations numerically using the Verlet’s algorithm in its
velocity form, with a time step of 2 fs. A simulation box with $N=1$ $024$
$000$ atoms is considered with periodic boundary conditions (PBC) in the three
directions of space. This system is first equilibrated in the liquid state
above the melting points at $T=2500$ K and quenched with a cooling rate of
$10^{12}$ K/s at ambient pressure down to the nose of the time-temperature-
transformation (TTT) curve at $T=1250$ K as shown in Fig. 1. The homogeneous
nucleation process is then observed along this isotherm.

When the onset of nucleation is observed, typically at the nucleation time, a
$N$-atom configuration is extracted and its inherent structure is used to
build a training set. The inherent structures [25] are obtained by minimizing
the energy by means of a conjugate gradient algorithm, implemented in lammps,
to bring the system to its local minimum in the potential energy surface to
suppress the thermal noise.

Figure 1: Time-temperature-transformation curve with 1 024 000 atoms in a
temperature range near the nose, which was obtained over $5$ independent runs
for each temperature . Inset: evolution of the potential energy as a function
of time along a nucleation process. See more details in [21].

###  II.2 Data preparation

Local atomic structures in the configuration are defined through a
coordination sphere, consisting of a central particle and its surrounding
environment up to a chosen cut-off radius. The cut-off radius is defined
thanks to the radial distribution function $g(r)$, which gives the probability
to find a particle at a distance $r$ to a particle at the origin and is
defined by

$g(r)=\frac{N}{V}\frac{n(r)}{4\pi r^{2}\Delta r},$ (1)

where $n(r)$ stands for the mean number of particles in a spherical shell of
radius $r$ and thickness $\Delta r$ centered on the origin. A typical
representation of the radial distribution function $g(r)$ for the liquid state
of Zr is depicted in Fig. 2(a). Each minimum of the radial distribution
function can reasonably be used as a cut-off radius to obtain local
environments that represent consecutive neighbor shells of each central
particle. The first minimum (beyond the first maximum) defines the end of the
first neighbor shell ($4.43$ Å from the central particle in Fig. 2(a)), the
second minimum defines the end of the second neighbor shell ($7.28$ Å from the
central particle) and so on.

Among the million particles in the configuration, local atomic structures are
subsampled to construct a training set for the learning. To be representative
of the configuration while satisfying statistical independence, the
subsampling should pave the whole simulation box with respect to the PBC. For
a given cut-off radius, the subsample is one among the (almost) maximal sets
such that the central particles of the local atomic structures are separated
by at least twice this cut-off radius. More precisely, such sets are computed
through an iterative process summarized on the flowchart Fig. 2(b). The choice
of the cut-off radius (between first shell, second shell, etc.) as well as the
number of structures in the training set will be justified in Section II.3.

Figure 2: (a) Schematic illustration of the first and second neighbor shells
of a central particle (in red) with cut-off radii determined by the radial
distribution function g(r) of undercooled Zr at 1250 K. (b) Flowchart
illustrating how to define a training set with independent structures. First,
all particle positions in the simulation box are randomly ordered and set in a
list. Then, the first particle of this list is set as a reference, and all
particles which are closer than twice the selected cut-off radius $r_{c}$ from
this reference particle are removed. Changing the reference particle, this
process is iteratively repeated until all central particles meet the distance
criterion based on $r_{c}$.

Once the central particles are identified, the Python package pyscal [26] is
used on the full configuration to efficiently extract the particles
coordinates of the local atomic structures.

###  II.3 Persistent homology as a local atomic structure descriptor

To encode topological information of the local atomic structures, we use PH, a
popular TDA tool [27, 28] that detects relevant topological features from a
point cloud.

Given a set of points $X$, one can construct a collection of simplicial
complexes [29], called Vietoris-Rips complexes as follows. Given $t\geq 0$, we
consider the simplicial complex $\mathcal{X}_{t}$ which is the union of
$k$-simplices for $k\in\mathbb{N}$, where 0-simplices are all the points of
$X$, 1-simplices are segments with extremities in $X$ of length smaller than
$t$, 2-simplices are full triangles with vertices in $X$ and distant of at
most $t$ from one another, etc. In the end, we get an increasing sequence
(called filtration) $\mathcal{X}=(\mathcal{X}_{t})_{t\geq 0}$ of simplicial
complexes. A schematic representation of this filtration process is depicted
in Fig. 3. Actually, since the number of points is finite, changes in these
complexes only appear at a finite sequence of steps $t_{1}<\cdots<t_{n}$ which
gives a finite filtration of simplicial complexes
$\mathcal{X}_{0}=\mathcal{X}_{t_{0}}\subset\mathcal{X}_{t_{1}}\subset\dots\subset\mathcal{X}_{t_{n}}$.
For each space $\mathcal{X}_{t}$ of the filtration, we compute the simplicial
homology, which is the sequence of their homology groups
$(H_{k}(\mathcal{X}_{t}))_{k\geq 0}$, further denoted $(H_{k})_{k\geq 0}$ when
the space is understood or not specified for general consideration. The
dimension of these homology groups gives insights on the dataset. For
instance, the dimension of $H_{0}$ gives the number of connected components,
the dimension of $H_{1}$ the number of holes, and the dimension of $H_{2}$ the
number of cavities inside the simplicial complex. Given a filtration of
simplicial complexes $(\mathcal{X}_{t})_{t\geq 0}$, one can keep track of
topological features (encoded by elements in the homology groups) and how
”persistent” they are, i.e., at which $t$ they appear (the birth) and at which
$t$ they disappear (the death). A persistence diagram (PD) summarizes this
information by plotting the pair (birth, death) of each topological feature of
the filtration on a graph, as illustrated in Fig. 3 or Fig. 4.

Figure 3: Schematic snapshots of the Vietoris-Rips complexes filtration
$(\mathcal{X}_{t})_{t\geq 0}$ applied to a local structure with two neighbor
shells. The _barcode_ below encodes the persistence of each topological
feature by a line starting at its birth (the radius where its appear) and
finishing at its death (the radius where it disappears). The colors correspond
to the homological dimension (red, blue and green for $H_{0}$, $H_{1}$, and
$H_{2}$ respectively) of each topological feature. The insert represents the
corresponding PD where each point corresponds to a line, and thus its
associated topological feature with coordinates (birth, death).

Using Python packages gudhi [30] and ripser.py [31], the individual PDs of the
local atomic structures from the previously extracted training set are
computed for each homological dimensions $H_{0}$, $H_{1}$, and $H_{2}$. To
remove topological noise when considering $H_{1}$ and $H_{2}$, we use a
subsampling approach as introduced in [32].

While PDs are usually compared using the Bottleneck distance, their space
equipped with this metric cannot be embedded into a Euclidean space nor even a
normed vector space [33]. To tackle this problem, several mapping into vector
spaces have been proposed in the literature. In this paper, we use a method
developed in [34] and classically used to study 3D shapes. Each coordinate of
the topological vector is associated to a pair of points $(x,y)$ in a PD $D$
for a fixed level of homology, except for the infinite point, and is
calculated by

$m_{D}(x,y)=\min\\{\|x-y\|_{\infty},d_{\Delta}(x),d_{\Delta}(y)\\},$ (2)

where $d_{\Delta}(\cdot)$ denotes the $\ell^{\infty}$ distance to the
diagonal, and those coordinates are sorted by decreasing order. Remark that
the dimension of each topological vector depends on the local atomic
structure, so we fill in each vector with a non-informative value to reach the
maximal dimension of the descriptor space. Here we decided to fill the vectors
with the value $-1$ instead of $0$ as proposed in [34], as distances between
points in PDs are always nonnegative, and a zero, i.e., two points of the PD
with the same birth and death, corresponds to relevant information in our
context.

To highlight the importance of these zeros in our topological vectors, let’s
consider Zr atoms on a body-centered cubic (bcc) lattice up to the second
neighbor shell ($7.16$ Å from the central particle) and Zr atoms on a face-
centered cubic (fcc) lattice in the same shell as shown in Fig. 4. Even though
these structures are fundamentally different, it can be noted that the only
thing that differs between the topological vectors obtained from each
structure is the dimension of the vector (64+120 [0] values for the bcc
lattice, 140+52+68 [0] values for the fcc lattice). Thus, if we fill the
smaller with zeros, we will obtain two identical vectors, i.e., the same point
in the descriptor space.

Figure 4: (a) Zr atoms on a bcc lattice and the corresponding PD. (b) Zr atoms
on a fcc lattice and the corresponding PD. The homological dimensions $H_{0}$,
$H_{1}$, and $H_{2}$ are depicted in red, blue, and green respectively. As
these structures are on periodic lattices, all topological features of the
filtration from a specific homological dimension have the same pair (birth,
death).

Finally, to illustrate the relevance of our topological signature to study the
local atomic structures, we confront it on our Zr simulation to widely used
classical physical descriptors such as the Bond Angle Analysis (BAA) [35], the
Common Neighbor Analysis (CNA) [36], and the Bond-Orientational Order Analysis
(BOOA) [37], in its averaged definition [3]. The latter computes the order
parameters $\bar{q}_{l}$ given by

$\bar{q}_{l}(i)=\sqrt{\frac{4\pi}{2l+1}\sum_{m=-l}^{l}\lvert\frac{1}{\tilde{N}_{b}(i)}\sum_{k=0}^{\tilde{N}_{b}(i)}q_{lm}(k)\rvert^{2}}$
(3)

with the complex vector

$q_{lm}(i)=\frac{1}{N_{b}(i)}\sum_{j=1}^{N_{b}(i)}Y_{lm}(\mathbf{r}_{ij}),$
(4)

in which $\tilde{N}_{b}(i)$ is accounting for all the neighbors of a central
particle $i$ plus $i$ itself, $N_{b}(i)$ for the sole number of nearest
neighbors and $Y_{lm}(\mathbf{r}_{ij})$ are the spherical harmonics functions
of the vector $\mathbf{r}_{ij}$ from particle $i$ to $j$, with $l$ an integer
ranging from 2 to 12 and $m\in[-l,+l]$. A direct comparison of the CNA and
BOOA, both being mostly used, can be found in Ref. [39]. These classical
descriptors characterize the radial and/or angular distribution of bonded
pairs of atoms in the coordination sphere. The averaged version of the BOOA is
preferred over the regular definition as a considerable improvement of the
accuracy by taking into account information from the neighboring particles
beyond the first neighbor shell. Fig. 5(a) gives an empirical correlation
matrix between the coordinates of our TDA descriptor, BAA, BOOA, and CNA. This
matrix is constructed on a sample of approximately $25$ $000$ local structures
extracted as described in Section II.2 with a cut-off radius corresponding to
the first minimum of the radial distribution function. The latter choice was
made as a descriptor, as the CNA cannot be applied beyond the first shell of
neighbors. The empirical correlation matrix shows that the topological
signature is highly correlated with all the other descriptors. Some
specificities of the data are highlighted: as all local neighborhood contain
at least $10$ particles, the first $10$ components of $H_{0}$ are highly
correlated; they are correlated with high values of $H_{1}$ too because most
of the local structures have few $H_{1}$ components ($50$ % of the population
has less than 5 components); $H_{2}$ components are not depicted here as there
is none in local atomic structures with only the first neighbor shell; the
correlation matrix associated to the CNA is sparse, because only a few bonds,
which are the first ones in the matrix (using the indexing of Faken and
Jónsson [40]: $[666]$ and $[444]$ for the bcc ordering and $[555]$, $[544]$
and $[433]$ for the icosahedral ordering) are present in most of the
structures. It also explains the low correlation between the last bonds in the
matrix of CNA with the other descriptors.

Figure 5: (a) Empirical correlation matrix (absolute value) between several
descriptor families: BAA, CNA, TDA, and BOOA. Low values are white, high
values are red. We distinguish between the descriptors by black lines, and
between levels of homology with a blue line ($H_{0}$ and $H_{1}$). (b) PDs of
a local atomic structure with one and two neighbor shells. The yellow line
corresponds to a threshold computed with the subsampling approach to remove
noises.

This example on the first neighbor shell illustrates how relevant signatures
from $H_{0}$ and $H_{1}$ are. To increase the number of $H_{0}$ and $H_{1}$
components and also to capture information on $H_{2}$ (which appears when
considering more than just one neighbor shell), we built a training set of $5$
$314$ local atomic structures with two neighbor shells to construct our model,
approaching the intermediate order. This leads here to $68$
$H_{0}$-components, $100$ $H_{1}$-components, and $22$ $H_{2}$-components,
instead of only $17$ $H_{0}$-components and $18$ $H_{1}$-components with one
neighbor shell. Fig. 5(b) shows a representative example of the differences in
the PDs of a structure with one neighbor shell and the same structure with two
neighbor shells. The second neighbor shell plays also a crucial role in the
structural description [41, 42].

Albeit considering higher-order shells is possible with our persistent
homological description of local structures, the advantage we gain in
topological information is balanced by a too high spatial resolution of the
local structures, leading to a loss of information in the GMM clustering.
Thus, we have restricted ourselves to the second neighbor shell, through a
trade-off between the structural information and a coarsening of the spatial
resolution, already observed for the averaged BOOA [3]. Note that the reason
for using the radial distribution function $g(r)$ to determine the second
shell of neighbors is to fulfill this compromise; adaptive methods which are
more general, such as SANN [43] or RAD [44] algorithms being limited to the
first shell of neighbors. It is worth mentioning that the clustering remains
essentially unchanged by varying the cut-off radius by $\pm 5$ % about the
second minimum, which is compatible with the volume change during the crystal
nucleation.

###  II.4 Gaussian mixture model clustering

A Gaussian Mixture Model (GMM) groups data points into clusters within an
unsupervised way through a mixture of $M$ Gaussian distributions
$(\phi(\;\cdot\;;\boldsymbol{\mu}_{m},\Sigma_{m}))_{1\leq m\leq M}$ of weights
$(\alpha_{m})_{1\leq m\leq M}$ as

$\sum_{m=1}^{M}\alpha_{m}\phi(\;\cdot\;;\boldsymbol{\mu}_{m},\Sigma_{m}),$ (5)

where $\boldsymbol{\mu}_{m}$ is the mean and $\Sigma_{m}$ the covariance
matrix of the $m$th Gaussian distribution. After standardizing the descriptor
space, an Expectation-Maximization (EM) algorithm [45] is used as an iterative
method to estimate the unknown parameters
$(\alpha_{m},\mu_{m},\Sigma_{m})_{1\leq m\leq M}$. We use the implementation
in the Python package scikit-learn [46] with full covariance matrices and $3$
$000$ $k$-means [47] initializations.

To select the number of Gaussian components, i.e. the number of clusters, the
integrated completed likelihood (ICL) [48] criterion has been used:

ICL
$\displaystyle=-2\ln(\hat{L})+D\ln(n)-2\sum_{k=1}^{n}\sum_{m=1}^{M}\hat{\tau}_{m,k}ln(\hat{\tau}_{m,k})=\text{BIC}-2\sum_{k=1}^{n}\sum_{m=1}^{M}\hat{\tau}_{m,k}ln(\hat{\tau}_{m,k}),$
(6)

where $\hat{L}$ denotes the likelihood evaluated on the estimator, $D$ the
number of parameters to be estimated ($D=M(d+d^{2}+1)$, with $d$ the dimension
of the topological vectors), and $\hat{\tau}_{m,k}$ denotes the probability to
belong to the $m$th component conditionally to the observation $x_{k}$. The
ICL generalizes the widely used Bayesian information criterion (BIC) [49] to
clustering methods by adding an entropic penalty computed from the posterior
probabilities of the data to be assigned to each Gaussian component. The
number of clusters in the model is set to the number achieving the minimum of
ICL.

##  III Application of the unsupervised approach on elemental zirconium
simulations

We focus in this section on the application of our method, called TDA-GMM, to
elemental Zr simulations. It demonstrates that this protocol provides relevant
structural information in a context as challenging as homogeneous nucleation.

###  III.1 Learning the model

As mentioned in Sec. II.2, configuration is chosen in course of the nucleation
process, where the supercooled liquid coexists with crystalline nuclei, to
capture all structural atomic events of interest. Fig. 6 depicts the ICL
computed on the extracted training set as described in Section II.3, leading
to select $7$ clusters [50]. Each cluster $C_{i}$ is represented by the local
structure that is the closest to the mean of the corresponding Gaussian
component, as illustrated in Fig. 7(a). Note that here, an embedded data
representation of the space of interest is not straightforward to interrogate
visually the width and shape of the different Gaussians [51], but an analysis
of the eigenvalues of the covariance matrices shows different elliptical
shapes. This proves the necessity of GMM over simpler unsupervised algorithms
such as $k$-means which would only be suitable to fit hyperspheres.

Figure 6: ICL criterion on the training set. Its minimum corresponds to the
maximum likelihood, which leads to the optimal number of clusters to build the
model.

Once the model is learned, it can be applied to new uncategorized local atomic
structures through their topological signatures. The robustness of our model
is checked by considering a test set consisting of all the structures in the
configuration but the one in the training set. Most of the structures in this
test set (99.998 %) have probabilities beyond 0.999 to be in one cluster of
the model, and the remaining ones still have probabilities higher than 0.5.
Further analysis of the fit of the model to the test set is performed using
the Mahalanobis distance (MaD) [52], which computes the distance between an
empirical Gaussian distribution and a dataset. Using hard-assignment to affect
each structure to a cluster, we compute the MaD for each cluster between the
test set and the multivariate Gaussian distribution. Outliers are then
considered based on the interquartile range rule. Table 1 presents the
percentage of structures in each cluster viewed as outliers according to this
criterion: more than 96 % of the structures in the test set have a MaD in
agreement with their assigned Gaussian distribution of the model.

Clusters | $C_{1}$ | $C_{2}$ | $C_{3}$ | $C_{4}$ | $C_{5}$ | $C_{6}$ | $C_{7}$ | Total
---|---|---|---|---|---|---|---|---
Proportion (%) | 11.70 | 7.22 | 9.28 | 24.45 | 33.46 | 11.03 | 2.86 | 100
Outliers (%) | 2.47 | 2.12 | 6.29 | 3.25 | 3.67 | 3.88 | 3.06 | 3.56
Table 1: Proportions of structures and percentage of outliers in each cluster
of the test set based on the interquartile range rule on the MaD between the
test set and their assigned multivariate Gaussian distribution of the model in
the descriptor space.

The final clustering is shown in Fig. 7(b) on the entire configuration
depicted in the real space of the simulation box with the help of ovito [53].
In the inset, the classification found by a classical adaptative CNA algorithm
as implemented in ovito is represented for comparison. In contrast to TDA-GMM,
only two types of structures that correspond to strictly crystalline or
icosahedral structures were successfully identified.

Figure 7: (a) Representation in the real space of the local structures
assigned to each cluster $C_{i}$. (b) Snapshot of the one million atoms in the
simulation box after the clustering in the descriptor space. In the inset, the
classification obtained from an adaptative CNA method is depicted.

###  III.2 Description of the bulk crystal and liquid

We have first applied our model to a configuration with $N=128$ $000$ atoms in
the bulk crystal and liquid at $T=1250$ K. Fig. 8 and Table 2 shows the
results obtained in both cases through the central atoms of the local
structures associated with each cluster.

Figure 8: Snapshots at $T=1250$ K of Zr bulk crystal (a) and liquid (b) brought in their inherent local structures. Atoms are colored according to the cluster they belong to (see Fig. 7). Bulk | Crystal | Liquid
---|---|---
$C_{1}$ (%) | 100 | 0.00
$C_{2}$ (%) | 0 | 0.03
$C_{3}$ (%) | 0 | 6.49
$C_{4}$ (%) | 0 | 29.83
$C_{5}$ (%) | 0 | 45.22
$C_{6}$ (%) | 0 | 14.70
$C_{7}$ (%) | 0 | 3.72
Table 2: Proportion of each cluster $C_{i}$ in the bulk crystal and liquid.

The structures of the bulk crystal are classified as local structures in
cluster $C_{1}$. Referring to Fig. 7(a), the atoms of the mean local structure
associated with $C_{1}$ are indeed on a periodic lattice. In the bulk liquid,
more than 99.96 % of the identified local structures belong to clusters
$C_{3},C_{4},C_{5},C_{6}$ and $C_{7}$. From this simple application, and
without any information about the physical nature of the atomic structures,
the trained TDA-GMM model distinguishes solid-like from liquid-like clusters.

###  III.3 Description of the homogeneous crystal nucleation

During the process of homogeneous nucleation along an isotherm, the internal
energy of the system undergoes a sharp drop. It follows the growth of the
nuclei that carry the crystalline periodicity, until the state of least energy
is reached. To autonomously detect these nuclei among the clusters obtained by
the TDA-GMM method, five configurations along an isotherm at the nose of the
TTT at $T=1250$ K with $N=1$ $024$ $000$ atoms are selected. They range from
the onset of nucleation at $160$ ps from the quench, up to a configuration in
the poly-crystalline bulk at $1560$ ps. The model built in Section III.1
corresponds to a configuration in the course of nucleation at $360$ ps. Table
3 shows the evolution of the proportion of each cluster in these
configurations.

Time (ps) | 160 | 200 | 240 | 280 | 310(M) | 1560(S)
---|---|---|---|---|---|---
$C_{1}$ (%) | 0.25 | 0.95 | 2.85 | 6.83 | 12.66 | 67.21
$C_{2}$ (%) | 0.67 | 1.37 | 2.76 | 4.96 | 6.85 | 22.65
$C_{3}$ (%) | 7.37 | 7.68 | 8.18 | 8.75 | 8.51 | 6.99
$C_{4}$ (%) | 29.99 | 29.47 | 28.41 | 26.52 | 24.22 | 1.58
$C_{5}$ (%) | 43.69 | 42.87 | 40.99 | 37.42 | 33.78 | 0.90
$C_{6}$ (%) | 14.44 | 14.14 | 13.48 | 12.36 | 11.16 | 0.35
$C_{7}$ (%) | 3.59 | 3.53 | 3.32 | 3.15 | 2.82 | 0.31
Table 3: Proportion of each cluster $C_{i}$ at different times of the
nucleation process. Superscripts (M) and (S) correspond respectively to the
configuration used to train the model and the solidified configuration.

Local structures from clusters $C_{1}$ and $C_{2}$ follow a fast growth until
they become the majority in the final bulk. On the opposite, the other
clusters undergo a more or less strong decrease in their proportions over
time. As mentioned in Section III.2, local structures from $C_{1}$ correspond
to a crystalline ordering whereas atoms of the mean local structure associated
with $C_{2}$ are on a partial periodic lattice. Fig. 9 represents the
snapshots of these clusters in these different configurations. The other
clusters, classified as liquid, show that apart from the nuclei, there exist
complex structural heterogeneities in the supercooled liquid, which we
describe in Section III.5, and which are successfully highlighted by the TDA-
GMM method. Such an assessment has recently been brought out as well in
supercooled liquid mixtures by another unsupervised method based on an
information-theoretic approach [12].

Figure 9: Snapshots of Zr’s configurations along nucleation at the nose of the
TTT at $T=1250$ K. The local structures in the clusters $C_{1}$ and $C_{2}$
identified as carrying the crystalline periodicity are depicted through their
central particle in red and pink, respectively.

The independent nuclei formed by local structures from $C_{1}$ and/or $C_{2}$
are extracted using a distance-based neighboring criterion as implemented in
ovito. This criterion is set as the first minimum of $g(r)$ in such a manner
that two particles are connected if they are at a distance less than or equal
to this threshold. The size distributions of these independent nuclei are then
obtained by accounting for all atoms of the overlapping local atomic
structures of two neighbor shells defined by each central particle belonging
to $C_{1}$ and/or $C_{2}$. The critical size is estimated by the smallest
cluster that persists over the several configurations of nucleation, with at
least the same number of atoms. In such a high undercooling regime, this leads
to a critical nucleus size of 70-80 atoms.

###  III.4 Translational and orientational orderings

Two physical orderings which have been shown to be important parameters
driving nucleation are the translational and orientational orderings [54, 55,
56].

With the nuclei previously identified, their centers of mass are easily
inferred. Based on these positions, the spatial evolution of the translational
ordering from the nuclei to the liquid can be quantified through density. The
density is indeed the dedicated indicator to measure the fluctuations in the
relative positions of neighboring atoms: for each cluster $C_{i}$ and its
relative number of atoms $K_{i}(r_{n})$ in a spherical shell of volume $V_{s}$
extending from the center of a nucleus to a distance $r_{n}$ to
$(r_{n}+\epsilon)$ with $\epsilon=1$ Å, the radial partial atomic density
profiles $\rho_{i}(r_{n})$ are computed as:

$\displaystyle\rho_{i}(r_{n})$
$\displaystyle=\frac{K_{i}(r_{n})}{V_{s}},\hskip 28.45274pt\text{with
}V_{s}=\frac{4}{3}\pi[(r_{n}+\epsilon)^{3}-r_{n}^{3}].$ (7)

On the other hand, the analysis of the orientational ordering provides
information on the fluctuations in relative geometric bonding between
neighboring atoms. Such a characterization can be performed with help of the
classical methods mentioned in Section II.3. For instance, the BOOA has shown
to be successful to analyze the crystallization of colloidal suspensions [57,
58]. Here, following the indexing of Faken and Jónsson, the CNA was performed
in each cluster as a representative measure [59]. It was verified, as
indicated by the high correlations in Fig. 5(a), that these two methods give
similar information about the orientational ordering [39].

This both orderings analysis shows us a general behavior, even for the
precritical nuclei, which is depicted in Fig. 10. Fig. 10(a) reveals that all
the nuclei emerge with a translational ordering which corresponds to the
density of the bulk crystal. Despite the small density window between the bulk
crystal and liquid at $T=1250$ K, it is clearly spotted that the sum of the
respective densities $\rho_{i}$ of the clusters $C_{i}$ associated with the
nuclei or liquid leads to the density values of the bulk crystal and liquid,
respectively. Fig. 10(b) confirms the crystalline structures in $C_{1}$ and
$C_{2}$ which exhibit bcc bonds ($[666]$: 8, $[444]$: 6) and thus a concurrent
emergence of the orientational ordering in the nucleation process. The slight
percentage of icosahedral bonds $[555]$ in $C_{2}$ tends to explain the
partial periodicity observed in Fig. 7(a). On the other hand, the local
structures from clusters $C_{4}$ to $C_{7}$ share strong five-fold symmetries
bonds characteristics of the undercooled liquid, although they also possess
non-negligible relative bcc bonds. About $4$ $\%$ of $[666]$ in these clusters
are in excess to form bcc structures and have the potential to combine with
$[555]$ bonds to form Z14, Z15, and Z16 Frank-Kasper polyhedra [60]. This can
lead to a geometrical frustration able to slow down the onset of the
nucleation in the liquid part [61, 62]. At last, $C_{3}$ is a more peculiar
case, with a strong distribution at the border of nuclei while also being
present in the deeper liquid. Its presence in the frontier region is explained
by the spatial resolution already mentioned in Section II.3, while its
presence in the liquid can be seen as a precursor for the emergence of the
embryos. This is in agreement with the latest view of a heterogeneous scenario
in which precursors of the crystalline ordering arise from structural
heterogeneities in the liquid [63, 64, 65] and which we develop in Section
III.5.

Figure 10: (a) Representation on a nucleus of the typical translational
ordering of the nuclei and their surrounding liquid through the radial density
profile of the $C_{i}$ clusters. The corresponding densities of the bulk
crystal and undercooled liquid have been determined at $T=1250$ K and constant
ambient pressure. This behavior of the translational ordering is
representative of all nuclei in the configurations, even the precritical ones.
(b) Typical orientational ordering of the $C_{i}$ clusters in the
configurations through CNA. $[555]$, $[544]$ and $[433]$ stand for five-fold
symmetries bonds and, $[666]$ and $[444]$ for bonds relative to a bcc crystal
ordering.

There is no consensus in the literature for a general pathway of the
nucleation process of materials, which seems to be system dependent: hard
spheres systems undergo concurrent translational and orientational ordering
[56], whereas a decoupling of these two symmetries has been observed in
colloidal systems [54]. For a pure metal like Zr, our unsupervised analysis
brings us to the conclusion that the translational and orientational orderings
arise synchronously.

###  III.5 Liquid heterogeneities

To assess the scenario mentioned in Section III.4, heterogeneity at an early
nucleation stage is evaluated by a comparison of the distribution of atoms
belonging to each cluster against the uniform distribution. As the comparison
of distributions is particularly complex in three dimensions since the
empirical cumulative distribution function is not defined, a Kolmogorov-
Smirnov (KS) test [66] was performed on the projection of atomic positions in
the three directions of space. For a level $0.01$, if the uniform distribution
is rejected on at least one projection, it is rejected for the entire
multivariate dataset. Table 4 depicts the results through the $p$-values of
the test obtained from the configuration at $160$ ps. For all the clusters the
test is always rejected, meaning that their distributions are not uniform, in
other words, indicating structural heterogeneities in the simulation box.

$p$-value | $x$ | $y$ | $z$
---|---|---|---
$C_{1}$ | 0.000 | 0.000 | 0.000
$C_{2}$ | 0.000 | 0.000 | 0.000
$C_{3}$ | 0.000 | 0.000 | 0.000
$C_{4}$ | 0.000 | 0.000 | 0.007
$C_{5}$ | 0.001 | 0.000 | 0.000
$C_{6}$ | 0.007 | 0.331 | 0.001
$C_{7}$ | 0.192 | 0.156 | 0.005
Table 4: $p$-values computed on the projection of atomic positions at an early
nucleation stage ($160$ ps) onto each direction of the simulation box from a
KS test against the uniform distribution.

##  IV Conclusions

We used PH, a main computational tool in TDA, to build descriptors of local
atomic structures. A GMM allowed us to autonomously identify clusters of
similar local structures in the space of these topological descriptors. We
applied this unsupervised approach to the analysis of several MD
configurations of Zr in the challenging context of homogeneous nucleation. We
highlighted some interesting properties of the phenomenon such as a concurrent
emergence of the translational and orientational orderings and structural
heterogeneities in the undercooled liquid. A closer look at the nucleation
process of other pure metals has been investigated through this TDA-GMM method
and the results will be published in future work. An extension of the analysis
on multicomponent systems where chemical orderings arise would be especially
relevant for further inspection and control of such mechanisms to enhance
materials design. More generally, this unsupervised learning methodology opens
the door to further study of structure-dependent phenomena at the atomic
scale.

## Acknowledgement

We acknowledge the CINES and IDRIS under Project No. INP2227/72914, as well as
CIMENT/GRICAD for computational resources. This work was performed within the
framework of the Centre of Excellence of Multifunctional Architectured
Materials “CEMAM”ANR-10-LABX-44-01 funded by the “Investments for the Future”
Program. This work has been partially supported by MIAI@Grenoble Alpes
(ANR-19-P3IA-0003). Fruitful discussions within the French collaborative
network in high-temperature thermodynamics GDR CNRS3584 (TherMatHT) are also
acknowledged.

## References

  * [1] G. C. Sosso, J. Chen, S. J. Cox, M. Fitzner, P. Pedevilla, A. Zen, and A. Michaelides, Chem. Rev. 116, 7078 (2016).
  * [2] D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications. 2nd Ed, Vol. 50 (1996).
  * [3] T. Hey, S. Tansley, and K. Tolle, The Fourth Paradigm: Data-Intensive Scientific Discovery (Microsoft Research, 2009).
  * [4] P. Geiger and C. Dellago, The Journal of Chemical Physics 139, 164105 (2013).
  * [5] C. Dietz, T. Kretz, and M. H. Thoma, Phys. Rev. E 96, 011301 (2017).
  * [6] E. Boattini, M. Ram, F. Smallenburg, and L. Filion, Molecular Physics 116, 3066 (2018).
  * [7] R. S. DeFever, C. Targonski, S. W. Hall, M. C. Smith, and S. Sarupria, Chem. Sci. 10, 7503 (2019).
  * [8] W. F. Reinhart, A. W. Long, M. P. Howard, A. L. Ferguson, and A. Z. Panagiotopoulos, Soft Matter 13, 4733 (2017).
  * [9] W. F. Reinhart and A. Z. Panagiotopoulos, Soft Matter 14, 6083 (2018).
  * [10] M. Ceriotti, J. Chem. Phys. 150, 150901 (2019).
  * [11] E. Boattini, M. Dijkstra, and L. Filion, J. Chem. Phys. 151, 154901 (2019).
  * [12] J. Paret, R. L. Jack, and D. Coslovich, J. Chem. Phys. 152, 144502 (2020).
  * [13] M. Spellings and S. C. Glotzer, AIChE Journal 64, 2198 (2018).
  * [14] F. C. Motta, in Advances in Nonlinear Geosciences, edited by A. A. Tsonis (Springer International Publishing, Cham, 2018), pp. 369–391.
  * [15] H. Tanaka, H. Tong, R. Shi, and J. Russo, Nat Rev Phys 1, 333 (2019).
  * [16] T. Nakamura, Y. Hiraoka, A. Hirata, E. G. Escolar, and Y. Nishiura, Nanotechnology 26, 304001 (2015).
  * [17] Y. Hiraoka, T. Nakamura, A. Hirata, E. G. Escolar, K. Matsue, and Y. Nishiura, Proc Natl Acad Sci USA 113, 7035 (2016).
  * [18] A. Hirata, T. Wada, I. Obayashi, and Y. Hiraoka, Commun Mater 1, 98 (2020).
  * [19] S. Hong and D. Kim, J. Phys.: Condens. Matter 31, 455403 (2019).
  * [20] K. Sasaki, R. Okajima, and T. Yamashita, Liquid Structures Characterized by a Combination of the Persistent Homology Analysis and Molecular Dynamics Simulation (Thessaloniki, Greece, 2018), p. 020015.
  * [21] S. Becker, E. Devijver, R. Molinier, and N. Jakse, Phys. Rev. B 102, 104205 (2020).
  * [22] S. J. Plimpton, J. Comput. Phys.117, 1 (1995).
  * [23] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids: Second Edition, Oxford Science Publication (Oxford University Press, Oxford, 2017).
  * [24] M. I. Baskes, Phys. Rev. B 46, 2727 (1992).
  * [25] F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 (1982).
  * [26] S. Menon, G. Leines, and J. Rogal, JOSS 4, 1824 (2019).
  * [27] H. Edelsbrunner, D. Letscher, and A. Zomorodian, Discrete Comput Geom 28, 511 (2002).
  * [28] A. Zomorodian and G. Carlsson, Discrete Comput Geom 33, 249 (2005).
  * [29] R. Ghrist, Bull. Amer. Math. Soc. 45, 61 (2007).
  * [30] C. Maria, J.-D. Boissonnat, M. Glisse, M. Yvinec. ”The Gudhi Library: Simplicial Complexes and Persistent Homology”, Mathematical Software, ICMS, 167-174, 2014.
  * [31] C. Tralie, N. Saul, and R. Bar-On, JOSS 3, 925 (2018).
  * [32] B. T. Fasy, F. Lecci, A. Rinaldo, L. Wasserman, S. Balakrishnan, and A. Singh, Ann. Statist. 42, (2014).
  * [33] M. Carrière and U. Bauer, 35th International Symposium on Computational Geometry, 21, (2019).
  * [34] M. Carrière, S. Y. Oudot, and M. Ovsjanikov, Computer Graphics Forum 34, 1 (2015).
  * [35] G. J. Ackland and A. P. Jones, Phys. Rev. B 73, 054104 (2006).
  * [36] J. D. Honeycutt and H. C. Andersen, J. Phys. Chem. 91, 4950 (1987).
  * [37] P. J. Steinhardt, D. R. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983).
  * [38] W. Lechner and C. Dellago, The Journal of Chemical Physics 129, 114707 (2008).
  * [39] See Supplemental Material at [URL will be inserted by publisher] for a comparison of the CNA and the BOOA.
  * [40] D. Faken and H. Jónsson, Computational Materials Science 2, 279 (1994).
  * [41] A. K. Soper and M. A. Ricci, Phys. Rev. Lett. 84, 2881 (2000).
  * [42] Z. Yan, S. V. Buldyrev, P. Kumar, N. Giovambattista, P. G. Debenedetti, and H. E. Stanley, Phys. Rev. E 76, 051201 (2007).
  * [43] J. A. van Meel, L. Filion, C. Valeriani, and D. Frenkel, J. Chem. Phys. 136, 234107 (2012).
  * [44] J. Higham, and R. H. Henchman, J. Chem. Phys. 145, 084108 (2016).
  * [45] A. P. Dempster, N. M. Laird, and D. B. Rubin, Journal of the Royal Statistical Society: Series B (Methodological) 39, 1 (1977).
  * [46] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, and D. Cournapeau, JMLR 12, 2825 (2011).
  * [47] S. Lloyd, IEEE Trans. Inform. Theory 28, 129 (1982).
  * [48] C. Biernacki, G. Celeux, and G. Govaert, IEEE Trans. Pattern Anal. Machine Intell. 22, 719 (2000).
  * [49] G. Schwarz, The Annals of Statistics 6, 461 (1978).
  * [50] See Supplemental Material at [URL will be inserted by publisher] for boxplots of the topological vectors in the homological dimensions $H_{0}$, $H_{1}$, and $H_{2}$ for each cluster $C_{i}$.
  * [51] See Supplemental Material at [URL will be inserted by publisher] for a principal component analysis on the first two axes and a t-SNE representation of the trained model.
  * [52] P. C. Mahalanobis, Proceedings of the National Institute of Sciences (Calcutta) 2, 49 (1936).
  * [53] A. Stukowski, Modelling Simul. Mater. Sci. Eng. 18, 015012 (2010).
  * [54] J. Russo and H. Tanaka, Sci Rep 2, 505 (2012).
  * [55] J. Russo and H. Tanaka, J. Chem. Phys. 18 (2016).
  * [56] J. T. Berryman, M. Anwar, S. Dorosz, and T. Schilling, The Journal of Chemical Physics 145, 211901 (2016).
  * [57] U. Gasser, E. R. Weeks, A. Schofield, P. N. Pusey and D. A. Weitz, Science 292, 258 (2001).
  * [58] U. Gasser, A. Schofield and D. A. Weitz, J. Phys.: Condens. Matter 15, 375 (2003).
  * [59] N. Jakse and A. Pasturel, Phys. Rev. Lett. 91, 195501 (2003).
  * [60] D. R. Nelson and F. Spaepen, Solid State Physics, Vol. 42 (Elsevier, 1989), pp. 1–90.
  * [61] H. Tanaka, Eur. Phys. J. E 35, 113 (2012).
  * [62] J. Russo, F. Romano, and H. Tanaka, Phys. Rev. X 8, 021040 (2018).
  * [63] A. Pasturel and N. Jakse, Npj Comput Mater 3, 33 (2017).
  * [64] Y. Shibuta, S. Sakane, E. Miyoshi, S. Okita, T. Takaki, and M. Ohno, Nat Commun 8, 10 (2017).
  * [65] G. Jug, A. Loidl, and H. Tanaka, EPL 133, 56002 (2021).
  * [66] A. Kolmogoroff, Giornale dell’Istituto Italiano degli Attuari 4, (1933).

Supplemental Material

## Boxplots of the descriptors for each cluster obtained by the TDA-GMM method
on elemental Zr

We show in the Figure S11 to Figure S31 below boxplots of the topological
vectors in the homological dimensions $H_{0}$, $H_{1}$, and $H_{2}$ for each
cluster $C_{i}$.

### Cluster 1

Figure S11: Boxplots of the topological vectors in homological dimension
$H_{0}$ for data in $C_{1}$. Figure S12: Boxplots of the topological vectors
in homological dimension $H_{1}$ for data in $C_{1}$. Figure S13: Boxplots of
the topological vectors in homological dimension $H_{2}$ for data in $C_{1}$.

### Cluster 2

Figure S14: Boxplots of the topological vectors in homological dimension
$H_{0}$ for data in $C_{2}$. Figure S15: Boxplots of the topological vectors
in homological dimension $H_{1}$ for data in $C_{2}$. Figure S16: Boxplots of
the topological vectors in homological dimension $H_{2}$ for data in $C_{2}$.

### Cluster 3

Figure S17: Boxplots of the topological vectors in homological dimension
$H_{0}$ for data in $C_{3}$. Figure S18: Boxplots of the topological vectors
in homological dimension $H_{1}$ for data in $C_{3}$. Figure S19: Boxplots of
the topological vectors in homological dimension $H_{2}$ for data in $C_{3}$.

### Cluster 4

Figure S20: Boxplots of the topological vectors in homological dimension
$H_{0}$ for data in $C_{4}$. Figure S21: Boxplots of the topological vectors
in homological dimension $H_{1}$ for data in $C_{4}$. Figure S22: Boxplots of
the topological vectors in homological dimension $H_{2}$ for data in $C_{4}$.

### Cluster 5

Figure S23: Boxplots of the topological vectors in homological dimension
$H_{0}$ for data in $C_{5}$. Figure S24: Boxplots of the topological vectors
in homological dimension $H_{1}$ for data in $C_{5}$. Figure S25: Boxplots of
the topological vectors in homological dimension $H_{2}$ for data in $C_{5}$.

### Cluster 6

Figure S26: Boxplots of the topological vectors in homological dimension
$H_{0}$ for data in $C_{6}$. Figure S27: Boxplots of the topological vectors
in homological dimension $H_{1}$ for data in $C_{6}$. Figure S28: Boxplots of
the topological vectors in homological dimension $H_{2}$ for data in $C_{6}$.

### Cluster 7

Figure S29: Boxplots of the topological vectors in homological dimension
$H_{0}$ for data in $C_{7}$. Figure S30: Boxplots of the topological vectors
in homological dimension $H_{1}$ for data in $C_{7}$. Figure S31: Boxplots of
the topological vectors in homological dimension $H_{2}$ for data in $C_{7}$.

## Embedded data representations

Figure S32: (a) PCA and (b) t-SNE first components of the trained model in the
descriptor space.

The most interesting in-depth that we were able to get to interrogate visually
the data was obtained by a classical Principal Component Analysis (PCA) [1]
represented on Figure S32(a). The dissociation between solid and liquid-like
structure is obvious on the first axis. However, there is no clear frontier
between all the clusters due to the dimensionality reduction into such a plane
(whereas our clusters are well separated, as given by the posterior
probabilities). A similar graph is obtained on Figure S32(b) by a t-SNE
algorithm [2] until convergence (around 20000 iterations) with a perplexity
value of 75. The representation onto two axis is not perfect (66 % of
information is covered by PCA) but we can see that the width is different for
each cluster ($C_{1}$ and $C_{2}$ being more diverse in the first axis) and
also different volumes ($C_{7}$ seems to have a high variance for few points).

## Bond-Orientational Order Analysis (BOOA)

Among the set of parameters based on spherical harmonics which can be computed
from the Bond-Orientational Order Analysis (BOOA), some of them are widely
used in the literature such as the $\bar{q}_{2}$-$\bar{q}_{6}$-plane for
identification of cristallinity, the $\bar{q}_{4}$-$\bar{q}_{6}$-plane for
distinction of crystal structure or the $\bar{q}_{6}$ distribution for
identification of solidity [3]. Figure S33 shows these representations on the
trained model. The distinction between the structures associated with the
liquid and solid clusters is highlighted. Figure S33(c) is particularly
informative on the fact that the structures in $C_{1}$ present a well-defined
crystalline order, while the distribution in $C_{2}$ is more spread out.
$C_{3}$ shows a boundary distribution between the clusters classified as
liquid and the cluster $C_{2}$ mainly located at the border of the nuclei.
These conclusions are in agreement and equivalent to our results obtained with
the CNA. However, although the liquid clusters are identifiable with BOOA, it
is not clear that the associated structures also possess a relative
crystalline order. On the contrary, the CNA clearly reveals the existence of
$[666]$ and $[444]$ bcc geometrical bonds in the structures of these liquid
clusters. This feature of the CNA allowed us to investigate the heterogeneous
scenario of nucleation as described in the manuscript.

Figure S33: (a) $\bar{q}_{2}$-$\bar{q}_{6}$-plane, (b)
$\bar{q}_{4}$-$\bar{q}_{6}$-plane and (c) probability distributions of
$\bar{q}_{6}$, for the structures in the trained model.

## References

  * [1] I. T. Jolliffe, Principal Component Analysis, Second Edition, Springer Series in Statistics (New York, 2002).
  * [2] L. van der Maaten and G. Hinton, J Mach Learn Res 9, 2579 (2008).
  * [3] W. Lechner and C. Dellago, The Journal of Chemical Physics 129, 114707 (2008).


# The Effect of the Parametric Decay Instability on the Morphology of Coronal
Type III Radio Bursts

Chaitanya Prasad Sishtla Department of Physics, University of Helsinki,
Helsinki, Finland  Immanuel Christopher Jebaraj Department of Physics &
Astronomy, University of Turku, Turku, Finland  Jens Pomoell Department of
Physics, University of Helsinki, Helsinki, Finland  Norbert Magyar Centre for
mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven,
Leuven, Belgium  Marc Pulupa Space Sciences Laboratory, University of
California, Berkeley, CA 94720-7450, USA  Emilia Kilpua Department of Physics,
University of Helsinki, Helsinki, Finland  Stuart D. Bale Space Sciences
Laboratory, University of California, Berkeley, CA 94720-7450, USA  Physics
Department, University of California, Berkeley, CA 94720-7300, USA

###### Abstract

The nonlinear evolution of Alfvén waves in the solar corona leads to the
generation of Alfvénic turbulence. This description of the Alfvén waves
involves parametric instabilities where the parent wave decays into slow mode
waves giving rise to density fluctuations. These density fluctuations, in
turn, play a crucial role in the modulation of the dynamic spectrum of type
III radio bursts, which are observed at the fundamental of local plasma
frequency and are sensitive to the local density. During observations of such
radio bursts, fine structures are detected across different temporal ranges.
In this study, we examine density fluctuations generated through the
parametric decay instability (PDI) of Alfvén waves as a mechanism to generate
striations in the dynamic spectrum of type III radio bursts using
magnetohydrodynamic simulations of the solar corona. An Alfvén wave is
injected into the quiet solar wind by perturbing the transverse magnetic field
and velocity components which subsequently undergo the PDI instability. The
type III burst is modelled as a fast-moving radiation source that samples the
background solar wind as it propagates to emit radio waves. We find the
simulated dynamic spectrum to contain striations directly affected by the
multi-scale density fluctuations in the wind.

Alfvén waves, PDI, Type III radio bursts, magnetohydrodynamics (MHD)

##  1 Introduction

The expanding solar wind contains fluctuations in the velocity and magnetic
field that are primarily transverse to the mean field and have characteristics
of Kolmogorov (1991) turbulence. In-situ studies in the heliosphere show that
most of the wave power in solar wind fluctuations result from Alfvénic waves,
while non-Alfvénic waves have a minor contribution (Higdon, 1984; Bruno &
Carbone, 2013; Chen, 2016). This allows modelling the ubiquitous solar wind
fluctuations as Alfvén waves (Belcher & Davis, 1971; D’Amicis & Bruno, 2015).
Subsequently, the heating and acceleration of the solar wind is explained via
the cascade of large-wavelength Alfvén waves leading to Alfvénic turbulence
(Coleman, 1968; Goldreich & Sridhar, 1995). The Alfvénic fluctuations are
generated, e.g., by the photospheric convective motions where the solar
magnetic field footpoints are anchored (Cranmer & van Ballegooijen, 2005). As
these waves propagate outward from the Sun, they can decay into slow modes
through modulational instabilities (Lashmore-Davies, 1976; Sakai & Sonnerup,
1983), and generate density variations through nonlinear wave-wave
interactions between fast, slow, and Alfvén waves (Nakariakov et al., 2000;
Chandran, 2005). In particular, Alfvén waves can couple with background
gradients in the background magnetic field to decay into a slow wave and a
reflected (sunward) Alfvén wave through the parametric decay instability (PDI)
(Galeev & Oraevskii, 1963; Sagdeev & Galeev, 1969; Derby, 1978; Hoshino &
Goldstein, 1989). The PDI instability has been observed both in the solar wind
(Bowen et al., 2018) and in the laboratory (Dorfman & Carter, 2016), and its
importance in generating counter-propagating waves leading to Alfvénic
turbulence (Usmanov et al., 2000; Suzuki & Inutsuka, 2006) is a topic of
continued interest (Chandran, 2018; Shoda et al., 2018a; Sishtla et al.,
2022). The relevance of the PDI for Alfvén wave dynamics is of particular
significance in the solar corona due to the low-beta plasma present in this
region. Hinode observations (De Pontieu et al., 2007) characterise the Alfvén
waves injected into the solar corona to have wave periods between $100-500$ s,
and amplitudes between $20-50$ km s-1. As the waves further propagate in the
solar corona, the low plasma beta environment (Iwai et al., 2014; Bourdin,
2017) increases substantially the decay rates of the pump Alfvén wave due to
the lower threshold of the PDI in a low beta plasma (Goldstein, 1978). In this
context, the PDI can contribute to coronal heating and the observed spectrum
of cross-helicity in the solar wind (Inhester, 1990; Del Zanna et al., 2001)
by generating counter-propagating Alfvén waves. The present study underscores
the significance of PDI by demonstrating that it has a direct effect on other
phenomena as well, in particular the formation of striations and other fine
structures within Type III solar radio bursts.

Figure 1: Example dynamic radio spectra produced from the measurements of
FIELDS/RFS onboard PSP during encounter 9 (22/11/2021). The probe was
approximately 17 $R_{\odot}$ during the time intervals at which these
measurements were made. A few examples of the “structured” type III radio
bursts are indicated in the top panel. The local plasma frequency at the
spacecraft ($f_{pe}$) is also marked in both spectra.

Type III radio bursts are signatures of suprathermal electron beams
propagating along open magnetic field lines (e.g., McLean & Labrum, 1985; Pick
& Vilmer, 2008), with typical beam speeds for $\approx 30$ keV electrons
between $0.1c-0.4c$ (where $c$ is the speed of light in vacuum) (Zhang et al.,
2018). Due to their near-relativistic velocities, electromagnetic emissions
from these beams can be used to probe the nature of the background plasma
through which the beams propagate. Type III radio bursts are the most commonly
observed and the most intense radio bursts of solar origin. Two example
intervals of the Radio Frequency Spectrometer (RFS; Pulupa et al., 2017) which
is part of the FIELDS instrument suite (Bale et al., 2016) onboard the Parker
Solar Probe (PSP; Fox et al., 2016) is shown in Figure 1. The measurements
were made during PSP’s 9th close encounter and show a number of type III radio
bursts. Type III bursts as the ones recorded in Figure 1, appear as fast
drifting bursts over a wide range of frequencies (e.g., Suzuki & Dulk, 1985).
The most widely accepted generation mechanism is coherent emission and was
first proposed by Ginzburg & Zhelezniakov (1958). This so-called “plasma
emission mechanism” is a two step process. The first step involves an electron
beam which becomes unstable due to time-of-flight effects and generates
Langmuir waves. In the second step, the Langmuir waves linearly or non-
linearly convert part of their energy into electromagnetic waves at radio
frequencies. As a natural consequence of the generation mechanism, emission is
predicted at the fundamental plasma frequency $f_{pe}$, and its harmonics
$nf_{pe}$ (where n = 2,3,… Robinson & Cairns, 1998a, b, c). Melrose (1980)
expanded upon the original idea of non-linear processes which included
scattering off ion-sound waves for the generation of fundamental radiation.
However, later models demonstrated that a linear conversion was preferable in
a randomly inhomogeneous plasma (Hinkel-Lipsker et al., 1992; Kim et al.,
2007; Voshchepynets et al., 2015; Voshchepynets & Krasnoselskikh, 2015;
Krasnoselskikh et al., 2019). Other competing mechanisms for the generation of
electromagnetic radiation through beam-driven Langmuir include, electrostatic
decay (Cairns & Melrose, 1985; Cairns & Fung, 1988; Robinson et al., 1994),
electromagnetic decay (Cairns, 1987; Robinson et al., 1993), quasi-mode
processes (Yoon et al., 1994), non-zero pitch angle electron beam driven
radiation (Tsiklauri, 2011; Schmitz & Tsiklauri, 2013), and antenna radiation
(Papadopoulos & Freund, 1978; Malaspina et al., 2010, 2012). Simulating the
local evolution of the beam-plasma system and the emission of electromagnetic
radiation self-consistently requires kinetic (Thurgood & Tsiklauri, 2015;
Krafft & Savoini, 2022) or other approaches used to treat Langmuir turbulence,
such as using the Zakharov-Hamiltonian set of equations (Krafft et al., 2013).
Taking these concerns into account, magnetohydrodynamic (MHD) simulations
cannot simulate such beam-plasma systems. Therefore, in this study we model
self-consistently the large scale inhomogenieties in the background plasma,
and use an emission proxy model to simulate the type III burst. Using this
approach we will demonstrate that inhomogeneities naturally give rise to the
fine-structures (striations) within fundamental emission of type III radio
bursts that resemble the fine-structure detected in the latest spacecraft
observation (Figure 1). Since it is the electrostatic non-linear interactions
which ultimately give rise to the harmonic emission, this process is not
modelled self-consistently here.

Previous remote sensing studies have suggested the structuring of type III
bursts by fast mode waves (Goddard et al., 2016; Kontar et al., 2017; Kolotkov
et al., 2018; Kaneda et al., 2018) which are formed by the dispersive
evolution of MHD waves propagating along a non-uniform plasma (Nakariakov et
al., 2004). However, studies made using high frequency and time resolution
radio observations have insisted that the fine structures are formed due to
coronal density fluctuations and are a multi-scale phenomena (Mugundhan et
al., 2017). The density fluctuations at large scales were reported to be
responsible for the interplanetary type III fine structures by Jebaraj et al.
(2023b) who also showed that the scale size of density fluctuations increase
with increasing distance from the Sun. Here, the scale sizes generating the
fine structures were prescribed to be 100-1000 Debye lengths, i.e., the length
corresponding to the relaxation of the electron beam.

Recent studies using the FIELDS/RFS measurements have made a number of
important discoveries regarding type III bursts (Pulupa et al., 2020; Jebaraj
et al., 2023a). Of particular importance is the distinct morphology of the
fundamental emission (highlighted as “structured” in Figure 1), which is
strongly structured and was previously seen as an outlier among hectometer
type III bursts (Jebaraj et al., 2023b). This finding, together with the
mechanism for their generation proposed in (Jebaraj et al., 2023b) would
indicate that structured emissions are a natural consequence when density
fluctuations are abundant. This leads to a qualitative assessment that density
fluctuations must be present in abundance at heights corresponding to hecto-
kilometer wavelengths.

In this study, we use MHD simulations along with an emission proxy to study
striations in the fundamental emission band of type III radio bursts as a
consequence of Alfvén wave propagation and evolution in the solar wind.
Previously, Kolotkov et al. (2018) modelled such striations by fitting a
modulated Newkirk density model (Newkirk, 1961) onto an observed type III
burst. They utilised a fitting approach to achieve agreement between the
observed intensity and modelled emission intensity variations to explain the
striations in type III bursts through density modulations. In this study, we
adopt a separate approach of modelling the radio burst striations through
self-consistent density fluctuations that are present in the solar wind. The
observed emission fluxes of such a radio burst might then similarly be a means
of probing the underlying density fluctuations inside a flux tube which is the
source of the observed type III burst. To perform this study, we first
initialise a solar wind with a continual injection of Alfvén waves which
undergo PDI to generate density variations. Then, a suprathermal radiation
source propagating in this quasi-steady wind is modeled to obtain the dynamic
spectrum. The study finds that striations in the type III burst are generated
by the density variations in the quasi-steady wind. The manuscript describes
the simulation setup in section 2, and the modelled type III burst in section
3.

##  2 Solar wind model with PDI-driven density fluctuations

The solar wind is modelled within a narrow, open magnetic flux-tube which is
centered on a radial magnetic field line. The geometry of the flux-tube is
parametrised by specifying the cross-sectional area ($a$) to be proportional
to the flux tube expansion ($f$) (Kopp & Holzer, 1976),

$a=a_{0}\left(\frac{r}{r_{0}}\right)^{2}f$ (1)

where $a_{0}$ is the cross-sectional area at the reference height $r_{0}$. By
flux conservation the magnetic field component along the flux-tube satisfies

$B_{r}=B_{r}(r_{0})\frac{a_{0}}{a}.$ (2)

The functional form providing the expansion factor $f$ is chosen as

$f=\frac{f_{\mathrm{max}}\exp\left(\left(r-R_{1}\right)/\sigma_{1}\right)+f_{1}}{\exp\left(\left(r-R_{1}\right)/\sigma_{1}\right)+1}.$
(3)

Here $f_{1}$ is chosen so that $f(r_{0})=1$, and $f_{\mathrm{max}}=3$,
$R_{1}=1.3~{}R_{\odot}$, and $\sigma_{1}=0.5~{}R_{\odot}$ are constants with
$R_{\odot}$ denoting the solar radius. Similar flux-tube geometries have been
used in previous studies to investigate solar wind properties in open magnetic
field regions (Chandran et al., 2011; Shoda et al., 2018b). In particular,
Sishtla et al. (2022) discusses the suppression of the PDI for lower Alfvén
wave frequencies by higher values of $f_{\mathrm{max}}$.

###  2.1 Forming the solar wind

The solar wind evolves in the flux-tube through the dynamical evolution of the
plasma by considering a one-dimensional MHD description including the relevant
physical processes of gravity and ad-hoc coronal heating,

$\frac{\partial}{\partial t}\rho+\frac{1}{a}\frac{\partial}{\partial
r}\left(a\rho v_{r}\right)=0$ (4) $\displaystyle\frac{\partial}{\partial
t}\left(\rho v_{r}\right)+\frac{1}{a}\frac{\partial}{\partial
r}\left[a\left(\rho
v_{r}^{2}+p+\frac{\mathbf{B^{2}_{\perp}}}{2\mu_{0}}\right)\right]=\left(p+\frac{\rho\mathbf{v^{2}_{\perp}}}{2}\right)\frac{1}{a}\frac{\partial}{\partial
r}a-\rho\frac{GM_{\odot}}{R_{\odot}^{2}}$ (5)
$\displaystyle\frac{\partial}{\partial
t}\left(\rho\mathbf{v_{\perp}}\right)+\frac{1}{a}\frac{\partial}{\partial
r}\left[a\left(\rho
v_{r}\mathbf{v_{\perp}}-\frac{B_{r}\mathbf{B_{\perp}}}{\mu_{0}}\right)\right]=-\frac{1}{2a}\frac{\partial
a}{\partial r}\left(\rho
v_{r}\mathbf{v_{\perp}}-\frac{B_{r}\mathbf{B_{\perp}}}{\mu_{0}}\right)$ (6)
$\frac{1}{a}\frac{\partial}{\partial r}\left(aB_{r}\right)=0$ (7)
$\displaystyle\frac{\partial}{\partial
t}\mathbf{B_{\perp}}+\frac{1}{a}\frac{\partial}{\partial
r}\left[a(\mathbf{B_{\perp}}v_{r}-B_{r}\mathbf{v_{\perp}})\right]=\frac{1}{2a}\frac{\partial
a}{\partial r}\left(\mathbf{B_{\perp}}v_{r}-B_{r}\mathbf{v_{\perp}}\right)$
(8) $\displaystyle\frac{\partial}{\partial
t}\mathcal{E}+\frac{1}{a}\frac{\partial}{\partial
r}\left[a\left(v_{r}\left\\{\mathcal{E}+p+\frac{\mathbf{B^{2}}}{2\mu_{0}}-\frac{B_{r}^{2}}{\mu_{0}}\right\\}-B_{r}\frac{\mathbf{B_{\perp}\cdot
v_{\perp}}}{\mu_{0}}\right)\right]=-\rho gv_{r}+S$ (9)

where,

$\mathcal{E}=\frac{1}{2}\rho v^{2}+\frac{P}{\gamma-1}+\frac{B^{2}}{2\mu_{0}}$
(10) $S=S_{0}\mathrm{exp}\left(-\frac{r}{L}\right).$ (11)

The quantities $\rho$, $\mathbf{v}$, $\mathbf{B}$, $\mathcal{E}$, and $p$
correspond to the mass density, plasma bulk velocity, magnetic field, total
energy density, and thermal pressure, respectively. The directions along and
transverse to the flux tube are denoted by $r$ and $\perp$, respectively. We
use an ideal-gas equation of state and set the adiabatic index to be
$\gamma=5/3$. To obtain a steady-state solar wind that approximates a Parker-
like outflow, we incorporate an additional energy source term in Equation 11
(Pomoell et al., 2015; Mikić et al., 2018) with $S_{0}=0.5\times 10^{-6}$ Wm-3
and $L=0.4~{}R_{\odot}$ m.

The MHD equations 4-9 are evolved in time using the strong stability
preserving (SSP) Runge-Kutta method to advance the semi discretised equations
in time. The simulation domain is discretized using $3000$ cells spaced
logarithmically from the solar corona (lower boundary). The method employs the
Harten–Lax–van Leer (HLL) approximate Riemann solver supplied by piece-wise,
linear slope limited interface states. The numerical method also employs the
constrained transport technique to ensure the magnetic field solution is
divergence free (Kissmann & Pomoell, 2012). These methods have been used
applied in previous studies of the solar corona (Pomoell & Vainio, 2012;
Sishtla et al., 2023). The simulation domain extends radially from
$r=1.03~{}R_{\odot}$ to $r=30~{}R_{\odot}$, with the $\theta,~{}\phi$
directions being transverse to the background magnetic field. At the low-
coronal inner boundary, the mass density $\rho$, temperature $T$, and radial
magnetic field $B_{r}$ are fixed to the values $\rho_{0}=8.5\times
10^{-13}~{}\mathrm{kg}~{}\mathrm{m}^{-3},~{}T_{0}=8\times
10^{5}~{}\mathrm{K},~{}B_{0}=2.5\,\mathrm{G}$.

###  2.2 Development of PDI in the solar wind

After a quasi-steady solar wind configuration is obtained by integrating the
MHD equations in time, we continually introduce monochromatic Alfvénic
perturbations at the lower boundary ($r=1.03~{}R_{\odot}$) by perturbing the
transverse components of the magnetic field. These Alfvénic fluctuations are
described using Elsässer variables (Elsasser, 1950), and denote correlations
between fluctuating velocity and magnetic field as
$\mathbf{z_{\theta,\phi}}^{\pm}=\mathbf{v_{\theta,\phi}}^{\pm}\pm{\mathbf{B_{\theta,\phi}}}/{\sqrt{\mu_{0}\rho}}$.
In particular, at the inner radial boundary we impose

$\mathbf{z}^{-}=\mathrm{Z_{0}}\sqrt{2}\sin(2\pi
f_{0}t)\mathbf{\hat{\theta}}+\mathrm{Z_{0}}\sqrt{2}\cos(2\pi
f_{0}t)\mathbf{\hat{\phi}}$ (12)

where $\mathrm{Z_{0}}=32$ km s-1 and $f_{0}=1~{}$mHz. This corresponds to
injecting an anti-sunward (outward) circularly polarised Alfvén wave into the
magnetic flux-tube. It is known that the MHD equations permit the existence of
finite-amplitude Alfvén waves with circular polarisation which do not generate
variations in the total magnetic field and thermal pressure as they propagate
(Barnes & Hollweg, 1974; Hollweg, 1974). In low $\beta$ plasma the circularly
polarised wave can couple with small random fluctuations in the magnetic
field, such as when the background field has gradients, to trigger the
parametric decay instability (Goldstein, 1978). This instability occurs in our
simulation as the solar wind obtained through Section 2.1 contains
$\beta<0.45$ in the simulation domain prior to the injection of the Alfvén
wave. The specific choice of $\mathrm{Z_{0}}$ and $f_{0}$ ensures that the
injected Alfvén wave does not significantly drive the solar wind as the
injected energy is much smaller than that required to sustain the wind in an
open field region (Withbroe & Noyes, 1977), while still ensuring a sufficient
growth rate for the PDI. The PDI results into decay of the anti-sunward
propagating Alfvén wave decays into a reflected sunward propagating Alfvén
wave, and anti-sunward propagating MHD slow wave. Therefore, the PDI satisfies
the conditions $\mathbf{k}_{1}=\mathbf{k}_{2}+\mathbf{k}_{3}$ and
$\omega_{1}=\omega_{2}+\omega_{3}$ where the subscripts 1-3 refer to the anti-
sunward and sunward Alfvén waves, and the slow wave, respectively (Chandran,
2018). The development of the PDI instability depends on the frequency of the
injected Alfvén wave (Shoda et al., 2018a). In addition, the expansion factor
of the flux-tube suppresses the development of PDI for lower frequency Alfvén
waves (Sishtla et al., 2022). In this study we have chosen appropriate values
of $f_{\mathrm{max}}$, and $f_{0}$ which excite the PDI instability, with
Sishtla et al. (2022) providing a complete description of the Alfvén wave
dynamics inside a flux-tube as prescribed in this simulation.

Figure 2: Dynamical evolution of the solar wind plasma along the flux-tube. In
every panel, each column represents the quantity (indicated with the colorbar)
at the given time along the flux tube. The panels show the parameters (a)
$v_{\phi}$, (b) $v_{\theta}$, (c) $v_{r}$, and (d) $v_{r}/v_{a}$ for the case
where a $1~{}$mHz Alfvén wave is continually injected starting at $t=0$. In
addition, the density fluctuations defined as $\delta\rho/\rho_{0}$ where
$\rho_{0}$ is the time averaged density are presented for different averaging
intervals of (e) $0.5~{}$mins, (f) $1.25~{}$mins, (g) $2.5~{}$mins, and (h)
$5~{}$mins.

In Figure 2(a)-(d) the dynamical evolution of the solar wind is presented
after the start of the injection of the Alfvén wave. Panels (a) and (b) show
the transverse ${v_{\phi}}$ and ${v_{\theta}}$ velocity components that are
perturbed due to the introduction of the circularly polarised wave. The start
of the wave injection is at $t=0$ mins ($x$-axis) and the waves are seen to
propagate through the corona towards the outer boundary. From $t\approx
225~{}$ mins onwards, a significant change in the dynamics originating at
approximately between $r\sim 6R_{\odot}$ and $r\sim 8R_{\odot}$ occurs. At
this time, the creation of two different space-time paths in $v_{\phi,\theta}$
is visible, with one path propagating sunward and the other anti-sunward. At
this same instance, the development of fluctuations in $v_{r}$ (panel (c)),
along with variations in the Alfvén critical point (panel (d)) is evident.

In panels (e)-(h) the evolution of the level of density fluctuations is shown,
defined as $\delta\rho/\rho_{0}$ where $\rho_{0}$ is the time-averaged
density, and $\delta\rho=\rho-\rho_{0}$ is the corresponding density
fluctuation around $\rho_{0}$. Panels (a)-(d) present the density fluctuations
computed using different averaging windows, which allows us to study the
fluctuations at different frequencies. This figure shows that the density
fluctuations start to develop at various scales at $t\approx 225~{}$ min, with
higher absolute values of $\delta\rho/\rho_{0}$ in panel (d) where the
averaging interval captures frequencies closer to the pump wave. It can then
be seen through the emerging density fluctuations in panels (e)-(h) along with
the decay of the injected wave at the same time in panels (a)-(d) that it is
the parametric decay instability that is excited at $t\approx 225~{}$ min.

The different tracks in panels (a), (b) correspond to the waves propagating
with phase speeds $v_{r}\pm v_{a}$ (where $v_{a}$ is the Alfvén speed) in the
anti-sunward and sunward directions (Verdini et al., 2009). This modified
$\mathbf{v_{\perp}}$ then drive variations in $v_{r}$ through Equation 5 which
we observe in panel (c). Moreover, the $v_{\phi,\theta}$ tracks indicate a
quasi-periodic structure of the solar wind as the Alfvén waves are continually
injected. This structure indicates that the density fluctuations would be
similarly quasi-period which is seen in panels (e)-(h).

##  3 Modelling the effect of PDI on the morphology of type III radio bursts

Figure 1 shows two example intervals during PSP’s close encounter 9. The most
striking feature in the figure is the strongly structured morphology of the
type III bursts that are uniquely observed within the fundamental emission.
Assuming that the fundamental emission is generated via linear mechanisms
(Voshchepynets et al., 2015; Voshchepynets & Krasnoselskikh, 2015;
Krasnoselskikh et al., 2019; Jebaraj et al., 2023b), the generation of
morphological features within a type III should be facilitated by density
inhomogeneities at different spatial scales present in the medium. To this
degree, we model the generation of the fundamental emission as a near-
relativistic propagating source that interacts with the ambient plasma to
generate electromagnetic (EM) waves based on the local electron plasma
frequency $f_{\mathrm{pe}}$,

$f_{\mathrm{pe}}\,\mathrm{[MHz]}\approx 9\times
10^{-3}\sqrt{n\,\mathrm{[cm^{-3}]}}$ (13)

where $n=\rho/m_{p}$ is the number density and $m_{p}$ is the proton mass.
Since it is not possible to simulate beam-plasma interactions within a MHD
plasma description, we construct an emission proxy approximating this process
for a given speed $v_{e}$ and time of injection $t_{i}$ of the radiation
source. The radiation source would represent a sampling of the suprathermal
electron beam propagating along the open magnetic field lines. If we assume a
constant velocity for the radiation source, it arrives at the outer boundary
at time $t=d/v_{e}$, with $d$ representing the length of the flux tube in the
MHD simulation. Finally, the location of the radiation source $r_{e}$ at any
intermediate time when it is in the simulation domain can be expressed as
$r_{e}(t)=r_{0}+v_{e}t$, where $r_{0}$ is the location where the radiation
source originates. Thus, we can sample the location of the radiation source as
it propagates and use the background density from the MHD simulation to
calculate the EM emission at the location $r_{e}$ using Equation 13. In other
words, we make a snapshot of the background plasma which represents the
variations in the background electron plasma frequency in the modelled flux
tube.

Figure 3: The spatio-temporal dependence of the (a) local plasma frequency in
the solar wind, and (b) the associated simulated dynamic spectrum for a
radiation source injected at $t=1000~{}$mins is shown. The path of the
radiation source is annotated in (a). Subsequently, panels (c) and (d) present
the local plasma frequency annotated with the path of a beamlet and the
associated dynamic spectrum. In panel (d) the dynamic spectrum is coloured
based on the $5~{}$mins averaged density fluctuations encountered by the
beamlet. Similar to panels (a) and (b), the panels (e) and (f) present the
modelled emission for a radiation source but in a solar wind injected with
Alfvén waves with frequencies $1,~{}2,~{}$and $3~{}$mHz. The panels (g) and
(h) present the modelled emissions for a similar beamlet as in panels (c) and
(d) but in a solar wind injected with the Alfvén waves of three different
frequencies.

This process of simulating the Type III radio burst is illustrated in Figure
3. In panel (a), we show the solar wind’s local electron plasma frequency,
which contains quasi-periodic variations as the PDI generates density
fluctuations. We then inject a fast moving radiation source at $t=1000~{}$min
with a velocity of $0.05c$, which propagates through the plasma along the
trajectory annotated in panel (a). We sample the radiation source’s location
as it propagates and plot the associated emission frequency in panel (b),
representing the snapshot of the plasma along the radiation source propagation
path. This simulated dynamic spectrum in panel (b) is plotted for emission
frequencies between $1$ MHz and $100$ MHz, after the radiation source
propagates away from the low solar corona which is associated with
significantly higher frequencies due to the increased density (Equation 13).
The striations of various scales and amplitudes caused by density fluctuations
in panel (b) is characteristic of a Type III radio burst. As this methodology
does not involve simulating beam-plasma interactions, the features of the
simulated dynamic spectrum are not affected by the exact value of the near-
relativistic velocity of the radiation source.

We can now extend this methodology to model a beamlet, a collection of
radiation sources with an associated beam angle, as illustrated in panels (c)
and (d). The modelled beamlet consists of multiple radiation sources injected
at the same time but with slightly varying constant speeds along the radial
magnetic field. We define a beam angle in phase-space, thus denoting the angle
between the fastest and slowest radiation source. In panels (c) and (d) we
launch $100$ radiation sources at $t=1000~{}$mins with the fastest speed being
$0.05c$ with a beam angle of $5\times 10^{-6}$ degrees. The local electron
plasma frequency as calculated from Equation 13 along with the path of the
beamlet is presented in panel (c). The simulated dynamic spectrum due to this
beamlet is shown in panel (d) color-coded by the $5~{}$min averaged density
fluctuations $\delta\rho/\rho_{0}$ at the location of emission. These
$5~{}$min averaged fluctuations are the larger scale density perturbations
present in the simulated solar wind. The density fluctuations at such scales
have greater $\left|\delta\rho/\rho_{0}\right|$ as they are closer to the
injected Alfvén wave frequency due to the decay process involving the pump
wave parametrically generating lower frequency wave distributions (see e.g.
Chandran, 2018). The dynamic spectrum in panel (d) is seen to rise with an
accumulations of plasma ($\delta\rho>0$) and falls due to a deficit
($\delta\rho<0$) as the emission frequency of the electron beamlet is higher
for an increased density (Equation 13). Therefore, the striations in the
dynamic spectrum are a consequence of a time of flight effect of the electrons
propagating in the wind.

###  3.1 Spectral striations are the natural response to beam propagation in
inhomogeneous plasma

Consider a radiation source propagating through the coronal plasma, leading to
emissions at the plasma frequency,

$2\pi f_{\mathrm{pe}}=\sqrt{\frac{ne^{2}}{\epsilon_{0}m_{e}}}$ (14)

where $n$ is the electron number density, $e$ is the elementary charge,
$m_{e}$ is the mass of the electron, and $\epsilon_{0}$ is the permittivity of
vacuum. Then, for any perturbations in $n$, the corresponding perturbation in
emission frequency is

$2\pi\delta
f_{\mathrm{pe}}=\frac{1}{2}\sqrt{\frac{e^{2}}{n\epsilon_{0}m_{e}}}\delta
n=\frac{\delta n}{2n}\sqrt{\frac{ne^{2}}{\epsilon_{0}m_{e}}}=\frac{\delta
n}{2n}2\pi f_{\mathrm{pe}}$ (15)

Therefore,

$\frac{\delta f_{\mathrm{pe}}}{f_{\mathrm{pe}}}=\frac{1}{2}\frac{\delta
n}{n}=\frac{1}{2}\frac{\delta\rho}{\rho}$ (16)

For cases when $-1\leq\delta\rho/\rho\leq 1$, then
$-f_{\mathrm{pe}}/2\leq\delta f_{\mathrm{pe}}\leq f_{\mathrm{pe}}/2$. Thus,
increased $\left|\delta\rho/\rho\right|$ will directly correlate with greater
striations in the dynamic spectrum. To illustrate this further we present
Figures 3(e), 3(f) and Figures 3(g), 3(h), which show the simulated dynamic
spectrum for both a radiation source and a beamlet propagating in a solar wind
which is subject to the simultaneous injection of Alfvén waves at the three
frequencies $1$ mHz, $2$ mHz, and $3$ mHz. This is in contrast to the
preceding sections where $\delta\rho/\rho_{0}$ variations were the consequence
of a single frequency wave injection. By introducing additional wave modes we
expect more wave-wave interactions leading to greater number of striations.

Panel (e) presents the dynamic spectrum for a radiation source injected at
$t~{}=~{}1000~{}$mins. As compared to the case of injecting a single wave mode
in panel (a), panel (e) shows finer fluctuations in the local plasma
frequency, and this is captured along the electron path in panel (f). If we
further launch $100$ radiation sources at $t~{}=~{}1000~{}$mins with a beam
angle of $5\times 10^{-6}$ degrees, the simulated dynamic spectrum in panels
(g), (h) similarly shows greater striations than panels (c), (d). Our model
therefore suggests that the striations in the dynamic spectrum are a
consequence of a time of flight effect of the electrons propagating in the
wind that embeds density perturbations. Thus, in a more general solar wind
plasma containing a well-developed spectrum of density fluctuations, the time
of flight effect of the propagating electrons would lead to structured
striations in the dynamic spectra of type III bursts such as those shown in
Figure 1.

##  4 Conclusion

In this study we simulated the effect of coronal density fluctuations on the
fundamental band of a type III radio burst. The fluctuations were self-
consistently generated by injecting an Alfvén wave into the low corona, which
subsequently propagates into the higher corona and triggers the PDI
instability via which density fluctuations are generated in the solar wind.
Such instabilities and the associated density fluctuations are expected to be
present in the solar corona where the low plasma beta conditions result in a
sufficient growth rate for PDI. The synthetic type III burst produced in this
study was affected by the multi-scale density fluctuations in the solar wind.
These scales correspond to hundreds to thousands of Debye length with a
maximum of the relaxation length of the beam. In this regard, the synthetic
type III represent a snapshot of the fluctuations of the same level in the
background density along a flux tube. This is relevant considering recent PSP
observations indicating that most type III bursts recorded during close
encounters are structured at various spatio-temporal scales (Pulupa et al.,
2020; Jebaraj et al., 2023a). Furthermore, the presence of dynamic wave
activity in the inner heliosphere has also been widely recognized (Malaspina
et al., 2020; Mozer et al., 2023). It is possible that the consistent
appearance of structured type III radio bursts observed when the probe was
close to the Sun indicates the presence of a fundamental process which
mediates this. The results presented in this study indicate the possibility of
radio bursts affected by Alfvénic turbulence in the solar corona which might
lead them to exhibit variations in the emission frequency based on the
spectral properties of turbulent density fluctuations. Conversely, examining
the emissions in the radio burst might provide a means of probing the density
fluctuations inside the flux-tube where the type III burst originated. The
exact spatial scales of the fluctuations which affect the beam-plasma system
are highly dependent upon the characteristics of the beam such as its density.
As previously mentioned, this relaxation length of the electron beam will then
form the upper limit for the observed striations (Jebaraj et al., 2023b).

While our simulations self-consistently generate density fluctuations via the
PDI, there are also other processes close to the Sun that may generate density
fluctuations. Inhomogenous plasma conditions in the direction transverse to
the propagating wave vector lead to phase mixing which is an alternate source
of generating fast waves, and by extension compressive fluctuations
(Nakariakov et al., 2000). This phase mixing can independently generate
turbulent structures in the solar wind (Magyar et al., 2017). More recently,
it has been seen that small-scale reconnection on the solar surface can drive
jetlets and switchbacks which contribute to the compressive nature of the wind
(Raouafi & Stenborg, 2014; Raouafi et al., 2023). This is well complemented by
Bale et al. (2019, 2023) who have found a strong correspondence between the
solar wind arising from coronal hole boundaries and reconnection driven
processes. Therefore, compressive fluctuations driven by recconection may also
contribute to the turbulent environment which can affect the morphology of the
type III bursts. Thus, density fluctuations that are responsible for causing
variations in the emission frequencies of radio sources may arise from a
multitude of non-linear processes in the solar corona and further out in
interplanetary space. Nevertheless, our results show that even a simple
scenario where PDI is responsible for producing density variations in the
solar wind is able to produce synthetic radio spectra consistent with the
observed emission. Subsequently, analysing type III bursts in a high temporal
resolution, such as through PSP, can provide indirect information on the
relative density fluctuations (Krupar et al., 2020) present in the wind. The
scenario presented in this study proves as a first step in utilizing remote
sensing observations such as the ones from radio wavelengths, e.g. using a
forward modelling approach, to constrain the evolution of plasma in regions
not generally amenable to direct measurements.

## Acknowledgements

Parker Solar Probe was designed, built, and is now operated by the Johns
Hopkins Applied Physics Laboratory as part of NASA’s Living with a Star (LWS)
program (contract NNN06AA01C). Support from the LWS management and technical
team has played a critical role in the success of the Parker Solar Probe
mission. The work has been supported by the Finnish Centre of Excellence in
Research on Sustainable Space (FORESAIL; grant no. 336807). This is a project
under the Academy of Finland, and this research has been supported by the
European Research Council (SolMAG; grant no. 724391) as well as Academy of
Finland project SWATCH (343581). Open access is funded by Helsinki University
Library. I.C.J, M.P, and S.D.B acknowledge support from the International
Space Science Institute (ISSI) in Bern, through ISSI International Team
project No.557, “Beam-Plasma Interaction in the Solar Wind and the Generation
of Type III Radio Bursts”. I.C.J also acknowledges support from the ISSI
visiting scientist program. The authors would also like to thank Dr. V. V.
Krasnoselskikh for the active discussions on the topic of parametric decay
instability and radio generation mechanisms.

## References

  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, Space Sci. Rev., 204, 49, doi: 10.1007/s11214-016-0244-5
  * Bale et al. (2019) Bale, S. D., Badman, S. T., Bonnell, J. W., et al. 2019, Nature, 576, 237, doi: 10.1038/s41586-019-1818-7
  * Bale et al. (2023) Bale, S. D., Drake, J. F., McManus, M. D., et al. 2023, Nature, 618, 252, doi: 10.1038/s41586-023-05955-3
  * Barnes & Hollweg (1974) Barnes, A., & Hollweg, J. V. 1974, J. Geophys. Res., 79, 2302, doi: 10.1029/JA079i016p02302
  * Belcher & Davis (1971) Belcher, J. W., & Davis, Leverett, J. 1971, J. Geophys. Res., 76, 3534, doi: 10.1029/JA076i016p03534
  * Bourdin (2017) Bourdin, P. A. 2017, ApJ, 850, L29, doi: 10.3847/2041-8213/aa9988
  * Bowen et al. (2018) Bowen, T. A., Badman, S., Hellinger, P., & Bale, S. D. 2018, ApJ, 854, L33, doi: 10.3847/2041-8213/aaabbe
  * Bruno & Carbone (2013) Bruno, R., & Carbone, V. 2013, Living Reviews in Solar Physics, 10, 2, doi: 10.12942/lrsp-2013-2
  * Cairns (1987) Cairns, I. H. 1987, Journal of Plasma Physics, 38, 169, doi: 10.1017/S0022377800012496
  * Cairns & Fung (1988) Cairns, I. H., & Fung, S. F. 1988, J. Geophys. Res., 93, 7307, doi: 10.1029/JA093iA07p07307
  * Cairns & Melrose (1985) Cairns, I. H., & Melrose, D. B. 1985, J. Geophys. Res., 90, 6637, doi: 10.1029/JA090iA07p06637
  * Chandran (2005) Chandran, B. D. G. 2005, Phys. Rev. Lett., 95, 265004, doi: 10.1103/PhysRevLett.95.265004
  * Chandran (2018) Chandran, B. D. G. 2018, Journal of Plasma Physics, 84, 905840106, doi: 10.1017/S0022377818000016
  * Chandran et al. (2011) Chandran, B. D. G., Dennis, T. J., Quataert, E., & Bale, S. D. 2011, ApJ, 743, 197, doi: 10.1088/0004-637X/743/2/197
  * Chen (2016) Chen, C. H. K. 2016, Journal of Plasma Physics, 82, 535820602, doi: 10.1017/S0022377816001124
  * Coleman (1968) Coleman, Paul J., J. 1968, ApJ, 153, 371, doi: 10.1086/149674
  * Cranmer & van Ballegooijen (2005) Cranmer, S. R., & van Ballegooijen, A. A. 2005, ApJS, 156, 265, doi: 10.1086/426507
  * D’Amicis & Bruno (2015) D’Amicis, R., & Bruno, R. 2015, ApJ, 805, 84, doi: 10.1088/0004-637X/805/1/84
  * De Pontieu et al. (2007) De Pontieu, B., McIntosh, S. W., Carlsson, M., et al. 2007, Science, 318, 1574, doi: 10.1126/science.1151747
  * Del Zanna et al. (2001) Del Zanna, L., Velli, M., & Londrillo, P. 2001, A&A, 367, 705, doi: 10.1051/0004-6361:20000455
  * Derby (1978) Derby, N. F., J. 1978, ApJ, 224, 1013, doi: 10.1086/156451
  * Dorfman & Carter (2016) Dorfman, S., & Carter, T. A. 2016, Phys. Rev. Lett., 116, 195002, doi: 10.1103/PhysRevLett.116.195002
  * Elsasser (1950) Elsasser, W. M. 1950, Physical Review, 79, 183, doi: 10.1103/PhysRev.79.183
  * Fox et al. (2016) Fox, N. J., Velli, M. C., Bale, S. D., et al. 2016, Space Sci. Rev., 204, 7, doi: 10.1007/s11214-015-0211-6
  * Galeev & Oraevskii (1963) Galeev, A., & Oraevskii, V. 1963, in Soviet Physics Doklady, Vol. 7, 988
  * Ginzburg & Zhelezniakov (1958) Ginzburg, V. L., & Zhelezniakov, V. V. 1958, Soviet Ast., 2, 653
  * Goddard et al. (2016) Goddard, C. R., Nisticò, G., Nakariakov, V. M., Zimovets, I. V., & White, S. M. 2016, A&A, 594, A96, doi: 10.1051/0004-6361/201628478
  * Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763, doi: 10.1086/175121
  * Goldstein (1978) Goldstein, M. L. 1978, ApJ, 219, 700, doi: 10.1086/155829
  * Higdon (1984) Higdon, J. C. 1984, ApJ, 285, 109, doi: 10.1086/162481
  * Hinkel-Lipsker et al. (1992) Hinkel-Lipsker, D. E., Fried, B. D., & Morales, G. J. 1992, Physics of Fluids B, 4, 559, doi: 10.1063/1.860255
  * Hollweg (1974) Hollweg, J. V. 1974, J. Geophys. Res., 79, 1539, doi: 10.1029/JA079i010p01539
  * Hoshino & Goldstein (1989) Hoshino, M., & Goldstein, M. L. 1989, Physics of Fluids B, 1, 1405, doi: 10.1063/1.858971
  * Inhester (1990) Inhester, B. 1990, J. Geophys. Res., 95, 10525, doi: 10.1029/JA095iA07p10525
  * Iwai et al. (2014) Iwai, K., Shibasaki, K., Nozawa, S., et al. 2014, Earth, Planets and Space, 66, 149, doi: 10.1186/s40623-014-0149-z
  * Jebaraj et al. (2023a) Jebaraj, I. C., Krasnoselskikh, V., Pulupa, M., Magdalenic, J., & Bale, S. D. 2023a, ApJ, 955, L20, doi: 10.3847/2041-8213/acf857
  * Jebaraj et al. (2023b) Jebaraj, I. C., Magdalenic, J., Krasnoselskikh, V., Krupar, V., & Poedts, S. 2023b, A&A, 670, A20, doi: 10.1051/0004-6361/202243494
  * Kaneda et al. (2018) Kaneda, K., Misawa, H., Iwai, K., et al. 2018, ApJ, 855, L29, doi: 10.3847/2041-8213/aab2a5
  * Kim et al. (2007) Kim, E.-H., Cairns, I. H., & Robinson, P. A. 2007, Phys. Rev. Lett., 99, 015003, doi: 10.1103/PhysRevLett.99.015003
  * Kissmann & Pomoell (2012) Kissmann, R., & Pomoell, J. 2012, SIAM Journal on Scientific Computing, 34, A763, doi: 10.1137/110834329
  * Kolmogorov (1991) Kolmogorov, A. N. 1991, Proceedings of the Royal Society of London Series A, 434, 9, doi: 10.1098/rspa.1991.0075
  * Kolotkov et al. (2018) Kolotkov, D. Y., Nakariakov, V. M., & Kontar, E. P. 2018, ApJ, 861, 33, doi: 10.3847/1538-4357/aac77e
  * Kontar et al. (2017) Kontar, E. P., Yu, S., Kuznetsov, A. A., et al. 2017, Nature Communications, 8, 1515, doi: 10.1038/s41467-017-01307-8
  * Kopp & Holzer (1976) Kopp, R. A., & Holzer, T. E. 1976, Sol. Phys., 49, 43, doi: 10.1007/BF00221484
  * Krafft & Savoini (2022) Krafft, C., & Savoini, P. 2022, ApJ, 924, L24, doi: 10.3847/2041-8213/ac46a7
  * Krafft et al. (2013) Krafft, C., Volokitin, A. S., & Krasnoselskikh, V. V. 2013, ApJ, 778, 111, doi: 10.1088/0004-637X/778/2/111
  * Krasnoselskikh et al. (2019) Krasnoselskikh, V., Voshchepynets, A., & Maksimovic, M. 2019, ApJ, 879, 51, doi: 10.3847/1538-4357/ab22bf
  * Krupar et al. (2020) Krupar, V., Szabo, A., Maksimovic, M., et al. 2020, ApJS, 246, 57, doi: 10.3847/1538-4365/ab65bd
  * Lashmore-Davies (1976) Lashmore-Davies, C. N. 1976, Physics of Fluids, 19, 587, doi: 10.1063/1.861493
  * Magyar et al. (2017) Magyar, N., Van Doorsselaere, T., & Goossens, M. 2017, Scientific Reports, 7, 14820, doi: 10.1038/s41598-017-13660-1
  * Malaspina et al. (2010) Malaspina, D. M., Cairns, I. H., & Ergun, R. E. 2010, Journal of Geophysical Research (Space Physics), 115, A01101, doi: 10.1029/2009JA014609
  * Malaspina et al. (2012) —. 2012, ApJ, 755, 45, doi: 10.1088/0004-637X/755/1/45
  * Malaspina et al. (2020) Malaspina, D. M., Halekas, J., Berčič, L., et al. 2020, The Astrophysical Journal Supplement Series, 246, 21, doi: 10.3847/1538-4365/ab4c3b
  * McLean & Labrum (1985) McLean, D. J., & Labrum, N. R. 1985, Solar radiophysics : studies of emission from the sun at metre wavelengths
  * Melrose (1980) Melrose, D. B. 1980, Space Sci. Rev., 26, 3, doi: 10.1007/BF00212597
  * Mikić et al. (2018) Mikić, Z., Downs, C., Linker, J. A., et al. 2018, Nature Astronomy, 2, 913, doi: 10.1038/s41550-018-0562-5
  * Mozer et al. (2023) Mozer, F., Bale, S., Kellogg, P., et al. 2023, Physics of Plasmas, 30, 062111, doi: 10.1063/5.0151423
  * Mugundhan et al. (2017) Mugundhan, V., Hariharan, K., & Ramesh, R. 2017, Sol. Phys., 292, 155, doi: 10.1007/s11207-017-1181-5
  * Nakariakov et al. (2004) Nakariakov, V. M., Arber, T. D., Ault, C. E., et al. 2004, MNRAS, 349, 705, doi: 10.1111/j.1365-2966.2004.07537.x
  * Nakariakov et al. (2000) Nakariakov, V. M., Ofman, L., & Arber, T. D. 2000, A&A, 353, 741
  * Newkirk (1961) Newkirk, Gordon, J. 1961, ApJ, 133, 983, doi: 10.1086/147104
  * Papadopoulos & Freund (1978) Papadopoulos, K., & Freund, H. P. 1978, Geophys. Res. Lett., 5, 881, doi: 10.1029/GL005i010p00881
  * Pick & Vilmer (2008) Pick, M., & Vilmer, N. 2008, A&A Rev., 16, 1, doi: 10.1007/s00159-008-0013-x
  * Pomoell et al. (2015) Pomoell, J., Aran, A., Jacobs, C., et al. 2015, Journal of Space Weather and Space Climate, 5, A12, doi: 10.1051/swsc/2015015
  * Pomoell & Vainio (2012) Pomoell, J., & Vainio, R. 2012, ApJ, 745, 151, doi: 10.1088/0004-637X/745/2/151
  * Pulupa et al. (2017) Pulupa, M., Bale, S. D., Bonnell, J. W., et al. 2017, Journal of Geophysical Research (Space Physics), 122, 2836, doi: 10.1002/2016JA023345
  * Pulupa et al. (2020) Pulupa, M., Bale, S. D., Badman, S. T., et al. 2020, ApJS, 246, 49, doi: 10.3847/1538-4365/ab5dc0
  * Raouafi & Stenborg (2014) Raouafi, N. E., & Stenborg, G. 2014, ApJ, 787, 118, doi: 10.1088/0004-637X/787/2/118
  * Raouafi et al. (2023) Raouafi, N. E., Stenborg, G., Seaton, D. B., et al. 2023, ApJ, 945, 28, doi: 10.3847/1538-4357/acaf6c
  * Robinson & Cairns (1998a) Robinson, P. A., & Cairns, I. H. 1998a, Sol. Phys., 181, 363, doi: 10.1023/A:1005018918391
  * Robinson & Cairns (1998b) —. 1998b, Sol. Phys., 181, 395, doi: 10.1023/A:1005033015723
  * Robinson & Cairns (1998c) —. 1998c, Sol. Phys., 181, 429, doi: 10.1023/A:1005023002461
  * Robinson et al. (1993) Robinson, P. A., Cairns, I. H., & Gurnett, D. A. 1993, ApJ, 407, 790, doi: 10.1086/172560
  * Robinson et al. (1994) Robinson, P. A., Cairns, I. H., & Willes, A. J. 1994, ApJ, 422, 870, doi: 10.1086/173779
  * Sagdeev & Galeev (1969) Sagdeev, R. Z., & Galeev, A. A. 1969, Nonlinear Plasma Theory
  * Sakai & Sonnerup (1983) Sakai, J. I., & Sonnerup, U. O. 1983, J. Geophys. Res., 88, 9069, doi: 10.1029/JA088iA11p09069
  * Schmitz & Tsiklauri (2013) Schmitz, H., & Tsiklauri, D. 2013, Physics of Plasmas, 20, 062903, doi: 10.1063/1.4812453
  * Shoda et al. (2018a) Shoda, M., Yokoyama, T., & Suzuki, T. K. 2018a, ApJ, 860, 17, doi: 10.3847/1538-4357/aac218
  * Shoda et al. (2018b) —. 2018b, ApJ, 853, 190, doi: 10.3847/1538-4357/aaa3e1
  * Sishtla et al. (2022) Sishtla, C. P., Pomoell, J., Kilpua, E., et al. 2022, A&A, 661, A58, doi: 10.1051/0004-6361/202142999
  * Sishtla et al. (2023) Sishtla, C. P., Pomoell, J., Vainio, R., Kilpua, E., & Good, S. 2023, A&A, 679, A54, doi: 10.1051/0004-6361/202347250
  * Suzuki & Dulk (1985) Suzuki, S., & Dulk, G. A. 1985, in Solar Radiophysics: Studies of Emission from the Sun at Metre Wavelengths, ed. D. J. McLean & N. R. Labrum, 289–332
  * Suzuki & Inutsuka (2006) Suzuki, T. K., & Inutsuka, S.-I. 2006, Journal of Geophysical Research (Space Physics), 111, A06101, doi: 10.1029/2005JA011502
  * Thurgood & Tsiklauri (2015) Thurgood, J. O., & Tsiklauri, D. 2015, A&A, 584, A83, doi: 10.1051/0004-6361/201527079
  * Tsiklauri (2011) Tsiklauri, D. 2011, Physics of Plasmas, 18, 052903, doi: 10.1063/1.3590928
  * Usmanov et al. (2000) Usmanov, A. V., Goldstein, M. L., Besser, B. P., & Fritzer, J. M. 2000, J. Geophys. Res., 105, 12675, doi: 10.1029/1999JA000233
  * Verdini et al. (2009) Verdini, A., Velli, M., & Buchlin, E. 2009, ApJ, 700, L39, doi: 10.1088/0004-637X/700/1/L39
  * Voshchepynets & Krasnoselskikh (2015) Voshchepynets, A., & Krasnoselskikh, V. 2015, Journal of Geophysical Research (Space Physics), 120, 10,139, doi: 10.1002/2015JA021705
  * Voshchepynets et al. (2015) Voshchepynets, A., Krasnoselskikh, V., Artemyev, A., & Volokitin, A. 2015, ApJ, 807, 38, doi: 10.1088/0004-637X/807/1/38
  * Withbroe & Noyes (1977) Withbroe, G. L., & Noyes, R. W. 1977, ARA&A, 15, 363, doi: 10.1146/annurev.aa.15.090177.002051
  * Yoon et al. (1994) Yoon, P. H., Wu, C. S., Vinas, A. F., et al. 1994, J. Geophys. Res., 99, 23,481, doi: 10.1029/94JA02489
  * Zhang et al. (2018) Zhang, P. J., Wang, C. B., & Ye, L. 2018, A&A, 618, A165, doi: 10.1051/0004-6361/201833260