3.1. Adjoint relativistic Fokker–Planck equation

This section describes the implementation of the adjoint to the relativistic Fokker–Planck equation. While a direct evaluation of the particle distribution via a PINN is possible (see McDevitt & Tang (Reference McDevitt and Tang2024) for example), obtaining an accurate description of the runaway electron population is particularly challenging due to the very small number of electrons expected to run away. This electron population, typically orders of magnitude less than the bulk electron population, will provide a negligible contribution to the loss function defined by (2.1) unless special care is taken. An alternative means of describing REs is via the solution of the adjoint to the relativistic Fokker–Planck equation. This approach has been pursued in a variety of contexts, including wave-driven currents in magnetized fusion plasmas (Antonsen & Chu Reference Antonsen and Chu1982; Taguchi Reference Taguchi1983; Fisch Reference Fisch1986; Karney & Fisch Reference Karney and Fisch1986) and RE formation (Karney & Fisch Reference Karney and Fisch1986; Liu et al. Reference Liu, Brennan, Bhattacharjee and Boozer2016, Reference Liu, Brennan, Boozer and Bhattacharjee2017; Zhang & del Castillo-Negrete Reference Zhang and del Castillo-Negrete2017; McDevitt Reference McDevitt2023). As will be discussed below, the adjoint problem treated in this analysis describes the probability $P$ of an electron at an initial momentum space location running away at a later time. Thus, its range will be $P\in (0,1)$ and provides a quantity well suited for the optimization problem posed in (2.1). In the absence of synchrotron radiation, the adjoint to the steady-state relativistic Fokker–Planck equation can be written as

(3.1a)\begin{equation} \left[{-}E_\Vert \xi - C_F \right]\frac{\partial P }{\partial p} - \left( 1-\xi^2\right)\frac{E_\Vert}{p} \frac{\partial P }{\partial \xi} ={-}\frac{\nu_D}{2} \frac{\partial}{\partial \xi} \left[ \left( 1-\xi^2\right) \frac{\partial P }{\partial \xi}\right] , \end{equation}

where the collisional coefficients are taken to be

(3.1b,c)\begin{equation} \nu_D = \left( 1+Z_{{\rm eff}} \right)\frac{\gamma}{p^3}, \quad C_F = \frac{1 + p^2}{p^2} . \end{equation}

Here, the relativistic momentum $p$ is normalized to the electron mass $m_e$ and speed of light $c$ $p\to p/ ( m_e c)$, the electron's pitch is defined by $\xi \equiv p_\Vert / p$, time is normalized as $t \to t / \tau _c$, where $\tau _c \equiv 4{\rm \pi} \epsilon ^2_0 m^2_e c^3/ ( e^4 n_e \ln \varLambda )$ is the collision time of a relativistic electron for a given free electron density $n_e$ and Coulomb logarithm $\ln\varLambda$ with the physical constants of the permittivity of free space $\epsilon_0$ and elementary charge $e$, the collisional coefficients $\nu _D$ and $C_F$ are normalized to $\tau _c$ and the parallel electric field is normalized to the Connor–Hastie electric field $E_\Vert \to E_\Vert /E_c$, where $E_c \equiv m_e c / ( e \tau _c )$ (Connor & Hastie Reference Connor and Hastie1975). Energy diffusion has been neglected due to this term being exceptionally small for conditions typical of a tokamak disruption. In particular, noting that the energy diffusivity scales with the electron temperature $T_e/( m_e c^2 )$ for the high energies characteristic of REs, we expect this approximation to be well satisfied for a low temperature post-thermal quench plasma. Furthermore, the collision frequencies used in (3.1) assume the limit $v>v_{{\rm Te}}$, where $v_{{\rm Te}}$ is the electron thermal velocity. For the parameters of interest, the critical speed for an electron to run away will be much larger than $v_{{\rm Te}}$, hence we anticipate that this approximation will not substantially impact our results.

Concerning the boundary conditions, we will enforce $P=0$ at $p=p_{\min }$ and $P=1$ along the upper boundary, where the energy flux $U_p \equiv -E_\Vert \xi - C_F$ is positive. Specifically, the high energy boundary condition enforces that $P$ is unity for values of the pitch $\xi$ where electric field acceleration exceeds collisional drag such that the electron will be accelerated out of the simulation domain. Further noting that for values of the pitch near $| \xi | \approx 1$, (3.1a) is nearly hyperbolic, we anticipate that electrons accelerated through the high energy boundary will be accelerated to arbitrarily high energies, and thus have a zero probability of returning to the simulation domain. Hence, any electron located at the high energy boundary $p_{\max }$ with $U_p > 0$, is treated as a RE. With these boundary conditions, the quantity $P ( p, \xi )$ indicates the probability that an electron initially located at $( p, \xi )$ will run away at a later time, and is often referred to as the runaway probability function (RPF) (Karney & Fisch Reference Karney and Fisch1986).

