2.1. Tritium beta decay

The kinetic source rate due to tritium decay rate $(\mathrm {d} f/\mathrm {d} t)_{\rm T}$ can be derived using the expression for the energy spectrum of electrons emitted through beta decay (Krane Reference Krane1988)

(2.1)\begin{equation} f_\beta(W) \propto F(p,2)p \mathcal{W} (W_{\rm max}-W)^2\quad \text{for }W\le W_{\rm max}, \text{ and }f_\beta(W)=0\text{ above} \;W_{\rm max}. \end{equation}

Here $\mathcal {W}={m_{\rm e}} c^2\gamma$ is the total energy, $W={m_{\rm e}} c^2(\gamma - 1)$ is the kinetic energy, $\gamma =\sqrt {p^2+1}$ is the Lorentz factor and $p$ is the momentum normalized to ${m_{\rm e}} c$. Since the maximum energy of the emitted electrons is $W_{\rm max}=18.6\,{\rm keV}\ll {m_{\rm e}} c^2$, we may take the non-relativistic limit of the Fermi function (Fermi Reference Fermi1934),

(2.2) \begin{equation} F(p, Z_{\rm f})= \frac{2{\rm \pi}\alpha Z_{\rm f}/\beta} {1-\exp({-}2{\rm \pi}\alpha Z_{\rm f}/\beta)}, \end{equation}

with the normalized speed $\beta =p/\gamma$, the charge of the final state nucleus $Z_{\rm f}$ (for tritium decay $Z_{\rm f}=2$) and the fine structure constant $\alpha \approx 1/137$.

Since the source is isotropic, we may write

(2.3)\begin{equation} 4 {\rm \pi}p^2 \,\mathrm{d} p \left(\frac{\mathrm{d} f}{\mathrm{d} t}\right)_{\rm T} = \mathrm{d} W f_\beta(W). \end{equation}

Using $\mathrm {d} W / \mathrm {d} p={m_{\rm e}} c^2\beta$, we find

(2.4)\begin{equation} \left(\frac{\mathrm{d} f}{\mathrm{d} t}\right)_{\rm T} = \frac{1}{4{\rm \pi} p^2} f_\beta(W){m_{\rm e}} c^2\beta \propto \frac{1}{p^2} \frac{p \gamma (\gamma_{\rm max}-\gamma)^2} {1-\exp({-}4{\rm \pi}\alpha /\beta )}. \end{equation}

The requirement that the source integrated over the entire momentum space yields the tritium decay rate, $(\ln 2)\, n_{\rm T}/\tau _{\rm T}$ – where $\tau _{\rm T}\approx 4500$ days is the half-life for tritium – determines the proportionality factor in (2.4) such that

(2.5)\begin{equation} \left(\frac{\mathrm{d} f}{\mathrm{d} t}\right)_{\rm T}\approx C\frac{\ln2}{4{\rm \pi}} \frac{n_{\rm T}}{\tau_{\rm T}}\frac{1}{p^2} \frac{p \gamma \left(\gamma_{\rm max}-\gamma\right)^2} {1-\exp({-}4{\rm \pi}\alpha/\beta)}, \end{equation}

with the constant

(2.6)\begin{equation} C = 1\left/ \int_0^{p_{\rm max}}\frac{1}{p^2} \frac{p \gamma \left(\gamma_{\rm max}-\gamma\right)^2} {1-\exp({-}4{\rm \pi}\alpha/\beta)} p^2\,\mathrm{d} p\right.\approx 31\,800, \end{equation}

where $p_{\rm max}$ is the momentum corresponding to $W_{\rm max}$.

2.2. Compton scattering rate

The kinetic Compton scattering rate can be calculated using the Compton cross-section and the gamma energy flux spectrum emitted by a given wall material. When using the kinematic relation between the photon energy $W_\gamma$, and the energy $W$ and deflection angle $\theta$ of the scattered electron, we can write

