# Ising model on a _restricted_ scale-free network

R. A. Dumer<EMAIL_ADDRESS>Instituto de Física - Universidade
Federal de Mato Grosso, 78060-900, Cuiabá, Mato Grosso, Brazil.  M. Godoy
<EMAIL_ADDRESS>Instituto de Física - Universidade Federal de Mato
Grosso, 78060-900, Cuiabá, Mato Grosso, Brazil.

###### Abstract

The Ising model on a _restricted_ scale-free network (SFN) has been studied
employing Monte Carlo simulations. This network is described by a power-law
degree distribution in the form $P(k)\sim k^{-\alpha}$, and is called
restricted, because independently of the network size, we always have fixed
the maximum $k_{m}$ and a minimum $k_{0}$ degree on distribution, being that
for it, we only limit the minimum network size of the system. We calculated
the thermodynamic quantities of the system, such as, the magnetization per
spin $\textrm{m}_{\textrm{L}}$, the magnetic susceptibility
$\textrm{$\chi$}_{\textrm{L}}$, and the reduced fourth-order Binder cumulant
$\textrm{U}_{\textrm{L}}$, as a function of temperature $T$ for several values
of lattice size $N$ and exponent $1\leq\alpha\leq 5$. For the values of
$\alpha$, we have obtained the finite critical points due to we also have
finite second and fourth moments in the degree distribution, and the phase
diagram was constructed for the equilibrium states of the model in the plane
$T$ versus $k_{0}$, $k_{m}$, and $\alpha$, showing a transition between the
ferromagnetic $F$ to paramagnetic $P$ phases. Using the finite-size scaling
(FSS) theory, we also have obtained the critical exponents for the system, and
a mean-field critical behavior is observed.

##  I Introduction

Hyperlinks pointing from one web page to another (World Wide Web), computers
physically linked (internet), actors that have acted in a movie together,
scientists that have an article together, and proteins that bind together
experimentally, are some of the cases that, when analyzed in terms of nodes
and edges of a network, they are part of a broad group of real systems in
which the degree distribution has a power-law tail [1, 2]. This degree
distribution has the form $P(k)\sim k^{-\alpha}$, representing the probability
of a site in the network to have a degree $k$, i.e., $k$ edges linked to it,
with exponent $\alpha$. Barabási-Albert [3] proposed a growing process of
network creation, in which the degree distribution resultant also has the form
of power-law. That growing process is based only on two generic mechanisms:
(i) networks expand continuously by the addition of new vertices, and (ii) new
vertices attach preferentially to sites that are already well connected. With
these mechanisms, most connected sites are most likely to receive new
connections, and a “rich-get-richer”, and self-organization phenomenon, as in
real networks, is observed.

Due to the applicability of the networks with a power-law degree distribution,
also called scale-free networks, it was been implemented in many physical
problems [4, 5, 6, 7, 8, 9, 10, 11]. Highlighting the critical phenomena
arising from the Ising model, from an analytical way, Dorogovtsev, Goltsev,
and Mendes showed that its critical behavior is very dependent on the
distribution of connections [12, 13]. From this dependence, when $\alpha>5$ in
$P(k)$, its fourth moment $\left\langle k^{4}\right\rangle$ is convergent and
a mean-field critical behavior is obtained. When, $3<\alpha\leq 5$, anomalous
behavior of the thermodynamics quantities is observed, due to divergence in
$\left\langle k^{4}\right\rangle$, but when $\alpha\leq 3$, the divergence is
in the second moment, $\left\langle k^{2}\right\rangle$, and the criticality
vary with the size of the system, being an infinity order phase transition in
the thermodynamic limit. In the network proposed in the Barabási-Albert (BA)
model [3], the exponent is limited by $\alpha=3$, and in previews of
approximate and numerical results, the infinity order phase transition is
verified [14, 15]. In addition to $\alpha=3$, the cases where $\left\langle
k^{4}\right\rangle$ and $\left\langle k^{2}\right\rangle$ are convergent and
divergent, by Monte Carlo simulations [16] and replica method [17], have been
proven the non-trivial critical exponents, beyond the and size-dependent
critical temperature.

In these Monte Carlo simulations, the SFN, constructed with a selected value
of $\alpha$, is absent from the two growing mechanisms in the BA model and
created distributing degrees for the vertices based on the exact value of
$P(k)=Ak^{-\alpha}$ [16]. For this exact value is predefined the minimum
$k_{0}$ and maximum $k_{m}$ degrees of the network, and calculated the
normalization constant of the distribution,
$A=\sum_{k=k_{0}}^{k_{m}}k^{\alpha}$. In this sense, supposing that the number
of sites with the degree $k_{m}$ is $N_{k_{m}}=1$, the network size with these
characteristics is given by $N=k_{m}^{\alpha}/A$. From this way, fixing
$k_{0}$ and varying the network sizes $N$, Herrero [16] could reproduce the
main critical phenomena seen analytically [12, 13], by the FSS analysis of
networks with also non-fixed $k_{m}$.

In this work, we investigated the Ising model on a _restricted_ SFN, where
each site of the network is occupied by a spin variable $\sigma=1/2$ that can
assume values $\pm 1$. Our network was built for various integer values of the
exponent $\alpha$, and divided into two sublattices, each connection
distributed by $P(k)$, a power-law degree distribution, which should connect
these sublattices. Besides that, in a similar way to what was proposed by
Herrero [16] to construct its random uncorrelated network, we also preview
define who will be $k_{0}$ and $k_{m}$ in our system. These values are kept
fixed while we vary the size $N$ of the network. Thereunto, the minimal
network size that we can use is defined by $N_{0}=k_{m}^{\alpha}/A$ and is not
always obtained $N_{k_{m}}=1$. However, as we have the same degrees on all
network sizes, we always have a convergent $\left\langle k^{2}\right\rangle$
and $\left\langle k^{4}\right\rangle$, and consequently finite transitions
point, for ferromagnetic to paramagnetic phases, in the whole values of
$\alpha$. Thus here, through Monte Carlo simulations, we have built phases
diagrams of critical temperature $T_{c}$ as a function of $k_{0}$ and $k_{m}$
for the studied values of $\alpha$, and using the FSS theory, we have obtained
the critical exponents for the system.

This article is organized as follows: In Section II, we describe the network
used and the Hamiltonian model of the system. In Section III, we present the
Monte Carlo simulation method, some details concerning the simulation
procedures and the thermodynamic quantities of the system also necessary for
the application of FSS analysis. The behavior of thermodynamic quantities,
phases diagrams, and critical exponents are described in Section IV. Finally,
in Section V, we present our conclusions.

##  II Model

Figure 1: Schematic representation of the _restricted_ SFN. The red circles
indicate the sites on one of the sublattices, the green circles are the sites
on the other sublattice, and the black solid lines are the connections between
the two sublattices. The size of the circles is proportional to the sites
degree, varying from $k_{0}=2$ to $k_{m}=8$, with $\alpha=3$, and $N=10^{2}$.

The Ising model studied in this work have $N=L^{2}$ spins $\sigma_{i}=\pm 1$,
on a _restricted_ SFN and a ferromagnetic interaction of strength $J_{ik}$. To
distribute the connections on this network, we first defined $k_{0}$, $k_{m}$
and $\alpha$, i.e., the minimum and maximum degree that the network is
required to have and its exponent on the distribution, respectively. Next, we
calculated the normalization constant of the distribution,
$A=\sum_{k=k_{0}}^{k_{m}}k^{\alpha}$, and found the smallest network that we
will use in the system, $N_{0}=k_{m}^{\alpha}/A$. With these values, we create
a set of site numbers, $\\{N_{k}\\}$, that will have the respective degrees
$k$, $N_{k}=AN/k^{\alpha}$, and distribute them by the network. For this
distribution of connections, we have divided the network into two sublattices,
where one sublattice plays the role of central spins, while the other
sublattice contains the spins in which the central spins can connect. Thus,
starting with the lowest degree, we select one site $i$ on the network, and
its sublattice will be the sublattice of central spins, then, from the other
sublattice, we select one site $j$ at random, that has not yet received their
respective connections. With this, we couple $j$ to the neighbors of $i$, and
for the site $j$ we couple the site $i$ to their neighbors. This process is
done until $i$ has their $k$ connections and visited all the set
$\\{N_{k}\\}$. It is valid to say that, for this network, there is not need to
create these sublattices and the connections can be created completely
randomly between free sites, but here, this implementation was done as a way
to prepare the system for future works in non-equilibrium systems without loss
of generality.

Figure 2: Log-log plot of the degree distribution on _restricted_ SFN for some
selected values of $\alpha$, as shown in the figure. All curves refer to the
network size $N=256^{2}$, minimum degree $k_{0}=4$ and maximum degree
$k_{m}=10$. The dotted lines have the exact expected slopes, $\alpha$, and the
errors of the slopes are in order of $10^{-4}$and the error bars are the same
size or smaller than the symbols.

In Fig. 1, we displayed an example of the network, with $\alpha=3$, $k_{0}=2$,
$k_{m}=8$ and $N=10^{2}$. The sites in the middle of the figure are the more
connected, while the peripheral sites are the less connected, and sites from
the sublattice green are only connected with sites from the sublattice red.

Based on this construction, throughout work, we have selected the integer
values of $1\leq\alpha\leq 5$, and network sizes $(32)^{2}\leq N\leq(256)^{2}$
to study the Ising model. In Fig. 2, we shown the degree distribution in the
largest network size of the system for these selected values of $\alpha$. With
these distributions in the log-log plot, we can see that the construction
method used here, guarantees the power-law form, as predicted in SFN.

The ferromagnetic Ising spin energy is described by the Hamiltonian on the
form

$\mathcal{H}=-\sum_{\left\langle i,j\right\rangle}J_{ij}\sigma_{i}\sigma_{j}$
(1)

where the sum is over all pair of spins, and $J_{ij}$ is the ferromagnetic
interaction, and assuming the value of unity if sites $i$ and $j$ are
connected by a link.

##  III Monte Carlo simulations

In the simulation of the system specified by the Hamiltonian in Eq. (1) and
with a _restricted_ SFN, we have chosen the initial state of the system with
all spins in the random states, and a new configuration is generated by a
Markov process. In this process, for a given temperature $T$, exponents
$\alpha$ of the degree distribution, network size $N$, and minimum $k_{0}$ and
maximum $k_{m}$ degree, we choose at random a spin $\sigma_{i}$ on the network
and change its state by the one-spin flip mechanism with a transition rate
given by the following Metropolis prescription

