A physics-constrained deep learning treatment of runaway electron dynamics Open Access

The treatment of runaway electron (RE) dynamics is a critical component in the integrated description of tokamak disruptions^1,2 and startup scenarios.^3,4 While first-principle kinetic solvers are available to evaluate the transient RE dynamics during these phases of a tokamak discharge,^5–9 such approaches are computationally intensive, thus, sharply limiting their application. This has led to the development of several reduced models of RE formation that may be conveniently coupled to magnetohydrodynamic (MHD) codes.^10–12 These reduced models, however, are often based on the assumption of a slowly evolving background plasma and are, thus, unable to accurately capture highly transient phases of RE generation or decay. An exception to this trend includes reduced models of hot tail RE seed formation,^13–15 where RE formation during a rapid thermal quench is estimated. Such models, however, have to date been developed specifically for the thermal quench, and, thus, are not applicable to other transient phases of a disruption, and do not self-consistently include the presence of large-angle collisions. The amplification of the RE seed by large-angle collisions, the so-called RE avalanche,¹⁶ is instead incorporated by a separate module^17–19 that does not directly account for the transiently evolving momentum space distribution of REs.

The present work, together with Paper II,²⁰ aims to address these limitations by employing an adjoint-deep learning framework that will allow the transient dynamics of REs to be evolved with an online execution time that is orders of magnitude less than a traditional relativistic Fokker–Planck solver. To achieve this, we will first focus on the solution to an adjoint problem,^14,21–23 rather than directly solving the relativistic Fokker–Planck equation. A distinguishing feature of this approach is that a single solution of the adjoint problem can be used to evolve the RE density forward in time for an arbitrary initial momentum space distribution. Thus, once a solution to the adjoint problem has been identified, the impact of distinct initial momentum space distributions on RE density evolution can be rapidly evaluated for a given set of physics parameters. The second component of the framework will be the use of a physics-informed neural network (PINN)²⁴ to identify a parametric solution to the adjoint of the relativistic Fokker–Planck equation. Physics-informed neural networks are a form of deep learning, whereby physical constraints are embedded in the training of a neural network (NN). While this approach has gained prominence due to its ability to integrate physical constraints with data,²⁵ or treat inverse problems,²⁶ for cases where a closed system of equations is available, PINNs can be used to learn the solution to a PDE across a broad range of physics parameters in the absence of data.^27,28 Further noting that while the offline training time of a PINN is long, its online inference time is orders of magnitude faster than a traditional solver. The PINN will thus provide a rapid surrogate of the physical system of interest. Combining the above two components, the resulting adjoint-deep learning framework will be able to generate time histories of the RE number density at a fraction of the cost of a traditional solver, for an arbitrary initial momentum space distribution of REs and across a range of plasma conditions. The present paper, along with Paper II,²⁰ will focus on developing the fundamental components of this adjoint-deep learning framework, while future works will focus on employing it for describing tokamak disruptions.

As specific example problems, we consider the decay of a primary RE distribution, while Paper II²⁰ will treat the avalanche amplification of an initial RE seed. While large-angle collisions are essential to providing a description of RE generation during the current quench phase of a disruption, an accurate treatment of RE decay is essential when describing the decay of plasma current²⁹ during the RE plateau. This decay rate will be inferred as a function of electric field strength, effective charge, and the strength of synchrotron radiation. The RE decay rate is shown to have a strong nonlinear dependence on the strength of the electric field, in contrast to commonly employed analytic expressions for RE growth and decay such as Ref. 17, which indicate a linear dependence. Furthermore, we find that characterizing the decay of REs using a simple decay rate, as is often done in current applications, is often inaccurate either due to the long time required for the RE distribution to settle into a well defined exponential decay or due to the RE density exhibiting temporal dynamics that cannot be fit by an exponential, even in the long time limit. While the present paper does not include large-angle collisions, Paper II²⁰ finds that the decay of the overall RE distribution is well approximated by the decay of the primary population neglecting large-angle collisions. The present estimate, thus, provides an efficient means of evaluating the decay of a RE population when below threshold. As a further application, time histories of the RE density evolution are computed for distinct initial momentum space distributions of REs. It is shown that the adjoint-deep learning framework is able to accurately resolve these time histories for RE distributions with sharply different initial pitch distributions. The specific time histories predicted by the adjoint-deep learning framework are compared with those predicted by the particle-based RE solver RAMc,⁷ with good agreement evident across a range of initial momentum space distributions.