(2.7)\begin{equation} \cos\theta=1-\frac{{m_{\rm e}} c^2}{W_\gamma}\frac{W}{W_\gamma'}, \end{equation}

where $W_\gamma ' = W_\gamma - W$ is the scattered photon energy. The Compton scattering process is described by the Klein–Nishina differential cross-section

(2.8)\begin{equation} \frac{\mathrm{d} \sigma}{\mathrm{d} \varOmega}= \frac{r_{\rm e}^2}{2} \frac{W_\gamma'^2}{W_\gamma^2} \left[\frac{W_\gamma}{W_\gamma'}+\frac{W_\gamma'}{W_\gamma} -\sin^2(\theta)\right], \end{equation}

with $r_{\rm e}=e^2/4{\rm \pi} \varepsilon _0m_{\rm e}c^2$ being the classical electron radius.

Thus, in a plasma with the total electron density $n_{{\rm e},{\rm tot}}$, i.e. the electron population of both free and bound electrons, and isotropic photon distribution, we can write

(2.9)\begin{equation} 4 {\rm \pi}p^2 \,\mathrm{d} p \left(\frac{\mathrm{d} f}{\mathrm{d} t}\right)_{\rm C} = n_{{\rm e},{\rm tot}} \int_{W_{\gamma0}}^\infty \,\mathrm{d} W_\gamma \left[\varGamma_\gamma(W_\gamma) \left(2{\rm \pi}\sin\theta \frac{\mathrm{d} \sigma}{\mathrm{d} \varOmega} \right) \,\mathrm{d}\theta(W,W_\gamma)\right], \end{equation}

where $\varGamma _\gamma (W_\gamma )$ is the gamma energy flux spectrum. The lower integration limit $W_{\gamma 0}$ can be derived from (2.7) as $W_{\gamma 0}=(p+\gamma -1)/2$.

Martín-Solís et al. (Reference Martín-Solís, Loarte and Lehnen2017) estimated the photon flux energy spectrum in ITER by $\varGamma _\gamma (W_\gamma )=\varGamma _0\exp [-\exp (-z)-z+1]$, where ${z = [\ln (W_\gamma /10^6)+1.2]/0.8}$ and

(2.10)\begin{equation} \varGamma_0=\varGamma_{\rm flux}\displaystyle \left/ \int_0^\infty \exp[-\exp({-}z)-z+1]\, \mathrm{d} W_\gamma\right. =\varGamma_{\rm flux}/5.8844, \end{equation}

with the total flux ${\varGamma _{\rm flux}=10^{18}}\,({\rm m}^2\,{\rm s})^{-1}$. This applies for an H-mode discharge at 15 MA and 500 MW fusion power. For an L-mode case, when the plasma energy is lower, the total flux would be approximately 25 % of that of the H-mode (Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017). Part of the flux is caused by the time-integrated activation of the walls, and the typical ratio between the prompt gamma flux from fusion neutrons and the gamma flux after the fusion reactions cease is $1000$. These calculations were performed for a beryllium first wall. ITER currently plans to have a first wall made of tungsten, which might change the photon flux and spectrum.

The runaway production rate from Compton scattering can be obtained as

(2.11a)\begin{align} \left(\frac{\mathrm{d} f}{\mathrm{d} t} \right)_{\rm C} & =\dfrac{n_{{\rm e},{\rm tot}}}{2}\dfrac{1}{p^2}\int_{W_{\gamma0}}^\infty\,\mathrm{d} W_\gamma \varGamma_\gamma(W_\gamma) \sin\theta \frac{\mathrm{d} \sigma}{\mathrm{d} \varOmega}\frac{\mathrm{d} \theta}{\mathrm{d} p} \end{align}

(2.11b)\begin{align} & =\dfrac{n_{{\rm e},{\rm tot}}}{2}\dfrac{1}{p^2}\int_{W_{\gamma0}}^\infty\,\mathrm{d} W_\gamma \varGamma_\gamma(W_\gamma)\frac{\mathrm{d} \sigma}{\mathrm{d} \varOmega} \dfrac{\beta}{\left(\dfrac{W_\gamma}{{m_{\rm e}} c^2}+1-\gamma\right)^2}, \end{align}

where we have used

(2.12) \begin{equation} \frac{\mathrm{d} \theta}{\mathrm{d} p}= \frac{\mathrm{d} }{\mathrm{d} p} \arccos\left(1 - \dfrac{{m_{\rm e}} c^2/W_\gamma}{\dfrac{W_\gamma}{{m_{\rm e}} c^2(\gamma-1)}-1} \right) =\dfrac{p / \gamma} {\sqrt{1-\cos^2\theta}\left(\gamma-1\right)^2 \left(\dfrac{W_\gamma/{m_{\rm e}} c^2}{\gamma-1}-1\right)^2}. \end{equation}

Both the kinetic tritium and Compton sources are consistent with the fluid runaway generation rate (Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017)

(2.13)\begin{equation} \left(\frac{\mathrm{d} n_{\rm re}}{\mathrm{d} t}\right)_{{\rm T}/{\rm C}}= \int_{p>p_{\rm c}}\,\mathrm{d}^3\boldsymbol{p}\left(\frac{\mathrm{d} f}{\mathrm{d} t}\right)_{{\rm T}/{\rm C}}, \end{equation}

when the runaway region is set to be above the critical momentum for RE generation $p_{c}$.

3.1. Simulation set-up

The simulations are performed using the numerical tool Dream (Disruption Runaway Electron Analysis Model). Dream simulates the plasma evolution of tokamak disruptions self-consistently, including the evolution of the plasma temperature, current densities, electron and ion densities, ion charge states, and electric field. The simulations are performed with a multi-fluid plasma model, with the possibility of simulating some, or all, subsets of the total electron population kinetically. Dream divides the electrons into up to three populations.

Both in the fluid and the kinetic simulations performed here, the bulk and RE populations are treated as distinct fluids. The bulk electrons are characterized by their density ${n_{\rm e}}$, temperature ${T_{\rm e}}$ and the Ohmic current density ${j_\Omega }$. The REs are described by their density ${n_{\rm re}}$, which also yields the runaway current density ${j_{\rm re}}=ec{n_{\rm re}}$, because they are assumed to move with the speed of light parallel to the magnetic field. In reactor-scale tokamak disruptions, this approximation, i.e. $v_{\rm re}=c$, is typically valid, as supported by Buchholz (Reference Buchholz2023). The total current density ${j_\Omega } + {j_{\rm re}}$ is governed by the evolution of the poloidal flux $\psi (r)$ (Hoppe et al. Reference Hoppe, Embréus and Fülöp2021a). For the electrical conductivity $\sigma$ which enters in this evolution, we use the model described by Redl et al. (Reference Redl, Angioni, Belli and Sauter2021) that is valid for arbitrary plasma shaping and collisionality, and takes into account the effects of trapping. Neutral collisions – which are only relevant at very low post-TQ temperatures that yield strongly unfavourable outcomes – are not accounted for in the resistivity, as they are not relevant in the parameter ranges of interest (Vallhagen et al. Reference Vallhagen, Embréus, Pusztai, Hesslow and Fülöp2020).

In the kinetic simulations, there is an additional population of superthermal electrons that are represented by the distribution function ${f_{\rm hot}}$, which is evolved using a pitch-angle averaged Fokker–Planck equation (Hoppe et al. Reference Hoppe, Embréus and Fülöp2021a). This division of the hot and cold populations is especially suitable when there are two distinct Maxwellian electron populations present in the plasma, as is the case when cold materials are injected into a hot plasma during massive material injection.

The temperature of the cold population is determined by a time dependent energy balance equation, accounting for Ohmic heating, line and recombination radiation, and bremsstrahlung, as well as a radial heat transport, according to (43) of Hoppe et al. (Reference Hoppe, Embréus and Fülöp2021a). The ionization, recombination and radiation rates are taken from the ADAS database for neon and the AMJUEL databaseFootnote ¹ for deuterium. AMJUEL contains coefficients which account for opacity to Lyman radiation, which is important at high deuterium densities (Vallhagen et al. Reference Vallhagen, Pusztai, Hoppe, Newton and Fülöp2022).Footnote ² The cold electron density is constrained by the condition of quasi-neutrality.

In the kinetic model, the Dreicer and hot-tail generation mechanisms are naturally present as velocity space particle fluxes, while in the fluid model, they are both implemented as generation rates. The Dreicer generation rate is computed by a neural network (Hesslow et al. Reference Hesslow, Unnerfelt, Vallhagen, Embréus, Hoppe, Papp and Fülöp2019b) if the input parameters are within the training interval. However, for low electric fields, outside the training interval ($E/E_{\rm D}<0.01$), the values given by the neural network are inaccurate. In such cases, an analytic extrapolation is performed to recover a similar asymptotic behaviour as to the Connor–Hastie formula (Connor & Hastie Reference Connor and Hastie1975), i.e. $\mathrm {d} {n_{\rm re}} / \mathrm {d} t \propto \exp (-A (E_{\rm D}/E)^b)\rightarrow 0$ as $E_{\rm D}/E\rightarrow \infty$. Here, $A$ and $b$ are positive parameters fitted to the dependence of the neural network on the electric field for $E/E_{\rm D} > 0.01$, with the density and temperature fixed at the values of the current time step. They are fitted under the constraint that $(\mathrm {d} {n_{\rm re}} / \mathrm {d} t)_{\rm Dreicer}(E/E_{\rm D})$ should be once differentiable for all values of $E/E_{\rm D}$.

The hot-tail generation rate is evaluated from the model distribution function of Smith & Verwichte (Reference Smith and Verwichte2008) and using the formula derived by Svenningsson (Reference Svenningsson2020) – it is presented in appendix C.4 of Hoppe et al. (Reference Hoppe, Embréus and Fülöp2021a). In both fluid and kinetic models, the avalanche generation is determined by the model by Hesslow et al. (Reference Hesslow, Embréus, Vallhagen and Fülöp2019a).

In the fluid simulations, the activated generation sources from tritium decay and Compton scattering are evaluated according to Vallhagen et al. (Reference Vallhagen, Embréus, Pusztai, Hesslow and Fülöp2020); in the kinetic simulations, we use the sources derived in § 2. During the TQ, the photon flux is set to ${\varGamma _{\rm flux}=10^{18}}\,({\rm m}^2\,{\rm s})^{-1}$ to be consistent with the prediction of the photon flux at ITER for an H-mode discharge at 15 MA and 500 MW fusion power (Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017). However, with the drop in plasma temperature, the fusion power ceases. For this reason, the photon flux after the TQ is decreased by a factor of $10^3$, as was done by Vallhagen et al. (Reference Vallhagen, Hanebring, Artola, Lehnen, Nardon, Fülöp, Hoppe, Newton and Pusztai2024).

The disruption simulations are performed in an ITER-like tokamak scenario with either a deuterium plasma of density $10^{20}\,{\rm m}^{-3}$ for the non-activated scenarios or a deuterium-tritium plasma with the same total ion density and 50–50 % isotope concentrations. In both scenarios, the fuel density is spatially uniform. Discharges without tritium are simulated for 200 ms, while activated discharges are simulated for 400 ms, to capture the interesting runaway current dynamics in both scenarios.

The magnetic field ${B_0=5.3\,{\rm T}}$ on axis, the resistive wall time $\tau _{\rm w}=0.5\,{\rm s}$, the major radius $R_0=6\,{\rm m}$ and the plasma minor radius $a=2\,{\rm m}$. Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023) evaluated the effective radius of the first toroidally closed conducting structure by matching the poloidal magnetic energy to that in simulated ITER discharges, yielding $b=2.833\,{\rm m}$. However, during disruptions, the plasma can be vertically displaced, resulting in a continuous contraction of the confined region and a corresponding reduction of the poloidal magnetic field energy reservoir that contributes to runaway generation. In this case, a more tightly fit wall ($b=2.15\,{\rm m}$) might be more appropriate. Simulations are therefore performed with both of these values of the conducting wall radius to bracket a plausible range of magnetic to kinetic energy conversion. Furthermore, the geometry is determined by realistic and radially varying shaping parameters congruent with the Miller model equilibrium (Miller et al. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998), used also by Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023).