When near the threshold electric field for avalanche generation, synchrotron radiation substantially impacts the RPF. While adding synchrotron radiation leads to a modest modification of (3.1), it does complicate the physical interpretation of the RPF. In particular, with the inclusion of synchrotron radiation, the adjoint equation takes the form

(3.2)\begin{align} & \left[{-}E_\Vert \xi - C_F - \alpha \gamma p \left( 1 - \xi^2 \right) \right]\frac{\partial P }{\partial p} + \left( 1-\xi^2\right) \left[ -\frac{E_\Vert}{p} + \alpha \frac{\xi}{\gamma} \right] \frac{\partial P }{\partial \xi}\nonumber\\ & \quad ={-}\frac{\nu_D}{2} \frac{\partial}{\partial \xi} \left[ \left( 1-\xi^2\right) \frac{\partial P }{\partial \xi}\right] , \end{align}

where the strength of synchrotron radiation is set by the parameter $\alpha \equiv \tau _c/\tau _s$, where $\tau _s \equiv 6 {\rm \pi}\epsilon _0 m^3_e c^3 / ( e^4 B^2 )$. Here, (3.2) is solved similarly as before with a boundary condition of $P=0$ at $p=p_{\min }$ along with the condition that $P=1$ on the high energy boundary $p=p_{\max }$ when the energy flux, now defined by $U^{(\alpha )}_p \equiv -E_\Vert \xi - C_F - \alpha \gamma p ( 1 - \xi ^2 )$, is positive. While this problem formulation is directly analogous to the case that neglects synchrotron radiation, the physical interpretation of the RPF is slightly modified. Specifically, as shown in Andersson, Helander & Eriksson (Reference Andersson, Helander and Eriksson2001), Decker et al. (Reference Decker, Hirvijoki, Embreus, Peysson, Stahl, Pusztai and Fülöp2016), Guo, McDevitt & Tang (Reference Guo, McDevitt and Tang2017) and McDevitt, Guo & Tang (Reference McDevitt, Guo and Tang2018), synchrotron radiation damping and pitch-angle scattering results in electrons obtaining a saturated energy, achieved via the formation of a circulation pattern in momentum space centred about an O-point. Thus, an electron accelerated through the high energy boundary will have a non-zero probability of returning to the simulation domain after a finite time, rather than being accelerated to arbitrarily high energy, as was the case when synchrotron radiation was neglected. For the case with synchrotron radiation, the RPF should thus be given the narrower interpretation as the probability that an electron reaches the high energy boundary of the simulation domain before slowing down to the low energy bulk, rather than the probability of an electron being accelerated to arbitrarily high energy. For this analysis the high energy boundary was chosen to be 5 MeV, which is a sufficiently high energy to capture the critical energy for an electron to run away in all cases except when asymptotically close to threshold.

3.2. Embedding physical constraints into the PINN

Our aim in this section will be to develop a PINN framework capable of robustly representing solutions to (3.2) across the three-dimensional parameter space $( E_\Vert, Z_{{\rm eff}}, \alpha )$. Here, the inputs of the PINN will be $(p,\xi,E_\Vert,Z_{{\rm eff}},\alpha )$ and the output will be the RPF $P$. A key component of our description is enforcing a subset of physical properties of the RPF as hard constraints. In particular, we will (i) enforce the low energy boundary condition $P=0$ at $p=p_{\min }$, (ii) constrain the RPF to have a range between zero and one and (iii) ensure the RPF vanishes when $| E_\Vert | < 1$. These three constraints are enforced by introducing a customized output layer to the neural network of the form

(3.3a)\begin{gather} P^\prime \left( p, \xi \right) \equiv \frac{1}{2} \left[ 1 + \tanh \left( \frac{E_\Vert - 1}{\Delta E}\right)\right] \left( \frac{p-p_{\min}}{p_{\max}-p_{\min}} \right) P_{{\rm NN}} \left( p, \xi \right) , \end{gather}

(3.3b)\begin{gather}P \left( p, \xi \right) \equiv \tanh\left( {P^\prime}^2 \left( p, \xi \right) \right) . \end{gather}