The remainder of this paper is organized as follows. Section II discusses the adjoint formulation of the relativistic Fokker–Planck equation. A brief description of PINNs is provided in Sec. III, along with a discussion of the tailored NN architecture employed in this paper. The temporal evolution of the number density of primary electrons for distinct initial momentum space distributions of REs is given in Sec. IV, along with evaluating RE decay for a range of plasma parameters. A brief discussion along with conclusions are given in Sec. V.

The physical interpretation of the RPF can be made clear by considering the orbits of several electrons initialized at different pitches as shown in Fig. 1. These orbits have been computed using characteristic equations derived from Eq. (1), where pitch-angle scattering was neglected leading to deterministic orbits. The parameters (⁠ $E∥=2$ and $α=0.1$⁠) were selected such that the system is slightly above the avalanche threshold, where effects from synchrotron radiation are most pronounced. After integrating forward in time until $t=τc$ [see Fig. 1(a)], all of the electron orbits have slowed down, except an orbit with an initial pitch of $ξ=−0.99$⁠. This last orbit was initialized with a pitch that put it in a region with $Up>0$ (below the dashed black curve) allowing it to be accelerated. It is useful to note that the rate at which electrons in the $Up<0$ region decelerate depends sensitively on their pitch, with electrons near $ξ≈0$ decelerating the most rapidly due to synchrotron radiation being maximal in this region, where the presence of synchrotron radiation results in the region with $Up>0$ contracting as p increases. At a later time [ $t=2τc$ Fig. 1(b)], the majority of electron orbits continue to slow down, but are also tilted toward the negative pitch region. This is due to the electric field pushing the electrons toward the $ξ=−1$ axis. Finally, by $t=10τc$⁠, all of the electrons have either been accelerated through the $p=10$ surface, or decelerated to the bulk (i.e., $p=0$ since we have assumed a cold bulk plasma). Those that have been accelerated through the $p=10$ surface [which could be taken to be the upper momentum boundary] will be treated as REs according to Eq. (10). In contrast, two orbits in Fig. 1(c), the purple and red orbits initialized with $ξ>0.8$ have fallen back to the bulk, and are, thus, deemed to not run away. Thus, for the case where pitch-angle scattering is neglected, a separatrix will form delineating runaway and non-runaway orbits.

The primary impact of pitch-angle scattering on the RPF will be to broaden the sharp transition region and shift its location to higher energy. This is illustrated in Fig. 2, where we have computed the RPF from the steady state adjoint equation, but artificially multiplied the pitch-angle scattering term by factors of 0.01 [Fig. 2(a)], 0.1 [Fig. 2(b)], and 1 [Fig. 2(c)] for a plasma with $Zeff=1$⁠, while keeping all other terms fixed, including the collisional drag. It is evident that for Fig. 2(a) (where the pitch-angle scattering term is multiplied by 0.01), a sharp transition region is present, where the region with $P≈1$ closely aligns with the RE orbits shown in Fig. 1(c). However, as the strength of pitch-angle scattering is increased in Figs. 2(b) and 2(c), the transition region broadens, and the location of the $P=0.5$ contour shifts to higher energy, implying that RE generation will be less efficient as pitch-angle scattering is increased.

The adjoint problem discussed in Sec. II allows the RE density to be evolved in time beginning with an arbitrary initial momentum space distribution of REs. Our aim in this section will be to extend the generalizability of the adjoint framework by utilizing a PINN to evaluate the parametric solution to Eq. (5). In this way, the resulting adjoint-deep learning framework will be able to evaluate the RE density for both an arbitrary initial electron distribution, and a range of plasma conditions. Once trained, the resulting framework will allow for nearly instantaneous projections of RE density incorporating fully kinetic physics.