At $t=0$, the plasma current ${I_{\rm p}}=15\,{\rm MA}$ with current density profile $\hat {j}(r)=[1-(r/a)^2]^{0.41}$, the plasma temperature is parabolic with $T(r=0)=20\,{\rm keV}$ and the plasma is fully ionized. The deuterium and neon of the MMI are assumed to consist of cold atom populations, deposited instantaneously and homogeneously, as was done by Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023).

To account for the effect of the stochastic magnetic field during the TQ, we employ spatially and temporally constant magnetic perturbations of ${\delta B/B}=0.5\,\%$. Such a value of ${\delta B/B}$ is consistent with representative transport levels in MHD simulations of ITER disruptions (Hu et al. Reference Hu, Nardon, Hoelzl, Wieschollek, Lehnen, Huijsmans, van Vugt and Kim2021). Although we do not consider different values of ${\delta B/B}$ here, we note that lower values tend to reduce conducted heat losses as well as the runaway seed, thereby increasing the region in parameter space that yields tolerable disruptions, while the location of the optimum is less affected, consistent with the findings of Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023). Following the Rechester–Rosenbluth model (Reference Rechester and Rosenbluth1978), this magnetic perturbation leads to a transport of hot electrons, REs and electron heat, with a diffusion coefficient $D\propto {\left |v_\parallel \right |} R_0({\delta B/B})^2$. Here, ${\left |v_\parallel \right |}$ is replaced by ${v_{\rm te}=\sqrt {2{T_{\rm e}}/m_{\rm e}}}$ (local electron thermal speed) for the heat transport and by the speed of light $c$ for the RE transport, where we assume a parallel streaming of the REs along the perturbed field lines. As noted by Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023), this approach gives an upper bound on the effect of the RE transport as it does not account for the energy and angular dependence of the RE distribution, as well as finite Larmor radius and orbit width effects. In the kinetic simulations, the parallel velocity dependence for the hot electrons is determined by ${f_{\rm hot}}$.

The transport of superthermal and runaway electrons is switched off during the current quench to signify that the magnetic flux surfaces are reformed. We still use a small but finite heat transport corresponding to ${\delta B/B}=0.04\,\%$ to avoid the development of non-physical narrow hot Ohmic channels (Putvinski et al. Reference Putvinski, Fujisawa, Post, Putvinskaya, Rosenbluth and Wesley1997; Fehér et al. Reference Fehér, Smith, Fülöp and Gál2011). As long as such channels do not try to form, this low transport level is subdominant to radiative heat losses at post-TQ temperatures. In the absence of a clear, widely accepted definition of the end of the TQ, in the simulations, we define the end of the thermal quench by the mean temperature falling below 20 eV. We discard simulations where the TQ condition has not been fulfilled within 20 ms, which is a factor of $\sim 10$ greater than the anticipated TQ time of 1–2 ms according to Hollmann et al. (Reference Hollmann, Aleynikov, Fülöp, Humphreys, Izzo, Lehnen, Lukash, Papp, Pautasso and Saint-Laurent2015). In such cases, we can assume that the TQ was not successfully completed and other complications might arise, such as reheating and unsuccessful current quench.

We note that the criterion used to switch the energetic electron transport off and reduce the heat transport is a source of uncertainty, as the time it takes for flux surfaces to reform might be delayed somewhat after the temperature has reached its post-TQ value. Clearly, longer windows of runaway loss are expected to reduce the final RE current. However, the presence of a relatively small reformed region may be sufficient to support significant RE currents (Tinguely et al. Reference Tinguely, Pusztai, Izzo, Särkimäki, Fülöp, Garnier, Granetz, Hoppe, Paz-Soldan and Sundström2023b).

3.2. Cost function

The cost function of an optimization problem is critical for the optimization to yield meaningful results and, here, the cost function should reliably quantify the risks associated with tokamak disruptions, in the form of a single scalar. First, the runaway current should be accounted for in the cost function. We will consider two plausible choices for representing the runaway current contribution to the cost function – the maximum runaway current ${I_{\rm re}^{\rm max}}=\max _{t}{I_{\rm re}}(t)$ and the runaway current at the time when it makes up 95 % of the plasma current ${I_{\rm re}^{95\,\%}}={I_{\rm re}}(t_{{I_{\rm re}}=0.95{I_{p}}})$, which we will call the 95 % runaway current. To account for the advantages and disadvantages of both options (discussed in Appendix A), we choose to use the first occurrence of ${I_{\rm re}^{95\,\%}}$ unless ${I_{\rm re}^{\rm max}}$ occurs before ${I_{\rm re}^{95\,\%}}$ or if ${I_{\rm re}}(t)<0.95{I_{p}}(t)$ throughout the simulation. This will be called the representative runaway current ${I_{\rm re}^{\rm repr}}$.

