An efficient surrogate model of secondary electron formation and evolution

McDevitt, Christopher J.; Arnaud, Jonathan S.; Tang, Xian-Zhu

doi:10.1063/5.0253374

This work extends the adjoint-deep learning framework for runaway electron (RE) evolution, developed by McDevitt et al. [Phys. Plasmas 32, 042503 (2025)], to account for large-angle collisions. By incorporating large-angle collisions, the framework allows the avalanche of REs to be captured, an essential component of RE dynamics. This extension is accomplished by using a Rosenbluth–Putvinski approximation to estimate the distribution of secondary electrons generated by large-angle collisions. By evolving both the primary and multiple generations of secondary electrons, the present formulation can capture both the detailed temporal evolution of a RE population beginning from an arbitrary initial momentum space distribution, along with providing approximations to the saturated growth and decay rates of the RE population. Predictions of the adjoint-deep learning framework are verified against a traditional RE solver, with good agreement present across a broad range of parameters.

I. INTRODUCTION

The description of runaway electron (RE) formation and evolution in magnetized fusion devices poses a substantial scientific challenge. While reduced analytic,^2–6 semi-analytic,^7–9 or machine learning^10–13 models are often used to infer RE generation rates, threshold electric fields, or the number of seed REs, such models are unable to recover the detailed dynamics of RE evolution. Instead, the description of RE dynamics is often performed using continuum^14–18 or particle-based RE solvers.^19–22 In Paper I,¹ the decay rate of a population of primary electrons was evaluated using a parametric solution to the adjoint of the relativistic Fokker–Planck equation. By targeting the adjoint to the relativistic Fokker–Planck equation, a single solution to this equation for a given set of parameters allowed the RE density to be evolved forward in time beginning from an arbitrary initial momentum space distribution of electrons. The framework was further generalized by utilizing a physics-informed neural network (PINN) to identify a parametric solution to the adjoint of the relativistic Fokker–Planck equation. Hence, once trained, the combined adjoint-deep learning framework enabled the dynamics of the RE density to be inferred beginning from an arbitrary initial momentum space distribution for a range of physics parameters. While the offline training time of the PINN was substantially greater than the solution time of a traditional RE solver, its online inference time was orders of magnitude faster, suggesting its potential as a surrogate model for RE evolution in integrated descriptions of tokamak disruptions.

Reference 1, however, did not describe the exponential growth of the RE population due to the avalanche mechanism,²³ an essential component when describing the conversion of Ohmic to RE current in tokamak devices.^24,25 The present paper seeks to remove this limitation by extending the adjoint-deep learning framework to incorporate large-angle collisions, thus allowing for the amplification of an initial seed population. A second generalization of the framework will be achieved by training across a broader range of physics parameters compared to Ref. 1. This is accomplished by introducing a normalized time and energy coordinate that allows features of the runaway probability function (RPF) to be accurately tracked while using a comparable number of training points as the narrow parameter range employed in Ref. 1. The derived adjoint-deep learning approach will thus provide an efficient framework through which RE dynamics can be accurately described.

The present paper is not the first to evaluate the avalanche growth rate of REs using an adjoint approach,²⁶ nor the first to use a PINN for evaluating the avalanche growth rate.¹³ With regard to Ref. 13, that paper employed a PINN to solve the steady state adjoint of the relativistic Fokker–Planck equation. Using that solution, the saturated avalanche growth rate of the RE population was inferred. A limitation of Ref. 13, however, was that it computed the rate that secondary electrons are generated and run away, without accounting for the decay of the primary distribution. Thus, when below the threshold, Ref. 13 simply predicted a negligibly small avalanche growth, without accounting for the decay of the overall RE population. This resulted in a slight ambiguity in the determination of the RE avalanche threshold, together with the inability to correctly describe RE decay when below the threshold. Furthermore, Ref. 13 focused on the saturated growth rate of REs, and thus did not treat the transient evolution of the RE population. Results from the present model will fill these gaps by providing an efficient means of evaluating the temporal evolution of the RE density both above and below the threshold. The near and below threshold limits are of particular importance during the RE plateau, where accounting for the decay of the RE distribution is essential for accurately predicting the dissipation of the RE beam. While analytic approximations to this rate are available, see for example, Ref. 4, such formulas are highly inaccurate when below the threshold.¹ When carrying out these extensions, we will consider the limit of a completely screened plasma for a constant electric field. A more comprehensive treatment of REs during a disruption will require including a temporally varying electric field and partial screening corrections,²⁷ which will be left to a future work.

The remainder of this paper is organized as follows. Section II discusses an extension of the adjoint-deep learning formulation of the relativistic Fokker–Planck equation to include a secondary source term. Physics-informed neural networks are briefly described in Sec. III. The temporal evolution of the number density of an avalanching RE population both above and below the threshold is given in Sec. IV. Section V evaluates the avalanche growth and decay rates across a broad range of plasma conditions and compares with a traditional RE solver. A brief discussion along with conclusions are given in Sec. VI.

II. ADJOINT FRAMEWORK FOR RELATIVISTIC ELECTRON EVOLUTION

A. The adjoint of the relativistic Fokker–Planck equation

This section will provide a brief description of an adjoint formulation for evolving the RE density beginning from an arbitrary initial momentum space distribution. A more detailed description is given in Ref. 1. The adjoint of the relativistic Fokker–Planck equation can be written as^11,26,28,29

(1)

\frac{\partial P}{\partial t} - E^{*} (P, t) - C^{*} (P, t) - R^{*} (P, t) = 0

(1a)

with the adjoint operators defined by

E^{*} (P, t) = E_{| |} ξ \frac{\partial P}{\partial p} + (\frac{1 - ξ^{2}}{p}) E_{| |} \frac{\partial P}{\partial ξ},

(1b)

C^{*} (P, t) = C_{F} \frac{\partial P}{\partial p} - \frac{C_{B}}{p^{2}} \frac{\partial}{\partial ξ} [(1 - ξ^{2}) \frac{\partial P}{\partial ξ}],

(1c)

R^{*} (P, t) = α p γ (1 - ξ^{2}) \frac{\partial P}{\partial p} - α \frac{ξ (1 - ξ^{2})}{γ} \frac{\partial P}{\partial ξ} .

(1d)

Here, we have defined the electron's pitch by

ξ \equiv p_{| |} / p \in [- 1, 1]

⁠, time is normalized as

t \to t / τ_{c}

⁠, with

τ_{c} \equiv 4 π ϵ_{0}^{2} m_{e}^{2} c^{3} / (e^{4} n_{e} ln Λ)

and Coulomb logarithm

ln Λ

⁠, momentum as

p \to p / (m_{e} c)

⁠, the electric field as

E_{| |} \to E_{| |} / E_{c}

⁠, where

E_{c} \equiv m_{e} c / (e τ_{c})

is the Connor–Hastie electric field,³ the Lorentz factor is defined as

γ = \sqrt{1 + p^{2}}

⁠,

α \equiv τ_{c} / τ_{s}

⁠, and

τ_{s} \equiv 6 π ϵ_{0} m_{e}^{3} c^{3} / (e^{4} B^{2})

is the timescale associated with synchrotron radiation, and since Eq. does not include toroidal geometry, magnetic trapping effects on RE generation are not included. The coefficients for the collisional drag

C_{F}

and pitch-angle scattering

C_{B}

in Eq. are defined by

C_{F} \equiv \frac{γ^{2}}{p^{2}},

(2a)