Here, $P_{{\rm NN}}$ is the output of the hidden layers of the neural network, and $\Delta E$ is a hyperparameter whose value should satisfy $\Delta E < 1$, where for all cases in this paper it is taken to be $\Delta E = 0.1$. From (3.3) it can be verified that, regardless of the value of $P_{{\rm NN}}$, the predicted RPF $P ( p, \xi )$: (i) vanishes at the low energy boundary, (ii) has a range of zero to one and (iii) tends to zero for $E_\Vert < 1$ and $\Delta E \ll 1$ (only positive electric fields are considered when training the RPF PINN). We note in passing that the output layer defined by (3.3) results in both the value and the first derivative of $P$ vanishing at $p=p_{\min }$. While this latter condition is not required when defining the RPF, it will nevertheless be satisfied as long as the value of $p_{\min }$ chosen is well below the momentum $p_{{\rm RE}} \equiv 1/\sqrt {E_\Vert -1}$ where collisional drag exceeds electric field acceleration. Specifically, since energy diffusion is neglected in the present analysis, electrons located below $p_{{\rm RE}}$ will have zero probability of running away. Thus, for $p < p_{{\rm RE}}$ the solution of the RPF will be a constant with a magnitude of zero, and thus is consistent with pure Dirichlet and Neumann boundary conditions at $p=p_{\min }$.

The loss function employed is taken to have the form

(3.4)\begin{equation} \text{Loss} = \frac{1}{N_{{\rm PDE}}} \sum^{N_{{\rm PDE}}}_i \left[\left( \frac{p^2_i}{1+p^2_i}\right)\mathcal{R} \left( p_i, \xi_i ; \boldsymbol{\lambda}_i \right) \right]^2 + \frac{1}{N_{{\rm bdy}}} \sum^{N_{{\rm bdy}}}_i \left[ 1 - P \left( p_i, \xi_i ; \boldsymbol{\lambda}_i \right)\right]^2 , \end{equation}

where $\mathcal {R} ( p_i, \xi _i ; \boldsymbol {\lambda }_i )$ is the residual of (3.2) and $\boldsymbol {\lambda }$ represents the physics parameters $( E_\Vert, Z_{{\rm eff}}, \alpha )$. Here, the first term in (3.4) penalizes deviations from the governing PDE given by (3.2), where the prefactor $p^2_i/( 1+p^2_i)$ removes the low energy divergence of the pitch-angle scattering operator. The second term in the loss function defined by (3.4) enforces the high energy boundary condition, i.e. $P=1$ at $p_{\max }$ when $U^{(\alpha )}_p > 0$. In particular, the boundary points $N_{{\rm bdy}}$, will only be applied at locations that satisfy both $U^{(\alpha )}_p > 0$ and $p=p_{\max }$.

3.3. Parametric dependence of the runaway probability function

In this section, we will seek to obtain solutions to the PINN in the five-dimensional space defined by the two independent coordinates $( p, \xi )$ and the three physics parameters $( E_\Vert, Z_{{\rm eff}}, \alpha )$. The loss history of the PINN is shown in figure 1. Here, after roughly 100 000 epochs, the loss associated with the PDE drops below $10^{-7}$, along with the boundary loss dropping to ${\approx }10^{-9}$, with the loss dropping more slowly for the remaining ${\sim }400\,000$ epochs. The test loss of the PDE reaches ${\approx }10^{-9}$ by the end of the training, indicating that an accurate solution was found. The periodic spikes in the training data are due to additional training points sampled after every 100 000 epochs at locations where the residual is maximal. The test loss is only updated after every 100 000 epochs when using the L-BFGS optimizer, leading to the sharp variations evident in the dashed curves.

Figure 1. Loss history (blue lines for the PDE and red lines for the boundary condition (BC)) for a feedforward neural network with six hidden layers each with a width of 64 neurons, along with roughly 1 000 000 training points. Here, 15 000 epochs were performed with the ADAM optimizer, with the remaining steps performed using L-BFGS. The model was trained across $E_\Vert \in (1,10)$, $Z_{{\rm eff}} \in (1,10)$ and $\alpha \in (0, 0.2)$. The range of $p$ was chosen such that the low energy boundary was $10\ \text {keV}$ and the high energy boundary was $5\ \text {MeV}$.