Additionally, any remaining Ohmic current at the end of the simulation has the potential of being converted into runaway current, as well as indicating if the termination of the discharge was successful. Consequently, it should also be included in the cost function. Based on Lehnen & the ITER DMS task force (Reference Lehnen2021), the upper safe operational limit of the runaway current is 150 kA, which will also be used as the upper safe operational limit for the Ohmic current.

Another indication of the success of the disruption mitigation system is the duration of the CQ (CQ time). In this work, the CQ time is evaluated through extrapolation, as ${\tau _{\rm CQ}}=(t_{{I_\Omega }=3\,{\rm MA}}-t_{{I_\Omega }=12\,{\rm MA}})/0.6$ if the Ohmic current drops below 3 MA during the simulation and, otherwise, ${\tau _{\rm CQ}}=(t_{\rm final}- t_{{I_\Omega }=12\,{\rm MA}})/(0.8-{I_\Omega ^{\rm final}}/I_{\rm p}^{t=0})$. Here $t_{{I_\Omega }=3\,{\rm MA}}$ ($t_{{I_\Omega }=12\,{\rm MA}}$) is the time when ${I_\Omega }$ is 20 % (80 %) of the initial plasma current. As stated by Hollmann et al. (Reference Hollmann, Aleynikov, Fülöp, Humphreys, Izzo, Lehnen, Lukash, Papp, Pautasso and Saint-Laurent2015), its safe operational interval is $[50, 150]\,{\rm ms}$ for 15 MA discharges.

Finally, the last figure of merit considered is the transported heat loss fraction ${\eta _{\rm cond}}$, meaning the fraction of the initial plasma kinetic energy which has been lost from the plasma due to energy transport. It is accounted for using the Rechester–Rosenbluth transport in the collisionless limit, which means that the radial heat diffusion is caused by parallel streaming of electrons along perturbed field lines (Rechester & Rosenbluth Reference Rechester and Rosenbluth1978), as detailed in § 3.1. For a safe disruption, the transported heat load fraction should be below 10 %, according to Hollmann et al. (Reference Hollmann, Aleynikov, Fülöp, Humphreys, Izzo, Lehnen, Lukash, Papp, Pautasso and Saint-Laurent2015).

Our cost function ${\mathcal {L}}({I_{\rm re}^{\rm repr}}, {I_\Omega ^{\rm final}}, {\tau _{\rm CQ}}, {\eta _{\rm cond}})$ is designed to yield a value below one if all components are within their intervals of safe operation. To achieve this, ${I_{\rm re}^{\rm repr}}$, ${I_\Omega ^{\rm final}}$ and ${\eta _{\rm cond}}$ are normalized to their respective safe operational limits. The CQ time is shifted and normalized according to ${\left |{\tau _{\rm CQ}}-100\,{\rm ms}\right |}/50\,{\rm ms}$ to achieve the same behaviour for both the lower and upper limit of safe operation. The four figures of merit can be combined using different orders of $p$-norms and here we have chosen to use the Euclidean norm ($p=2$). In Appendix A, we discuss the impact of different choices for $p$.

The final consideration in the design of our cost function is the relative importance of the components. With what we have discussed so far, a disruption with ${I_{\rm re}^{\rm repr}}=120\,{\rm kA}$ and ${\tau _{\rm CQ}}=100\,{\rm ms}$ would yield the same value of ${\mathcal {L}}$ as ${I_{\rm re}^{\rm repr}}=0\,{\rm kA}$ and ${\tau _{\rm CQ}}=60\,{\rm ms}$ (if ${I_\Omega ^{\rm final}}=50\,{\rm kA}$ and ${\eta _{\rm cond}}=3\,\%$), whereas the latter would be preferable in an actual disruption. More specifically, the relevant consideration for the CQ time is that it is between 50 ms and 150 ms (and preferably with some margin for uncertainties), not necessarily that it is as close to the middle of this safe interval (i.e. 100 ms) as possible. For the RE current however, it is preferable that it is as close to zero as possible, even if RE currents up to 150 kA are tolerable. Furthermore, minimizing ${I_{\rm re}^{\rm repr}}$ was, in this work, deemed more important than minimizing ${\eta _{\rm cond}}$. To achieve this order of significance when the figures of merit are within their safety intervals, we introduced nonlinear (but smooth) rescaling functions $f_i$, which take the form

(3.1)\begin{equation} f_i= \begin{cases} (x_i)^{k_i} & \text{if } x_i\leq1,\\ k_i(x_i-1)+1 & \text{if } x_i\geq1, \end{cases} \end{equation}

for the currents $k_I=1$ and $x_I=I/150\,{\rm kA}$, for the conducted heat load $k_{\eta _{\rm cond}}=3$ and $x_{\eta _{\rm cond}}={\eta _{\rm cond}}/10\,\%$, and for the CQ time $k_{\tau _{\rm CQ}}=6$ and $x_{\tau _{\rm CQ}}={\left |{\tau _{\rm CQ}}-100\,{\rm ms}\right |}/50\,{\rm ms}$. Note that this implementation of $f_i$ is differentiable.

The final expression for the cost function can thus be written as

(3.2) \begin{equation} {\mathcal{L}}=\sqrt{(C_{I_{\rm re}^{\rm repr}} f_{I_{\rm re}^{\rm repr}})^2+(C_{I_\Omega^{\rm final}} f_{I_\Omega^{\rm final}})^2 +(C_{\eta_{\rm cond}} f_{\eta_{\rm cond}})^2+(C_{\tau_{\rm CQ}} f_{\tau_{\rm CQ}})^2}. \end{equation}

Each of the components are weighted equally, with weight $C=0.5$, to ensure values below one for scenarios of safe operation.

For some simulations, the cost function cannot be evaluated – mainly due to lack of convergence or incomplete thermal quench within 20 ms. In these cases, the cost function was set to a high value (${\mathcal {L}}=75 \gg 1$) to decrease the probability of the optimizer exploring the surrounding area further.

3.3. Optimization set-up

Bayesian optimization is an optimization strategy which uses Bayesian inference to estimate a probability distribution function for the entire objective (or cost) function based on a prior set of samples (input-output pairs). This is achieved by assuming that all unknowns of the system to be inferred are random variables.