For the present study, we will consider the limit where no data are available, such that we will rely exclusively on physical constraints for training the PINN. A key component when developing a robust PINN will be to enforce several properties of the solution as hard constraints in the NN architecture itself, rather than as soft constraints in the loss function. This prevents the PINN from making unphysical predictions while simultaneously, and perhaps more importantly, severely restricting the solutions that the optimizer searches across, leading to more robust convergence of the PINN. These hard constraints will be enforced by a customized “physics layer,” detailed in Sec. III B below.

As a specific application of the above-mentioned framework, we will treat the decay of an initial primary RE distribution. Here, an initial population of REs will be initialized with a given momentum space distribution and the adjoint-deep learning framework described in Secs. II and III will be used to project the RE density forward in time. We will begin by describing representative solutions to the adjoint of the relativistic Fokker–Planck equation (Sec. IV A) and subsequently describe the decay of a primary RE population for a range of initial conditions (Sec. IV B). Finally, the predictions of the adjoint-deep learning framework will be compared with a traditional RE solver in Sec. IV C.

Using the NN architecture and loss function described in Sec. III B above, example solutions to the adjoint problem are shown in Figs. 5 and 7 with the test and training loss history shown in Fig. 4. Here, we have trained the PINN for electric fields in the range $E∥∈(0,3)$⁠, effective charges $Zeff∈(1,2)$⁠, and synchrotron radiation $α∈(0,0.1)$⁠. A much broader range of physics parameters will be considered in Ref. 20. We have also chosen a domain with a minimum energy of $10 keV$⁠, a maximum energy of $5 MeV$⁠, and $tfinal=20$⁠. The energy defining which electrons are considered to be REs is taken to be $pRE=pmax/4$⁠, which corresponds to an energy of roughly 1 MeV. The loss history of the PINN is shown in Fig. 4. Here, the first 15 000 epochs are performed using the ADAM optimizer, while L-BFGS is used thereafter. After periods of 50 000 epochs of training with L-BFGS, the training point distribution is resampled using the residual based adaptive distribution (RAD) method described in Ref. 41. Such a sampling routine redistributes roughly half of the training points using a uniform random distribution, with the remaining points selected with a probability proportional to the magnitude of the residual. In this way, the entire training domain is sampled, but with regions that have relatively large residuals receiving more training points [the final training point distribution is shown in Fig. 4(b)]. This training point resampling results in spikes in the training loss every 50 000 epochs as evident from Fig. 7(a). Despite these spikes in the training loss of the PDE, the test loss of the PDE monotonically decays to a value of $≈10−7$⁠, with the boundary loss term dropping to a lower value, suggesting that the PDE is well satisfied across the range of parameters under consideration.

Two example solutions are shown in Figs. 5 and 7 for electric fields below and above the avalanche threshold (the avalanche threshold is $E∥≈1.8$ for these parameters). Considering the subthreshold case first (⁠ $E∥=1.5$⁠, Fig. 5), the $P=0.5$ contour (black contour in the top row of Fig. 5) shifts to higher momentum and smaller pitch when moving from the $t=20$ to $t=0$ time slice. Physically, this implies that an electron requires an energy near 5 MeV and a pitch near $ξ=−1$ at $t=0$ to have a better than fifty percent chance of remaining in the RE region or being accelerated through the high momentum boundary, by $t=tfinal$⁠. This is due to the combination of drag and synchrotron losses exceeding electric field acceleration throughout the majority of momentum space, with only a narrow channel where $Up>0$ (below the dashed black contour in the bottom row of Fig. 5). Even inside this narrow channel, the presence of pitch-angle scattering will result in electrons initially inside the $Up>0$ region being scattered to a larger pitch-angle where $Up<0$⁠, such that they are slowed down^42–45 resulting in the RPF taking on a small value across nearly all of momentum space. More insight into the accuracy of the solution in regions where $P≪1$ can be gained by considering the $log10$ of both the RPF and the residual (see Fig. 6). Here, it is evident that the magnitude of the residual of the RPF is small throughout the momentum space, with the magnitude of the residual normalized to the RPF [Fig. 6(c)] being small across regions with significant values of the RPF. We note that at low energies this ratio can exceed one, however, the RPF has already decayed to a value of the order of $10−5$ at these energies, such that the contribution to the RE density will nearly always be negligible from this region.