Four example predictions of the RPF are shown in figure 2. Here, the RPF vanishes at low energies, where drag exceeds electric field acceleration, but increases at higher energy due to the $( 1+p^2)/p^2$ drop in the collisional drag. In particular, for the parameters indicated in figure 2(a), the $P=0.5$ contour is located at an approximate energy of $200\ \text {keV}$ for $\xi =-1$. Considering a case with a large $Z_{{\rm eff}}$ (see figure 2(b) with $Z_{{\rm eff}}=10$), the $P=0.5$ contour shifts to higher energy, with a more gradual transition between the $P\approx 0$ and $P\approx 1$ regions. The impact of synchrotron radiation on the RPF is shown in figure 2(c), where the magnitude of synchrotron radiation was taken to be $\alpha = 0.2$. Compared with an otherwise identical case, but without synchrotron radiation (figure 2a), it is evident that the location of the $P =0.5$ contour at $\xi =-1$ has only shifted slightly. This is due to synchrotron radiation vanishing for $\xi =-1$ and having a modest magnitude at low energies. In contrast, the RPF at high energies is more substantially impacted, where the region with $P \approx 1$ is now largely localized to negative values of pitch for the energy range considered. Finally, increasing the electric field to $E_\Vert = 10$, with $Z_{{\rm eff}}$ = 1 and $\alpha = 0$, results in a drop in the location of the $P=0.5$ contour (compare figure 2a,d), due to the electric field being able to overcome collisional drag for a larger range of energies. The residual of (3.2), multiplied by the prefactor $p^2/( 1+p^2)$ contained in the loss defined in (3.4), is shown in figure 3 for the four cases described above. Here, it is evident that the residuals have been reduced to a small value with a maximum of ${\approx }0.03$ for the largest electric field case (see figure 3d). A direct assessment of the accuracy of the PINN's prediction of the avalanche growth rate will be made in § 4.3, where the predictions of the PINN are compared with first principle Monte Carlo simulations.

Figure 2. Runaway probability functions for different values of the physics parameters $( E_\Vert, Z_{{\rm eff}}, \alpha )$. Panel (a) is for $( E_\Vert = 3, Z_{{\rm eff}} = 1, \alpha = 0)$, panel (b) is for $( E_\Vert = 3, Z_{{\rm eff}} = 10, \alpha = 0)$, panel (c) is for $( E_\Vert = 3, Z_{{\rm eff}} = 1, \alpha = 0.2)$ and panel (d) is for $( E_\Vert = 10, Z_{{\rm eff}} = 1, \alpha = 0)$.

Figure 3. Residual of (3.2) (multiplied by the prefactor $p^2/( 1+p^2)$) for different values of the physics parameters $( E_\Vert, Z_{{\rm eff}}, \alpha )$. Panel (a) is for $( E_\Vert = 3, Z_{{\rm eff}} = 1, \alpha = 0)$, panel (b) is for $( E_\Vert = 3, Z_{{\rm eff}} = 10, \alpha = 0)$, panel (c) is for $( E_\Vert = 3, Z_{{\rm eff}} = 1, \alpha = 0.2)$ and panel (d) is for $( E_\Vert = 10, Z_{{\rm eff}} = 1, \alpha = 0)$.

3.4. Critical energy to run away

The location of the $P = 0.5$ contour $( p_{{\rm crit}}, \xi _{{\rm crit}} )$ provides a useful reference point for identifying the critical energy and pitch above which electrons are likely to run away. The parametric solution of the PDE given by the PINN thus provides $( p_{{\rm crit}},\xi _{{\rm crit}} )$ across the entire domain by simply extracting the $P = 0.5$ surface. To show the dependence of ($E_\Vert, Z_{{\rm eff}}, \alpha$) on the critical energy to run away, we evaluate the energy at which $P = 0.5$ and $\xi = -1$ across the entire training region, which is shown in figure 4. Here, the $\log _{10}$ of the critical energy in units of eV is plotted, where the empty regions represent scenarios where $E_\Vert$ is below the avalanche threshold $E_{{\rm av}}$ (see § 4.1 below). The dependence of $E_\Vert$ and $Z_{{\rm eff}}$ on the critical energy can be seen in figure 4(a) for $\alpha = 0$, where increasing the electric field from $E_\Vert = 2$ to $E_\Vert = 10$ decreases the critical energy from $\sim$1 MeV to $\sim$36 keV. At a given electric field (e.g. $E_\Vert = 4$), the critical energy increases by almost an order of magnitude as $Z_{{\rm eff}}$ is increased from one to ten, indicating that higher $Z$ elements in the plasma increase the critical energy for REs. Moreover, the addition of synchrotron radiation (compare figure 4b,c), shows that the region where $E_\Vert < E_{{\rm av}}$ is increased (larger white region). As $\alpha$ is increased further (compare figure 4a,c), however, we also see that the shape and range of the critical energy varies modestly, indicating the critical energy is modestly impacted by $\alpha$ when $E_\Vert \gg 1$.

Figure 4. Value of $\log _{10} E_{{\rm crit}}$, where $E_{{\rm crit}}$ is the critical energy to run away in eV. Panel (a) is for $\alpha = 0$, panel (b) is for $\alpha = 0.1$ and panel (c) is for $\alpha = 0.2$.

4.1. Secondary source of runaway electrons