To make viable predictions for the objective function, an assumption of the nature of its probability distribution is needed. For this, Bayesian optimization uses stochastic processes, which are infinite collections of random variables (Garnett Reference Garnett2023). A common stochastic process used for Bayesian optimization is the Gaussian process (GP), in which the random variables are distributed according to multivariate Gaussian distributions. In the context of Bayesian optimization, a GP is specified by a mean function $\mu (x)$, which determines the expected objective function value at any point $x$ of the search space, and a covariance function $\phi (x, x')$, which measures the correlation between points $x$ and $x'$ (Garnett Reference Garnett2023). Consequently, this optimization method not only yields an optimum, but also an estimate of the objective function itself in the form of the mean function $\mu (x)$ as well as its uncertainty.

During the optimization, the optimizer advances by taking samples of the objective function and successively updating its sample set and subsequently its GP estimate of the probability distribution function. The strategy for choosing which point to sample next is defined by the acquisition function, which uses the GP estimate to determine which point has the highest potential of improving the current optimum.

In this work, the Bayesian optimization was performed using the Python package by Nogueira (Reference Nogueira2014) which uses GPs. The Matérn kernel with smoothness parameter $\nu =3/2$ was used for the GPs and the expected improvement acquisition function was used as the optimization strategy (Garnett Reference Garnett2023).

4.1. Optimization of non-activated discharges

The non-activated scenarios are characterized by only having the Dreicer and hot-tail generation mechanisms as seed generation sources,Footnote ⁴ which tend to be more susceptible to runaway mitigation measures, as we shall see. For the non-activated scenarios, the optimization was performed on the MMI densities within ${\log ({n_{\rm D}}/[1\,{{\rm m}^{-3}}])\in [20, 22]}$ and ${\log ({n_{\rm Ne}}/[1\,{{\rm m}^{-3}}])\in [15, 20]}$ for the fluid model. For the kinetic model, the intervals were set as ${\log ({n_{\rm D}}/[1\,{{\rm m}^{-3}}])\in [21, 22]}$ and $\log ({n_{\rm Ne}}/[1\,{{\rm m}^{-3}}])\in [17.5, 20]$. The total injected number of atoms of $\sim 10^{25}$ corresponding to the upper end of the deuterium density ranges, $10^{22}\,{\rm m}^{-3}$, is an approximate upper bound on material assimilation using shattered pellet injection at ITER (Vallhagen et al. Reference Vallhagen, Hanebring, Artola, Lehnen, Nardon, Fülöp, Hoppe, Newton and Pusztai2024). When using the fluid (kinetic) model, the optimization was performed using 250 (200) disruption simulations.

Figures 1(a) and 1(b) show the results from the optimization of the non-activated scenario for both the fluid and kinetic models, respectively. The most prominent difference between the results of the fluid and kinetic models is that there exists a region of safe operation (as indicated by blue shades in the figure) for the simulations done with the kinetic model, but not with the fluid model. With the kinetic model, the band of safe operation reaches from $({n_{\rm D}}, {n_{\rm Ne}})=(10^{21.5}, 10^{19.5})\,{\rm m}^{-3}$ to $({n_{\rm D}}, {n_{\rm Ne}})=(10^{22}, 10^{19})\,{\rm m}^{-3}$. Furthermore, the kinetic model yields lower values of the cost function overall compared with the fluid model for ${\log ({n_{\rm D}}/[1\,{{\rm m}^{-3}}])\in [21, 22]}$ and ${\log ({n_{\rm Ne}}/[1\,{{\rm m}^{-3}}])\in [17.5, 20]}$.

Figure 1. Logarithmic contour plots of the cost function estimate $\mu$ for the non-activated scenario, generated using (a) the fluid and (b) the kinetic model in Dream. Note that the colour mapping is adapted such that blue shades represent regions of safe operation. The black star (diamond) indicates the optimal samples found using the fluid (kinetic) model, while the black dots indicate all optimization samples. The grey area covers the region of incomplete TQ.

We find that the main difference between the two models is the runaway current. This is illustrated in figure 2, where the safe regions of each of the cost function components are indicated by the coloured areas. For the fluid model, the safe region of ${I_{\rm re}^{\rm repr}}$ only covers the lower part of the ${\log ({n_{\rm D}}/[1\,{{\rm m}^{-3}}])\in [21, 22]}$ and ${\log ({n_{\rm Ne}}/[1\,{{\rm m}^{-3}}])\in [17.5, 20]}$ square, but for the kinetic simulations, it covers all but the upper right corner. Thus, with the kinetic model, we find that the band of safe operation that can be observed in figure 1(b) corresponds to the overlap of ${I_{\rm re}^{\rm repr}}<150\,{\rm kA}$ and ${\eta _{\rm cond}}<10\,\%$, which are the two competing components. Similarly, the optimum found using the fluid model is located at a point between the region of safety for ${I_{\rm re}^{\rm repr}}$ and ${\eta _{\rm cond}}$, since here there is no overlap. It is worth noting, however, that such a region of overlap of ${I_{\rm re}^{\rm repr}}<150\,{\rm kA}$ and ${\eta _{\rm cond}}<10\,\%$ might not be possible with a tungsten wall, as the additional high-$Z$ impurity can have a more significant negative impact concerning the runaway current (Walkowiak et al. Reference Walkowiak, Hoppe, Ekmark, Jardin, Bielecki, Król, Savoye-Peysson, Mazon, Dworak and Scholz2024) than the improvement the extra radiative losses represent for ${\eta _{\rm cond}}$.

Figure 2. Regions of safe operation (shaded) for the non-activated case with regards to ${I_{\rm re}^{\rm repr}}$ (red), ${I_\Omega ^{\rm final}}$ (blue), ${\tau _{\rm CQ}}$ (green) and ${\eta _{\rm cond}}$ (yellow). Panel (a) uses fluid simulations and panel (b) uses kinetic simulations. The optimal samples are indicated by a star in panel (a) and a diamond in panel (b). Note that only in the kinetic case, there is a finite intersection of all the regions of safe operation.

The shapes of the region of safety for both ${I_\Omega ^{\rm final}}$ and ${\tau _{\rm CQ}}$ are also noticeably different, though not to the same degree as for ${I_{\rm re}^{\rm repr}}$. This is simply a consequence of the difference in runaway current. In cases when less of the Ohmic current is transformed into runaway current, the Ohmic current is higher at the end of the simulation and the CQ time is therefore longer. This causes both of the regions of ${I_\Omega ^{\rm final}}<150\,{\rm kA}$ and $50\,{\rm ms}<{\tau _{\rm CQ}}<150\,{\rm ms}$ to shift upwards when the ${I_{\rm re}^{\rm repr}}<150\,{\rm kA}$ region expands. The choice of plasma model does not have a significant impact on the conducted heat load – for both models, the safe region covers the upper right corner of the region explored with the kinetic model.

The results shown in figures 1 and 2 correspond to simulations with a wall radius of $2.833\,{\rm m}$, but the results are qualitatively similar for a wall radius of $2.15\,{\rm m}$ – the corresponding landscapes as presented in figure 2 are much the same. The main difference is that with the smaller wall radius, the currents are in general lower and the CQ times are shorter. In figure 3, the plasma and runaway currents are presented for the optimum found using the kinetic model, as well as fluid simulations for the same parameters using $a=2.833\,{\rm m}$ and $a=2.15\,{\rm m}$. With the larger wall radius, ${I_{\rm re}^{\rm repr}}\approx 4\,{\rm MA}$, while the smaller wall radius (which also uses a resistive wall time of $0.5\,{\rm s}$) yields ${I_{\rm re}^{\rm repr}}\approx 400\,{\rm kA}$. This dramatic reduction in the representative runaway current caused by placing the conducting wall near the plasma, thereby reducing the magnetic energy reservoir for runaway acceleration, is consistent with recent findings reported by Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023) and Vallhagen et al. (Reference Vallhagen, Hanebring, Artola, Lehnen, Nardon, Fülöp, Hoppe, Newton and Pusztai2024). The reason to consider such a reduced wall radius – compared with the nominal one matching the total available magnetic energy content in ITER – is to emulate the reduction of the magnetic energy that is converted to runaway acceleration in a vertical displacement event, as the confined region is gradually scraped off from the plasma channel. This approach can give a lower bound on the representative runaway current.

