2.1. Simulation set-up

The disruption simulations assume an initially ($t<0$) pure fully ionized deuterium–tritium (D–T) plasma with $50$–$50\,\%$ isotope concentrations. Specifically, the initial electron density is spatially constant $10^{20}\,{\rm m}^{-3}$, the temperature is parabolic with $20\,{\rm keV}$ on axis and the total plasma current is $15\,{\rm MA}$. The simulations use an ITER-like magnetic geometry with major radius $R_0=6\,{\rm m}$, minor radius $a=2\,{\rm m}$, wall radius $b=2.833\, {\rm m}$, on-axis toroidal magnetic field $B(r=0)=5.3\,{\rm T}$ and a resistive wall time of $\tau _w=0.5\,{\rm s}$, as well as a Miller model equilibrium (Miller et al. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998) with realistic, radially varying shaping parameters; further information is given in Appendix A.

The simulations are performed by the Dream (Disruption Runaway Electron Analysis Model) code that captures the particle acceleration and energy dissipation processes following a disruption (Hoppe, Embreus & Fülöp Reference Hoppe, Embreus and Fülöp2021). It solves a set of coupled transport equations describing the evolution of temperature, ion charge state densities, current density and electric field in arbitrary axisymmetric geometry. The temperature evolution includes ohmic heating, radiated power using atomic rate coefficients, collisional energy transfer from hot electrons and ions, as well as dilution cooling.

Dream allows modelling of the REs at different degrees of approximation ranging from fluid to fully kinetic. As we do not require kinetic outputs, we limit our modelling to the least computationally expensive, fluid treatment of the plasma. This means that the thermal bulk of cold electrons and the small runaway population are modelled as two separate fluid species. The former is characterized by a density $n_e$, a temperature $T_e$ as well as an ohmic current density $j_{\rm ohm}$, and the REs are described by their density $n_{\rm RE}$. It is assumed that the REs move with the speed of light parallel to the magnetic field, hence their associated current density is $j_{\rm RE} = ecn_{\rm RE}$. The simulations include Dreicer, hot-tail and avalanche sources, as well as REs generated by Compton scattering of $\gamma$ photons and tritium decay. These are modelled as quasi-stationary sources feeding electrons into the runaway population (Fülöp et al. Reference Fülöp, Helander, Vallhagen, Embréus, Hesslow, Svensson, Creely, Howard and Rodriguez-Fernandez2020). The runaway generation rates used in the simulations have been benchmarked with the corresponding kinetic results (Hoppe et al. Reference Hoppe, Embreus and Fülöp2021). Further details on the simulations are given in Appendix A.

Neutral neon and deuterium are introduced with zero temperature at the start of the simulation ($t=0$). At the same time an elevated transport of electron heat and energetic electrons is activated, using a Rechester–Rosenbluth-type model (Rechester & Rosenbluth Reference Rechester and Rosenbluth1978) with a radially constant normalized magnetic perturbation amplitude $\delta B/B$. This is done to emulate the break-up of flux surfaces during the TQ, and leads to heat losses, with a heat diffusivity proportional to $R_0 v_{te} (\delta B/B)^2$, where $v_{te}=\sqrt {2 T_e/m_e}$ is the local electron thermal speed. The full expression is given by (B.5) of Hoppe et al. (Reference Hoppe, Embreus and Fülöp2021). In four optimization runs $\delta B/B$ is scanned over the range 0.2 %–0.5 % that falls within the range of values observed in magnetohydrodynamic simulations of the TQ (Hu et al. Reference Hu, Nardon, Hoelzl, Wieschollek, Lehnen, Huijsmans, van Vugt and Kim2021).

During the TQ we also account for a diffusive transport of REs using a diffusion coefficient of similar form, but assuming a parallel streaming along the perturbed field lines at the speed of light, $D_{\rm RE}={\rm \pi} R_0 c (\delta B/B)^2$. This approach neglects the momentum space variation of the transport coefficients (Särkimäki et al. Reference Särkimäki, Embreus, Nardon and Fülöp2020), as well as the form of the RE distribution function, which would reduce the effect of runaway transport. Thus, using this expression provides an upper bound on the effect of runaway transport for a given magnetic perturbation amplitude (Svensson et al. Reference Svensson, Embreus, Newton, Särkimäki, Vallhagen and Fülöp2021). We employ here the same $\delta B/B$ as for electron heat transport for consistency.

The injected material is ionized by its interaction with the plasma, and cools it by radiation and dilution. When the average electron temperature falls below $10^{-3}$ times the maximum initial temperature (here $20\,{\rm eV}$), we assume that the TQ is completed and the flux surfaces reform. After the TQ the transport of energetic electrons is switched off, and a significantly reduced, but finite electron heat diffusivity is used ($\delta B/B=0.04\,\%$). This is to avoid the development of non-physical narrow hot ohmic channels during the CQ. Such ohmic channels are soliton-like solutions of the problem (Putvinski et al. Reference Putvinski, Fujisawa, Post, Putvinskaya, Rosenbluth and Wesley1997) without sufficient heat diffusivity. In a physical system the corresponding excessive temperature and current gradients would be expected to destabilize these formations well before they could fully form. Note, that the diffusive heat transport is subdominant compared with radiative heat losses at the low post-TQ temperatures, thus this heat transport has no effect besides not allowing hot channels to form.