C_{B} \equiv \frac{γ}{2} \frac{(Z_{eff} + 1)}{p} .

(2b)

An adjoint problem for the runaway probability function (RPF) P can be defined by enforcing the terminal condition and momentum boundary conditions,

(3)

P (t = t_{final}) = Θ (p - p_{R E}),

(3a)

P (p = p_{min}) = 0,

(3b)

P (p = p_{max}) = {\begin{matrix} 1, & U_{p} (p = p_{max}) > 0, \\ unconstrained, & U_{p} (p = p_{max}) < 0, \end{matrix}

(3c)

where

Θ (x)

is a Heaviside function, and

U_{p}

is defined by

U_{p} \equiv - E_{| |} ξ - C_{F} - α p γ (1 - ξ^{2}) .

(4)

Here,

p_{R E}

indicates the momentum above which an electron will be counted as a RE, whose specific value will be discussed below. Once the adjoint problem defined by Eq. with the terminal and boundary conditions defined by Eq. has been solved, the RE density at a time between

0 \leq t \leq t_{final}

is given by¹

n_{R E} (t) = \int d^{3} p f_{e}^{(init)} (p, ξ) P (p, ξ, τ; t_{final}),

(5)

where

τ \equiv t_{final} - t

⁠. This latter formulation will allow the time evolution of the density moment of an arbitrary initial momentum space distribution of REs

f_{e}^{(init)} (p, ξ)

to be projected to any time between 0 and

t_{final}

⁠. It does not, however, include large-angle collisions, which are essential when describing RE populations for the large electric fields expected in tokamak disruptions. The extension to include large-angle collisions is discussed in Sec. II B as follows.

B. Runaway electron avalanche

Once the time evolution of the RE density is inferred, this can be used to evaluate the rate that a given seed population of REs will be enhanced by the avalanche mechanism. From Eq. , the evolution of the RE seed is given by

n_{seed} (t) = \int d^{3} p f_{e}^{(seed)} (p, ξ) P (p, ξ, τ; t_{final}) .

(6)

This seed population will generate additional REs via large-angle collisions. By introducing a temporal grid from

t = 0

to

t = t_{final}

⁠, in increments of

Δ t_{a v} \equiv t_{final} / N_{a v}

⁠, where

N_{a v}

is the number of avalanche time steps, the distribution of secondary electrons after one time step can be approximated by

f_{sec}^{(1)} (p, ξ, t_{1}) = Δ t_{a v} S (p, ξ, t_{1}),

(7)

where

t_{i} = i Δ t_{a v}

⁠,

f_{sec}^{(1)}

indicates the first generation of secondary electrons, and

S (p, ξ)

is the source of secondary electrons defined by

S (p, ξ, t_{1}) = \int d^{3} p^{'} S_{0} (p^{'}, ξ^{'}, p, ξ) f_{e} (p^{'}, ξ^{'}, t_{1}) .

(8)

Here,

f_{e} (p^{'}, ξ^{'}, t_{1})

is the seed electron population evaluated after one time step, and

S_{0} (p^{'}, ξ^{'}, p, ξ)

is given by

S_{0} (p^{'}, ξ^{'}, p, ξ) = n_{e} c r_{e}^{2} \frac{v^{'}}{2 π p^{2}} \frac{d σ_{M} (p^{'}, p)}{d p} Π (p^{'}, ξ^{'}, p; ξ),

(9)

where

r_{e} = e^{2} / (4 π ϵ_{0} m_{e} c^{2})

is the classical electron radius,

d σ_{M} / d p

is the Møller cross section,^30,31

Π (p^{'}, ξ^{'}, p; ξ)

describes the pitch-angle dependence of secondary electron generation (see Ref. 32 for an explicit expression),

d^{3} p = 2 π p^{2} dpd ξ

⁠, all variables have been dedimensionalized according to

p^{'} \to p^{'} / m_{e} c

⁠,

v^{'} \to v^{'} / c

⁠, and

σ_{M} \to σ_{M} / r_{e}^{2}

⁠. Since the adjoint formulation does not immediately give us the time evolution of the seed electron distribution

f_{e} (p^{'}, ξ^{'}, t)

⁠, only its density moment, we will need to introduce a closure. The simplest closure available, introduced in Ref. 4, takes the limit whereby the energy of the primary electrons is assumed to be infinite

p^{'} \to \infty

⁠, and the pitch is taken to be aligned with the magnetic field

ξ^{'} = - 1

⁠. Taking these limits, Eq. reduces to

S (p, ξ, t_{1}) = n_{e} n_{R E} (t_{1}) c r_{e}^{2} \frac{v}{γ^{2} - 1} \frac{1}{{(γ - 1)}^{2}} δ (ξ - ξ_{1}),

(10)

where

ξ_{1}

is defined by

ξ_{1} = - \sqrt{\frac{γ - 1}{γ + 1}} .

(11)

Substituting Eq. into yields an explicit expression for the momentum space distribution of the first generation of secondaries, i.e.,

f_{sec}^{(1)} (p, ξ) = Δ t_{a v} n_{e} n_{R E} (t_{1}) c r_{e}^{2} \frac{v}{γ^{2} - 1} \frac{1}{{(γ - 1)}^{2}} δ (ξ - ξ_{1}),

(12)

where the RE density at

t = t_{1} = Δ t_{a v}

will be given by the number of seed electrons

n_{R E} (t_{1}) = n_{seed} (t_{1})

⁠, which can be evaluated from Eq. . The evolution of the number density of the first generation of secondary REs can then be expressed as

n_{sec}^{(1)} (t) = \int d^{3} p f_{(sec)}^{(1)} (p, ξ) P (p, ξ, τ + t_{1}; t_{final}),

(13)

where the RPF is now evaluated at a time

τ + t_{1} = t_{final} - (t - t_{1})

⁠, since the first generation of secondaries is born at

t = t_{1} = Δ t_{a v}

⁠. Substituting Eq. into yields an explicit expression for the time evolution of the number density of the first generation of secondary electrons, i.e.,

n_{sec}^{(1)} (t) = 2 π Δ t_{a v} n_{e} n_{R E} (t_{1}) c r_{e}^{2} \int d p \frac{v}{{(γ - 1)}^{2}} P (p, ξ_{1}, τ + t_{1}; t_{final}) .

(14)

Similarly, the momentum space distribution of the second generation of secondaries is given by

f_{sec}^{(2)} (p, ξ) = Δ t_{a v} n_{e} n_{R E} (t_{2}) c r_{e}^{2} \frac{v}{γ^{2} - 1} \frac{1}{{(γ - 1)}^{2}} δ (ξ - ξ_{1}),

(15)

where

n_{R E} (t_{2}) = n_{seed} (t_{2}) + n_{sec}^{(1)} (t_{2})

⁠. An explicit expression for the next generation of secondary electrons can then be written as

\begin{matrix} n_{sec}^{(2)} (t) = \int d^{3} p f_{(sec)}^{(2)} (p, ξ) P (p, ξ, τ + t_{2}) \\ = 2 π Δ t_{a v} n_{e} n_{R E} (t_{2}) c r_{e}^{2} \int d p \frac{v}{{(γ - 1)}^{2}} P (p, ξ_{1}, τ + t_{2}; t_{final}) . \end{matrix}

(16)

where we have used Eq. in the second line. The ith generation of secondaries is then given by

(17)

n_{sec}^{(i)} (t) = 2 π Δ t_{a v} n_{e} n_{R E} (t_{i}) c r_{e}^{2} \int d p \frac{v}{{(γ - 1)}^{2}} P (p, ξ_{1}, τ + t_{i}; t_{final}) .

(17a)

The density of REs is

n_{R E} (t_{i}) = n_{seed} (t_{i}) + \sum_{j = 1}^{i - 1} n_{sec}^{(j)} (t_{i}) .

(17b)

By evaluating each successive generation of secondaries, together with the initial evolution of the seed RE population, this will enable the time history of an avalanching RE distribution to be determined. We note that while the present adjoint formulation only predicts the number density of REs, the RPF incorporates fully kinetic physics into the evolution of the RE population. The main approximation is due to the use of a Rosenbluth–Putvinski source [Eq. ], which does not account for the energy or pitch distribution of the RE population.

III. PHYSICS-INFORMED NEURAL NETWORKS

While a range of approaches to the numerical solution of Eq. could be applied, we will utilize a PINN^33,34 in the present manuscript. In so doing, this will allow us to identify the parametric solution to Eq. such that once trained, we will be able to evaluate the RPF for a range of physics parameters. The specifics of our implementation of a PINN can be found in Ref. 1, however, here, we will provide the essential formulas. In order to enforce several properties of the solution as hard constraints, we will introduce a physics layer of the form

P^{'} (p, ξ, t; λ) = P_{term} (p) + (\frac{p - p_{min}}{p_{max} - p_{min}}) tanh (t_{final} - t) P_{N N} (p, ξ, t; λ),

(18)

where

λ

is a vector containing the physics parameters

(E_{| |}, Z_{eff}, α)

⁠, and the terminal condition is defined by

P_{term} (p) \equiv \frac{1}{Δ P} tanh (\frac{p - p_{R E}}{Δ p}),

(19)

where the final RPF is computed by passing Eq. through

P (p, ξ, t; λ) = \frac{1}{2} {1 + tanh [P^{'} (p, ξ, t; λ)]} .

(20)

Here,

P_{N N}

represents the output of the hidden layers of the neural network (NN),

Δ p / p_{max}

and

Δ P

are hyperparameters whose value should be set to be much less than one, and

λ

is a vector representing the physics parameters, which are given by

λ = (E_{| |}, Z_{eff}, α)

⁠. It can be verified that for small

Δ p / p_{max}

and

Δ P

⁠, the physics layer enforces (i) the terminal condition given by Eq. , (ii) lower momentum boundary condition [Eq. ], and (iii) constrains the RPF to have a range between zero and one. With this physics layer, only the residual of the PDE and the boundary condition on the upper momentum boundary will be included in the loss function, such that we will minimize a loss of the form

loss = \frac{w_{PDE}}{N_{PDE}} \sum_{i}^{N_{PDE}} {[G (p) (\frac{p_{i}^{2}}{1 + p_{i}^{2}}) R (p_{i}, ξ_{i}, t_{i}; λ_{i})]}^{2} + \frac{1}{N_{bdy}} \sum_{i}^{N_{bdy}} {[P_{i} - P (p_{i}, ξ_{i}, t_{i}; λ_{i})]}^{2},

(21)

where

G (p) = 1 - exp [- \frac{{(p_{max} - p)}^{2}}{Δ p_{max}^{2}}],

(22)

where

N_{bdy}

is the number of points on the upper momentum boundary that will be sampled, and

R

is the residual of Eq. . The role of the factors

G (p)

and

p_{i}^{2} / (1 + p_{i}^{2})

are to remove large contributions to the residual that emerge at either high or low momenta boundaries, as discussed further in Ref. 1. The Python script used for training the PINN is written using the DeepXDE library,³⁵ with TensorFlow³⁶ as the backend, and will be made available upon acceptance at https://github.com/cmcdevitt2/RunAwayPINNs.

IV. TEMPORAL EVOLUTION OF AN AVALANCHING RUNAWAY ELECTRON POPULATION

This section will utilize the algorithm described in Sec. II B to describe the temporal evolution of a RE population including large-angle collisions. We will assume a seed electron population defined by

f_{e}^{(init)} (p, ξ) \propto exp [- \frac{{(p - p_{0})}^{2}}{Δ p^{2}} - \frac{{(ξ - ξ_{0})}^{2}}{Δ ξ^{2}}],

(23)

where the evolution of the number density of seed electrons will be evolved using Eq. , with each generation of secondary electrons described by Eq. . The momentum above which an electron will be counted as a RE is taken to be

p_{R E} = p_{max} / 4

⁠, which corresponds to an energy of roughly 1 MeV. While the use of a Rosenbluth–Putvinski secondary source term will result in a weak dependence on the choice of

p_{R E}

⁠, such that modestly improved agreement with first-principle Monte Carlo simulations could be achieved by tuning this parameter, we have opted not to do this since such tuning may be parameter dependent. We have set the hyperparameters to

Δ P = 0.15

⁠,

Δ p = 0.1 p_{max}

⁠,

Δ p_{max} = 0.05 p_{max}

⁠, and

w_{PDE} = 10

when training the model. The resulting loss history for the PINN is shown in Fig. 1. Here, we have trained the PINN for electric fields in the range

E_{| |} \in (0, 3)

⁠, effective charges

Z_{eff} \in (1, 5)

⁠, and synchrotron radiation

α \in (0, 0.2)

⁠, a somewhat broader range of parameters compared to Ref. 1, though at the cost of a much larger number of training points. A residual-based adaptive sampling scheme is employed,³⁷ leading to periodic spikes in the training loss.

FIG. 1.

View large Download slide

Training and test loss history. The solid blue curve indicates the training loss of the PDE, the solid red curve indicates the training loss for the boundary condition, and the “×” markers indicate the test loss (the test and training losses are the same for the boundary term). 3 000 000 training and testing points are used. The training points initially obey a Hammersley sequence, with the test distribution being uniform random. The training points are resampled using the residual-based adaptive distribution described in Ref. 37, every 50 000 epochs of L-BFGS.

Two example time histories are shown in Fig. 2. Here, the time evolution of the total number of REs (solid blue curve), initial seed population (dashed black curve), and several generations of secondary REs (multicolored dashed curves) are shown. For panel (a) of Fig. 2, the electric field was chosen to be above the threshold (⁠ $E_{| |} = 3$ ⁠), with the total number of REs increasing exponentially (solid blue curve). Considering each generation, the seed electron population is normalized to have a value of one at $t = 0$ ⁠, where its value remains largely unchanged since $E_{| |} > E_{a v}$ (⁠ $E_{a v} \approx 1.8$ for these parameters). For each avalanche time step, a small number of secondary electrons are born, where each generation grows rapidly until plateauing later in time. The initial rapid increase in the secondary electron population is due to the energy and pitch distribution of secondaries immediately after they are born obeying Eq. (10), where the majority of these electrons will have energies below the 1 MeV threshold energy for an electron to be counted as a RE. For cases above the threshold, a large number of the secondary electrons will be quickly accelerated above 1 MeV, resulting in an increase in the total number of REs. This process is repeated for each generation of secondaries, with each generation having a larger magnitude due to the number of REs increasing according to Eq. (17b). Summing together the number of primary electrons with all generations of secondary electrons results in an exponentially increasing number of REs, as indicated by the solid blue curve in Fig. 2. Turning to Fig. 2(b), here, the electric field was chosen to be below the threshold (⁠ $E_{| |} = 1.5$ ⁠), with the total number of REs initially increasing slightly before decaying slowly in time. Considering first the seed runaway population, the number of seed electrons initially remains approximately constant, but then begins to decay at a later time. The finite number of seed electrons leads to the generation of secondary electrons, and hence an initial increase in the total number of REs. However, since each secondary generation decays exponentially in time, this does not lead to the long-term growth of the overall population. At late times, the total number of REs decays exponentially at a rate that is similar, though not identical, to the decay of the seed RE population.

FIG. 2.

View large Download slide

Runaway electron evolution over time for two different values of the electric field. The number of total REs is shown by the blue curves; the seed population is indicated by the black dashed curve; and each generation of secondaries are indicated by the multicolor dashed curves. Panel (a) is for $E_{| |} = 3$ (above marginality), whereas panel (b) is for $E_{| |} = 1.5$ (below marginality). The other parameters are given by $Z_{eff} = 1$ ⁠, $α = 0.1$ ⁠, and $ln Λ = 15$ ⁠, with 19 generations of secondaries included (i.e., 20 time steps). The seed electron distribution was taken to be Eq. (23) with $(p_{0} = 3 p_{max} / 4, Δ p = p_{max} / 10, ξ_{0} = - 1, Δ ξ = 0.1)$ ⁠.

The RE density histories for several different electric fields are shown in Fig. 3, where the PINN's predictions (blue curves) are compared with the RE solver RAMc (black curves). RAMc is a particle-based RE solver²¹ that evolves guiding center orbits together with small and large angle collisions, as well as synchrotron radiation. For the present case, the electrons were initialized on the magnetic axis for a large ITER-like device, such that neoclassical effects arising from magnetic trapping^4,38–40 along with spatial transport⁴¹ are negligible. This will allow for comparisons with the present formulation, which does not account for toroidal geometry or spatial transport. Furthermore, RAMc uses a Møller source [Eqs. (8) and (9)] to describe secondary electron generation, a higher physics fidelity treatment of large-angle collisions compared to the Rosenbluth–Putvinski source used in the adjoint-deep learning approach. The present comparison will thus test the accuracy of the adjoint-deep learning framework along with regimes where a Rosenbluth–Putvinski source provides an adequate approximation to the large-angle collision operator. In all cases, the electron distribution was initialized to be a block of particles at high energies with pitch near $ξ \approx - 1$ [see Fig. 3(f)] to be consistent with the initial RE distribution implemented in the RAMc code. From Figs. 3(a) and 3(b), it is evident that when well above the threshold, the two approaches are in good agreement. However, when approaching the threshold [panels (c) and (d) of Fig. 3], more substantial differences between the two approaches are evident, with the PINN systematically overpredicting the RE density. The cause of this discrepancy can be linked to the Rosenbluth–Putvinski secondary source, which assumes the primary electrons to have asymptotically large energy. While very high energy primary electrons will be present when the system is above marginality, no such electrons are expected to be present when at or below marginality. As a result, the Rosenbluth–Putvinski source will prove inaccurate when near the threshold. Comparing the time evolution of the primary electron population [dashed curves in Figs. 3(d) and 3(e)] with the result from RAMc, it is evident that the decay rate of the primary electron population is in better agreement with the results of RAMc. The reason for this is that RAMc employs a Møller secondary source term that accounts for the energy and pitch of the primary electrons. When below the threshold, the energy of the primary distribution will be modest, such that the secondary electrons generated by the Møller source will have little impact on the RE number density, with the majority of them being born with too low an energy to run away. Thus, as the system drops further below marginality, neglecting the Rosenbluth–Putvinski source in the PINN leads to an increasingly good approximation to the RE density evolution.

FIG. 3.

View large Download slide

Comparison of RE density evolution with the Monte Carlo code RAMc for electric fields $E_{| |} = 3$ [panel (a)], $E_{| |} = 2.5$ [panel (b)], $E_{| |} = 2$ [panel (c)], $E_{| |} = 1.5$ [panel (d)], and $E_{| |} = 1$ [panel (e)]. The initial electron distribution assumed for both RAMc and the PINN is shown in panel (f), where a rectangular block of electrons were inserted. The other parameters were taken to be $Z_{eff} = 1$ ⁠, $α = 0.1$ ⁠, and $ln Λ = 15$ ⁠.

V. PARAMETRIC VARIATION OF RUNAWAY ELECTRON AVALANCHE GROWTH AND DECAY RATE

Our aim in the present section will be to train the PINN across a broader range of physics parameters. The treatment of a broader range of parameters will pose a substantial challenge to training an accurate PINN. Specifically, as the size of the parameter space is increased, the number of training points required to densely sample the space will grow substantially, thus requiring a longer period of offline training along with increasing memory usage. In addition, the location of the transition region of the RPF (i.e., where the RPF transitions between zero and one) and the time to reach a steady state solution will traverse a broad range of energy and time scales, further complicating capturing the RPF with a single PINN. To accommodate both of these challenges, we will introduce a normalized time and energy coordinate that will allow the PINN to adapt the time period and the energy range that must be captured in order to accurately infer the RPF evolution. Such an adaptive approach will allow us to capture a broader range of physics parameters for a given number of training points. To normalize the time period over which the RPF is evolved, we will utilize the Rosenbluth–Putvinski growth rate as a rough estimate of the rate that the RE density changes, i.e.,⁴

τ_{c} γ_{R P} = \frac{1}{ln Λ} \sqrt{\frac{π}{3 (Z_{eff} + 5)}} (E_{| |} - 1),

(24)

where the parallel electric field is normalized to the Connor–Hastie threshold

E_{c}

⁠. The time period over which the RPF will be evolved, denoted by

{\bar{t}}_{final}

⁠, will be taken to have the form

{\bar{t}}_{final} = \frac{t_{final}}{\sqrt{1 + γ_{a v}^{2} t_{final}^{2} / N_{a v}^{2}}},

(25)

where

N_{a v}

is the number of avalanche timescales that will be simulated and

t_{final}

provides an upper bound on the time period of the simulation. It may be readily verified that near the threshold where

t_{final} γ_{a v} / N_{a v} ≪ 1

⁠, Eq. asymptotes to

{\bar{t}}_{final} \approx t_{final}

⁠, whereas when far from the threshold where

t_{final} γ_{a v} / N_{a v} ≫ 1

⁠,

{\bar{t}}_{final}

asymptotes to

{\bar{t}}_{final} \approx N_{a v} / γ_{a v}

⁠. Equation thus allows for relatively long simulations to be carried out when the system is near the threshold, but substantially shrinks the time domain for scenarios far above the threshold where the RPF quickly reaches a quasi-steady state (see Fig. 6 below).

The upper bound in the momentum domain will be chosen to be larger than the momentum, where the electric field acceleration and drag balance for an electron with

ξ = - 1

⁠, which defines a critical momentum given by

p_{crit} = 1 / \sqrt{E_{| |} - 1}

⁠. To avoid the singularity at

E_{| |} = 1

⁠, we will approximate this critical momentum by

p_{crit} \approx 1 / \sqrt{E_{| |} + δ}

⁠, where

δ

is a small value that removes the divergence at

E_{| |} = 0

⁠, taken to be

δ = 0.01

in the present study. The upper momentum bound will then be taken to have the form

{\bar{p}}_{max} = \frac{p_{max}}{\sqrt{1 + p_{max}^{2} / {(N_{p} p_{crit})}^{2}}},

(26)

where

N_{p}

is a factor that sets how far above

p_{crit}

the high momentum boundary is set, and

p_{max}

limits the upper momentum boundary for the case of a weak electric field. For

p_{max}^{2} / {(N_{p} p_{crit})}^{2} ≪ 1

(i.e., a weak electric field), the upper momentum boundary is approximately given by

{\bar{p}}_{max} \approx p_{max}

⁠, whereas for a strong electric field

p_{max}^{2} / {(N_{p} p_{crit})}^{2} ≫ 1

⁠, Eq. asymptotes to

{\bar{p}}_{max} \approx N_{p} p_{crit}

⁠. An upper momentum bound given by Eq. thus allows a high momentum boundary of roughly

p_{max}

to be used when near or below marginality, but will substantially shrink the momentum domain for large electric fields where the transition region is located at low momenta.

Utilizing these adaptive time and momentum ranges, the loss history and final training point distribution are shown in Fig. 4. Here, we have trained the PINN to learn the RPF for electric fields in the range

E_{| |} \in (0, 10)

⁠,

Z_{eff} \in (1, 10)

⁠, and

α \in (0, 0.2)

for the hyperparameters

Δ P = 0.2

⁠,

Δ \bar{p} = 0.1 p_{max}

⁠,

Δ p_{max} = 0.05 {\bar{p}}_{max}

⁠, and

w_{PDE} = 10

⁠. The time and momentum domains are defined using Eqs. , with the values

t_{max} = 20

⁠,

N_{a v} = 1 / 2

⁠,

N_{p} = 10

⁠, and a

p_{max}

is chosen consistent with

10 MeV

⁠. The low momentum boundary is taken to be consistent with

10 keV

⁠, and

p_{R E}

is set to be

p_{R E} = {\bar{p}}_{max} / 4

⁠. We also slightly modified the loss function to

loss = \frac{w_{PDE}}{N_{PDE}} \sum_{i}^{N_{PDE}} {[G (p) \frac{E_{avg}}{1 + E_{| |}} (\frac{p_{i}^{2}}{1 + p_{i}^{2}}) R (p_{i}, ξ_{i}, t_{i}; λ_{i})]}^{2} + \frac{1}{N_{bdy}} \sum_{i}^{N_{bdy}} {[P_{i} - P (p_{i}, ξ_{i}, t_{i}; λ_{i})]}^{2},

(27)

FIG. 4.

View large Download slide

(a) Training and test loss history. The solid blue curve indicates the training loss of the PDE, the solid red curve indicates the training loss for the boundary condition, and the “×” markers indicate the test loss (the test and training losses are the same for the boundary term). 1 000 000 training and test points were used. (b) Final training point distribution.

where we have included the additional factor

E_{avg} / (1 + E_{| |})

in front of the residual, where

E_{avg} = (E_{min} + E_{max}) / 2

⁠, due to the magnitude of the residual of Eq. increasing with the electric field strength. Dividing by the electric field in the loss, this prevents the training of the PINN from being unduly biased toward cases with large electric fields. A residual-based adaptive training point sampling scheme³⁷ has been adopted, where from Fig. 4(b), it is evident that the density of training points is highest near

p / {\bar{p}}_{max} \approx 0.2

and at the high momentum boundary near

ξ \approx \pm 1

⁠. The concentration of training points near

p / {\bar{p}}_{max} \approx 0.2

is due to the value of

p_{R E}

being chosen to be

p_{R E} / {\bar{p}}_{max} = 0.25

⁠, whereas the concentration of points near the upper momentum boundary with

ξ \approx \pm 1

often contains the maximum value of the RPF when below the threshold (see the second column of Fig. 5 below).

FIG. 5.

Time slices of the RPF evolution (top row) along with associated residuals (bottom row) for E||=0.8 and pRE=p¯max/4. The bottom row indicates the residual to Eq. (1a) multiplied by p2/(1+p2), such that the low energy divergence is removed. The first column indicates the terminal condition at t≈19.36, the second column t≈17.43 and the last column t=0. The dashed white curves are the location of secondary injection ξ1 defined by Eq. (11). The other parameters are given by Zeff=5 and α=0.1, and we chose tfinal=20, Np=10, Nav=0.5, along with a pmax consistent with 10 MeV.

View large Download slide

Time slices of the RPF evolution (top row) along with associated residuals (bottom row) for $E_{| |} = 0.8$ and $p_{R E} = {\bar{p}}_{max} / 4$ ⁠. The bottom row indicates the residual to Eq. (1a) multiplied by $p^{2} / (1 + p^{2})$ ⁠, such that the low energy divergence is removed. The first column indicates the terminal condition at $t \approx 19.36$ ⁠, the second column $t \approx 17.43$ and the last column $t = 0$ ⁠. The dashed white curves are the location of secondary injection $ξ_{1}$ defined by Eq. (11). The other parameters are given by $Z_{eff} = 5$ and $α = 0.1$ ⁠, and we chose $t_{final} = 20$ ⁠, $N_{p} = 10$ ⁠, $N_{a v} = 0.5$ ⁠, along with a $p_{max}$ consistent with $10 MeV$ ⁠.

Two example RPFs for both weak and large electric fields are shown in Figs. 5 and 6, respectively. For the case of a weak electric field (⁠

E_{| |} = 0.8

⁠,

Z_{eff} = 5

⁠,

α = 0.1

⁠), a large energy domain is used, with an upper bound of

\approx 4.5 MeV

along with an integration time of nearly twenty

τ_{c}

⁠, such that the relatively slow evolution of the RPF can be captured. Here, the terminal condition is chosen such that

p_{R E} = {\bar{p}}_{max} / 4

⁠, which for the weak electric field case corresponds to an energy slightly below

1 MeV

⁠. As time evolves, the transition region shifts to higher energy, with maxima near

ξ \approx \pm 1

at

t \approx 17.43

(middle column in Fig. 5). Here, these maxima at high energy are due to synchrotron radiation preferentially slowing down electrons with large quantities of perpendicular energy (i.e., with

ξ \approx 0

⁠). By

t = 0

⁠, the RPF is nearly zero everywhere, except near the low energy boundary where it has a magnitude of

\sim 10^{- 6}

⁠. The reason the RPF does not vanish identically at the low energy boundary is due to the physics layer defined by Eqs. not precisely vanishing at

p = p_{min}

⁠. In particular, at

p = p_{min}

⁠, the RPF reduces to [see Eq. ]

P (p_{min}) = 0.5 {1 + tanh [P_{term} (p_{min})]},

FIG. 6.

Time slices of the RPF evolution (top row) along with the associated residual (bottom row) for E||=8 and pRE=p¯max/4. The bottom row indicates the residual to Eq. (1a) multiplied by p2/(1+p2), such that the low energy divergence is removed. The first column indicates the terminal condition at t≈2.19, the second column t≈1.65, and the last column t=0. The dashed white curves are the location of secondary injection ξ1 defined by Eq. (11). The other parameters are given by Zeff=5 and α=0.1, and we chose tfinal=20, Np=10, Nav=0.5, along with a pmax consistent with 10 MeV.

View large Download slide

Time slices of the RPF evolution (top row) along with the associated residual (bottom row) for $E_{| |} = 8$ and $p_{R E} = {\bar{p}}_{max} / 4$ ⁠. The bottom row indicates the residual to Eq. (1a) multiplied by $p^{2} / (1 + p^{2})$ ⁠, such that the low energy divergence is removed. The first column indicates the terminal condition at $t \approx 2.19$ ⁠, the second column $t \approx 1.65$ ⁠, and the last column $t = 0$ ⁠. The dashed white curves are the location of secondary injection $ξ_{1}$ defined by Eq. (11). The other parameters are given by $Z_{eff} = 5$ and $α = 0.1$ ⁠, and we chose $t_{final} = 20$ ⁠, $N_{p} = 10$ ⁠, $N_{a v} = 0.5$ ⁠, along with a $p_{max}$ consistent with $10 MeV$ ⁠.

which for

Δ P = 0.2

and

Δ p = 0.1 {\bar{p}}_{max}

does not precisely vanish. This small inaccuracy will have a negligible impact for parameters far above the threshold, and will only impact predictions on the PINN at late times when well below marginality. Furthermore, since the spurious contribution to the RPF is located at the low energy boundary, it can be removed via a post processing routine as discussed below.

Considering now a case well above marginality (⁠ $E_{| |} = 8$ ⁠, $Z_{eff} = 5$ ⁠, $α = 0.1$ ⁠, see Fig. 6), the RPF quickly reaches a steady state, with the RPF at $t \approx 1.65$ having nearly the same form as at $t = 0$ ⁠, with the total time of the simulation being 2.19. The steady state solution has a contour of roughly one half at an energy of $60 keV$ ⁠, suggesting that the maximum energy used in this case of $\approx 1 MeV$ is sufficient to capture the critical portion of the RPF.

Using the PINN described in Figs. 4–6, we will be interested in comparing predictions of the RE avalanche growth or decay rate with Monte Carlo simulations carried out with RAMc. Before carrying out this comparison, we note that the use of a Rosenbluth–Putvinski source to describe large-angle collisions introduces an ambiguity in the evaluation of the avalanche growth rate. In particular, since the Rosenbluth–Putvinski source only depends on the number of REs, and not on the energy distribution, it is necessary to introduce a criterion for classifying which electrons should be counted as REs. In the present formulation, the quantity $p_{R E}$ defined by Eq. (3) acts as the effective definition of which electrons are counted as REs. In Fig. 7, this value is roughly $2.5 p_{crit}$ for the hyperparameters selected, where $p_{crit}$ is the critical energy to run away given by $p_{crit} \approx 1 / E_{| |}$ ⁠. This definition ensures that only electrons located substantially above the critical energy to run away are included in the large-angle collision operator.

FIG. 7.

View large Download slide

Comparison of avalanche growth rates between the Monte Carlo code RAMc with a Møller secondary source and predictions of the PINN. The solid blue curve includes large-angle collisions, whereas the solid red curve includes large-angle collisions when above the threshold, but uses the primary decay rate when below the threshold, and the black “×” markers are values from RAMc. The other parameters were $Z_{eff} = 5$ ⁠, $α = 0.1$ ⁠, and $ln Λ = 16.15$ ⁠. The initial electron distribution was defined by Eq. (23) with $(p_{0} = 3 {\bar{p}}_{max} / 4, Δ p = {\bar{p}}_{max} / 10, ξ_{0} = - 1, Δ ξ = 0.1)$ ⁠.

An additional subtlety that arises when using an adaptive energy region is that, for the case of a large electric field, the modest value of the upper momentum boundary results in the integral used to evaluate the number of secondary electrons [Eq. (17a)] being prematurely cutoff. This can result in underestimating the number of secondary electrons generated. To account for this, we extended the energy integral in Eq. (17a) to infinity, where we take the value of the RPF to be $P ({\bar{p}}_{max}, ξ_{1})$ for the extended region of the integral. For cases well above threshold, we will have $P ({\bar{p}}_{max}, ξ_{1}) \approx 1$ (see the dashed white curves in Fig. 6, for example), and we would expect this quantity to remain near unity for $p > {\bar{p}}_{max}$ ⁠. In contrast, when near or below the threshold, $P ({\bar{p}}_{max}, ξ_{1})$ will no longer be unity, and will likely vary significantly at higher momentum. For this latter case, however, the high energy boundary will be located at several MeV (see Fig. 5, for example), where at such high energies, the Rosenbluth–Putvinski source has a negligibly small value. In this limit, the contribution from the high energy extension of the integral defined by Eq. (17a) will make a negligible contribution to the rate of RE growth or decay, such that our extension of the integral defined in Eq. (17a) past the ${\bar{p}}_{max}$ boundary will have little impact on the predicted number of secondary electrons.

Finally, as evident in Fig. 5(c), when well below the threshold, a spurious contribution to the RPF can emerge at the low energy boundary. To remove this contribution, we will note two momenta. The first is the critical momentum for an electron to overcome drag, defined by $p_{c} \equiv 1 / \sqrt{E_{| |} - 1}$ ⁠, for $E_{| |} \leq 1$ ⁠, or $p_{c} \to \infty$ otherwise. The second is a hybrid momentum formed by taking the weighted average of $p_{min}$ and $p_{R E}$ ⁠, i.e., $p_{ref} = (p_{R E} + 3 p_{min}) / 4$ ⁠. To remove the spurious contribution to the RPF near $p_{min}$ ⁠, we will multiply the integrand appearing in Eq. (17a) by a Heaviside function $Θ (p - p_{c})$ for $p_{c} < p_{ref}$ ⁠, thus removing electrons below the critical momentum, or $Θ (p - p_{ref})$ otherwise. In so doing, contributions from momenta below $p_{ref}$ will be removed when near or below the threshold, but when far above the threshold (i.e., when $p_{c} < p_{ref}$ ⁠), only momenta below $p_{c}$ are removed. Thus, in both limits, the spurious contribution to the RPF near $p = p_{min}$ is removed, but when well above the threshold, only very small values of momenta are impacted.

Using the above algorithm, the avalanche growth or decay of an initial seed RE population has been evaluated for cases both well above and below the threshold, with the result compared with RAMc simulations that employ a Møller secondary source (see Fig. 7). Good agreement is found both for a strongly decaying RE distribution along with a RE distribution that is rapidly growing. Considering a broader comparison between predictions from the PINN and RAMc, 50 values of the parameters $(E_{| |}, Z_{eff}, α)$ were randomly sampled over the range of parameters the PINN was trained across. Cases that decayed too rapidly for a good exponential fit to be identified, which corresponded to cases far below the threshold as discussed in Ref. 1, were discarded leaving 44 cases that will be compared with PINN predictions. The resulting comparison between the two approaches is shown in Fig. 8(a). Generally, good agreement is evident, suggesting that the PINN is able to describe the growth or decay of a RE population across a range of plasma conditions. We note that the rapid increase in the magnitude of the RE decay rate as the electric field drops below the threshold (i.e., the steep drop off evident in Fig. 7) resulted in only a few data points that were substantially below the threshold, but not so far below the threshold to prevent an exponential fit, being present in our randomly sampled dataset.

FIG. 8.

View large Download slide

Comparison of avalanche growth and decay rates between the Monte Carlo code RAMc and the PINN over a broad range of plasma conditions. Panel (a) includes a Rosenbluth–Putvinski secondary source term for all the PINN's predictions, whereas panel (b) uses the decay rate of the primary distribution when the PINN predicts RE decay. These latter points are highlighted in red. The electric field, effective charge $Z_{eff}$ ⁠, and synchrotron radiation $α$ were selected randomly over the region the PINN was trained. The Coulomb logarithm was taken to be $ln Λ = 9.998$ in this scan. The initial electron distribution was defined by Eq. (23) with $(p_{0} = 3 {\bar{p}}_{max} / 4, Δ p = {\bar{p}}_{max} / 10, ξ_{0} = - 1, Δ ξ = 0.1)$ ⁠.

The largest disagreement between the predictions of the PINN and RAMc occurs when below the threshold. As noted in Fig. 3, and the ensuing discussion, when below the threshold, the Rosenbluth–Putvinski source overpredicts the impact of secondary electron generation on the RE density, where a better approximation can be achieved by neglecting the secondary source and computing the decay rate of the primary distribution. This is done in Fig. 8(b), where the decay rate of the values (red “×” markers) that the PINN predicted were below the threshold, is now computed using the decay rate of the primary distribution, yielding improved agreement with predictions from RAMc, particularly for the strongest decaying cases.

VI. CONCLUSION

An adjoint-deep learning framework for evaluating the time evolution of the RE density incorporating large-angle collisions was derived. By utilizing a PINN to identify the parametric solution to the adjoint problem, this enables the time evolution of the RE density to be rapidly inferred across a broad range of parameters for an arbitrary initial momentum space distribution. This framework was applied to evaluate time histories of the RE density, along with saturated growth and decay rates across a broad range of physics parameters. Excellent agreement was found between predictions of the adjoint-deep learning framework and a traditional RE solver for the saturated avalanche growth rate. It was found that near and below marginality, the Rosenbluth–Putvinski approximation to the secondary source term lead to quantitatively inaccurate time histories of the RE density, though the qualitative behavior was captured. Furthermore, it was found that when below marginality, the decay rate of the RE distribution is well approximated by the decay rate of the primary distribution function (i.e., neglecting large-angle collisions). We anticipate that extending the present framework to include a more comprehensive set of physics, including time varying electric fields, temporally varying plasma compositions, and partial screening corrections to collisional coefficients, will allow the present approach to provide an efficient, yet accurate tool for evaluating saturated RE growth and decay rates, as well as describing RE dynamics during rapid variations in the background plasma. Finally, the present approach to describing RE generation can be extended to describe other RE generation mechanisms including nuclear seeds arising from tritium decay or Compton scattering,⁶ thus, providing a comprehensive treatment of RE formation across a range of plasma scenarios. Such extensions will be the topic of future work.

ACKNOWLEDGMENTS

This work was supported by the Department of Energy, Office of Fusion Energy Sciences, at the University of Florida under Award Nos. DE-SC0024649 and DE-SC0024634, and at Los Alamos National Laboratory (LANL) under Contract No. 89233218CNA000001. The authors acknowledge the University of Florida Research Computing for providing computational resources that have contributed to the research results reported in this publication. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a Department of Energy Office of Science User Facility, using NERSC Award No. FES-ERCAP0028155.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Christopher J. McDevitt: Conceptualization (lead); Data curation (equal); Formal analysis (lead); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Jonathan S. Arnaud: Conceptualization (equal); Data curation (lead); Software (equal); Visualization (equal); Writing – review & editing (equal). Xian-Zhu Tang: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Writing – review & editing (equal).

DATA AVAILABILITY

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

REFERENCES

1.

C. J.

McDevitt

et al, “

A physics-constrained deep learning treatment of runaway electron dynamics

,”

Phys. Plasmas

32

,

042503

(

2025

).

https://doi.org/10.1063/5.0253370

Google Scholar

2.

M.

Kruskal

and

I.

Bernstein

, Report No. Matt-Q-20 (

Princeton Plasma Physics Laboratory

,

1962

).

3.

J.

Connor

and

R.

Hastie

,

Nucl. Fusion

15

,

415

(

1975

).

https://doi.org/10.1088/0029-5515/15/3/007

Google Scholar

Crossref

4.

M.

Rosenbluth

and

S.

Putvinski

,

Nucl. Fusion

37

,

1355

(

1997

).

https://doi.org/10.1088/0029-5515/37/10/I03

Google Scholar

Crossref

5.

H. M.

Smith

and

E.

Verwichte

,

Phys. Plasmas

15

,

072502

(

2008

).

https://doi.org/10.1063/1.2949692

Google Scholar

Crossref

6.

J.

Martín-Solís

,

A.

Loarte

, and

M.

Lehnen

,

Nucl. Fusion

57

,

066025

(

2017

).

https://doi.org/10.1088/1741-4326/aa6939

Google Scholar

Crossref

7.

P.

Aleynikov

and

B. N.

Breizman

,

Phys. Rev. Lett.

114

,

155001

(

2015

).

https://doi.org/10.1103/PhysRevLett.114.155001

Google Scholar

Crossref

PubMed

8.

L.

Hesslow

,

O.

Embréus

,

O.

Vallhagen

, and

T.

Fülöp

,

Nucl. Fusion

59

,

084004

(

2019a

).

https://doi.org/10.1088/1741-4326/ab26c2

Google Scholar

Crossref

9.

C. J.

McDevitt

,

Z.

Guo

, and

X.-Z.

Tang

,

Plasma Phys. Controlled Fusion

60

,

024004

(

2018

).

https://doi.org/10.1088/1361-6587/aa9b3f

Google Scholar

Crossref

10.

L.

Hesslow

,

L.

Unnerfelt

,

O.

Vallhagen

,

O.

Embréus

,

M.

Hoppe

,

G.

Papp

, and

T.

Fülöp

,

J. Plasma Phys.

85

,

475850601

(

2019

).

https://doi.org/10.1017/S0022377819000874

Google Scholar

Crossref

11.

C. J.

McDevitt

,

Phys. Plasmas

30

,

092501

(

2023

).

https://doi.org/10.1063/5.0164712

Google Scholar

Crossref

12.

M.

Yang

,

P.

Wang

,

D.

del-Castillo-Negrete

,

Y.

Cao

, and

G.

Zhang

,

SIAM J. Sci. Comput.

46

,

C508

(

2024

).

https://doi.org/10.1137/23M1585635

Google Scholar

Crossref

13.

J. S.

Arnaud

,

T. B.

Mark

, and

C. J.

McDevitt

,

J. Plasma Phys.

90

,

905900409

(

2024

).

https://doi.org/10.1017/S0022377824000679

Google Scholar

Crossref

14.

R. W.

Harvey

,

V. S.

Chan

,

S. C.

Chiu

,

T. E.

Evans

,

M. N.

Rosenbluth

, and

D. G.

Whyte

,

Phys. Plasmas

7

,

4590

(

2000

).

https://doi.org/10.1063/1.1312816

Google Scholar

Crossref

15.

E.

Nilsson

,

J.

Decker

,

Y.

Peysson

,

R. S.

Granetz

,

F.

Saint-Laurent

, and

M.

Vlainic

,

Plasma Phys. Controlled Fusion

57

,

095006

(

2015

).

https://doi.org/10.1088/0741-3335/57/9/095006

Google Scholar

Crossref

16.

Z.

Guo

,

C.

Mcdevitt

, and

X.

Tang

,

Phys. Plasmas

26

,

082503

(

2019

).

https://doi.org/10.1063/1.5055874

Google Scholar

Crossref

17.

M.

Hoppe

,

O.

Embreus

, and

T.

Fülöp

,

Comput. Phys. Commun.

268

,

108098

(

2021

).

https://doi.org/10.1016/j.cpc.2021.108098

Google Scholar

Crossref

18.

J.

Rudi

,

M.

Heldman

,

E. M.

Constantinescu

,

Q.

Tang

, and

X.-Z.

Tang

,

J. Comput. Phys.

507

,

112954

(

2024

).

https://doi.org/10.1016/j.jcp.2024.112954

Google Scholar

Crossref

19.

L.-G.

Eriksson

and

P.

Helander

,

Comput. Phys. Commun.

154

,

175

(

2003

).

https://doi.org/10.1016/S0010-4655(03)00293-5

Google Scholar

Crossref

20.

C.

Sommariva

,

E.

Nardon

,

P.

Beyer

,

M.

Hoelzl

,

G. T. A.

Huijsmans

,

D.

van Vugt

, and

J.

Contributors

,

Nucl. Fusion

58

,

016043

(

2018

).

https://doi.org/10.1088/1741-4326/aa95cd

Google Scholar

Crossref

21.

C. J.

McDevitt

,

Z.

Guo

, and

X.-Z.

Tang

,

Plasma Phys. Controlled Fusion

61

,

054008

(

2019

).

https://doi.org/10.1088/1361-6587/ab0d6d

Google Scholar

Crossref

22.

M. T.

Beidler

,

D.

del-Castillo-Negrete

,

D.

Shiraki

,

L. R.

Baylor

,

E. M.

Hollmann

, and

C. J.

Lasnier

,

Nucl. Fusion

64

,

076038

(

2024

).

https://doi.org/10.1088/1741-4326/ad4c77

Google Scholar

Crossref

23.

I.

Sokolov

,

JETP Lett.

29

,

218

(

1979

).

Google Scholar

24.

T. C.

Hender

,

J. C.

Wesley

,

J.

Bialek

,

A.

Bondeson

,

A. H.

Boozer

,

R. J.

Buttery

,

A.

Garofalo

,

T. P.

Goodman

,

R. S.

Granetz

,

Y.

Gribov

et al,

Nucl. Fusion

47

,

S128

(

2007

).

https://doi.org/10.1088/0029-5515/47/6/S03

Google Scholar

Crossref

25.

B. N.

Breizman

,

P.

Aleynikov

,

E. M.

Hollmann

, and

M.

Lehnen

,

Nucl. Fusion

59

,

083001

(

2019

).

https://doi.org/10.1088/1741-4326/ab1822

Google Scholar

Crossref

26.

C.

Liu

,

D. P.

Brennan

,

A. H.

Boozer

, and

A.

Bhattacharjee

,

Plasma Phys. Controlled Fusion

59

,

024003

(

2017

).

https://doi.org/10.1088/1361-6587/59/2/024003

Google Scholar

Crossref

27.

L.

Hesslow

,

O.

Embréus

,

A.

Stahl

,

T. C.

DuBois

,

G.

Papp

,

S. L.

Newton

, and

T.

Fülöp

,

Phys. Rev. Lett.

118

,

255001

(

2017

).

https://doi.org/10.1103/PhysRevLett.118.255001

Google Scholar

Crossref

PubMed

28.

C. F.

Karney

and

N. J.

Fisch

,

Phys. Fluids

29

,

180

(

1986

).

https://doi.org/10.1063/1.865975

Google Scholar

Crossref

29.

G.

Zhang

and

D.

del-Castillo-Negrete

,

Phys. Plasmas

24

,

092511

(

2017

).

https://doi.org/10.1063/1.4986019

Google Scholar

Crossref

30.

C.

Møller

,

Ann. Phys.

406

,

531

(

1932

).

https://doi.org/10.1002/andp.19324060506

Google Scholar

Crossref

31.

A.

Ashkin

,

L. A.

Page

, and

W.

Woodward

,

Phys. Rev.

94

,

357

(

1954

).

https://doi.org/10.1103/PhysRev.94.357

Google Scholar

Crossref

32.

A. H.

Boozer

,

Phys. Plasmas

22

,

032504

(

2015

).

https://doi.org/10.1063/1.4913582

Google Scholar

Crossref

33.

M.

Raissi

,

P.

Perdikaris

, and

G. E.

Karniadakis

,

J. Comput. Phys.

378

,

686

(

2019

).

https://doi.org/10.1016/j.jcp.2018.10.045

Google Scholar

Crossref

34.

G. E.

Karniadakis

,

I. G.

Kevrekidis

,

L.

Lu

,

P.

Perdikaris

,

S.

Wang

, and

L.

Yang

,

Nat. Rev. Phys.

3

,

422

(

2021

).

https://doi.org/10.1038/s42254-021-00314-5

Google Scholar

Crossref

35.

L.

Lu

,

X.

Meng

,

Z.

Mao

, and

G. E.

Karniadakis

,

SIAM Rev.

63

,

208

(

2021

).

https://doi.org/10.1137/19M1274067

Google Scholar

Crossref

36.

M.

Abadi

,

P.

Barham

,

J.

Chen

,

Z.

Chen

,

A.

Davis

,

J.

Dean

,

M.

Devin

,

S.

Ghemawat

,

G.

Irving

,

M.

Isard

et al, in

12th USENIX Symposium on Operating Systems Design and Implementation (OSDI ‘16

) (

USENIX

,

2016

), pp.

265

–

283

.

37.

C.

Wu

,

M.

Zhu

,

Q.

Tan

,

Y.

Kartha

, and

L.

Lu

,

Comput. Methods Appl. Mech. Eng.

403

,

115671

(

2023

).

https://doi.org/10.1016/j.cma.2022.115671

Google Scholar

Crossref

38.

S.

Chiu

,

M.

Rosenbluth

,

R.

Harvey

, and

V.

Chan

,

Nucl. Fusion

38

,

1711

(

1998

).

https://doi.org/10.1088/0029-5515/38/11/309

Google Scholar

Crossref

39.

C.

McDevitt

and

X.-Z.

Tang

,

Europhys. Lett.

127

,

45001

(

2019

).

https://doi.org/10.1209/0295-5075/127/45001

Google Scholar

Crossref

40.

J. S.

Arnaud

and

C. J.

McDevitt

,

Phys. Plasmas

31

,

062509

(

2024

).

https://doi.org/10.1063/5.0198338

Google Scholar

Crossref

41.

C. J.

McDevitt

,

Z.

Guo

, and

X.-Z.

Tang

,

Plasma Phys. Controlled Fusion

61

,

024004

(

2019

).

https://doi.org/10.1088/1361-6587/aaf4d1

Google Scholar

Crossref

This work is licensed under a Creative Commons Attribution 4.0 International License .

An efficient surrogate model of secondary electron formation and evolution

I. INTRODUCTION

II. ADJOINT FRAMEWORK FOR RELATIVISTIC ELECTRON EVOLUTION

A. The adjoint of the relativistic Fokker–Planck equation

B. Runaway electron avalanche

III. PHYSICS-INFORMED NEURAL NETWORKS

IV. TEMPORAL EVOLUTION OF AN AVALANCHING RUNAWAY ELECTRON POPULATION

V. PARAMETRIC VARIATION OF RUNAWAY ELECTRON AVALANCHE GROWTH AND DECAY RATE

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

An efficient surrogate model of secondary electron formation and evolution Open Access

I. INTRODUCTION

II. ADJOINT FRAMEWORK FOR RELATIVISTIC ELECTRON EVOLUTION

A. The adjoint of the relativistic Fokker–Planck equation

B. Runaway electron avalanche

III. PHYSICS-INFORMED NEURAL NETWORKS

IV. TEMPORAL EVOLUTION OF AN AVALANCHING RUNAWAY ELECTRON POPULATION

V. PARAMETRIC VARIATION OF RUNAWAY ELECTRON AVALANCHE GROWTH AND DECAY RATE

VI. CONCLUSION

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

An efficient surrogate model of secondary electron formation and evolution