Figure 3. Total plasma current (solid) and runaway current (dashed) of the optimum found using the kinetic model, simulated with (a) the kinetic model and (b) the fluid model and $a=2.833\,{\rm m}$, as well as (c) the fluid model and $a=2.15\,{\rm m}$.

The difference in runaway current for the two models is dominated by differences in the modelling of the hot-tail RE-generation mechanism, which is highly dependent on the energy distribution of the superthermal electrons. In the kinetic simulations, ${f_{\rm hot}}$ describes the energy distribution of the superthermal electrons, which enables the hot-tail generation to be included in the velocity flux determined by the Fokker–Planck collision operator. This energy distribution of the superthermal electrons is not evolved in the fluid simulations and, therefore, a model for the energy distribution based on the simulated plasma parameters is needed. Such a model, with the purpose of being used to determine the hot-tail generation rate for fluid simulations, was derived by Smith & Verwichte (Reference Smith and Verwichte2008). They neglect the electric field and diffusion terms in the Fokker–Planck equation, and instead only retain collisional friction, which gives the solution

(4.1)\begin{equation} {f_{\rm hot}^{\rm fl}}(t,p) = \frac{n_0}{{\rm \pi}^{3/2}p_{\rm Te}^3} \exp\left[-\frac{\left(p^3 + 3\tau(t)\right)^{2/3}}{p_{\rm Te}^2}\right], \end{equation}

where $n_0$ is the initial free electron density, $p_{\rm Te} = \sqrt {2T_0/m_{\rm e}c^2}$ is the initial thermal momentum of the electrons normalized to ${m_{\rm e}} c$ (with $T_0$ being the initial electron temperature) and $\tau (t)$ is the time-integrated collision frequency. Using this model distribution function for the superthermal electrons, the hot-tail generation rate of the fluid model is derived by Svenningsson (Reference Svenningsson2020) as

(4.2) \begin{equation} \frac{\partial{n_{\rm re}}}{\partial{t}}={-}4{\rm \pi} {{p}_{\rm c}^{\rm fl}}^2\frac{\partial{{p}_{\rm c}^{\rm fl}}}{\partial{t}}{f_{\rm hot}^{\rm fl}}(t,{p}_{\rm c}^{\rm fl}), \end{equation}

where the critical momentum ${p}_{\rm c}^{\rm fl}$ is defined as the momentum at which the electric field acceleration and collisional friction balance each other; see (C.24) in Hoppe et al. (Reference Hoppe, Embréus and Fülöp2021a).Footnote ⁵

The optimum found using the kinetic model is a clear example of the fluid and kinetic models yielding significantly different RE-current dynamics, as shown in figure 3, due to the difference in hot-tail generation implementation in the two models. For the same MMI densities in a fluid simulation, we obtain a several MA runaway current, where the seed REs are exclusively generated by the hot-tail mechanism, while the kinetic simulation gives practically zero runaway current. In figure 4, the distribution function derived by Smith & Verwichte (Reference Smith and Verwichte2008) and the distribution function evolved in the kinetic simulation are plotted just before and at end of the TQ for $r=0.85\,{\rm m}$, which is the radial position with the highest hot-tail generation rate in the fluid simulation. The two distribution functions ${f_{\rm hot}}$ and ${f_{\rm hot}^{\rm fl}}$ are similar initially, as shown in figure 4(a). The discrepancy at $p\gtrsim 1$ can be fully explained by ${f_{\rm hot}}$ being a Maxwell–Jüttner distribution function, while ${f_{\rm hot}^{\rm fl}}$ is a non-relativistic Maxwell–Boltzmann distribution.

Figure 4. The distribution function for the optimal non-activated case, found using the kinetic model (a) just before and (b) at the end of the TQ. The distribution function ${f_{\rm hot}}$ (dashed) has been evolved in the kinetic simulation and ${f_{\rm hot}^{\rm fl}}$ (solid) is the model distribution used in the fluid simulation, given by (4.1). The latter is based on values of the plasma parameters evolved in the fluid simulation.

For this radius, the hot-tail generation peaks at the end of the TQ in the fluid simulation and then the two distribution functions are significantly different, as shown in figure 4(b). Here, ${f_{\rm hot}^{\rm fl}}\gg {f_{\rm hot}}$ around ${p}_{\rm c}^{\rm fl}$ (by many orders of magnitude; note the log-scale on a wide range), which can be explained by the kinetic simulations accounting for a radial diffusive transport of ${f_{\rm hot}}$ due to magnetic perturbations, while the evaluation of ${f_{\rm hot}^{\rm fl}}$ does not. A generalization of the fluid hot-tail rate, accounting for transport of the hot electrons, would not be straightforward, due to the velocity dependence of the Rechester–Rosenbluth transport. The neglect of this transport will lead to a significant overestimation of the runaway seed due to the factor of ${f_{\rm hot}^{\rm fl}}(t,{p}_{\rm c}^{\rm fl})$ in (4.2). This was confirmed by comparing the fluid and kinetic simulations with the radial transport removed.

As shown in figure 1, the optima found with both the fluid and kinetic models are located at high deuterium densities. As noted in § 3.1, at high deuterium densities, the effects of accounting for opacity to Lyman radiation become important. When not accounting for opacity to Lyman radiation at these optimum, the temperature decays faster, resulting in a significantly larger runaway current, which is in agreement with the conclusions of Vallhagen et al. (Reference Vallhagen, Pusztai, Hoppe, Newton and Fülöp2022).

In conclusion, the model distribution function used by the fluid model lacks electric field acceleration and radial transport effects, and it is the transport that dominates the differences observed between the two models. We note that previous work showed that if the losses due to magnetic perturbations are neglected, i.e. taking into account all the hot-tail electrons, the final runaway current is overestimated by a factor of approximately four in ASDEX Upgrade (Hoppe et al. Reference Hoppe, Hesslow, Embréus, Unnerfelt, Papp, Pusztai, Fülöp, Lexell, Lunt and Macusova2021b), which is in qualitative agreement with our conclusions here.