Our aim in this section will be to utilize the RPF PINN described in § 3.3 to estimate the avalanche growth rate. Denoting the source of secondary electrons as $S ( p, \xi )$, the rate that REs form due to the avalanche mechanism can be expressed as (Liu et al. Reference Liu, Brennan, Boozer and Bhattacharjee2017)

(4.1)\begin{equation} \left. \frac{{\rm d}n_{RE}}{{\rm d}t} \right|_{{\rm av}} = \int {\rm d}^3 \,p S \left( p, \xi \right) P \left( p, \xi \right) , \end{equation}

where $n_{{\rm RE}}$ is the density of REs. Here, $S ( p, \xi )$ indicates the rate and momentum space distribution of REs, whereas $P ( p, \xi )$ indicates the probability that an electron at a given momentum space location $( p, \xi )$ runs away. By integrating over momentum space, $S ( p, \xi ) P ( p, \xi )$ will thus indicate the expected rate that REs form due to the avalanche mechanism. The primary challenge with evaluating (4.1) is due to $S ( p, \xi )$ depending on the electron distribution $f_e$, i.e.

(4.2)\begin{equation} S \left( p, \xi \right) = \int {\rm d}^3 \,p^\prime S_0 \left(p^\prime , \xi^\prime , p , \xi \right) f_e \left( p^\prime, \xi^\prime \right) , \end{equation}

where $S_0 (p^\prime, \xi ^\prime, p , \xi )$ is defined by

(4.3)\begin{equation} S_0 \left(p^\prime , \xi^\prime , p , \xi \right) = n_e c r^2_e \frac{v^\prime}{2{\rm \pi} p^2} \frac{{\rm d} \sigma_M \left( p^\prime , p \right)}{{\rm d}p} \varPi \left( p^\prime, \xi^\prime , p ; \xi \right) , \end{equation}

where $r_e = e^2/( 4{\rm \pi} \epsilon _0 m_e c^2)$ is the classical electron radius, ${\rm d}\sigma _M/{\rm d}p$ is the Møller cross-section (Møller Reference Møller1932; Ashkin, Page & Woodward Reference Ashkin, Page and Woodward1954), $\varPi ( p^\prime, \xi ^\prime, p ; \xi )$ describes the pitch-angle dependence of secondary electron generation (see Boozer (Reference Boozer2015) for an explicit expression), ${\rm d}^3p = 2{\rm \pi} p^2\,{\rm d}p\,{\rm d}\xi$ and all variables have been non-dimensionalized $p \to p/m_ec$, $v^\prime \to v^\prime /c$, $\sigma _M \to \sigma _M / r^2_e$. Noting that the solution of the adjoint relativistic Fokker–Planck equation does not directly yield the RE distribution $f_e ( p^\prime, \xi ^\prime )$, a closure relation will need to be introduced to evaluate the rate of RE generation via avalanching. The simplest closure, introduced in Rosenbluth & Putvinski (Reference Rosenbluth and Putvinski1997), involves taking the limit where REs are assumed to have asymptotically high energies and a pitch of $\xi =-1$. While idealized, this closure has been shown to provide a good approximation to the full Møller source evaluated using a self-consistently computed RE distribution (McDevitt et al. Reference McDevitt, Guo and Tang2018). In the limit of $p^\prime \to \infty$ and $\xi ^\prime = -1$, (4.2) asymptotes to

(4.4)\begin{equation} S \left( p, \xi \right) = n_e n_{{\rm RE}} c r^2_e \frac{v}{\gamma^2-1} \frac{1}{\left( \gamma-1 \right)^2} \delta \left( \xi - \xi_1 \right) , \end{equation}

where we have introduced the Lorentz factor $\gamma \equiv \sqrt {1 + p^2}$ and $\xi _1$ is defined by

(4.5)\begin{equation} \xi_1 ={-}\sqrt{\frac{\gamma-1}{\gamma+1}} . \end{equation}

Using (4.4), (4.1) reduces to

(4.6)\begin{equation} \left. \frac{{\rm d} n_{RE}}{{\rm d}t} \right|_{av} = 2{\rm \pi} n_e n_{{\rm RE}} c r^2_e \int {\rm d}p\, p^2 \frac{v}{\left( \gamma^2-1\right) \left( \gamma - 1 \right)^2} P \left( p, \xi_1 \right) . \end{equation}

Noting that the right-hand side of (4.6) is directly proportional to the RE density $n_{{\rm RE}}$, this implies an exponentially growing solution with a growth rate given by

(4.7)\begin{equation} \tau_c \gamma_{{\rm av}} = \frac{1}{2\ln \varLambda} \int {\rm d}p \,\frac{v}{ \left( \gamma - 1 \right)^2} P \left( p, \xi_1 \right) , \end{equation}