2.2. Optimization

The optimization problem involves multiple objectives, i.e. multiple quantities need to be within certain limits simultaneously. The maximum value of the total RE current and the fraction of transported heat losses must be small, while the CQ time should be within certain limits. These quantities are normalized and combined into a single scalar cost function $\mathcal {L} \ge 0$ which is to be minimized. Denoting the control vector containing the parameters by $\boldsymbol {x}$, we wish to find the $\boldsymbol {x}^{\ast }$ that minimizes $\mathcal {L}$, where $\boldsymbol {x}$ resides in a specified volume $\mathcal {V}\subset \mathbb {R}^d$ of the control space (where $d$ is the dimensionality of the optimization).

We employ Bayesian optimization (Brochu, Cora & de Freitas Reference Brochu, Cora and de Freitas2010) using Gaussian process regression (Rasmussen & Williams Reference Rasmussen and Williams2005), using the Bayesian Optimization (Nogueira Reference Nogueira2014) Python package. A Gaussian process is fitted to the already sampled points $\{\boldsymbol {x}_i\}_{i=1}^n$, and the expected improvement acquisition strategy (described in Appendix B) is used to choose the next point to be sampled, $\boldsymbol {x}_{n+1}$. The Gaussian process contains information on both the expected value $\mu (\boldsymbol {x})$ and the uncertainty of the estimate of $\mathcal {L}$, quantified in terms of the covariance $k(\boldsymbol {x}, \boldsymbol {x}')$ between any two points $\boldsymbol {x}$ and $\boldsymbol {x}'$. In this process there is a balance between exploration and exploitation, i.e. search within regions with high uncertainty, as well as in regions that are most likely to host the global optimum.

The cost function we use is of the form

(2.1)\begin{equation} \mathcal{L}=\frac{I_{\rm RE}^{\rm max}}{I_{\rm RE}^{\rm tol}} + \frac{I_{\rm ohm}^{\rm fin}}{I_{\rm ohm}^{\rm tol}} + 10 \frac{\eta_{\rm cond}}{\eta_{\rm cond}^{\rm tol}} + 100\,\theta(t_{\rm CQ}) , \end{equation}

where $I_{\rm RE}^{\rm max}$ is the maximum RE current in the simulation, $I_{\rm RE}^{\rm tol}=150\,{\rm kA}$ represents the tolerable RE current in ITER and $I_{\rm ohm}^{\rm fin}$ is the ohmic current at the end of the ($150\,{\rm ms}$ long) simulation. A significant remnant ohmic current may be the sign of an incomplete TQ, and it can still potentially be converted to a RE current. Thus it is treated on equal footing with the RE current, so we also set $I_{\rm ohm}^{\rm tol}=150\,{\rm kA}$. Here, $\eta _{\rm cond}^{\rm tol}=0.1$ is the tolerable transported heat loss fraction. The prefactor $10$ in the $\eta _{\rm cond}$ term is used to get a penalty for non-tolerable transported heat losses comparable to typical penalties obtained for mega-ampere (${\rm MA}$) size currents. Finally, to penalize CQ times $t_{\rm CQ}$ below $t_{\rm L}=50\,{\rm ms}$ and above $t_{\rm U}=150\,{\rm ms}$ we use the penalty function

(2.2)\begin{equation} \theta(t_{\rm CQ})= \tilde{\varTheta}(t_{\rm L} - t_{\rm CQ}) + \tilde{\varTheta}(t_{\rm CQ} - t_{\rm U}), \end{equation}

where $\tilde {\varTheta }(t)= \tfrac {1}{2}[1+\tanh (t/\Delta t)]$ is a function similar to a step function but smooth with a transition width set by $\Delta t=3.3\,{\rm ms}$. Values of $t_{\rm CQ}$ outside the tolerable range yield a penalty as high as the maximum achievable penalty for any of the other terms in (2.1), due to the prefactor $100$ in front of $\theta$. We calculate the CQ time as $t_{\rm CQ}=[t(I_{\rm ohm}=0.2 I_p^{0})-t(I_{\rm ohm}=0.8 I_p^{0})]/0.6$ (Hender et al. Reference Hender, Wesley, Bialek, Bondeson, Boozer, Buttery, Garofalo, Goodman, Granetz and Gribov2007), where $I_{\rm ohm}(t)$ is the total ohmic current and $I_p^{0}$ is the initial plasma current.

In addition we set $\mathcal {L}=500$ for simulations where the TQ is not complete within $20\,{\rm ms}$, our condition for which is that the average temperature falls to below $10^{-3} T_e(r=0,t=0)$. Finally, as it is difficult to completely avoid simulations that fail due to numerical issues, we use $\mathcal {L}=500$ for these as well.

Independently of their dimensionality, the optimizations use $400$ samples, chosen according to the acquisition function through sequential function evaluations, following $20$ randomly selected initial samples. As the parameter space of injected quantities ranges across multiple orders of magnitude, the logarithm of the injected quantities is used as optimization parameters.

3.1. Optimization landscape with constant concentrations

First, we perform optimization in the two-dimensional (2-D) parameter space of radially constant injected deuterium and neon densities, $n_{{\rm D,inj}}$ and $n_{{\rm Ne,inj}}$. The ranges of injected densities we consider are $n_{{\rm D,inj}} \in [10^{18},\,3.16\times 10^{22}]\,{\rm m}^{-3}$, and $n_{{\rm Ne,inj}} \in [10^{16},\,10^{20}]\,{\rm m}^{-3}$.

Figure 1 shows the estimated mean of the cost function $\mu$ on a logarithmic contour plot for four different values of $\delta B/B$, with blue shades representing favourable and red shades unfavourable values. Each panel used $420$ samples, indicated by grey dots, while the optima are indicated with black stars. The area of favourable values (with blue shades) decreases with increasing $\delta B/B$, and this is mostly due to the increasing transported heat fraction, and to a lesser degree to an increasing RE current, to be discussed further in relation to figure 5. In general, the lower left corner of the plots is occupied by cases with an incomplete TQ. In this case the plasma tends to get reheated after the prescribed transport event, leading to long CQ times (i.e. $t_{\rm CQ}>150\,{\rm ms}$). With an increasing $\delta B/B$ the incomplete TQ region shrinks somewhat.

Figure 1. The estimated mean of the cost function of Bayesian optimizations in the $n_{{\rm D,inj}}$–$n_{{\rm Ne},{\rm inj}}$ space for various normalized magnetic perturbation amplitudes. The colour code varies from blue to red tones, representing favourable and unfavourable values of $\mu$; (a) $\delta B/B=0.2\,\%$, (b) $0.3\,\%$, (c) $0.4\,\%$, (d) $0.5\,\%$. Black stars indicate the locations of the optima. Grey dots show the samples taken; note that these are more numerous in the vicinity of the optima. Circles with case identifiers in panel (b) indicate the cases discussed in § 3.2.

Another general feature is a relatively narrow corridor of favourable parameters in the vicinity of $n_{{\rm D,inj}} = 10^{22}\,{\rm m}^{-3}$, extending from the lowest $n_{{\rm Ne,inj}}$ values plotted to a bit above $n_{{\rm Ne,inj}} = 10^{19}\,{\rm m}^{-3}$. The optima also reside in these corridors at $n_{{\rm Ne,inj}}$ values of a few times $10^{18}\,{\rm m}^{-3}$. Additionally, a wider corridor of moderate values of $\mathcal {L}$ extends to the left of the optima, in the range $n_{{\rm Ne,inj}}\approx 3\times 10^{18}$–$2\times 10^{19}\,{\rm m}^{-3}$, which is most pronounced in the $\delta B/B=0.3\,\%$ and $0.4\,\%$ cases. At $n_{{\rm Ne,inj}}$ values above that, $\mathcal {L}$ increases, and decreases again around the highest $n_{{\rm Ne,inj}}$ values included in this optimization. Before analysing the characteristic behaviours in these different regions of the current optimization landscape in § 3.2, we shall discuss the detailed dynamics at the optima.

We consider the behaviour of the representative optimum obtained in the $\delta B/B=0.3\,\%$ case (indicated by the black star in figure 1b), located at $n_{\rm D,inj}=9.4\times 10^{21}\,{\rm m}^{-3}$ and $n_{{\rm Ne,inj}}=2.9\times 10^{18}\,{\rm m}^{-3}$, shown in figure 2. Following the instantaneous material injection at $t=0$, the temperature profile drops by a factor of ${\approx }100$ within a $\mathrm {\mu }{\rm s}$, due to dilution. This is followed by an approximately exponential cooling with a characteristic time of $\tau \approx 1.5 \,{\rm ms}$. After the initial exponential cooling, a cold front starts to propagate radially inward from the edge. This inward propagating cooling is seen in the $t=2\,{\rm ms}$ curve in figure 2(b). This cooling proceeds until almost the entire plasma settles at around $5\,{\rm eV}$ (see the $t=10\,{\rm ms}$ curve), representing the equilibrium between ohmic heating and radiation corresponding to the ion composition and current density of the plasma. Then there is another inward propagating cooling happening over the next $50\,{\rm ms}$. As the ohmic current density drops during the CQ, the equilibrium temperature falls from ${\approx }5$ to $\approx 1.2\,{\rm eV}$. The temperature is radially uniform at this level at $60\,{\rm ms}$ (black curve), and remains there until the end of the simulation.

Figure 2. The best performing case for the optimization in the $n_{{\rm D,inj}}$–$n_{{\rm Ne},{\rm inj}}$ space, for $\delta B/B=0.3\,\%$. (a) The time evolution of the total plasma current (dashed), and its ohmic (solid) and RE (dash-dotted) components. (b–d) Radial profiles of quantities in a few time points, indicated by their respective figure legends; with increasing time corresponding to darker colours. (b) Electron temperature. (c) Ohmic (solid) and RE (dashed) current density. (d) Parallel electric field (solid). The effective critical electric field is also indicated for $t=60\,{\rm ms}$ (dotted); note that it does not vary appreciably over time.

The ohmic current density gradually decreases in the core and it drops rapidly across the cold front at the edge; compare the $10\,{\rm ms}$ curves in figure 2(b,c). This front propagates inward in the first $40\,{\rm ms}$, after which the ohmic current gets rapidly replaced by RE current; see the process in terms of total current components in figure 2(a), and RE current density at $60\,{\rm ms}$ in figure 2(c) (dashed line). The electric field exceeds the effective critical field $E_{c}^{{\rm eff}}$ – calculated as in appendix C2 of Hoppe et al. (Reference Hoppe, Embreus and Fülöp2021) – first in the edge, then it grows to an approximately radially constant value around $30 \,{\rm V}\,{\rm m}^{-1}$, approximately 4 times $E_{c}^{{\rm eff}}$, where it stays until the macroscopic RE conversion starts. Then it drops into the vicinity of $E_{c}^{{\rm eff}}$, such that in the core $E_\|$ is pinned to $E_{c}^{{\rm eff}}$, and it takes radially decreasing values at the edge; compare $E_\|$ (solid black curve) with $E_{c}^{{\rm eff}}$ (dotted) in figure 2(d). Then the electric field remains like that until most of the RE current dissipates away. Physically, the dissipation of the RE current, in the absence of transport losses, is caused by a collisional slowing down and thermalization of the REs when $E_\|< E_{c}^{{\rm eff}}$. In the employed fluid RE model it is technically accounted for by allowing the avalanche growth rate to become negative for $E_\|< E_{c}^{{\rm eff}}$ values. The corresponding decay of the RE current is quite pronounced in this case.

3.2. Characteristic cases with constant concentrations

In order to understand the typical dynamics in various regions in the $n_{{\rm D,inj}}$–$n_{{\rm Ne},{\rm inj}}$ space, we consider six representative cases in the $\delta B/B=0.3\,\%$ optimization, with case C1 being the optimum discussed above. The cases are indicated in figure 1(b) and corresponding injected quantities and figures of merit are listed in table 1. Cases C2 and C3 are taken in the high $n_{{\rm D,inj}}$ region of the space; C2 is located in the favourable channel at low $n_{{\rm Ne},{\rm inj}}$, and C3 at even higher $n_{{\rm D,inj}}$ than the optimum. Cases C4 to C6 are taken at a fixed $n_{{\rm D,inj}}=10^{20}\,{\rm m}^{-3}$, at respectively increasing value of $n_{{\rm Ne},{\rm inj}}$. We discuss C1–C3 and C4–C6 in the following subsections.

Table 1. Characteristic cases from the $n_{{\rm D,inj}}$–$n_{{\rm Ne},{\rm inj}}$ optimization landscape for $\delta B/B=0.3\,\%$, their four figures of merit and corresponding cost function values. The cases are marked in figure 1(b); C1 is the optimum.

3.2.1. Representative cases at high $n_{{\rm D,inj}}$

In the RE plateau the electric field tends to stay close to the effective field $E_{c}^{{\rm eff}}$, as expected (Breizman Reference Breizman2014). In particular, following the RE conversion, all $E_\|/E_{c}^{{\rm eff}}$ values (taken at mid-radius) settle around unity, as shown in figure 3(c). The high $n_{{\rm D,inj}}$ cases are all characterized by a significant decay rate of the RE current after it reaches its maximum value; see figure 3(a). This is consistent with their $E_\|/E_{c}^{{\rm eff}}$ being lower than unity towards the edge, as we have seen for C1 in figure 2(d). This is due to the relatively high value of $E_{c}^{{\rm eff}}$ typical at these high $n_{{\rm D,inj}}$ values.

Figure 3. Time evolution of quantities of interest for the high $n_{D,{\rm inj}}$ representative cases: C1–C3. Line colour darkens and dashes shorten with increasing case number, and case numbers are indicated with callouts. (a) Runaway electron current. (b) Electron temperature at mid-radius. (c) Electric field normalized to critical electric field at mid-radius. (Note the longer time range plotted in panel a.)

The dynamics of the RE current in C1 and C2 is fairly similar, as seen in figure 3(a). It may be surprising that $E_\|/E_{c}^{{\rm eff}}$ is almost all the time higher in C1 than C2 – shown in figure 3(c) – but the maximum RE current in C2 is still higher, and is reached a bit earlier. The reason for this is that the temperature drops to $1.09\,{\rm eV}$ in C2 already at $t\approx 10\,{\rm ms}$ (figure 3b), a temperature where $44\,\%$ of the hydrogenic species is recombined, thereby increasing the total-to-free electron density ratio and the avalanche growth rate in proportion. Meanwhile in C1 the temperature does not drop to this low temperature until the RE conversion is over.

The effect of the hydrogen recombination is even more pronounced in C3, where the temperature drops to $1.02\,{\rm eV}$ within a millisecond. The reason for this fast cooling is the very high dilution that brings the temperature down to a range where radiative losses are strong and can effectively (and rapidly) cool the plasma further. That the temperature drops immediately to its final value, without stopping at some higher, intermediate value, can be explained by how the temperature dependence of the total radiative losses ($P$) is affected by the very high hydrogen content. Depending on the hydrogen (including D and T) and neon densities, the curve $P(T_e)$ can exhibit a local minimum in the few ${\rm eV}$ range between a low $T_e$ peak caused by hydrogen and a higher $T_e$ peak from neon. The large hydrogen density in C3 leads to an elevated value of $P$ at this minimum, thereby effectively eliminating the bottleneck this minimum represents concerning the cooling. While $1.02\,{\rm eV}$ is just slightly cooler than the final temperature in C2, now $70\,\%$ of the hydrogenic species are recombined, which, in combination with the early high value of $E_\|/E_{c}^{{\rm eff}}$, leads to an extremely fast RE conversion and the highest RE current among these three cases.

In terms of figures of merit, C3 is not only problematic due to a high $I_{\rm RE}^{\rm max}$ value, but also because of the extremely short $t_{\rm CQ}\approx 5\,{\rm ms}$. While $I_{\rm RE}^{\rm max}$ is not too much higher in C2 than in C1, it has a $\eta _{\rm cond}\approx 44\,\%$, exceeding the tolerable $10\,\%$, unlike C1 and C3. This is due to the small neon content in C2.

The remarkably short cooling times, of the order of $2\,{\rm ms}$, observed at large deuterium injections, such as C3, may be partly due to our simplifying assumption of instantaneous deposition. However, in realistic material injection scenarios, the cooling at a given flux surface can be as rapid as observed here, even if the time scale needed for pellet shards flying at $500\,{\rm m}\,{\rm s}^{-1}$ to travel between the edge and the centre of an ITER plasma is longer (${\approx }4\,{\rm ms}$). As the local cooling time is the crucial factor to get a large hot-tail seed, and furthermore, the rapid avalanche rate depends on the final temperature, similar behaviour is also observed in shattered pellet injection simulations (Vallhagen et al. Reference Vallhagen, Pusztai, Hoppe, Newton and Fülöp2022). Ion convection time scales across the radius in a TQ can also be in the ${\rm ms}$ range. The excessive runaway generation is thus not an artefact of the instantaneous deposition, however, the detailed temperature evolution is expected to be different once the injection dynamics is resolved.

3.2.2. Representative cases at low $n_{{\rm D,inj}}$

The cases at $n_{D,{\rm inj}}=10^{20}\,{\rm m}^{-3}$ – C4 to C6 – are not affected by hydrogen recombination as their temperature never drops below $2\,{\rm eV}$. They reach much higher values of $E_\|/E_{c}^{{\rm eff}}$ than the high $n_{{\rm D,inj}}$ cases, as they have low $E_{c}^{{\rm eff}}$; compare figures 3(c) and 4(c). Their RE conversion timing and magnitude well correlates with when the peak of $E_\| /E_{c}^{{\rm eff}}$ is reached, and its magnitude. This, in turn, depends on the first equilibrium temperature reached, varying between approximately $5$ and $11\,{\rm eV}$, see figure 4(b). This temperature decreases monotonically with increasing injected neon quantity, while the magnitude of the final RE current increases, and the time of RE conversion shifts earlier. Once the conversion is complete, the temperature falls further into the $2$–$4\,{\rm eV}$ range. Note that at these low $n_{{\rm D,inj}}$ cases the dissipation rate of the RE current in the RE plateau is negligible during the simulation, due to the lower $E_{c}^{{\rm eff}}$ values.

Figure 4. Time evolution of quantities of interest for the low $n_{D,{\rm inj}}$ representative cases: C4–C6. Line colour darkens and dashes shorten with increasing case number, and case numbers are indicated with callouts. (a) Runaway electron current. (b) Electron temperature at mid-radius. (c) Electric field normalized to critical electric field at mid-radius. (Note the longer time range plotted in panel a.)

Only in C4 does $t_{\rm CQ}$ fall into the acceptable range, in the other two cases it is too short due to the early RE conversion. The reason for the non-monotonic dependence of $\mathcal {L}$ with increasing $n_{{\rm Ne},{\rm inj}}$, i.e. it is higher for C5 than for C6, is that it is caused by the reduction in the transported heat loss fraction from the $70\,\%$–$80\,\%$ range to $23\,\%$ (that is still not acceptably low though).

3.3. Radially varying material injection

Next, we relax the assumption of spatially homogeneous injection, and allow profile variations with a simple model for the injected densities, where the inward or outward peaking of the profile is set by a single parameter $c_i$ per species $i$

(3.1)\begin{equation} \tilde{n}_{i,{\rm inj}} \propto 1+\tanh \left[c_i\left(\frac{r}{a}-\frac{1}{2}\right)\right], \end{equation}

where the tilde indicates that $\tilde {n}_{i,{\rm inj}}$ is a radially varying quantity. The notation $n_{i,{\rm inj}}$ is reserved to the scalar parameter that appears in the optimization. The factor multiplying the expression in (3.1) is determined such that the total number of injected particles in the plasma is the same as in an injection of a constant density $n_{i,{\rm inj}}$. Negative/positive values of $c_i$ correspond to densities peaked in the plasma centre/edge, and in the optimization we allow values in the $[-10,\, 10]$ range.

Figure 5 shows the $I_{\rm RE}^{\rm max}$ and $\eta _{\rm cond}$ figures of merit, along with the cost function $\mathcal {L}$ at the optima found for different $\delta B/B$ values, when radially constant injection is employed (dotted line, referred to as 2D) and when profile variation is allowed (dashed, 4D). In the latter case, the additional degrees of freedom allow us to find optima with better properties. Since in all cases the remaining ohmic current is much smaller than $I_{\rm RE}^{\rm max}$ (in the $300$-$400\,{\rm kA}$ range), and $t_{\rm CQ}$ is also in the tolerable range, $\mathcal {L}$ is dominated by the two figures of merit plotted. In none of the cases considered is $I_{\rm RE}^{\rm max}$ tolerably small; it is around $4\,{\rm MA}$ independently of $\delta B/B$ in the 2-D optimization, and it reduces almost by a factor of 2 in four dimensions (without any clear trend with $\delta B/B$), as seen in figure 5(a).

Figure 5. Variation of (a) the maximum RE current, (b) the transported heat loss fraction and (c) the corresponding cost function in optimizations, for a range of $\delta B/B$ values, when optimizing only for injected densities (circle markers, blue short dashed curve) and when including profile variation as well in the optimization (squares, red long dashed). In panels (a,b), below the thin solid line the values are considered tolerable. In panel (a) simulations with the parameters corresponding to the 2-D optimum at $\delta B/B=0.3\,\%$ but without activated sources is indicated with a black rectangle marker, and a simulation with a reduced wall radius of $2.15\,{\rm m}$ is shown with a black asterisk.

There are two main reasons for obtaining such high values even in the optimal cases. We consider D–T plasmas, and we include RE seed sources relevant for activated operation, tritium decay and Compton scattering of $\gamma$ photons, in addition to Dreicer and hot-tail RE generation. The tritium decay and Compton sources can provide a significant RE seed even after the TQ, during which the transport due to magnetic perturbations decimates the initial hot-tail and Dreicer seed population. This circumstance also explains the weak sensitivity of $I_{\rm RE}^{\rm max}$ to $\delta B/B$ in the 2-D simulations. A simulation identical to the 2-D optimum at $\delta B/B=0.3\,\%$, but without activated seed sources (i.e. only Dreicer, hot-tail and avalanche sources active) yields a negligibly small $I_{\rm RE}^{\rm max}=4.1\,{\rm kA}$ instead of $4.2 \,{\rm MA}$.

A similarly important factor is the realistic radius of the conducting wall, which is chosen to match the energy in the poloidal magnetic field due to the plasma current within the conducting wall to that observed in JOREK simulations. If in the 2-D optimum at $\delta B/B=0.3\,\%$ we reduce the wall radius from $2.833$ to $2.15\,{\rm m}$, which was used in previous work, e.g. by Vallhagen et al. (Reference Vallhagen, Embreus, Pusztai, Hesslow and Fülöp2020), the RE current reduces to the – non-negligible, but still significantly lower – value of $1\,{\rm MA}$.

The fraction of transported heat losses, shown in figure 5(b), increases strongly with $\delta B/B$ in the 2-D cases, which is not surprising, since the heat transport during the TQ is then increasing, while the radiated losses are not directly impacted by $\delta B/B$. However, when profile variation is allowed $\eta _{\rm cond}$ is almost independent of $\delta B/B$; the reason for this will be explained in relation to figure 7.

3.4. Sensitivity of the optima

To gauge the sensitivity of the optima to the input parameters, we investigate the regions occupied by samples within some range of $\mathcal {L}$ above the optimal values. The location of the optima in the optimization space is marked in figure 6 ($\otimes$ markers). In the 2-D optimization study we also scatter plot all samples in the $10\,\%$ vicinity of the optimum, figure 6(a); this is such a narrow range in $\mathcal {L}$, that any point in this point cloud can be considered equally well performing as the optimum itself. In the 4-D optimization study we show points in the $25\,\%$ vicinity of the optima, figure 6(b–d). As the total number of samples is the same in both the 2-D and 4-D optimization studies, the higher dimensionality in four dimensions implies a sparser exploration in the vicinity of the optimum compared with two dimensions; hence the lower number of points in spite of the wider relative range included.

Figure 6. Scatter plot of input parameters for samples with the lowest $\mathcal {L}$ values in each optimization case. When (a) optimizing only for injected densities (two dimensions) they represent an additional $10\,\%$ range above the optimum, and when (b–d) including profile variation as well in the optimization (four dimensions), they represent a $25\,\%$ range. Darkening colour indicates increasing value of $\delta B/B$, as given in panel (b), and the optima are indicated by $\otimes$ markers. (a,b) Concentration space, (c,d) correlating concentration with profile parameter of an injected species. Note that in the 4-D $\delta B/B=0.3\,\%$ case there is no sample within the $25\,\%$ range above the optimum.

First, considering the 2-D optimization, figure 6(a), we find that the relative extent of the point clouds is significantly larger in the $n_{{\rm Ne},{\rm inj}}$ direction, than in the $n_{{\rm D,inj}}$ direction; for instance in the $\delta B/B=0.2\,\%$ case $n_{{\rm Ne},{\rm inj}}$ spans more than an order of magnitude, while $n_{{\rm D,inj}}$ spans only a bit more than a factor of two. In practice it translates to the need of a higher precision concerning the injected amount of deuterium than that of neon. The negative correlation between $n_{{\rm D,inj}}$ and $n_{{\rm Ne},{\rm inj}}$ seen from the arrangement of the point cloud indicates that there are similarities in the effects of these two injected species. These features are also reflected in the favourable valleys (blue tone regions) seen in figure 1. The favourable parameter range indicated by the point clouds shrinks with $\delta B/B$. Note that the region covered by the optima at different $\delta B/B$ values is even smaller than the smallest (black) point cloud; thus we should not read much into how the actual location of the optima varies with $\delta B/B$.

In the 4-D optimization, the resulting point clouds are more scattered, when projected into the $n_{{\rm D,inj}}$–$n_{{\rm Ne},{\rm inj}}$ subspace, see figure 6(b). If anything, there is still a weak anti-correlation between the injected quantities, but the poor statistics makes it less clear. Similarly to the 2-D optimization, the range covered in $n_{{\rm D,inj}}$ is smaller than that in $n_{{\rm Ne},{\rm inj}}$. We can also see that there are no cases within a relative range of $25\,\%$ of the optimum for $\delta B/B=0.3\,\%$. In addition, the optimum itself appears far in the parameter space from the other three overlapping clouds. Namely, it appears at the highest $n_{{\rm D,inj}}$ and lowest $n_{{\rm Ne},{\rm inj}}$ values. We omit this outlier case in the following discussion, but will return to it at the end of this section.

The point clouds occupy the relatively narrow $c_{\rm D}\in [-1.5,\ 1.2]$ range, as seen in figure 6(c), corresponding to modest profile variation. We find a positive correlation between $n_{{\rm D,inj}}$ and $c_{\rm D}$. It means that higher injected content corresponds to more edge-localized peaking. In particular, the injected densities at the plasma centre occupy a narrower range than at the edge (see table 2); apparently, the deuterium density value at the edge is less important. We also observe that lower $\delta B/B$ corresponds to higher $c_{\rm D}$ and $n_{{\rm D,inj}}$ values.

Table 2. Total hydrogenic (including the background) and neon densities at the plasma centre ($r=0$) and at the edge ($r=a$) in the 4-D optimization in the various $\delta B/B$ cases.

For the injected neon profiles, a strong outward peaking is preferred, with values of $c_{\rm Ne}\in [5,\ 10]$, as seen in figure 6(d). The total injected quantities are typically higher than those in the 2-D optimization, covering mostly the $n_{{\rm Ne},{\rm inj}}\in [10^{19}, 10^{20}]\, {\rm m}^{-3}$ range – an order of magnitude higher than in two dimensions. It is interesting to note that, similarly to deuterium, there is a positive correlation between $c_{\rm Ne}$ and $n_{{\rm Ne},{\rm inj}}$.

To understand why the optima in the 4-D optimization perform better than those of two dimensions, we compare the respective $\delta B/B=0.5\,\%$ cases, where the figures of merit are most disparate. The hydrogenic (blue curves) and neon (red) density profiles of the 2-D (dashed curves) and 4-D (solid) optima, are shown in figure 7(a). For deuterium, the 4-D optimization finds a moderate inward peaking ($c_{\rm D}=-0.92$), while the neon profile is strongly peaked at the edge ($c_{\rm Ne}=8.67$), covering a density range over three orders of magnitude.

Figure 7. Comparison of the optimal cases in the 2-D (dashed curves) and the 4-D (solid curves) optimization for $\delta B/B=0.5\,\%$. (a) Radial total hydrogenic density, $n_{{\rm D+T+D,inj}}$ (blue), and neon density, $n_{\rm Ne}$ (red). (b) The RE current density profiles taken at the time point when the total RE current takes its maximum, $t=42\,{\rm ms}$ ($50\,{\rm ms}$) in the 2-D (4-D) case. (c) Time evolution of the heat loss power in the first millisecond, when most of the thermal energy is lost from the plasmas (note the log scale). Blue curves represent the transported heat losses, red curves are the radiated losses.

The neon content has two major effects on our figures of merit. An increasing neon concentration corresponds to a lower quasi-equilibrium temperature during the RE conversion, typically leading to higher final RE currents. This is the same trend that we have witnessed moving from C4 to C6. At the same time, a higher neon concentration can help increase the radiated fraction of heat losses (this was also clear when comparing C5 with C6). However, the neon concentration affects the final RE current most strongly where the RE growth is strongest. This happens to be the plasma core in the parameter region of interest, without a radial variation of the neon density. In addition, to achieve a low $\eta _{\rm cond}$ value it is sufficient to have enough radiating impurities in the edge. Both requirements can be satisfied by an outward peaking neon concentration, which is indeed what the 4-D optima tend to develop.

We find that the 2-D optimum produces a centrally peaked RE current, as seen in figure 7(b) (dashed curve), while the 4-D optimum has a RE profile peaked off axis (solid curve), as expected for the low core concentration of neon. We note that, in this case, the runaway and the centrally peaked ohmic currents decay together after the RE current reaches its maximum, and only towards the end of the simulation ($\approx 140\,{\rm ms}$) does the total current profile become truly hollow.Footnote ¹

The time evolution of the volume-integrated heat losses is shown in figure 7(c) in the first millisecond. This is when the vast majority of the thermal energy content of the plasma is lost, while the fraction of magnetic-to-thermal energy conversion is still negligible. Again, the dashed curves correspond to the 2-D optimum; in this case, the transported loss (dashed blue) reaches comparable values to the radiated losses (dashed red). The entire energy loss process varies relatively smoothly over the plotted time range. In contrast, in the 4-D optimum case the transported heat losses (solid blue) are approximately two orders of magnitude lower than the radiated losses (solid red), and both of these channels have a strong peak at $t=0$, related to the ionization and equilibration of the injected material.

After having discussed the representative behaviour at the optima we return to the analysis of the outlier case, the 4-D optimum at $\delta B/B=0.3\,\%$. In this case the injected neon density is roughly three orders of magnitude lower than in the other three cases, and as such, it exhibits reheating following the TQ in the plasma centre. This reheated region supports a relatively slowly decaying ohmic current, hence the CQ time is on the long side $t_{\rm CQ}=123\,{\rm ms}$ (while still tolerable). The slowly decaying ohmic current and the high value of the effective critical electric field $E_{c}^{{\rm eff}}$, owing to the high $n_{{\rm D,inj}}$, lead to the RE growth stopping just before the RE current grows to macroscopic values. The strong dilution is able to rapidly reduce the temperature to sufficiently low values at the edge, so that even in the presence of a low neon content the cooling can continue to ${\approx }1 \,{\rm eV}$. As then most of the heat transported to the edge is radiated away by the recombined deuterium, the resulting transported heat loss also remains small in this case. This is a fragile case nevertheless; indeed there is no sample within $25\,\%$ of the $\mathcal {L}$ value reached by this optimum. Some parameter combinations in the vicinity of this optimum yield a behaviour reminiscent of C3, with an extremely rapid RE conversion and then a strongly decaying RE current. Thus, even though this optimum performs better than the other three cases in four dimensions, it should not be targeted in a experiment, due to the lack of robustness.

Finally, we comment on the numerical efficiency of the Bayesian approach. We estimate that to achieve a similar level of resolution in the regions that contain samples within $25\,\%$ of the optima would require more than $12\,000$ points in two dimensions and $800\,000$ points in four dimensions, should we decide to use equidistant scans over the entire search domains. These estimates are based on the average minimum distance between samples (in the search space mapped to the unit hyper-cube). As a reference, we use only 420 samples in both the 2-D and the 4-D optimizations. In uninteresting regions with high cost function values the resolution is much lower.

The Bayesian results can be confirmed with calculations on a uniform grid. In a detailed study of a similar problem presented by Bergström & Halldestam (Reference Bergström and Halldestam2022), it was shown that the mean function obtained by the Gaussian process regression accurately recovered the cost function calculated on a uniform grid in the vicinity of the optimum, and showed a good agreement even in regions with high cost values. In terms of finding the global optimum the Bayesian method outperformed Powell's method (Powell Reference Powell1964).