4.2. Optimization of activated discharges

In activated scenarios, the RE seed is expected to be dominated by the generation of REs from tritium beta decay and Compton scattering in cases when the hot-tail and Dreicer generation sources have been successfully mitigated. The same intervals and number of samples were used for the optimization of the activated scenarios as for the non-activated scenarios, except that the neon density interval was changed to $\log ({n_{\rm Ne}}/[1\,{{\rm m}^{-3}}])\in [15, 17.5]$ for the optimization with the kinetic model to include regions of lower runaway current.

Figure 5 shows the result from the optimizations of the activated scenarios. There are no regions of safe operation, neither when using the fluid nor the kinetic model – in fact, ${\mathcal {L}}>10$ (based on the mean prediction $\mu$) everywhere for both models and both models yield almost identical results. From the optimization using the fluid model, there are two local minima, which both are around ${\mathcal {L}}=15$. Since the cost function was designed to distinguish safe from unsafe disruptions, and this distinction happens around ${\mathcal {L}}\sim 1$, it is not as reliable at determining which of the two disruptions is better if their cost function values are similar and large while the individual figures of merit are different. Therefore, the global minimum (in figure 5a) is not necessarily safer than the other minimum.

Figure 5. Logarithmic contour plots of the cost function estimate $\mu$ for the activated scenario, generated using (a) the fluid and (b) the kinetic model in Dream. The black star (diamond) indicates the optimal samples found using the fluid (kinetic) model, while the black dots indicate all optimization samples.

The optimization using the kinetic model presented in figure 5(b) has not explored the region around the global minimum found with the fluid model, thus its minimum has the same location as the other local minimum found using the fluid model. However, the location of this global optimum did not change when performing a similar optimization using the kinetic model which covered regions surrounding both optima found using the fluid model.

That the two plasma models yield very similar results is also clear from figure 6, where the regions of safe operation for each figure of merit are plotted. The two models yield qualitatively the same results, since the shapes of the safe regions are remarkably alike to the point where it is difficult to distinguish differences due to plasma models from stochastic differences of the sampled locations.

Figure 6. Regions of safe operation (shaded) for the activated case with regards to ${I_{\rm re}^{\rm repr}}$ (red), ${I_\Omega ^{\rm final}}$ (blue), ${\tau _{\rm CQ}}$ (green) and ${\eta _{\rm cond}}$ (yellow). Panel (a) uses fluid simulations and panel (b) uses kinetic simulations. The optimal samples are indicated by a star in panel (a) and a diamond in panel (b). The other markers correspond to the cases in table 1.

In table 1, fluid and kinetic simulations are compared for different combinations of the deuterium and neon densities, including the two optima (the locations in ${n_{\rm D}}$–${n_{\rm Ne}}$ space are indicated by symbols in figure 6). The difference between the kinetic and fluid models are in most cases just a few percent, and the relative difference is larger when the figures of merit are small. The largest relative error occurs for the conducted heat load, since it is consistently a few percent larger in the kinetic case; this can be explained by the kinetic simulation capturing the transport of the hot electrons.Footnote ⁶

Table 1. Disruption figures of merits for fluid and kinetic simulations for parameters indicated in figure 6 with the same markers.

In absolute terms, only the runaway current varies significantly between the kinetic and fluid results – for large runaway current cases, the difference can be ${\sim }150\,{\rm kA}$, while the relative difference is small. As shown in table 1, the kinetic simulations have a higher runaway current in all cases but one. The reason for the kinetic model yielding larger runaway current is indirectly due to the activated sources. The direct generation of electrons above the critical momentum from the activated sources is practically identical for both the fluid and kinetic simulations. In the kinetic simulations however, the hot electrons generated by the activated sources with momenta below the critical momentum will cause an increase in RE generation due to flux in velocity space.

The dynamics of the RE seed generation in the activated scenarios are substantially different from those in the non-activated scenarios. In the non-activated case, the seed generation was dominated by the hot-tail mechanism. For an activated scenario, ${f_{\rm hot}}$ (which includes not only the hot-tail, but also superthermal electrons generated by the Dreicer, Compton scattering and tritium decay mechanisms) will not be depleted due to magnetic perturbations, because there will be a constant generation of hot electrons due to the activated sources. This means that the runaway seed generation will not necessarily be overestimated with the fluid model, as it was in the non-activated case.

In figure 7, the distribution functions ${f_{\rm hot}}$ and ${f_{\rm hot}^{\rm fl}}$ are plotted at the beginning and end of the TQ, at $r=1.05\,{\rm m}$, for the optimal activated scenario found using the fluid model. We see that, in contrast to the non-activated case, ${f_{\rm hot}^{\rm fl}}\ll {f_{\rm hot}}$ around ${p}_{\rm c}^{\rm fl}$ at the end of the TQ for the activated case (note the scale of the $y$-axis). The fluid simulation should still overestimate the runaway seed generation rate due to the neglect of superthermal electron transport in the fluid hot-tail model, which should lead to ${f_{\rm hot}^{\rm fl}}\gg {f_{\rm hot}}$ as in the non-activated cases. However, the fluid model does not capture the additional dynamics in the superthermal momentum range caused by the generation of hot electrons from the activated sources, as discussed above, and this explains why ${f_{\rm hot}^{\rm fl}}\ll {f_{\rm hot}}$ in the activated scenarios. Furthermore, for the activated cases, the inaccurate runaway seed generation will have a smaller impact on the runaway current evolution since the activated sources, and particularly the Compton source, will constitute a significant fraction of the seed. For both optima found in the activated scenarios, the Compton scattering generation mechanism dominates the runaway seed generation by several orders of magnitude.

Figure 7. The distribution function for the optimal activated case found using the fluid model (a) just before and (b) at the end of the TQ. The distribution function ${f_{\rm hot}}$ (dashed) has been evolved in the kinetic simulation and ${f_{\rm hot}^{\rm fl}}$ (solid) is the model distribution used in the fluid simulation, given by (4.1). The latter is based on values of the plasma parameters evolved in the fluid simulation.

The physical model and methodology of this study are similar to that of Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023), but the cost function landscapes and the found optima are significantly different from those presented in that work. The most significant difference compared with the simulations presented by Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023) is the cost function. Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023) used the following cost function:

(4.3a)\begin{gather} {\mathcal{L}}_{\rm IP}({I_{\rm re}^{\rm max}}, {I_\Omega^{\rm final}}, {\tau_{\rm CQ}},{\eta_{\rm cond}})=\frac{{I_{\rm re}^{\rm max}}}{150\,{\rm kA}} + \frac{{I_\Omega^{\rm final}}}{150\,{\rm kA}} + 100\,\theta({\tau_{\rm CQ}}) + 10 \frac{{\eta_{\rm cond}}}{10\,\%}, \end{gather}