where $\xi _1$ is defined by (4.5). Thus, once $P ( p, \xi )$ has been evaluated, the avalanche growth rate can be directly inferred from (4.7).

A caveat when evaluating (4.7) is that, while this integral formally extends to $p\to \infty$, the RPF is evaluated assuming a finite $p_{\max }$. The error induced by this approximation can be estimated by considering the magnitude of the integrand in (4.7) across the range of momenta used in this paper (see figure 5, where $p_{\min } \approx 0.2$ and $p_{\max } \approx 10.74$). Here, we have taken $P=1$ inside the integrand of (4.7), such that figure 5 provides an upper bound on the value of the integrand. Noting that the integrand has decayed to a value of ${\approx }3.5\times 10^{-4}$ at the upper boundary, this implies a small contribution to the avalanche growth rate for secondary electrons born with $p > p_{\max }$, particularly for parameters where the system is well above marginality.

Figure 5. The integrand $v / ( \gamma - 1)^2 /( 2 \ln \varLambda )$ of (4.7) when $P=1$. The low and high energy bounds are $p_{\min } \approx 0.2$ and $p_{\max } \approx 10.74$, respectively. The Coulomb logarithm was taken to be $\ln \varLambda = 15$.

4.2. Parametric dependence of RE avalanche

Using the RPF evaluated in § 3, (4.7) can be used to infer the avalanche growth rate across the parameter space $( E_\Vert, Z_{{\rm eff}}, \alpha )$ (see figure 6). Here, the avalanche growth rate increases approximately linearly with the electric field when $E_\Vert \gg 1$, with the slope and threshold sensitive to $Z_{{\rm eff}}$ and $\alpha$. Specifically, the value of $\alpha$ significantly impacts the avalanche threshold (i.e. where $\gamma _{{\rm av}} \approx 0$), but has a negligible impact at large electric fields. In contrast, $Z_{{\rm eff}}$ strongly impacts the avalanche growth rate for all values of the electric field. A feature of (4.7) is that since it indicates the number of REs generated via large-angle collisions, it is thus positive definite. As a result, avalanche growth rates predicted by (4.7) will not account for the decay of the RE population when $E_\Vert < E_{{\rm av}}$ and will instead asymptote to zero. While such behaviour is strictly correct when focusing solely on the avalanche growth mechanism, when developing a RE model appropriate for coupling with a self-consistent MHD solver it will be necessary to account for the decay of the RE distribution when $E_\Vert < E_{{\rm av}}$. The inclusion of this additional physics will be the subject of future work.

Figure 6. Avalanche growth rate vs electric field for different values of $Z_{{\rm eff}}$ and $\alpha$. Panel (a) is for $\alpha = 0$, panel (b) is for $\alpha = 0.1$ and panel (c) is for $\alpha = 0.2$. The blue points in panel (b) represent Monte Carlo simulations with $n_e = 5 \times 10^{21}\ {\rm m}^{-3}$, and (4.8) is shown as the dashed green curve. The Coulomb logarithm was taken to be $\ln \varLambda = 15$.

It will be of interest to compare the predictions of the PINN with available analytic theories. Considering the avalanche growth rate given by (Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997)

(4.8)\begin{equation} \tau_c \gamma_{{\rm RP}} = \frac{1}{\ln \varLambda} \sqrt{\frac{\rm \pi}{3\left( Z_{{\rm eff}}+5\right)}} \left( E_\Vert - 1 \right) , \end{equation}

a comparison between the PINN predictions and (4.8) is shown in figure 6(b). For reference, a small number of avalanche growth rates computed by a Monte Carlo code (McDevitt, Guo & Tang Reference McDevitt, Guo and Tang2019) using the complete Møller source are also shown (the blue points on figure 6b). Details on how this Monte Carlo data set was generated are given in § 4.3 below. It is apparent that both the PINN and (4.8) yield results in reasonable agreement with the Monte Carlo points, where the PINN yields more accurate results near threshold. In particular, while (4.8) implies an avalanche threshold field of $E_\Vert = 1$ regardless of the values of $( Z_{{\rm eff}}, \alpha )$, the PINN is able to account for the variation of the avalanche threshold for non-zero values of $\alpha$ and different $Z_{{\rm eff}}$. In addition, at larger electric fields a modest difference in the slope of the avalanche growth rate predicted by (4.8) compared with the Monte Carlo data is evident. The predictions of the PINN, in contrast, are able to recover the correct slope of the avalanche growth rate, albeit with a modest offset in the magnitude as will be discussed further in § 4.3 below.

When well above marginality the avalanche growth rate is most conveniently described by evaluating the amount of poloidal flux needed to increase the amplitude of a RE seed by one order of magnitude (Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997; Boozer Reference Boozer2018). This quantity can be evaluated by noting that, for $E_\Vert \gg 1$, the avalanche growth rate scales linearly with $E_\Vert$, i.e.