$W(\sigma_{i}\to\sigma_{i}^{\prime})=\left\\{\begin{array}[]{cccc}e^{\left(-\Delta
E/k_{B}T\right)}&\textrm{if}&\Delta E>0\\\ 1&\textrm{if}&\Delta E\leq
0&,\end{array}\right.$ (2)

where $\Delta E$ is the change in energy after flipping the spin,
$\sigma_{i}\to\sigma_{i}^{\prime}$, $k_{B}$ is the Boltzmann constant, and $T$
the temperature of the system. Therefore, in this scenario, the acceptance of
a new state is made if $\Delta E\leq 0$, but, in the case where $\Delta E>0$,
the acceptance is pondered by the probability $\exp\left(-\Delta
E/k_{B}T\right)$ and just it is accepted by choosing a random number
$0<\xi<1$, where $\xi\leq\exp\left(-\Delta E/k_{B}T\right)$. On the other
hand, if no one of these conditions was satisfied, we do not change the state
of the spin.

Repeating the Markov process $N$ times, we have a Monte Carlo Step (MCS). In
our simulations, we have waited $10^{4}$ MCS for the system reach the
equilibrium state, in all lattice sizes and value of the parameters. To
calculate the thermal averages of the interest quantities, we used more
$4\times 10^{4}$ MCS. The average over samples were done using $10$
independent samples for any configuration.

After reaching the equilibrium state, we have measured the following
thermodynamic quantities: magnetization per spin $\textrm{m}_{\textrm{L}}$,
magnetic susceptibility $\textrm{$\chi$}_{\textrm{L}}$ and reduced fourth-
order Binder cumulant $\textrm{U}_{\textrm{L}}$:

$\textrm{m}_{\textrm{L}}=\frac{1}{N}\left[\left\langle\sum_{i=1}^{N}\sigma_{i}\right\rangle\right],$
(3)

$\textrm{{\hbox{\chi}}}_{\textrm{L}}=\frac{N}{k_{B}T}\left[\left\langle\textrm{m}_{\textrm{L}}^{2}\right\rangle-\left\langle\textrm{m}_{\textrm{L}}\right\rangle^{2}\right],$
(4)

$\textrm{U}_{\textrm{L}}=1-\frac{\left[\left\langle\textrm{m}_{\textrm{L}}^{4}\right\rangle\right]}{3\left[\left\langle\textrm{m}_{\textrm{L}}^{2}\right\rangle^{2}\right]},$
(5)

where $\left[\ldots\right]$ represents the average over the samples, and
$\left\langle\ldots\right\rangle$ the thermal average over the MCS in the
equilibrium state. In the vicinity of the critical temperature $T_{c}$, the
above-defined quantities obey the following finite-size scaling relations:

$\textrm{m}_{\textrm{L}}=L^{-\beta/\nu}m_{0}(L^{1/\nu}\epsilon),$ (6)

$\textrm{{\hbox{\chi}}}_{\textrm{L}}=L^{\gamma/\nu}\chi_{0}(L^{1/\nu}\epsilon),$
(7)

$\textrm{U}_{\textrm{L}}^{\prime}=L^{1/\nu}\frac{U_{0}^{\prime}(L^{1/\nu}\epsilon)}{T_{c}},$
(8)

where $\epsilon=(T-T_{c})/T_{c}$, $m_{0}(L^{1/\nu}\epsilon)$,
$\chi_{0}(L^{1/\nu}\epsilon)$ and $U_{0}(L^{1/\nu}\epsilon)$ are the scaling
functions, and $\beta$, $\gamma$ and $\nu$ are the magnetization, magnetic
susceptibility and length correlation critical exponents, respectively.

Using the Eqs. (6), (7), (8) and the data from simulations for the network
sizes $32\leq L\leq 256$, we have obtained the critical exponents ratios,
$\beta/\nu$, $\gamma/\nu$ and $1/\nu$ from the slope of
$\textrm{m}_{\textrm{L}}(T_{c})$, $\textrm{$\chi$}_{\textrm{L}}(T_{c})$ and
$\textrm{U}_{\textrm{L}}^{\prime}(T_{c})$ as a function of $L$ in a log-log
plot. Aside from that, we also used data collapse from scaling functions to
estimate the critical exponent values.

##  IV Results

Figure 3: Thermodynamic quantities as a function of temperature $T$ for a
fixed value of $\alpha=1$, $k_{0}=4$ and $k_{m}=10$, and for different network
sizes, as present in the figure. (a) Magnetization $\textrm{m}_{\textrm{L}}$,
(b) reduced fourth-order Binder cumulant $\textrm{U}_{\textrm{L}}$, and (c)
susceptibility $\textrm{$\chi$}_{\textrm{L}}$.

The interesting results about the critical phenomena in complex networks, more
specifically on random uncorrelated networks, led us to understand the
importance of its degree distribution and respective moments. The Ising model
on the uncorrelated SFN, in the limit where $N\to\infty$, $k_{m}\to\infty$ and
the changing in its structure for $\alpha>3$, decreases the number of more
connected spins, admitting a finite-order phase transition until reaches the
standard mean-field critical behavior due to strong correlations of most
connected vertices in their neighborhood [12, 13]. However, here, by the
restriction of maximum $k_{m}$ and minimum $k_{0}$ degree, and minimum network
size, when $N\to\infty$, $k_{m}$ keeps finite, consequently changing the
number of more and less connected sites, their correlations, and the critical
phenomena as we can see in our results.

Figure 4: (a) Phase diagram of temperature $T$ as a function of $k_{0}$ for a
fixed value of $k_{m}=10$ and (b) as a function of $k_{m}$ for a fixed value
of $k_{0}=4$. Both figures representing transitions between $F$ to $P$ phases.
The solid lines are just guide for the eyes and the dashed lines represent the
analytical result obtained by Eq. (9).

To begin with, we have identified the point of transition between the
ferromagnetic $F$ to the paramagnetic $P$ phases through the curves of the
fourth-order Binder cumulant $\textrm{U}_{\textrm{L}}$ with different network
sizes [18, 19, 20, 21]. This critical point and second-order phase transition
can be identified by the crossing point of $\textrm{U}_{\textrm{L}}$ curves,
and an example is shown in Fig. 3. In this figure, we have presented as one of
the best results for thermodynamic quantities obtained by Eqs. (3), (4) and
(5), being that in Fig. 3(a), we can see the behavior of the magnetization
$\textrm{m}_{\textrm{L}}$ as a function of $T$, the reduced fourth-order
Binder cumulant $\textrm{U}_{\textrm{L}}$ (see Fig. 3(b)), and the magnetic
susceptibility $\textrm{$\chi$}_{\textrm{L}}$ (see Fig. 3(c)). In this case,
we have used a fixed value of $\alpha=1$, $k_{0}=4$ and $k_{m}=10$, and
different network sizes, as present in Fig. 3.

Figure 5: (a) Phase diagram of temperature $T$ as a function of $\alpha$, for
the transition between $F$ to $P$ phases, with fixed values of $k_{0}=4$ and
$k_{m}=10$. The dashed line is the theoretical result expected for random
uncorrelated networks and obtained using the Eq. (9). The red squares are the
simulation data points on the _restricted_ SFN, and calculated by the crossing
of $\textrm{U}_{\textrm{L}}$ curves. (b) Degree-degree correlations $r$ as a
function of $\alpha$ and obtined by the Eq. 10.

With the critical points obtained, the phase diagrams were built, see Fig. 4.
For these diagrams, the temperature as a function of $k_{0}$, $k_{m}$ and
$\alpha$ is studied, and for that, we have used $4\leq k_{0},k_{m}\leq 10$ and
integer values of $\alpha$ . The continuous transition lines can be seen in
Fig. 4(a) for temperature $T$ versus $k_{0}$, and for a fixed value of
$k_{m}=10$, and some selected values of $\alpha$, where we can verify that
have a finite critical point for all the values of $\alpha$. Because we have a
fixed value of $k_{m}$, as we increase $k_{0}$ the number of different degrees
on the network decrease until the whole system has the coordination number
$k_{0}=k_{m}=10$, and when that happens, the exponent of the degree
distribution is no longer important, returning us to a unique critical point
for all the exponents $\alpha$. On the other hand, in the case where
$k_{0}\neq k_{m}$, the degrees can be distributed on the network, and for
different values of $\alpha$ is showed the presence of different critical
points in the system, being that each value of $\alpha$ is responsible for a
distinct curve in Fig. 4(a). Instead, we can fix $k_{0}=4$ and varying
$k_{m}$, and as can seen in Fig. 4(b) for $T$ versus $k_{m}$. Thus, we can
observe the presence of the same qualitative behavior, like, continuous
transitions lines for the whole parameter values, distinct curves for distinct
exponents $\alpha$, and a unique critical point when $k_{m}=k_{0}$. Despite
these similarities, if we approximate these curves by a linear fit, the slopes
of curves in Fig. 4(a) have an increasing behavior according we increase
$\alpha$, while the slopes of curves in Fig. 4(b) decrease according to
increase $\alpha$. This is visually perceptible, but also is explicit in the
variables $\Theta_{k_{0}}$ and $\Theta_{k_{m}}$ presented in Tab. 1.

$\alpha$ | $\Theta_{k_{0}}$ | $\Theta_{k_{m}}$ | $-\beta/\nu$ | $\gamma/\nu$ | $1/\nu$
---|---|---|---|---|---
$1$ | $0.46\pm 0.04$ | $0.55\pm 0.04$ | $0.45\pm 0.04$ | $1.04\pm 0.05$ | $0.99\pm 0.05$
$2$ | $0.54\pm 0.02$ | $0.47\pm 0.01$ | $0.53\pm 0.04$ | $0.92\pm 0.06$ | $0.95\pm 0.04$
$3$ | $0.62\pm 0.04$ | $0.38\pm 0.03$ | $0.52\pm 0.04$ | $0.91\pm 0.05$ | $0.97\pm 0.06$
$4$ | $0.71\pm 0.05$ | $0.30\pm 0.03$ | $0.51\pm 0.06$ | $0.90\pm 0.07$ | $0.96\pm 0.05$
$5$ | $0.79\pm 0.08$ | $0.21\pm 0.06$ | $0.52\pm 0.07$ | $0.90\pm 0.08$ | $0.95\pm 0.04$
Table 1: Slopes of curves presented in Fig. 4(a), $\Theta_{k_{0}}$, and Fig.
4(b), $\Theta_{k_{m}}$, and the ratio between the critical exponents obtained
by the method presented in Fig. 6.

When $\alpha>3$, the analytical calculations predict well-defined finite
critical points in the $F$ to $P$ phase transition [12, 13, 17]. As a matter
of comparison, the analytical results for the random uncorrelated networks in
which is known its first and second moment of degree distribution, and
introduced by equation

$T_{c}=\frac{2}{\ln\left(\frac{\left\langle k^{2}\right\rangle}{\left\langle
k^{2}\right\rangle-2\left\langle k\right\rangle}\right)},$ (9)

also were plotted in Fig. 4 (see dashed lines). As we can see in this figure,
we also have plotted values of $\alpha<3$, and its possible because instead of
the usual integral approximation of the sum in $\left\langle k\right\rangle$
and $\left\langle k^{2}\right\rangle$, where result in $T_{c}\to 0$ and
$T_{c}\to\infty$, with $\alpha=2$ and $\alpha=3$, respectively, we have kept
the sum, once that $k_{0}$ and $k_{m}$ are restricted. The dashed lines on
Figs. 4(a) and 4(b) are described by Eq. (9), and apparently, we can observe
the same signature, both in the analytical result and simulations. However, in
the building process of our _restricted_ SFN, more connected sites are the
last to be chosen to add their connections, and their missing connections only
can be attached to sites that were not chosen yet, i.e., we have implicitly an
increase of the degree-degree correlations [22, 23], once in the last steps of
building the network, more connected sites only can connect with more
connected sites. These correlations can be identified in the higher values of
$T_{c}$ in the simulations if compared with analytical results, because as we
increase the difference between $k_{0}$ and $k_{m}$, we also increase the
number of different degrees on the network and consequently increase the
possibility of degree-degree correlations.

The degree-degree correlations were calculated here using the equation

$r=\frac{M^{-1}\sum_{i}u_{i}v_{i}-\left[M^{-1}\sum_{i}\frac{1}{2}\left(u_{i}+v_{i}\right)\right]^{2}}{M^{-1}\sum_{i}\frac{1}{2}\left(u_{i}^{2}+v_{i}^{2}\right)-\left[M^{-1}\sum_{i}\frac{1}{2}\left(u_{i}+v_{i}\right)\right]^{2}},$
(10)

for a network with $M$ edges connecting the pair of vertices, and their
respective $u_{i}$ and $v_{i}$ degrees, as defined in Ref. [22]. These
correlations are shown in Fig. 5(b) as a function of $\alpha$, and besides
being independent of network size, have a descending behavior according the
number of more connected sites increase, and ascending when this number become
small, but relevant yet. The peak of the correlation $r$ was observed in
$\alpha=3.5$ and has value $r=0.268$. Therefore, the network presents a
associative mixing, confirming that sites with high degree, on the network,
prefers to connect on more connected sites, which is the case for high values
of $\alpha$. It is interesting to mention that many social networks also have
significant associative mixing, and do not present some in complex network
models, like random graphs or growing network model (BA model [22]).

Now that we have established the critical behavior of the system as a function
of the degrees on the network, using the most distinct and distributed values
of degree, $k_{0}=4$ and $k_{m}=10$, we will effectively verify the critical
behavior of the _restricted_ SFN as a function of $\alpha$. As presented in
Fig. 5(a), this verification was first done by using some points in Fig. 4 by
simulation data, and Eq. (9). By Eq. (9) is explicit the lower values of
$T_{c}$ where we have a high degree-degree correlation, but, for small $r$,
i.e., in lower values of $\alpha$, the simulation results approach to
analytical result of a random uncorrelated network. That approach to
analytical results is also resultant of small influence of more connected
sites, and also observed for $\alpha>3.5$. By the fixed values of $k_{m}$ and
$k_{0}$, the minimum and maximum critical temperature values is limited to the
case where $k_{m}=k_{0}=4$ on Fig. 4(b) as we increase $\alpha$, because we
will never reach a network with coordination number $4$, without losing the
structure of degrees. The same is observable decreasing $\alpha$, because we
have the limit of $T_{c}$, in the case $k_{m}=k_{0}=10$ present in Fig. 4 (a)
and unreachable in our _restricted_ SFN. For some critical points in Fig. 5,
$1\leq\alpha\leq 5$, we have its explicit estimation and respective errors
presented in Tab. 2.

Analytical results for the random uncorrelated networks also extend to
critical exponents in the Ising model, being that is predicted for magnetic
critical exponent, a mean-field character, $\beta=1/2$ for $\alpha>5$, with
logarithmic corrections in $\alpha=5$, and $\beta=1/(\alpha-3)$, and for
$3<\alpha<5$. Interestingly, magnetic susceptibility critical exponent,
$\gamma$, has a mean-field universal character for $\alpha>3$, i.e.,
$\gamma=1$ for the whole values of the exponents of the degree distribution in
which a finite order phase transition is predicted, $\left\langle
k^{2}\right\rangle<\infty$. And in accordance with the scaling law relation
$\gamma/\nu=2-\eta$, is also expected a universal mean-field character for
correlation length, $\nu=1/2$, and Fisher exponent, $\eta=0$ [24]. With this
results, here, only for $1\leq\alpha\leq 5$ by limitation of the network size
and computational time, we have made a direct comparative with this analytical
results, using our _restricted_ SFN. Therefore, we have calculated the
critical exponents for the system, and this was made by two methods. One of
these methods, using the FSS relations, is based on the slope of its linear
fit, with the thermodynamic quantities near to critical point, and the second
one, is based on the data collapse of thermodynamic quantities in the form of
scaling functions [18, 19, 20].

Figure 6: Log-log plots of (a) the magnetization $\textrm{m}_{\textrm{L}}$,
(b) the susceptibility$\textrm{$\chi$}_{\textrm{L}}$, and (c) the derivative
of the cumulant$\textrm{U}_{\textrm{L}}^{\prime}$, at the critical point, as a
function of the effective length $L$. We have fixed the values of $k_{0}=4$,
$k_{m}=10$, and some selected values $\alpha$, as indicated in the figure. As
we are interest is only on the slope, the linear coefficients are changed for
better visualization of each of the curves. The slopes obtained here can be
seen in Tab. 1.

We can find the critical exponents of the system by using two methods. In the
first method, we have used the data of thermodynamic quantities near the
critical point and for different network sizes. For instance, in a log-log
plot in Fig. 6, we have fited the thermodynamic quantities as a function of
the effective length of the system $L$ and its slope returns us the ratios
between the critical exponents present in Eqs. (6), (7) and (8). For the ratio
$-\beta/\nu$ present in Eq. (6), the curves of magnetization for different
network sizes and values $\alpha$ is presented in Fig. 6(a), where we shown
its linear fits. In the same way, but now for the ratio $\gamma/\nu$, using
the Eq. (7), the linear fit of magnetic susceptibility curves can be seen in
Fig. 6(b). Only these two ratios does not give us all the critical exponents
of the system. Therefore, we also have used the forth-order Binder cumulant
curves and its derivative, Eq. (8), being that its linear fits is presented in
Fig. 6(c), and give us the ratio $1/\nu$ that is the inverse of correlation
length exponent, and only then obtain the exponents $\beta$ and $\gamma$. The
set of exponents obtained by this method can be found in Tab. 1. However, it
is necessary to emphasize that, because we are dealing with a random network
with modifiable degree distribution, it is required to modify the Eqs. (6),
(7) and (8) to the mean-field scaling relations, in which we change the
effective length $L$ by the total number of spins on network, $N$. This
modification in the scaling relations can be easily implemented in the
obtained critical exponents, dividing $\nu$ by two, once that $N=L^{2}$.

Figure 7: Data collapse near the critical point for the magnetization $\textrm{m}_{\textrm{L}}$ (a) and susceptibility $\textrm{$\chi$}_{\textrm{L}}$ (b) for different network sizes, as shown in the figure. Here, the minimum and maximum degree are fixed $k_{0}=4$, $k_{m}=10$, and $\alpha=1$. The log-log plots were used to obtain the slope $\Theta$ of scaling functions asymptotic behavior, i.e., distant of $\epsilon=0$. The straight-dashed lines represent the asymptotic behavior of the scaling functions, Eq. (6) and (7). $\alpha$ | $T_{c}$ | $\beta$ | $\gamma$ | $\nu_{\textrm{m}_{\textrm{L}}}$ | $\nu_{\textrm{$\chi$}_{\textrm{L}}}$
---|---|---|---|---|---
$1$ | $6.228\pm 0.004$ | $0.47\pm 0.03$ | $1.03\pm 0.03$ | $1.00\pm 0.06$ | $1.00\pm 0.04$
$2$ | $5.734\pm 0.008$ | $0.55\pm 0.02$ | $0.95\pm 0.03$ | $1.03\pm 0.03$ | $1.05\pm 0.04$
$3$ | $5.225\pm 0.009$ | $0.57\pm 0.03$ | $0.91\pm 0.05$ | $0.96\pm 0.05$ | $1.08\pm 0.03$
$4$ | $4.710\pm 0.020$ | $0.56\pm 0.03$ | $0.90\pm 0.04$ | $1.05\pm 0.05$ | $1.10\pm 0.03$
$5$ | $4.210\pm 0.020$ | $0.58\pm 0.03$ | $0.92\pm 0.05$ | $1.03\pm 0.04$ | $1.10\pm 0.05$
Table 2: As a function of $\alpha$, here is showed critical temperature
$T_{c}$ present in the phases diagram of Fig. 5, and critical exponents
obtained by the data collapse method, being that
$\nu_{\textrm{m}_{\textrm{L}}}$ and $\nu_{\textrm{$\chi$}_{\textrm{L}}}$ are
the correlation length exponent obtained by the magnetic and magnetic
susceptibility data collapsed curves, respectively.

In the second method, we can improve the values found for the critical
exponents, by collapsing the data points. Our main goal is to use the
thermodynamic quantities curves and different network sizes to obtain the form
of its scaling functions, as a function of $L^{1/\nu}\epsilon$. It is possible
because, in the proximity of the critical points, scaling functions in Eqs.
(6), (7) and (8) must be independently of network sizes, if on its, is
utilized the correct critical exponents of the system, i.e., in the proximity
of $T_{c}$ using the correct exponents on scaling relations we obtain a
collapsed curve in the form of scaling functions. In Fig. 7, we show some
examples of the best collapsed curves obtained. In this figure, what we have
done is adjust the exponents on isolated scaling functions, $m_{0}$ and
$\chi_{0}$, as a function of $L^{1/\nu}\varepsilon$, and when the
thermodynamic quantity with different network size, best collapses into a
single curve, that exponents were considered the critical exponents of the
system. Fig. 7(a) shown, for $\alpha=1$, the collapsed curves of magnetization
data, Eq. (6), in which allowed us to adjust and obtain the exponents $\beta$
and $\nu$, and besides that, in Fig. 7(b) also for $\alpha=1$, magnetic
susceptibility data, Eq. (7), enable us to obtain exponents $\gamma$ and
$\nu$. That estimated critical exponents for $1\leq\alpha\leq 5$ is presented
in Tab. 2, but, we also have to take into account the modification in scaling
functions for the mean-field character, and the real values of exponent $\nu$,
is only obtained dividing its results by two.

Figure 8: Static critical exponents $\beta$, $\gamma$, and $\nu$ as functions
of the exponent $\alpha$.

To the best comprehension of the critical exponents obtained for the
_restricted_ SFN, the average of its values, which is equivalent in both
methods utilized in this work, was plotted in Fig. 8 as a function of
$\alpha$. In lower values of $\alpha$, we tend to a system with all the
degrees having the same number of sites and the mean-field critical behavior
is present, but, when $\alpha$ increase, the number of more connected sites
decreases, and a slight deviation from this behavior is observed. With that
observations, we can say that degree-degree correlations, also causes a
deviation from the expected mean-field behavior on random networks with degree
distribution convergent moments [24, 25, 26].

##  V Conclusions

Here, we have employed Monte Carlo simulations to the study of thermodynamic
quantities and the critical behavior of the Ising model on a _restricted_ SFN.
When we fix the maximum $k_{m}$ and minimum $k_{0}$ number of degrees for the
whole network sizes, as we made on our _restricted_ SFN, we always have
convergent moments based in the degree distribution $P(k)$. We have used a
power-law degree distribution, and as the analytical results predict, the
convergent second and fourth moments in the arbitrary degree distributions
present a finite order phase transition [12, 13, 17, 24]. We have obtained the
critical points of second-order phase transitions and built phase diagram of
temperature $T$ as a function of $k_{0}$ and $k_{m}$, and in which when we
have $k_{0}=k_{m}$ a random uncorrelated network is obtained. Therefore, the
critical points are in agreement with analytical calculations, but, the
increase on difference between $k_{0}$ and $k_{m}$, also increase the degree-
degree correlations and causes a deviation from that calculations. The phase
diagram of temperature $T$ as a function of $\alpha$ also was built, where we
always have a finite critical temperature and a decrease in the critical point
as we increase $\alpha$, once that we also decrease the number of more
connected sites on network. With these critical points, we estimated the
critical exponents for the system as a function of $\alpha$, and different
from what is predicted by analytical results in a random uncorrelated SFN,
here, with lower values of $\alpha$, the increase in the number of more
connected sites reaches a mean-field critical behavior. Otherwise, when
$\alpha$ increase, a mean-field critical behavior is also observed, but how
are we dealing with correlated degrees on network, the decreasing in the
number of more connected sites causes a slight deviation on that exponents.
That mean-field behavior was predicted on networks with convergent moments on
its degree distribution and it was found on a diversity of complex networks
[24, 25, 26]. In our work, it was made on the Ising model with a power-law
degree distributed network and unexpected values of $\alpha$.

## References

  * [1] R. Albert and A.-L. Barabási. Rev, Mod. Phys., 74, 47 (2002);
  * [2] A.-L. Barabási. _Network Science_ , (Cambridge University Press, NewYork, EUA, 2016);
  * [3] A.-L. Barabási and R. Albert. Science, 286, 509 (1999);
  * [4] M. Krasnytska, B. Berche, Y. Holovatch and R. Kenna. Entropy, 23, 1175 (2021);
  * [5] D.-H. Kim, G. J. Rodgers, B. Kahng and D. Kim. Phys. Rev. E, 71, 056115 (2005);
  * [6] S. N. Dorogovtsev, A. V. Goltsev and J. F. F. Mendes. Eur. Phys. J. B, 38, 177 (2004);
  * [7] G. Bianconi and A.-L. Barabási. Phys. Rev. Lett., 86, 5632 (2001);
  * [8] C. P. Herrero. Phys. Rev. E, 91, 052812 (2015);
  * [9] S. Aparicio, J. Villazón-Terrazas and G. Álvarez. Entropy, 17, 5848 (2015);
  * [10] A. L. M. Vilela, B. J. Zubillaga, C. Wang, M. Wang, R. Du and H. E. Stanley. Sci. Rep., 10, 8255 (2020);
  * [11] T. Gradowski and A. Krawiecki. A. Phys. Pol. A, 127, 1540 (2015);
  * [12] S. N. Dorogovtsev, A. V. Goltsev and J. F. F. Mendes. Phys. Rev. E, 66, 016104 (2002);
  * [13] A. V. Goltsev, S. N. Dorogovtsev and J. F. F. Mendes. Phys. Rev. E, 67, 026123 (2003);
  * [14] A. Aleksiejuk, J. A. Holyst and D. Stauffer. Physica A, 310, 260 (2002);
  * [15] G. Bianconi. Phys. Lett. A, 303, 166 (2002);
  * [16] C. P. Herrero. Phys. Rev. E, 69, 067109 (2004);
  * [17] M. Leone, A. Vazquez, A. Vespignani and R. Zecchina. Eur. Phys. J. B, 28, 191 (2002);
  * [18] K. Binder and D. W. Heermann. Monte Carlo Simulation in Statistical Physics. An Introduction, 6rd ed. (Springer, Cham, Switzerland, 2019);
  * [19] K. Binder and D. P. Landau. A Guide to Monte Carlo Simulations in Statistical Physics, 4rd ed. (TJ International Ltd, Padstow, UK, 2015);
  * [20] L. Böttcher and H. J. Herrmann. Computational Statistical Physics, 1rd ed. (Cambridge University Press, NewYork, EUA, 2021);
  * [21] S.-H. Tsai and S. R. Salinas. Braz. J. Phys., 28, 1, (1998);
  * [22] M. E. J. Newman. Phys. Rev. Lett., 89, 208701 (2002);
  * [23] M. Catanzaro, M. Boguña and R. Pastor-Satorras. Phys. Rev. E, 71, 027103 (2005);
  * [24] S. N. Dorogovtsev, A. V. Goltsev and J. F. F. Mendes. Rev, Mod. Phys., 80, 1275 (2008);
  * [25] R. A. Dumer and M. Godoy. Eur. Phys. J. B 95, 159 (2022);
  * [26] A. Barrat and M. Weigt. Eur. Phys. J. B, 13, 547 (2000).


# A Physics-Informed Deep Learning Model of the Hot Tail Runaway Electron Seed

Christopher J. McDevitt  Nuclear Engineering Program, University of Florida

###### Abstract

A challenging aspect of the description of a tokamak disruption is evaluating
the hot tail runaway electron (RE) seed that emerges during the thermal
quench. This problem is made challenging due to the requirement of describing
a strongly non-thermal electron distribution, together with the need to
incorporate a diverse range of multiphysics processes including
magnetohydrodynamic instabilities, impurity transport, and radiative losses.
The present work develops a physics-informed neural network (PINN) tailored to
the solution of the hot tail seed during an idealized axisymmetric thermal
quench. Here, a PINN is developed to identify solutions to the adjoint
relativistic Fokker-Planck equation in the presence of a rapid quench of the
plasma’s thermal energy. It is shown that the PINN is able to accurately
predict the hot tail seed across a range of parameters including the thermal
quench time scale, initial plasma temperature, and local current density, in
the absence of experimental or simulation data. The hot tail PINN is verified
by comparison with a direct Monte Carlo solution, with excellent agreement
found across a broad range of thermal quench conditions.

##  I Introduction

Runaway electrons (REs) are formed due to the force exerted by a DC electric
field overcoming the Coulomb drag, thus enabling the acceleration of the
electron to relativistic energies Dwyer and Uman (2014); Breizman et al.
(2019). Mechanisms through which REs may be generated can be separated into
two broad categories. The first such category corresponds to a _nonlocal_
process whereby an initially low energy electron makes a discrete jump in
energy due to a large-angle collision with a pre-existing high energy
electron. If the energy exchanged during this collision is sufficient to
scatter the initially low energy electron to high enough energy that it is
accelerated by the DC electric field, this will enable the pre-existing RE
population to grow. This mechanism corresponds to the avalanche mechanism of
RE generation SOKOLOV (1979).

A fundamentally distinct means through which REs can be generated corresponds
to the _local_ acceleration of electrons from the bulk electron population via
a strong DC electric field. Such a local means of RE generation often provides
the ‘seed’ population of REs, which is subsequently amplified by the avalanche
mechanism described above. Perhaps the most well studied example of a local
process of RE generation corresponds to Dreicer acceleration Dreicer (1959).
Here, an electric field accelerates electrons on the tail of the electron
distribution with energy diffusion providing a means of filling in the regions
of momentum space vacated by the accelerated electrons, thus providing a
continuous source of REs. While the Dreicer mechanism has been extensively
studied Connor and Hastie (1975); Hesslow et al. (2019); McDevitt and Tang
(2019); Chacon et al. (2022), it is often subdominant for the high plasma
densities present during a tokamak disruption Martín-Solís et al. (2017);
Aleynikov and Breizman (2017). A potentially larger seed mechanism active
during a tokamak disruption corresponds to the hot tail mechanism of RE
generation. This mechanism arises during a rapid collapse of the electron
temperature, such as the thermal quench of a tokamak disruption. During such a
rapid quench of the electron’s thermal energy, which may be induced by the
injection of impurities into the plasma for example, bulk electrons rapidly
cool due to their relatively high collisionality. Energetic electrons,
however, will cool on a far slower timescale thus resulting in the formation
of a strongly non-thermal ‘hot tail’ distribution of electrons Helander et al.
(2004); Smith and Verwichte (2008); Stahl et al. (2016); Matsuyama and Yagi
(2017). As the bulk plasma cools, this leads to a rise in the plasma’s
resistivity, resulting in the sudden growth in the strength of the electric
field, and the acceleration of a fraction of the ‘hot tail’ electron
distribution to relativistic energies.

Despite the importance of hot tail generation to fusion plasmas, an accurate
treatment of this process has resisted theoretical description to date for two
fundamental reasons. The first is that hot tail RE seed generation depends
sensitively on the time history of the temperature profile, and thus cannot be
described as a function of local parameters. While direct numerical
simulations can be performed to simulate any specific temperature history, no
general theory for this process is presently available. A second fundamental
challenge is that part of what induces the loss of thermal energy during a
disruption is the rapid transport of heat along open magnetic field lines.
High energy electrons (i.e the hot tail seed) will also be susceptible to loss
along such open magnetic field lines. Hence, an accurate treatment of hot tail
generation requires incorporating both the specific time history of the
electron temperature, along with the detailed magnetohydrodynamic (MHD)
spectrum.

The present paper seeks to address the first of these challenges by providing
a description of hot tail formation in the limit of an axisymmetric plasma.
Since a fraction of the hot tail seed will likely be lost along open magnetic
field lines, the present analysis should be viewed as an upper bound on the
hot tail seed. Future work will seek to address losses via spatial transport.
To carry out this analysis we will employ a deep learning framework with the
aim of learning the parametric dependence of the hot tail seed. Noting that
even a very simple model of a thermal quench implies a high dimensional
parameter space, a purely data driven approach will require the use of a
prohibitively large quantity of data to densely sample the relevant regions of
parameter space. In order to overcome this limitation, we will seek to utilize
physical constraints to eliminate the need to generate data to train a neural
network (NN). Since such physical constraints can be sampled at an arbitrary
number of points in the parameter space of the model, this will allow for the
solution to be constrained even in a high dimensional parameter space.

The present work will develop a physics-informed neural network (PINN) Raissi
et al. (2019); Karniadakis et al. (2021) tailored to the solution of the
adjoint relativistic Fokker-Planck equation, in the presence of an evolving
electric field. While PINNs allow for the incorporation of both data and
physical constraints in the training of a NN, in the present analysis we will
focus on the zero data limit, such that we will rely exclusively on physical
constraints to identify solutions of the adjoint Fokker-Planck equation. Our
aim will thus be to learn a NN based representation of the hot tail mechanism,
which we will refer to as a hot tail PINN. The idealized model employed will
involve the relativistic adjoint Fokker-Planck equation coupled to an electric
field evaluated from Ohm’s law during a prescribed temperature history. The
hot tail seed is subsequently learned across a broad range of thermal quench
times $\Delta t_{TQ}$, initial plasma temperatures $T_{init}$, and current
densities $j_{0}$, yielding a rapid surrogate model of the hot tail seed. As a
means of verifying the accuracy of the hot tail PINN’s predictions, a direct
Monte Carlo solution for the runaway probability function (RPF) Karney and
Fisch (1986) will be carried out using the RE solver RAMc McDevitt et al.
(2019).

The remainder of this paper is organized as follows. In Sec. II the general
PINN framework employed in this work is briefly reviewed. Section III
describes the idealized thermal quench model employed, along with
customizations of the PINN framework to the specific case of the adjoint
relativistic Fokker-Planck equation. The Monte Carlo code used for
verification is described in Sec. IV. Section V provides an overview of the
predictions of the hot tail PINN across a broad range of thermal quench
scenarios, with conclusions provided in Sec. VI.

##  II Physics-Constrained Deep Learning

A primary aim of the present paper is the development of a PINN customized to
treat the adjoint relativistic Fokker-Planck equation. Before doing so, it
will be useful to briefly review fundamental aspects of the PINN framework,
which has emerged as a prominent example of physics-informed deep learning
methods Lagaris et al. (1998); Karpatne et al. (2017); Karniadakis et al.
(2021); Lusch et al. (2018); Wang et al. (2020). The present discussion will
only focus on the essential concepts, where the interested reader is referred
to Ref. Karniadakis et al. (2021) for a more detailed discussion. Here, the
underlying strategy is to use physical constraints (primarily PDEs) to
regularize the training of a neural network. Such an approach not only
provides a natural means of avoiding overfitting, thus providing a more robust
interpolative tool, but also opens up the possibility of greater
generalizability to unseen parameter regimes. A PINN in its simplest form can
be expressed as Raissi et al. (2019); Karniadakis et al. (2021):

$\text{Loss}=\frac{1}{N_{data}}\sum^{N_{data}}_{i}\left[P_{i}-P\left(x_{i},t_{i},\lambda_{i}\right)\right]^{2}+\frac{1}{N_{PDE}}\sum^{N_{PDE}}_{i}\mathcal{R}^{2}\left(x_{i},t_{i},\lambda_{i}\right),$
(1)

where $P_{i}$ represent data points for the quantity of interest (the RPF in
the present paper), $x_{i}$ is a phase space coordinate (pitch or momentum,
for example), $t_{i}$ is time, and $\lambda_{i}$ represent parameters of the
physical system such as the time scale of the thermal quench or current
density, $N_{data}$ is the number of data points used in the training of the
NN, $N_{PDE}$ is the number of points where the physical constraint is
sampled, and $\mathcal{R}\left(x_{i}\right)$ is the residual of the PDE, which
penalizes deviations from the imposed physical constraint of the system. The
first term in Eq. (1) represents loss against available data, whereas the
second term enforces the physical constraints. We note that while the number
of data points $N_{data}$ will be sharply limited for many challenging
problems related to tokamak disruptions, the value of $N_{PDE}$ may be taken
to be very large, where increasing its value will only impact the
computational cost of the offline training of the model.

For the present study we will be interested in the limit where no data is
available. Our motivation will be to demonstrate the ability of PINNs to
accurately learn solutions to the adjoint Fokker-Planck equation across a
broad range of parameters $\lambda$ in this data free limit. As described in
Sec. III.3 below, the use of a vanilla PINN such as Eq. (1) will fail to
accurately evaluate the RPF in the absence of data during the highly transient
conditions that are characteristic of a thermal quench. A primary motivation
of the present work will thus be to develop a tailored hot tail PINN, with a
loss function and neural network architecture that enables an accurate
solution of the RPF across a broad range of thermal quench conditions.

##  III Reduced Model of the Thermal Quench

###  III.1 Adjoint Fokker-Planck Equation

The momentum space evolution of a relativistic population of electrons is
described by the relativistic Fokker-Planck equation, i.e.

$\displaystyle\frac{\partial f_{e}}{\partial t}$
$\displaystyle+\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left[p^{2}\left(-E_{\|}\xi-\frac{1+p^{2}}{p^{2}}\right)f_{e}\right]-\frac{\partial}{\partial\xi}\left[\left(\frac{1-\xi^{2}}{p}\right)E_{\|}f_{e}\right]$
$\displaystyle=\frac{\nu_{D}}{2}\frac{\partial}{\partial\xi}\left[\left(1-\xi^{2}\right)\frac{\partial
f_{e}}{\partial\xi}\right]+\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left(p^{2}C_{A}\frac{\partial f_{e}}{\partial p}\right),$ (2)

where the collisional coefficients are taken to be

$\nu_{D}=\left(1+Z_{eff}\right)\frac{\gamma}{p^{3}},\quad
C_{A}=\left(\frac{T_{e}}{m_{e}c^{2}}\right)\frac{\gamma^{3}}{p^{3}}.$

Here, the relativistic momentum $p$ is normalized as $p\to
p/\left(m_{e}c\right)$, the electron’s pitch is defined by $\xi\equiv
p_{\|}/p$, time is normalized as $t\to t/\tau_{c}$, where $\tau_{c}\equiv
4\pi\epsilon^{2}_{0}m^{2}_{e}c^{3}/\left(e^{4}n_{e}\ln\Lambda\right)$ is the
collision time of a relativistic electron, the collisional coefficients
$\nu_{D}$ and $C_{A}$ are normalized to $\tau_{c}$, and we have normalized the
parallel electric field to the Connor-Hastie electric field $E_{\|}\to
E_{\|}/E_{c}$, where $E_{c}\equiv m_{e}c/\left(e\tau_{c}\right)$ Connor and
Hastie (1975). Applying a PINN directly to this equation is challenging due to
the magnitude of $f_{e}$ varying strongly with energy, where our interest will
be capturing the tail of the electron distribution which is typically orders
of magnitude smaller than the bulk plasma.

An alternate approach involves the solution of the adjoint equation Liu et al.
(2016); Zhang and del Castillo-Negrete (2017)

$\displaystyle\frac{\partial P\left(t\right)}{\partial t}$
$\displaystyle+\left[-E_{\|}\left(t\right)\xi-\frac{1+p^{2}}{p^{2}}\right]\frac{\partial
P\left(t\right)}{\partial
p}-\left(\frac{1-\xi^{2}}{p}\right)E_{\|}\left(t\right)\frac{\partial
P\left(t\right)}{\partial\xi}$
$\displaystyle=-\frac{\nu_{D}}{2}\frac{\partial}{\partial\xi}\left[\left(1-\xi^{2}\right)\frac{\partial
P\left(t\right)}{\partial\xi}\right]-\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left[p^{2}C_{A}\frac{\partial P\left(t\right)}{\partial p}\right],$ (3)

subject to the terminal condition
$P\left(t=t_{final}\right)=P_{RE}\left(p\right)$, where
$P_{RE}\left(p\right)=1$ for $p\geq p_{RE}$ and $P_{RE}\left(p\right)=0$
otherwise. The RPF $P\left(t\right)$ represents the probability that an
electron located at $\left(p,\xi\right)$ at time $t$ obtains an energy greater
than $p_{RE}$ on, or before $t_{final}$ Karney and Fisch (1986). Since $P$ is
a probability its value will range between zero and one, and thus will be far
easier to represent with a PINN compared to the electron distribution $f_{e}$.
Equation (3) can be rewritten as an initial value problem by making the
variable substitution $\tau\equiv t_{final}-t$, i.e.

$\displaystyle\frac{\partial P^{\prime}\left(\tau\right)}{\partial\tau}$
$\displaystyle-\left[-E_{\|}\left(t_{final}-\tau\right)\xi-\frac{1+p^{2}}{p^{2}}\right]\frac{\partial
P^{\prime}\left(t\right)}{\partial
p}+\left(\frac{1-\xi^{2}}{p}\right)E_{\|}\left(t_{final}-\tau\right)\frac{\partial
P^{\prime}\left(\tau\right)}{\partial\xi}$
$\displaystyle=\frac{\nu_{D}}{2}\frac{\partial}{\partial\xi}\left[\left(1-\xi^{2}\right)\frac{\partial
P^{\prime}\left(\tau\right)}{\partial\xi}\right]+\frac{1}{p^{2}}\frac{\partial}{\partial
p}\left[p^{2}C_{A}\frac{\partial P^{\prime}\left(\tau\right)}{\partial
p}\right],$ (4)

where $P^{\prime}\left(\tau\right)$ is required to satisfy the initial
condition $P^{\prime}\left(0\right)=P_{RE}\left(p\right)$. By making this
variable change it is equivalent to beginning with a distribution of electrons
above the RE energy $p_{RE}$, and then integrating the equations of motion
backward in time to identify the regions of momentum space with orbits that
lead to RE formation.

###  III.2 Axisymmetric model of the Thermal Quench

To identify the hot tail seed, we will solve Eq. (4) in the presence of an
evolving temperature and electric field. Specifically, the electric field will
be evolved by Ohm’s law, under the assumption of constant current density
$j_{0}$, i.e.

$E_{\|}=\eta j_{\|}=\eta j_{0},$ (5)

where we have used the normalizations $E_{\|}\to E_{\|}/E_{c}$, $j_{0}\to
a^{2}j_{0}/I_{A}$ and $\eta\to I_{A}\eta/\left(a^{2}E_{c}\right)$. For a
Spitzer resistivity $\eta\propto 1/T^{3/2}_{e}$, neglecting the logarithmic
dependence of the Coulomb logarithm, the electric field can then be linked to
the temperature history as

$E_{\|}=\left(\frac{T_{init}}{T_{e}}\right)^{3/2}\eta_{init}j_{0}=\left(\frac{T_{init}}{T_{e}}\right)^{3/2}E_{init},$
(6)

where $\eta_{init}$ is the resistivity evaluated at the initial temperature
$T_{init}$, and $E_{init}\equiv\eta_{0}j_{0}$ is the initial electric field.
The temperature history will be parameterized by

$T_{e}\left(t\right)=\left(T_{init}-T_{final}\right)\exp\left(\frac{-t}{\Delta
t_{TQ}}\right)+T_{final}.$ (7)

Here, $T_{final}$ corresponds to the post-thermal quench temperature and
$\Delta t_{TQ}$ indicates the thermal quench time scale normalized to
$\tau_{c}$. For simplicity, we will take the plasma to be a fully ionized
deuterium plasma with constant density $n_{D}=n_{e0}$. The Coulomb logarithm
used to define the time normalization $\tau_{c}$ will be evaluated using the
initial electron temperature $T_{init}$ and electron density $n_{e0}$.

The present model is characterized by four physics input parameters
$\left(\Delta t_{TQ},T_{init},T_{final},E_{init}\right)$. In addition to these
four physics parameters, we will also need to specify an end time $t_{final}$,
the energy range of the simulation (i.e. $p_{min}$ and $p_{max}$), and the RE
energy $p_{RE}$. The solution can be shown to be insensitive to these
parameters so long as they are chosen to be sufficiently large (or small for
the case of $p_{min}$).

###  III.3 Hot Tail PINN

A property of NNs that has led to their widespread use across a range of
applications is the impressive expressivity of NNs Hornik et al. (1989). When
applying NNs to scientific problems, however, it is often preferable to
restrict the range of functions that can represented by the NN to those that
are consistent with the physical problem of interest. In so doing, the space
of available solutions is drastically reduced leading to a far more robust
training process Lu et al. (2021a). Our aim in this subsection will be to
impose as hard constraints in the NN architecture several properties of the
RPF. Specifically, we will restrict the range of values of the RPF to be
between zero and one, enforce the initial condition to be
$P^{\prime}\left(t=0\right)=P_{RE}\left(p\right)$, and enforce the boundary
condition $P_{RE}\left(p=p_{min}\right)=0$ as hard constraints. In addition,
we will enforce the condition
$P^{\prime}\left(p=p_{max},\xi=\xi_{min}\right)=1$, where $\xi_{min}=-1$.
While this condition is not strictly true for all possible thermal quench
scenarios, will choose $p_{max}$ to be sufficiently large such that it will be
satisfied for all cases considered in this paper. In practice, forcing
$P^{\prime}$ to be one at this location helps ensure robust convergence of the
PINN.

The means through which we will enforce these conditions will be by
introducing additional layers to the NN. Specifically, defining the output of
the NN by $P_{NN}$, we will introduce a layer of the form

$\bar{P}^{\prime}\left(p,\xi,t;\mathbf{\lambda}\right)=P_{RE}\left(p\right)+\left(\frac{p-p_{min}}{p_{max}-p_{min}}\right)\left(\frac{t-t_{min}}{t_{max}-t_{min}}\right)\left(\frac{p_{max}-p\Xi\left(\xi\right)}{p_{max}-p_{min}}\right)P_{NN}\left(p,\xi,t;\mathbf{\lambda}\right).$
(8)

To ensure the RPF has a range between zero and one we will pass this the final
layer defined by

$P^{\prime}\left(p,\xi,t;\mathbf{\lambda}\right)=\tanh\left(\frac{{\bar{P}^{\prime}\left(p,\xi,t,\mathbf{\lambda}\right)}^{2}}{\Delta_{prob}^{2}}\right),$
(9)

where $\mathbf{\lambda}$ indicates the parameters of the hot tail model
$\left(\Delta t_{TQ},T_{init},T_{final},j_{0}\right)$, we will choose
$\Delta^{2}_{prob}\ll 1$, and we have defined the function

$\Xi\left(\xi\right)\equiv\exp\left(-\frac{\left(\xi-\xi_{min}\right)^{2}}{\Delta\xi^{2}}\right).$
(10)

These output layers can be verified to constrain $P^{\prime}$ to have a range
between zero and one, forces the solution to vanish at $p=p_{min}$, exactly
enforces the initial condition $P^{\prime}=P_{RE}$, and satisfies the property
$P^{\prime}\left(p=p_{max},\xi=\xi_{min}\right)=1$ which is true for
sufficiently large $p_{max}$.

Rather than employing a step function for the initial condition
$P_{RE}\left(p\right)$, we will introduce a modest smoothing of this initial
condition as follows

$P_{RE}\left(p\right)=\frac{1}{2}\left[1-\tanh\left(\frac{p_{RE}-p}{dp}\right)\right].$
(11)

Noting that the PINN will be learning deviations from $P_{RE}$ for $t>0$, by
using a smooth form of $P_{RE}$ the PINN does not need to learn a sharp
discontinuity.

For all of the cases considered below we will take the various parameters
characterizing these additional output layers to be fixed. Specifically, we
will take $\Delta\xi=0.5$, $dp=0.1p_{max}$, $\Delta_{prob}=0.25$ and
$p_{RE}=0.8p_{max}$. While the value for $dp$ may seem somewhat large, since
$P_{RE}$ will be squared in Eq. (11), and then evaluated inside a $\tanh$,
this choice of $dp$ can be verified to lead to a sharp transition between zero
and one. For all cases considered in this work we will take $p_{max}\approx
1.7$, or a maximum energy of $500\;\text{keV}$. With regard to the NN, we will
employ a simple feedforward fully connected NN, with six hidden layers and 64
neurons per hidden layer. Finally we note that we have made use of the DeepXDE
library Lu et al. (2021b) for the studies carried out below, where all source
code will be made available through GitHub after acceptance of this article.

##  IV Verification of Hot Tail PINN: Monte Carlo Runaway Electron Simulations

Figure 1: Example time history of the electron temperature [panel (a)] and
electric field [panel (b)] used in the RAMc simulations. The thermal quench
time was taken to be $\Delta t_{TQ}=0.5\;\text{ms}$,
$j_{0}=1\;\text{MA}/\text{m}^{2}$, and the initial and final temperatures are
given by $T_{init}=15\;\text{keV}$ and $T_{final}=30\;\text{eV}$,
respectively.

The accuracy of the hot tail PINN will be verified by performing direct
simulations of hot tail formation using the RunAway Monte carlo (RAMc) code.
This particle based RE solver evaluates the guiding center motion of
relativistic electrons including both small and large-angle collisions
McDevitt et al. (2019). The inductive electric field within RAMc is updated
self-consistently by a flux diffusion equation accounting for the current
carried by REs McDevitt and Tang (2022). In order to provide direct
comparisons with predictions of the hot tail PINN and RAMc, several
simplifications to the RAMc solver will be introduced. In particular, the
collisional coefficients and the flux diffusion equation used in RAMc have
been modified to be consistent with Secs. III.1 and III.2 above. Specifically,
the large-angle collision operator has been disabled such that we will only
describe seed RE formation. Also, the coupling of the RE current back to the
flux diffusion equation has also been disabled to better approximate the
simplified Ohm’s law used in the hot tail PINN model [i.e. Eq. (6)]. While the
plasma current will evolve slightly during the RAMc simulations, we will use a
relatively high post thermal quench plasma temperature
($T_{e}=30\;\text{eV}$), which implies a modest resistivity, to ensure the
plasma current density undergoes little change during the short timescale of
the thermal quench. Furthermore, the marker particles are all initialized at
$r\approx 0$ such that the impact of tokamak geometry on RE formation will be
negligible. This enables comparisons between this toroidal solver and the slab
RE model described by Eq. (4). A large ITER-like device with a minor radius of
$a=200\;\text{cm}$, inverse aspect ratio of $\epsilon=1/3$, and an on-axis
magnetic field of $B_{0}=5.3\;\text{T}$ is selected when performing the RAMc
simulations to ensure electrons remain well confined to their initial flux
surfaces. In addition, since low energy particles have a vanishingly small
probability of running away, an absorbing boundary will be implemented at an
energy of $100\;\text{eV}$.

Figure 2: Particle distributions in pitch and energy at $t=0$ [panel (a)],
$t\approx 3.80\;\text{ms}$ [panel (b)], and $t\approx 4.94\;\text{ms}$ [panel
(c)]. The RPF is shown in panel (d). The thermal quench time was taken to be
$\Delta t_{TQ}=0.5\;\text{ms}$, $j_{0}=1\;\text{MA}/\text{m}^{2}$, and the
initial and final temperatures are given by $T_{init}=15\;\text{keV}$ and
$T_{final}=30\;\text{eV}$, respectively.

With the above simplifications to the RAMc code, direct Monte Carlo
calculations of the RPF will be employed to verify the predictions of the hot
tail PINN. This is done by initializing marker particles according to a
uniform random distribution in momentum space at $t=0$. All electrons will
have an initial spatial location taken to be $r\approx 0$. By then evaluating
which electrons are able to achieve an energy of $p_{RE}$ or greater after
$t=t_{final}$, this allows $P^{\prime}\left(\tau=t_{final}\right)$ to be
directly evaluated from the particle data. An example solution is shown in
Figs. 1 and 2. Here, 100,000 marker particles are used to evaluate the RPF,
where for this example the parameters were assumed to be $\Delta
t_{TQ}=0.5\;\text{ms}$, $T_{init}=15\;\text{keV}$, $T_{final}=30\;\text{eV}$,
and $E_{0}$ is evaluated for a fully ionized deuterium plasma with
$n_{D}=10^{14}\;\text{cm}^{-3}$ and $j_{0}=1\;\text{MA/m}^{2}$. The time
history of the on-axis electric field and temperature are shown in Fig. 1,
with the resulting RPF and particle distribution shown in Fig. 2. Here, the
rapid drop in the temperature leads to the formation of a large electric field
with a maximum value of roughly $60$ times the Connor-Hastie threshold. As
evident from the particle distributions, after roughly ten thermal quench
times ($\approx 5\;\text{ms}$), all of the particles have either been absorbed
by the low energy boundary, or been accelerated past the RE energy $p_{RE}$ of
$\approx 300\;\text{keV}$. Noting that the particles that run away have an
energy well in excess of $p_{RE}$, implies that there is little sensitivity to
the specific choice of $p_{RE}$, as long as the RE energy is below an MeV.
From the form of the RPF [Fig. 2(d)], it is apparent that electrons whose
initial energy is less than $\sim 100\;\text{keV}$ have a negligible chance of
running away, implying that the dominant contribution to the RE seed will be
from electrons with energies far in excess of the thermal energy. The
magnitude of the hot tail seed can be evaluated from this RPF by evaluating
the integral

$n_{RE}=\int d^{3}pf^{MJ}_{e}P\left(p,\xi\right),$ (12)

where $P\left(p,\xi\right)$ is the RPF indicated in Fig. 2(d), and
$f^{MJ}_{e}$ is a Maxwell-Jüttner distribution evaluated using the initial
electron temperature $T_{init}$ and electron density $n_{e0}$. Physically,
$P\left(p,\xi\right)$ indicates the probability that an electron initially at
a given momentum $\left(p,\xi\right)$ is able to run away, thus by integrating
with the initial electron distribution we will be able to compute the number
of REs. While Monte Carlo solutions are computationally intensive, a small
number of RPFs of the form shown in Fig. 2(d) will be used to verify the hot
tail PINN solution. We emphasize that this data set will be used for
verification, and not included in the training of the hot tail PINN.

##  V Predictions of the Hot Tail PINN

###  V.1 Runaway Probability Function Evolution During the Thermal Quench

Figure 3: Example time history of the electron temperature [panel (a)] and
electric field [panel (b)]. The thermal quench time was taken to be $\Delta
t_{TQ}=0.5\;\text{ms}$, $j_{0}=1\;\text{MA}/\text{m}^{2}$, and the initial and
final temperatures are given by $T_{init}=15\;\text{keV}$ and
$T_{final}=30\;\text{eV}$, respectively.

We will be interested in evaluating the RPF from Eq. (4) for the time
dependent temperature and electric fields described in Sec. III.2. For this
first case, we will take $\left(T_{init},T_{final},E_{0}\right)$ to be fixed,
and seek to learn the solution as a function of $\Delta t_{TQ}$. Specifically,
we will set $T_{init}=15\;\text{keV}$ and $T_{final}=30\;\text{eV}$, with
$E_{0}$ evaluated for a quasi-neutral deuterium plasma with
$n_{D}=10^{14}\;\text{cm}^{-3}$ and $j_{0}=1\;\text{MA}/\text{m}^{2}$. The hot
tail PINN will be trained for thermal quench times ranging from
$0.25\;\text{ms}$ to $2\;\text{ms}$. The total time of each simulation is
taken to be ten times the thermal quench time. The time evolution of both the
temperature and electric field are shown in Fig. 3 for $\Delta
t_{TQ}=0.5\;\text{ms}$. Comparing with the time history used in the RAMc
simulation for identical parameters (Fig. 1), the two sets of fields are in
good agreement, indicating that we will be able to perform a direct comparison
between predictions of the hot tail PINN and the Monte Carlo simulations.

Figure 4: Loss history [panel (a)] and RPF at $\xi=-1$ [panel (b)]. 10,000
iterations were performed with the ADAM optimizer, followed by over a million
iterations performed with the L-BFGS optimizer. The thermal quench time was
taken to be $\Delta t_{TQ}=1\;\text{ms}$, $j_{0}=1\;\text{MA}/\text{m}^{2}$,
and the initial and final temperatures are given by $T_{init}=15\;\text{keV}$
and $T_{final}=30\;\text{eV}$, respectively.

The loss history of the RPF trained with these evolving temperature and
electric fields is shown in Fig. 4(a). At increments of twenty-five thousand
steps additional training points are added at locations where the residual of
the PDE is maximal, causing the periodic spikes in the training loss evident
in Fig. 4(a). Due to the sharp variation of the RPF when the electric field is
maximal [i.e. for $\tau\equiv t_{final}-t\approx 0$, see Fig. 4(b)],
additional training points have also been added to the initial $10$% of the
simulation, i.e. from $\tau=0$ to $\tau=0.5\;\text{ms}$ for the case shown in
Fig. 4(b). These two differences in the training versus test point
distribution lead to a significant deviation in the training versus test loss,
with the training loss dropping to a value of approximately $1.3\times
10^{-6}$ after $\sim 1,000,000$ iterations, with a test loss of roughly
$1.7\times 10^{-7}$. The approximate seven orders of magnitude drop of the
test loss indicates that the hot tail PINN was able to identify a well
converged solution.

Figure 5: RPF at $\tau=0$ [panel (a)], $\tau=0.15\;\text{ms}$ [panel (b)],
$\tau=1.5\;\text{ms}$ [panel (c)] and $\tau=5\;\text{ms}$ [panel (d)]. The
thermal quench time was taken to be $\Delta t_{TQ}=0.5\;\text{ms}$, with the
initial and final temperatures given by $T_{init}=15\;\text{keV}$ and
$T_{final}=30\;\text{eV}$, respectively.

Time slices of the RPF computed from the hot tail PINN are indicated in Fig.
5. The first time slice [Fig. 5(a)] indicates the terminal RPF, where
$P^{\prime}\left(\tau=0\right)=P_{RE}$. As indicated in Sec. III.3 the
terminal RPF is taken to be a modestly smoothed step function. At $\tau=0$
(i.e. $t=t_{final}$), the electric field is at its maximum value, which noting
the backward nature of Eq. (4), implies electrons with $\xi=-1$ will move to
lower energy, whereas those with $\xi=1$ will increase their energy. This
results in the RPF front for electrons shifting downward in energy for
$\xi=-1$ and moving to higher energy for $\xi=1$ [compare Fig. 5(a) with
5(b)]. By $\tau=1.5\;\text{ms}$, the RPF takes on a form similar to the steady
state solution Liu et al. (2017). As the system evolves further the electric
field decreases in magnitude, resulting in the RPF front shifting to higher
energy. At $\tau=t_{final}$ the RPF front has shifted to approximately
$E\approx 150\;\text{keV}$. This last RPF distribution indicates the
probability of an electron running away at or before $t_{final}$, and is thus
the appropriate value to use when evaluating the hot tail seed. Comparing this
RPF with the RPF evaluated by RAMc shown in Fig. 2(d) above, good agreement is
evident between these two results indicating that the hot tail PINN was able
to accurately learn the RPF.

###  V.2 Hot Tail Seed

####  V.2.1 Dependence on the Timescale of the Thermal Quench

Figure 6: RPF at $\tau=t_{final}$ [panel (a)], log scale plot of the RPF at
$\tau=t_{final}$ [panel (b)], integrand of Eq. (13) [panel (c)], filtered
integrand of Eq. (15) [panel (d)]. The thermal quench time was taken to be
$\Delta t_{TQ}=0.5\;\text{ms}$, with the initial and final temperatures given
by $T_{init}=15\;\text{keV}$ and $T_{final}=30\;\text{eV}$, respectively.

Our first task will be to use the RPF evaluated in Sec. V.1 to evaluate the
hot tail seed. Taking the electron distribution before the thermal quench to
be a Maxwell-Jüttner distribution $f^{MJ}_{e}$ with a temperature of
$T_{init}$ and density $n_{e0}$, the number of hot tail electrons will then be
given by:

$n_{RE}=\int d^{3}pf^{MJ}_{e}P^{\prime}\left(\tau=t_{final}\right),$ (13)

where $P^{\prime}\left(\tau=t_{final}\right)$ will be evaluated using the hot
tail PINN. A complication in performing this analysis arises due to the
Maxwell-Jüttner distribution being peaked at low energies, and then decaying
exponentially with energy, such that a small error in the RPF at low energies
could lead to a spurious contribution to $n_{RE}$. This complication is
illustrated in Fig. 6 for $\Delta t_{TQ}=1\;\text{ms}$. The RPF at
$\tau=t_{final}$ is shown in Fig. 6(a), where the RE transition region is
centered around $E\approx 260\;\text{keV}$, with the RPF having a small value
for energies less than $\approx 200\;\text{keV}$. However, considering the
integrand that enters Eq. (6) [Fig. 6(c)] two maximums are evident. The first
is centered at an energy of approximately $190\;\text{keV}$, and is peaked
near the $\xi=-1$ axis. This corresponds to the peak associated with hot tail
generation. In particular, while the RPF has a modest value at this energy,
noting that a Maxwell-Jüttner distribution decays exponentially with energy,
the dominant contribution to the hot tail seed is expected to occur at the
transition front of the RPF. We will need to capture this maximum of the
integrand to evaluate the hot tail seed.

A second maximum of the hot tail integrand is also evident from Fig. 6(c).
This maximum spans energies between $\sim 8\;\text{keV}$ up to values modestly
under $\sim 40\;\text{keV}$. There is no physical mechanism that would
generate REs from this region of momentum space. The origin of this peak can
instead be tracked to a small inaccuracy of the RPF at low energy. Considering
a log scale plot of the RPF [Fig. 6(b)], the shape of the RPF for energies
less than $100\;\text{keV}$ is slightly distorted, with a maximum centered at
$\sim 20\;\text{keV}$ evident. This unphysical feature is due to the RPF being
exceptionally small in this region, with the artificial maximum having a
magnitude of $\sim 10^{-4}$, such that the residual of Eq. (4) does not
heavily weigh this region.

In order to remove this low energy unphysical contribution to the hot tail
integrand, we will choose a weighting scheme for the loss function that
enables regions with $P^{\prime}\ll 1$ to be strongly weighted, along with
post-process the predictions of the hot tail PINN. With regard to the former
approach, to increase the weight placed on regions where the RPF is small, we
have defined the loss function to be of the form:

$\text{Loss}=\frac{1}{N_{PDE}}\sum^{N_{PDE}}_{i}\left[\left(\frac{1}{\delta+P^{\prime}\left(p_{i},\xi_{i},\tau_{i},\mathbf{\lambda}_{i}\right)}\right)\left(\frac{p^{2}_{i}}{1+p^{2}_{i}}\right)\mathcal{R}_{i}\right]^{2}.$
(14)

Here, $\mathcal{R}_{i}$ corresponds to the residual of Eq. (4), the factor
$p^{2}_{i}/\left(1+p^{2}_{i}\right)$ removes the divergence of the pitch-angle
scattering operator at low energy, and the factor
$1/\left(\delta+P^{\prime}\right)$, increases the weight of regions with small
values of $P^{\prime}$. A loss function of this form was used when training
the RPF shown in Fig. 6, with $\delta$ chosen to be $\delta=0.05$, which
enabled the hot tail PINN to learn the RPF with sufficient accuracy such that
a clear separation between the physical and unphysical contributions to the
hot tail integrand is evident. Smaller values of $\delta$ could be used to
further separate the physical and unphysical contributions to the hot tail
integrand, however, the $\delta$ used in the present example was sufficient to
enable a simple post processing algorithm to remove the spurious contribution
from the integrand.

The first step of the post processing algorithm will be to limit the
integration region used in Eq. (13) to energies three times greater than the
thermal energy. In so doing the very low energy region, where spurious
features are most likely, will be removed. For the example indicated in Fig. 6
this filter will remove the contribution to the integrand for electrons with
energies less than $45\;\text{keV}$ and is thus sufficient to nearly
completely remove the unphysical contribution. By limiting ourselves to
electrons whose initial energy is greater than three times the thermal energy
the maximum hot tail seed that the model will be able to predict is given by
$n_{RE}/n_{e0}\approx 0.12$. However, noting that a relativistic electron
population with $\xi=-1$ and a magnitude of $n_{RE}/n_{e0}\sim 10^{-4}$ is
sufficient to carry most, if not all, of the plasma current, the prediction of
a RE seed comparable to $n_{RE}/n_{e0}\sim 10^{-1}$ would be unphysical. In
particular, for cases where particularly large RE seeds are present, this
would lead to a modification of the Ohm’s law defined by Eq. (6).

Noting that the maximum contribution to the hot tail integrand increases as
$\Delta t_{TQ}$ is increased, adjusting the lower bound on the integration
will not be sufficient for very slow thermal quenches. An additional filter
will be introduced by removing regions from the integrand where the RPF is
smaller than a critical value, indicated by $\Delta P\ll 1$. In so doing, the
predicted hot tail seed will be computed only using regions where the RPF has
a value greater than $\Delta P$, which can be chosen to exclude unphysical
contributions such as those indicated in Fig. 6(b). As a final filter we will
use the residual of the PDE to remove regions where the relative accuracy of
the solution is low. This is accomplished by only including contributions
where the residual of the PDE $R_{PDE}$ normalized to $P^{\prime}$ is greater
than a threshold value $\Delta R_{thres}$. Specifically, defining the
normalized residual $R_{norm}\equiv R_{PDE}/P^{\prime}$, we will require the
integrand of Eq. (13) to vanish when $\left|R_{norm}\right|>\Delta R_{thres}$.
These filters will be implemented by defining the function

$\chi\left(P^{\prime}\right)=\left[1-\exp\left(-\frac{{P^{\prime}}^{10}}{\Delta
P^{10}}\right)\right]\exp\left(-\frac{{R}^{10}_{norm}}{\Delta
R^{10}_{thres}}\right),$

such that the filtered integrand for Eq. (13) becomes

$n^{filter}_{RE}=\int
d^{3}pf^{MJ}_{e}P^{\prime}\left(\tau=t_{final}\right)\chi\left(P^{\prime}\right).$
(15)

The integrand of Eq. (15) is shown in Fig. 6(d) for the parameters $\Delta
P=10^{-4}$ and $\Delta R_{thres}=1$. It is evident that the spurious maximum
in the integrand has been removed, whereas the maximum associated with hot
tail generation is unaffected.

Figure 7: Predicted hot tail seed for different values of $\Delta P$ [panel
(a)] and different values of $R_{thres}$ [panel (b)]. For panel (a),
$R_{thres}$ was set to $1$ whereas for panel (b) $\Delta P$ was set to
$10^{-4}$. The integration region was taken to be between
$E_{min}=3T_{init}=45\;\text{keV}$ and $E_{max}=500\;\text{keV}$.

To determine the sensitivity of the predicted hot tail seed to the choice of
$\Delta P$ and $\Delta R_{thres}$, the predicted hot tail seed versus the
thermal quench time is shown in Fig. 7 for several values of $\Delta P$ and
$\Delta R_{thres}$. It is evident that for rapid thermal quenches ($\Delta
t_{TQ}\lesssim 1\;\text{ms}$) all cases are in approximate agreement with each
other, implying that the only filtering required is to set the low energy
limit of the hot tail integrand to $3T_{init}$. However, for slower thermal
quenches, the unfiltered result deviates from the filtered cases. This is due
to the transition region of the RPF shifting to higher energies for larger
$\Delta t_{TQ}$, allowing the integrand to be more easily polluted by small
inaccuracies in the RPF. The number of predicted hot tail REs is weakly
sensitive to the values of $\Delta P$ and $\Delta R_{thres}$ used, where the
percent difference is largest for slow thermal quenches. For the remainder of
this paper we will use $\Delta P=10^{-4}$ and $\Delta R_{thres}=1$. We note
that for all the cases treated in this paper, a filter based on excluding
contributions to the integrand with a value of the RPF below $10^{-4}$ is
sufficient to remove spurious contributions to the integrand. We have
introduced the two other approaches, limiting the range of integrand and
residual based filtering, with an eye toward future studies where higher
dimensional hot tail PINNs will be employed, and a more aggressive post-
processing procedure may be needed.

####  V.2.2 General Solution

Figure 8: Loss history for a 6-D model. The model was trained for $\Delta
t_{TQ}$ between $0.25\;\text{ms}$ and $2\;\text{ms}$, $T_{init}$ between
$10\;\text{keV}$ and $20\;\text{keV}$, and current densities between
$j_{0}=0.5\;\text{MA}$ and $j_{0}=1.5\;\text{MA}$. The other parameters of the
model were taken to be $T_{final}=30\;\text{eV}$ and $Z_{eff}=1$.

In this section we will explore the broader parameter space of the idealized
thermal quench model described above. In particular, we will train the model
for a range of values for the parameters $\left(\Delta
t_{TQ},T_{init},E_{init}\right)$ for fixed $T_{final}$. This latter parameter
will be taken to be $T_{final}=30\;\text{eV}$, with the other three parameters
trained within ranges of $\Delta t_{TQ}$ between $0.25\;\text{ms}$ and
$2\;\text{ms}$, $T_{init}$ between $10\;\text{keV}$ and $20\;\text{keV}$, and
the range of values of $E_{init}$ are chosen such that a plasma at a density
of $n_{e}=10^{14}\;\text{cm}^{-3}$ will have plasma currents that range
between $0.5\;\text{MA}/\text{m}^{2}$ and $1.5\;\text{MA}/\text{m}^{2}$. The
loss history for the model is shown in Fig. 8, with a minimum training loss of
roughly $\approx 10^{-6}$ and a test loss of $\approx 2.4\times 10^{-7}$. The
hot tail seed evaluated as a function of $\Delta t_{TQ}$ for different initial
temperatures and current densities is shown in Fig. 9. Here, we have only
included values of the hot tail seed less than $n_{RE}/n_{e0}=5\times 10^{-4}$
since for RE seeds larger than this the number of REs would be large enough to
carry most, if not all, of the plasma current, thus violating the Ohm’s law
employed in Eq. (6). From Fig. 9 it is apparent that the magnitude of the hot
tail seed is strongly sensitive to the initial temperature of the plasma
$T_{init}$, with a somewhat weaker dependence on the current density $j_{0}$.
The strong dependence on $T_{init}$ is largely due to the number of energetic
electrons in the initial Maxwell-Jüttner distribution depending exponentially
on $T_{init}$. In addition, $T_{init}$ impacts the electric field evolution
through the $T_{e}$ dependence in Ohm’s law [Eq. (6)], though this dependence
will have a far more modest impact on the number of hot tail REs generated.
The current density $j_{0}$ (which enters via $E_{0}$) has a direct impact on
the magnitude of the electric field. This influences the number of hot tail
electrons via two distinct mechanisms. The first is that the magnitude of the
electric field sets the critical energy to run away, which directly impacts
the number of electrons that will be accelerated rather than slowed down. A
second impact is that the electric field will also affect the form of the
electron distribution function. This dependence will introduce additional
anisotropy in the RPF, and thus have a more subtle impact on the hot tail
mechanism.

Figure 9: The hot tail seed evaluated from a 6-D model as a function of the
thermal quench time scale $\Delta t_{TQ}$ for different initial plasma
temperatures $T_{init}$ and current densities $j_{0}$. Panel (a) is for
$j_{0}=0.5\;\text{MA}/\text{m}^{2}$ and panel (b) is for
$j_{0}=1\;\text{MA}/\text{m}^{2}$. Panel (c) shows a comparison with Eqs. (23)
and (22) of Ref. Smith and Verwichte (2008), for $T_{init}=15\;\text{keV}$ and
$j_{0}=1\;\text{MA}/\text{m}^{2}$. The ‘x’ markers indicate values found using
the Monte Carlo code RAMc McDevitt et al. (2019). The other parameters were
taken to be $T_{final}=30\;\text{eV}$ and $Z_{eff}=1$.

A comparison with Ref. Smith and Verwichte (2008) and first principle
simulations of the hot tail mechanism using the RAMc code McDevitt et al.
(2019) is shown in Fig. 9(d). Here, the cyan curve indicates the approximate
analytic theory defined by Eqs. (23) and (22) of Ref. Smith and Verwichte
(2008). From Fig. 9(d) the hot tail PINN is in excellent agreement with the
predictions from Monte Carlo simulations performed with RAMc, with the
approximate analytic theory of Ref. Smith and Verwichte (2008) overestimating
the magnitude of the hot tail seed. A comparison with Monte Carlo predictions
across a broader range of parameters is shown in Fig. 10. Here, we have
considered a random sampling of the parameters $\left(\Delta
t_{TQ},T_{init},j_{0}\right)$ across the range over which the hot tail PINN
was trained. Excellent agreement is apparent between the direct Monte Carlo
simulations and the predictions of the hot tail PINN, suggesting the hot tail
PINN was able to accurately learn the solution across the three-dimensional
parameter space.

Figure 10: Comparison between the hot tail PINN and direct Monte Carlo
simulations over a broad parameter space. $\left(\Delta
t_{TQ},T_{init},j_{0}\right)$ were sampled randomly, with
$T_{final}=30\;\text{eV}$ for all cases.

In addition to evaluating the number of seed electrons, the hot tail PINN can
be used to infer additional properties of the electrons that make up the hot
tail seed. Two properties that we will focus on here will be the average
initial pitch and energy of the electrons that eventually run away. These can
be computed by evaluating the moments

$\left\langle\xi\right\rangle=\frac{1}{n_{RE}}\int d^{3}p\xi
f^{MJ}_{e}P^{\prime}\left(\tau=t_{final}\right)\chi\left(P^{\prime}\right),$
(16) $\left\langle E\right\rangle=\frac{1}{n_{RE}}\int
d^{3}pm_{e}c^{2}\left(\gamma-1\right)f^{MJ}_{e}P^{\prime}\left(\tau=t_{final}\right)\chi\left(P^{\prime}\right).$
(17)

Considering the average pitch first, this quantity is shown as a function of
initial temperature and thermal quench time in Fig. 11(a) for an initial
current density of $j_{0}=1\;\text{MA}$. Here it is apparent that the
electrons responsible for running away are on average moving in the negative
pitch direction as expected, where this anisotropy becomes most acute for slow
thermal quenches (long $\Delta t_{TQ}$) and low initial electron temperatures.
Noting that several models of the hot tail neglect the anisotropy of the
electron distribution, it is apparent that this approximation is never
precisely satisfied, but breaks down quite severely for the case of a slow
thermal quench. Turning to the average initial energy, here it is apparent
that the typical initial energy of a hot tail RE ranges from slightly less
than $100\;\text{keV}$, but can reach magnitudes up to roughly
$200\;\text{keV}$, where the timescale of the thermal quench is the most
important driver. At such high energies, a relativistic collision operator is
required.

Figure 11: Average initial pitch [panel (a)] and energy [panel (b)] of the hot
tail electrons as a function of the initial temperature and thermal quench
timescale. The plasma current was taken to be $j_{0}=1\;\text{MA/m}^{2}$.

##  VI Discussion and Conclusions

In the above, a PINN representation of the adjoint relativistic Fokker-Planck
equation was developed in the presence of a rapidly evolving electron
temperature and electric field. This hot tail PINN was used to evaluate the RE
seed that emerges during an idealized axisymmetric thermal quench. The
predictions from the hot tail PINN were shown to be in excellent agreement
with direct solutions from the Monte Carlo code RAMc across a broad range of
parameters, thus providing an effective surrogate model of the hot tail. Aside
from the number of electrons predicted to run away, the PINN framework enabled
additional properties such as the average initial pitch and energy of the
electrons that would eventually run away to be extracted. This latter
characteristic represents a strong point of the PINN framework, where rather
than output a single quantity such as the runaway density, the entire RPF
solution is encoded. This not only enables additional insight into properties
of the solution, but also allows for greater interpretability of the
predictions made by the hot tail PINN, and thus provides additional guidance
as to whether a prediction of the surrogate model can be trusted.

While the present study provides a proof-of-principle demonstration of the
potential of the hot tail PINN, several idealizations would need to be removed
before a truly predictive model of the hot tail seed for a realistic tokamak
disruption can be developed. Some of these idealizations would be conceptually
simple to relax, such as treating a nonhydrogenic plasma or a more realistic
temperature history, others would require significant extensions of the hot
tail PINN developed above. Among the most challenging would be the inclusion
of transport due to the partial or complete stochastization of the magnetic
field arising from the presence of 3D MHD instabilities. While an ad-hoc
spatial diffusivity due to MHD transport could be straightforwardly included
in the above model, such an approach is unlikely to provide sufficient
accuracy for a predictive hot tail model. Future work will instead seek to
take advantage of the ease through which inverse problems may be treated
within the PINN framework Cai et al. (2021); Mathews et al. (2022), with the
aim of identifying a representative description of electron transport by a
given spectrum of MHD modes. Once inferred, this transport model would be
included in the hot tail PINN, where the properties of the MHD spectrum would
correspond to additional parameters over which the hot tail PINN will be
trained over.

###### Acknowledgements.

This work was supported by the SciDAC project Tokamak Disruption Simulation
(TDS) by Office of Fusion Energy Science and Office of Advanced Scientific
Computing, and by the University of Florida Informatics Institute seed award.
The authors acknowledge the University of Florida Research Computing for
providing computational resources that have contributed to the research
results reported in this publication.

## References

  * Dwyer and Uman (2014) J. R. Dwyer and M. A. Uman, Physics Reports 534, 147 (2014).
  * Breizman et al. (2019) B. N. Breizman, P. Aleynikov, E. M. Hollmann, and M. Lehnen, Nuclear Fusion 59, 083001 (2019).
  * SOKOLOV (1979) I. SOKOLOV, JETP Letters 29, 218 (1979).
  * Dreicer (1959) H. Dreicer, Physical Review 115, 238 (1959).
  * Connor and Hastie (1975) J. Connor and R. Hastie, Nuclear fusion 15, 415 (1975).
  * Hesslow et al. (2019) L. Hesslow, L. Unnerfelt, O. Vallhagen, O. Embréus, M. Hoppe, G. Papp, and T. Fülöp, Journal of Plasma Physics 85, 475850601 (2019).
  * McDevitt and Tang (2019) C. McDevitt and X.-Z. Tang, EPL (Europhysics Letters) 127, 45001 (2019).
  * Chacon et al. (2022) L. Chacon, D. Daniel, and W. T. Taitano, Journal of Computational Physics 449, 110772 (2022).
  * Martín-Solís et al. (2017) J. Martín-Solís, A. Loarte, and M. Lehnen, Nuclear Fusion 57, 066025 (2017).
  * Aleynikov and Breizman (2017) P. Aleynikov and B. N. Breizman, Nuclear Fusion 57, 046009 (2017).
  * Helander et al. (2004) P. Helander, H. Smith, T. Fülöp, and L.-G. Eriksson, Physics of Plasmas 11, 5704 (2004).
  * Smith and Verwichte (2008) H. M. Smith and E. Verwichte, Physics of plasmas 15, 072502 (2008).
  * Stahl et al. (2016) A. Stahl, O. Embréus, G. Papp, M. Landreman, and T. Fülöp, Nuclear Fusion 56, 112009 (2016).
  * Matsuyama and Yagi (2017) A. Matsuyama and M. Yagi, Plasma and Fusion Research 12, 1403032 (2017).
  * Raissi et al. (2019) M. Raissi, P. Perdikaris, and G. E. Karniadakis, Journal of Computational physics 378, 686 (2019).
  * Karniadakis et al. (2021) G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, Nature Reviews Physics 3, 422 (2021).
  * Karney and Fisch (1986) C. F. Karney and N. J. Fisch, The Physics of fluids 29, 180 (1986).
  * McDevitt et al. (2019) C. J. McDevitt, Z. Guo, and X. Tang, Plasma Physics and Controlled Fusion 61, 054008 (2019).
  * Lagaris et al. (1998) I. E. Lagaris, A. Likas, and D. I. Fotiadis, IEEE transactions on neural networks 9, 987 (1998).
  * Karpatne et al. (2017) A. Karpatne, G. Atluri, J. H. Faghmous, M. Steinbach, A. Banerjee, A. Ganguly, S. Shekhar, N. Samatova, and V. Kumar, IEEE Transactions on knowledge and data engineering 29, 2318 (2017).
  * Lusch et al. (2018) B. Lusch, J. N. Kutz, and S. L. Brunton, Nature communications 9, 1 (2018).
  * Wang et al. (2020) R. Wang, K. Kashinath, M. Mustafa, A. Albert, and R. Yu, in _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_ (2020), pp. 1457–1466.
  * Liu et al. (2016) C. Liu, D. P. Brennan, A. Bhattacharjee, and A. H. Boozer, Physics of Plasmas 23, 010702 (2016).
  * Zhang and del Castillo-Negrete (2017) G. Zhang and D. del Castillo-Negrete, Physics of Plasmas 24, 092511 (2017).
  * Hornik et al. (1989) K. Hornik, M. Stinchcombe, and H. White, Neural networks 2, 359 (1989).
  * Lu et al. (2021a) L. Lu, R. Pestourie, W. Yao, Z. Wang, F. Verdugo, and S. G. Johnson, SIAM Journal on Scientific Computing 43, B1105 (2021a).
  * Lu et al. (2021b) L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, SIAM Review 63, 208 (2021b).
  * McDevitt and Tang (2022) C. J. McDevitt and X.-Z. Tang, arXiv preprint arXiv:2211.02160 (2022).
  * Liu et al. (2017) C. Liu, D. P. Brennan, A. H. Boozer, and A. Bhattacharjee, Plasma Physics and Controlled Fusion 59, 024003 (2017), URL http://stacks.iop.org/0741-3335/59/i=2/a=024003.
  * Cai et al. (2021) S. Cai, Z. Mao, Z. Wang, M. Yin, and G. E. Karniadakis, Acta Mechanica Sinica 37, 1727 (2021).
  * Mathews et al. (2022) A. Mathews, J. W. Hughes, J. L. Terry, and S.-G. Baek, Physical Review Letters 129, 235002 (2022).