While the magnitude of the residual in Fig. 5 is small for the three time slices shown, suggesting an accurate solution, the final time slice [ $t=0$⁠, Fig. 7(c)] suggests a limitation of the PINN approach employed. In particular, the high momentum boundary condition employed is $P(p=pmax)=1$ when $Up>0$ [see Eq. (9c)]. This is not satisfied at $t=0$⁠. Specifically, as evident from Fig. 5(f), a narrow region near $ξ=−1$ is present where $Up>0$⁠. However, the RPF is not unity at the high momentum boundary in this narrow channel. The reason for this inconsistency is due to pitch-angle scattering quickly diffusing electrons out of the $Up>0$ channel where they are subsequently slowed down. Thus, the magnitude of the RPF will quickly drop below unity for $p<pmax$ even when $Up>0$⁠. Since the present PINN implements the high momentum boundary condition as a soft constraint, deviation from the prescribed boundary condition in this narrow boundary region only leads to a modest penalty in the boundary loss term. We also note that the region where the RPF deviates from the boundary condition (⁠ $p=pmax$⁠, $ξ≈−1$ and $t=tfinal$⁠) represents a corner in the training regime, it is, thus, likely that it is not directly sampled by the random distribution of training points. While this is indicative of a slight inconsistency in the PINN solution, we do not expect it to strongly impact predictions of the RE density since the PINN correctly captures the value of the RPF throughout the majority of the simulation domain, with only a small inconsistency present in the corner of the simulation domain when near marginality. This will be verified in Sec. IV B, where predictions of the PINN are compared with solutions from a traditional RE solver (see Fig. 13).

Considering now a case above threshold (⁠ $E∥=2$⁠, Fig. 7), the channel with $Up>0$ is substantially broadened, facilitating the acceleration of electrons. In particular, the location of the $P=0.5$ contour near $ξ=−1$ shifts to lower momenta as time evolves, indicating more electrons will be accelerated by the electric field. Since the present case is only modestly above marginality, pitch-angle scattering introduces a substantial broadening of the transition region. In addition, the RPF reaches a near steady state condition after ten $τc$⁠, such that there is only a modest difference between the RPF at $t=10$ compared to $t=20$⁠.

The time evolution of the RE density can be used to estimate the decay rate of a primary electron population when below threshold. Using a RE decay rate to characterize RE evolution comes with two significant caveats. The first is that while the RE density exhibits an exponential decay when slightly below threshold in the long time limit,⁴⁵ its temporal evolution cannot be described by an exponential when the system is well below threshold. As a concrete example, for an electric field $E∥<1$⁠, the entire population of electrons will slow down due to collisional drag exceeding electric field acceleration at all energies and pitches. This will result in the primary population of electrons slowing down to the bulk in a finite amount of time, a process that cannot be described by a decaying exponential. This physics is illustrated in Fig. 9(a), where the temporal evolution of the RE density for electric fields in the vicinity of $E∥=1$ for $Zeff=1$ and $α=0$ is shown. It is evident that for $E∥<1$ (green curve) an exponential curve cannot be fit to the RE density, as evidenced by the poor fit between the dashed and solid green curves. Furthermore, when at threshold (⁠ $E∥=1$⁠, red curves) the exponential fit curve only approximates the latter portion of the time history of the RE density. Hence, a second caveat is that it often takes a substantial amount of time, tens to hundreds of $τc$⁠,⁴⁵ for the REs to settle into a distribution with a clear exponential decay rate. Further noting that for a plasma with $ne=1020 m−3$ and a Coulomb logarithm of $ln Λ=15$⁠, the relativistic collision time is given by $τc≈22 ms$⁠, implying that the RE density often will not achieve an exponential decay for tens to hundreds of milliseconds even for $E∥≳1$⁠. Considering a case with $Zeff=2$ and $α=0.1$ [Fig. 9(b)], and, thus, a larger avalanche threshold of $Eav≈1.8$⁠, a good exponential fit is not present for $E∥=1$ over twenty $τc$⁠, with only a marginally better fit present for $E∥=1.5$⁠.