(4.9)\begin{equation} \gamma_{{\rm av}} = \gamma_0 \left( \frac{E_\Vert}{E_c} - 1\right) \approx \gamma_0 \frac{E_\Vert}{E_c} \approx \frac{\gamma_0}{E_c} \frac{1}{R_0} \frac{\partial \psi}{\partial t} , \end{equation}

where $\gamma _0$ is a constant that depends on $( Z_{{\rm eff}}, \alpha, \ln \varLambda )$, $R_0$ is the major radius and $\psi$ is the poloidal flux function. After integrating (4.9) over the time interval $t_f - t_i$, the number of exponentials of the RE population can be written in terms of the change of poloidal flux, and a constant $\psi _{\text {exp}}$ that defines the efficiency of the avalanche. In particular, integrating (4.9) over the time interval $t_f - t_i$, yields

(4.10)\begin{equation} N_{\exp} = \int^{t_f}_{t_i} {\rm d}t \,\gamma_{{\rm av}} \approx \frac{\gamma_0}{E_c} \frac{1}{R_0} \int^{t_f}_{t_i} {\rm d}t \,\frac{\partial \psi}{\partial t} = \frac{\Delta \psi}{ \psi_{\text{exp}}} , \end{equation}

where $\Delta \psi \equiv \psi ( t_f) - \psi ( t_i)$ and we have defined $\psi _{\exp } \equiv R_0 E_c/ \gamma _0$, which is related to the amount of poloidal flux required to effect one exponential amplification of the RE population. Thus, once $\gamma _0$ is inferred, the efficiency through which the decay of the poloidal flux leads to an amplification of the seed RE population can be evaluated. Further defining the quantity $\psi _{10}\equiv \ln 10 \psi _{\exp }$ (which delineates base ten amplifications), the efficiency of the avalanche growth rate for a broad range of parameters is shown in figure 7(a). Here, the avalanche growth rate is most efficient for low values of $Z_{{\rm eff}}$ and increases by nearly a factor of two for $Z_{{\rm eff}}=10$, agreeing with previous results (McDevitt et al. Reference McDevitt, Guo and Tang2019). We have not indicated the dependence of $\psi _{10}$ on $\alpha$ since, for large electric fields, the avalanche growth rate will be independent of $\alpha$.

Figure 7. (a) The value of $2{\rm \pi} \psi _{10}/\mu _0 R_0$ as a function of $Z_{{\rm eff}}$, with a Coulomb logarithm of $\ln \varLambda = 15$. (b) Avalanche threshold $E_{{\rm av}}$ as a function of the synchrotron radiation strength $\alpha$. The solid lines represent the PINN predictions, the dashed lines and ‘o’ markers represent (B 15) and Monte Carlo results from McDevitt et al. (Reference McDevitt, Guo and Tang2019), respectively, and WebPlotDigitizer was used to extract the Monte Carlo values. The blue values represent $Z_{{\rm eff}} = 1$ and the red values represent $Z_{{\rm eff}} = 5$. The Coulomb logarithm was taken to be $\ln \varLambda = 20$.

In addition, the avalanche threshold is strongly impacted by the parameters $( Z_{{\rm eff}}, \alpha )$. In particular, as $Z_{{\rm eff}}$ or $\alpha$ are increased, the threshold electric field $E_{{\rm av}}$, where the avalanche growth rate is zero, is increased. The dependence of the avalanche threshold on $( Z_{{\rm eff}}, \alpha )$ is shown in figure 7(b), where the PINN predictions are the solid lines. The Monte Carlo results are the ‘o’ markers, and the dashed lines represent empirical fit formula given by (B 15) of McDevitt et al. (Reference McDevitt, Guo and Tang2019). Here, since the avalanche growth rates predicted by the PINN are positive definite, we define the avalanche threshold as the value of the electric field where $\tau _c \gamma _{{\rm av}} = 2\times 10^{-3}$. From figure 7(b) it is apparent that the PINN predicts $E_{{\rm av}}$ particularly well across $\alpha$ for $Z_{{\rm eff}} = 1$, with somewhat larger deviations evident for ${Z_{{\rm eff}}=5}$. The systematic under-prediction of the avalanche threshold is due to the assumption in (4.7) that primary electrons have infinite energy. Such an assumption overestimates the number of secondary electrons generated and, perhaps more importantly, neglects that the primary electron population itself will slowly decay to the thermal bulk when $E \approx E_{{\rm av}}$.

4.3. Verification of the surrogate model