(4.3b)\begin{gather}\theta({\tau_{\rm CQ}})= \tilde{\varTheta}(50\,{\rm ms} - {\tau_{\rm CQ}}) + \tilde{\varTheta}({\tau_{\rm CQ}} - 150\,{\rm ms}), \end{gather}

where $\tilde {\varTheta }(\tau )= [1+\tanh (\tau /3.3\,{\rm ms})]/2$. In (4.3a), ${I_{\rm re}^{\rm max}}$ is used instead of ${I_{\rm re}^{\rm repr}}$, but this does not have a significant impact on the results. Rather, differences between how ${\tau _{\rm CQ}}$ and ${\eta _{\rm cond}}$ are treated are responsible for most of the difference in the results.

First, the step function behaviour of $\theta ({\tau _{\rm CQ}})$ limits the optimization of the other parameters to only the region of safe operation for ${\tau _{\rm CQ}}$. All samples with ${\tau _{\rm CQ}}\notin [50, 150]\,{\rm ms}$ will have ${\mathcal {L}}_{\rm IP}>100$, which is an order of magnitude larger than the optimum found by Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023), making it difficult to distinguish the impact of the other components on the cost function value.

Second, the transported heat fraction is weighted one order of magnitude higher than the runaway current and final Ohmic current. This means that ${\eta _{\rm cond}}=10\,\%$ will have the same impact on ${\mathcal {L}}_{\rm IP}$ as ${I_{\rm re}^{\rm max}}=1.5\,{\rm MA}$, but while ${\eta _{\rm cond}}=10\,\%$ is marginally tolerable, ${I_{\rm re}^{\rm max}}=1.5\,{\rm MA}$ is not.

The cost function differences explain the difference in landscape of ${\mathcal {L}}_{\rm IP}$ compared with ${\mathcal {L}}$ here. Furthermore, it explains why no case with ${I_{\rm re}^{\rm max}}<4\,{\rm MA}$ is found by Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023) – the cases in table 1 with ${I_{\rm re}^{\rm max}}<4\,{\rm MA}$ have ${\tau _{\rm CQ}}>150\,{\rm ms}$ and ${\eta _{\rm cond}}\approx 80\,\%$. This yields ${\mathcal {L}}_{\rm IP}\sim 200$ which is a rather high value of the cost function according to the scale in figure 1 of Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023).

Apart from the difference in cost function landscape, the optimum found using the fluid model is shifted to higher values of ${n_{\rm Ne}}$ and lower values of ${n_{\rm D}}$ compared with the optimum found by Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023). The main reason for this shift is the decrease of the photon flux from ${\varGamma _{\rm flux}=10^{18}}\,({\rm m}^2\,{\rm s})^{-1}$ to ${\varGamma _{\rm flux}=10^{15}}\,({\rm m}^2\,{\rm s})^{-1}$ after the TQ, which was not accounted for by Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023). This change in the photon flux leads to a delay in the runaway current ramp up, resulting in a longer decay time for the Ohmic current. Therefore, ${\tau _{\rm CQ}}$ is increased from 49 ms to 67 ms, and thus would have a large impact also on ${\mathcal {L}}_{\rm IP}$.

An important similarity between our result and those of Pusztai et al. (Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023) is that we do not find any ${n_{\rm D}}$–${n_{\rm Ne}}$ combinations yielding tolerable or close to tolerable figures of merit in the activated case. For the full explored region, ${\mathcal {L}}\gg 1$, and as we mentioned earlier, the cost function was not designed to be particularly informative at these values. To study correlations and trade-offs between different figures of merit, another approach is needed.

We can use that each simulation which we have performed corresponds to a point in the four dimensional parameter space spanned by our four figures of merit, i.e. ${I_{\rm re}^{\rm repr}}, {I_\Omega ^{\rm final}}, {\tau _{\rm CQ}}, {\eta _{\rm cond}}$. In figure 8, these points have been projected on two-dimensional subspaces spanned by pairs of the figures of merit, which illustrates potential trade-offs between each pair. Figure 8(a–c) clarify that for scenarios with ${I_{\rm re}^{\rm repr}}<1\,{\rm MA}$, we get ${\tau _{\rm CQ}}>400\,{\rm ms}$ and ${\eta _{\rm cond}}>75\,\%$. More specifically, the correlation between ${I_{\rm re}^{\rm repr}}$ and ${\eta _{\rm cond}}$ suggests that it is not possible to achieve moderate values for the runaway current and the transported heat fraction simultaneously – when ${I_{\rm re}^{\rm repr}}<4\,{\rm MA}$, then ${\eta _{\rm cond}}>75\,\%$ and vice versa. This is related to the almost rectangular envelope bounding the point cloud from the side of low values in the ${I_{\rm re}^{\rm repr}}$–${\eta _{\rm cond}}$ panel.

Figure 8. Projections of the simulation dataset to all the two-dimensional subspaces of the figure of merit space $({I_{\rm re}^{\rm repr}}, {I_\Omega ^{\rm final}}, {\tau _{\rm CQ}}, {\eta _{\rm cond}})$, corresponding to the optimizations of the activated scenario. Both kinetic (blue) and fluid (red) simulations are shown. The intervals of safe operation for each cost function component is indicated by the black lines. This figure illustrates the trade-off between the different cost function components. The optimal samples found during the optimizations using the fluid (kinetic) model is indicated by a black star (diamond), while the other cases in table 1 are indicated as black markers of different shapes. Samples located outside of the plotted domains are placed at the edges.

Disregarding the runaway current, there are some points of safety overlap for the other three figures of merit; see figure 8(d–f). All simulations with ${\tau _{\rm CQ}}\in [50,150]\,{\rm ms}$ have ${I_\Omega ^{\rm final}}<150\,{\rm kA}$, and since there is a slight overlap between ${\tau _{\rm CQ}}\in [50,150]\,{\rm ms}$ and ${\eta _{\rm cond}}<10\,\%$, here ${I_\Omega ^{\rm final}}<150\,{\rm kA}$ as well. This can be related back to figure 6, where we see a small region of overlap for safe values of ${\tau _{\rm CQ}}$ and ${\eta _{\rm cond}}$, as well as the safe region for ${\tau _{\rm CQ}}$ being almost fully covered by the safe region for ${I_\Omega ^{\rm final}}$.

Being able to visualize correlations between the cost function components and which regions of the ${n_{\rm D}}$–${n_{\rm Ne}}$ space correspond to safe regions of the figures of merits is among the benefits of using a stochastic and exploratory optimization algorithm such as Bayesian optimization. This feature has enabled the exploration of how kinetic and fluid simulations compare for a large variety of scenarios, while simultaneously yielding information on the disruption evolution not just for the most optimal case, but for the full region of interest in the ${n_{\rm D}}$–${n_{\rm Ne}}$ space. Such comprehensive comparison or thorough exploration of MMI of deuterium and neon would not have been possible with a more traditional optimization algorithm, such as gradient descent or the comparison of just a few scenarios.