The time required for the RE distribution to settle into a well defined exponential decay depends sensitively on the initial momentum space distribution. Specifically, while the inferred RE decay rate should be independent of the initial momentum space distribution, assuming the initial electrons are at sufficiently high energy, there will inevitably be a transient period while the RE momentum space distribution settles into its saturated form. To assess this uncertainty, a scan of the initial momentum and pitch of the primary electron distributions is shown in Figs. 10(a) and 10(b), respectively. Here, it is evident that the early time evolution of the RE population varies substantially depending on both the initial momentum and pitch of the primary distribution, with resulting fitted exponentially implying a similar, though not identical, decay rate. The modest deviation in inferred decay rates is due to the finite time of the simulations, such that for initial conditions far from those that characterize the decaying primary distribution, twenty $τc$ is not sufficient for the REs to settle into their saturated decay rate. A strength of the present approach is that by predicting the specific time history of the RE density, rather than only the decay rate, the accuracy of the exponential fit can be readily verified, and the specific time history can be used if required.

Insight into the accuracy of the exponential decay rate near threshold can be gained by comparing with RE decay rates inferred in Ref. 45. Here, the RE solver RAMc⁴⁶ was used to evolve the primary RE distribution for $250τc$⁠, leading to well converged exponential decay rates. A comparison between the PINN and results from Ref. 45 are shown by the blue “×” markers in Figs. 11(a) and 12. The solid blue curve in both figures indicates results from a PINN trained with $tfinal=20τc$⁠, whereas the dashed blue curve in Fig. 12 being an analogous PINN, but trained over a time period of $tfinal=40τc$⁠. Twice as many training points were used in this latter case (2 000 000) resulting in an increase in offline training time. It is evident that both results are in good agreement with the benchmark results, though the PINN with the larger $tfinal$ exhibits modestly improved agreement.

As a means of verifying predictions of the adjoint-deep learning framework we will compute several time histories of the RE density using the Monte Carlo code RAMc⁴⁶ for several different initial momentum space distributions and physical parameters. While a comparison of inferred RE decay rates could readily be performed, due to the ambiguities that often emerge when inferring the RE decay rate from the RE time history (as described in Sec. IV B 2), a more stringent test will be to compare detailed time histories of the RE density. A comparison of several temporal histories of RE density with the Monte Carlo code RAMc⁴⁶ is illustrated in Fig. 13. Here, we have used the initial momentum space distributions illustrated by the second row of Fig. 13, and the values of $α$ and $Zeff$ are the same as the RPFs shown in Fig. 5 above. The form of these momentum space distributions were chosen to be consistent with the particle initialization scheme used in RAMc. Good agreement is evident for the three electric fields and three initial conditions considered. For the practically important case of an initial RE distribution approximately aligned with the magnetic field (⁠ $ξ=−1$⁠) the adjoint-deep learning predictions are in excellent agreement with the RAMc code. For cases with electrons further from $ξ=−1$⁠, progressively larger deviations are present, though good overall agreement is evident with the largest deviation appearing for the lowest electric field with electrons initialized with $ξ=1$ [see Fig. 13(c) (red curves)]. The origin of this deviation is due to the limited accuracy of the PINN's prediction of the RPF in regions where the RPF is nearly vanishing. In particular, the small number of REs evident in the solid red curve in Fig. 13(c) results from integrating the product of the RPF shown in Fig. 5(c) with the initial momentum space distribution shown in Fig. 13(f). Since the RPF is nearly zero in the region where the initial electron distribution is finite (i.e., at $ξ≈1$ and $E≈3.5 MeV$⁠), this, thus, requires the RPF to be computed very accurately to resolve this near zero value. As evident from Fig. 13(c), the PINN is able to accurately track the initial density decay, but exhibits a quantitative discrepancy once only a trace RE population remains.