In this section we will verify the PINN's predictions of the avalanche growth rate against first principle Monte Carlo simulations across the training region ($E_\Vert, Z_{{\rm eff}}, \alpha$). The RunAway Monte Carlo (RAMc) code is employed (see McDevitt et al. (Reference McDevitt, Guo and Tang2019) for a detailed description), which evolves the guiding centre motion of relativistic electrons and includes effects from small-angle collisions, large-angle collisions and synchrotron radiation. Large-angle collisions are evaluated using the full Møller source (Møller Reference Møller1932). In order to avoid toroidal corrections to the avalanche growth rate (McDevitt & Tang Reference McDevitt and Tang2019; Arnaud & McDevitt Reference Arnaud and McDevitt2024), all REs are initialized near $r=0$, and a large tokamak device was chosen with a minor radius of $a=2\ \text {m}$ and major radius $R_0=6\ \text {m}$, in order to render spatial transport negligible. A geometry with circular flux surfaces was selected, with a safety factor profile taken to be $q ( r ) = 2.1 + 2( r/a)^2$.

The Monte Carlo avalanche simulations are set up by initializing a small population of electrons (8 for this analysis) at high momentum ($p \in (10,20)$) and strongly aligned with the magnetic field ($\xi \in (-0.9,-1.0)$). They are then allowed to exponentially grow in time until a saturated growth rate can be identified, generally resulting in the initial seed RE population growing by several orders of magnitude. The plasma parameters chosen were a density of $n_e = 10^{21}\ {\rm m}^{-3}$, a temperature of $T_e = 10$ eV and a toroidal magnetic field strength that was varied to give the appropriate value of $\alpha$. One caveat is that at a fixed density and temperature $\alpha \propto B^2$, which in turn can lead to orbit drifts at the low boundary of synchrotron radiation due to the weak magnetic field. These orbit drifts are thus ensured to be negligible by enforcing a minimum toroidal magnetic field of $B = 2$ T, which corresponds to a lower bound ($\alpha \approx 2.8 \times 10^{-3}$) in the training region for synchrotron radiation. Fifty randomly selected samples of ($E_\Vert, Z_{{\rm eff}}, \alpha$) are chosen as input parameters for the Monte Carlo solver and the avalanche PINN, where cases below the avalanche threshold are discarded, leaving forty-one samples. A comparison between the predictions of the avalanche PINN and Monte Carlo solver is shown in figure 8(a). Here, the grey dashed line represents equality ($y$-axis equal to $x$-axis) between the PINN and Monte Carlo predictions. The coefficient of determination (Guilford & Guilford Reference Guilford and Guilford1936) between the PINN predictions and the Monte Carlo avalanche growth rate was evaluated to be $R^2 \approx 0.9849$, indicating that the PINN was able to accurately predict the avalanche growth rate across a broad range of parameters. One noticeable feature present in figure 8(a) is the systematic over-prediction of the avalanche growth rate by the PINN. This feature is expected, and is due to the use of the simplified Rosenbluth–Putvinski secondary source term described by (4.4), which assumes all primary electrons to have $p \to \infty$ and $\xi =-1$. We note that this feature can be mitigated by introducing a simple discrepancy model, where for the parameter regime considered in this study, multiplying the predictions of the PINN by a factor of 0.94886 significantly reduced the discrepancy between the PINN and Monte Carlo predictions. The correlation between the avalanche PINN predictions including the discrepancy model and the Monte Carlo predictions are shown in figure 8(b), where the coefficient of determination increased to $R^2$ = 0.9988. Improving the RE closure used when approximating the Møller source will be the subject of future work.

Figure 8. (a) Avalanche growth rate comparison between the PINN and Monte Carlo solver. The grey dashed line represents the coefficient of determination of $R^2 = 1$. (b) Same as panel (a), but with the predictions of the avalanche PINN multiplied by the factor $0.94886$. (c) Avalanche growth rate comparison between (4.8) and the Monte Carlo solver. The avalanche growth rates were evaluated across $E_\Vert \in (1,10)$, $Z_{{\rm eff}} \in (1,10)$ and $\alpha \in (2.8 \times 10^{-3},0.2)$. The other parameters were chosen to be $T_e = 10\ {\rm eV}$ and $n_e = 10^{21}\ {\rm m}^{-3}$.

The PINN is also shown to robustly provide better avalanche growth rate predictions than the analytical expression provided by (4.8) (compare figure 8a,c). Here, the avalanche growth rate provided by (4.8) compared with the Monte Carlo solver has a weaker correlation than that between the PINN and Monte Carlo solver, where the coefficient of determination for figure 8(c) was $R^2 \approx 0.9644$. The avalanche PINN thus improves on the accuracy of the avalanche growth rate predictions, even in the absence of the 0.94886 factor.