Prediction of pellet mass thresholds for ELM triggering in low-collisionality, ITER-like discharges

In ITER, pellets are calculated to require more than 8 times the mass than currently planned to reliably trigger edge-localized modes (ELMs). Unmitigated heat flux impulses from ELMs are intolerable in ITER at full power and current. Therefore, ITER operation relies on multiple approaches to control ELM heat fluxes. One method is pellet ELM pacing to instigate small rapid ELMs with low heat flux. Predicting the performance of pellet pacing is critical for ITER, which is expected to operate in a regime with a low-collisionality, peeling-limited pedestal. However, to trigger ELMs the local pressure increase in the expanding pellet cloud pushes the equilibrium over the ballooning stability limit. In this work, linear and nonlinear M3D-C1 simulations are used to predict pellet mass thresholds in DIII-D discharges and ITER scenarios with peeling-limited pedestals. It is found that the distance of the equilibrium’s operational point from the ballooning branch of the pedestal stability boundary strongly changes thresholds. Linear M3D-C1 simulations find a strong dependence of the pellet mass threshold on the poloidal injection location for ITER’s 15 MA, Q = 10 scenario. The required pellet mass at the planned injection locations is 8 to 17 times larger than currently considered. However, such linear simulations do not include pellet ablation physics or time evolution of density and temperature. A new scheme of 2D nonlinear simulations, coupled with linear stability analysis at various steps throughout the nonlinear time evolution, was developed to include such physics and improve on the linear results. These new nonlinear-to-linear simulations confirm previous findings. This result suggests that pellet ELM triggering in ITER could require pellets much larger than those currently planned, which makes ELM-pacing operationally challenging. On the other hand, fueling pellets injected from the high-field side will likely not unintentionally trigger ELMs in an otherwise ELM-stable plasma.


Introduction
Mitigating heat fluxes from edge-localized modes (ELMs) is crucial in future tokamak fusion devices with current higher than most contemporary devices to ensure a reasonable lifetime of plasma facing components.Pellet pacing is an ELM control method designed to instigate small rapid ELMs with low heat flux.Experiments in present day devices with medium to high pedestal collisionality have demonstrated that small pellets can reliably pace ELMs at frequencies much larger than the natural ELM frequency without significantly increasing the edge density [1].Previous studies on DIII-D [2,3], ASDEX Upgrade [4] and JET [5,6] have explored pellet ELM pacing in regimes with high collisionality [7] in the pedestal region, ν * ped > 2, due to the difficulty of achieving low density edge plasma conditions while injecting pellets.It is the general understanding that pellets trigger ELMs by providing a pressure perturbation in the plasma [1,8].Therefore, it is beneficial for pellet ELM triggering to operate close to the ballooning stability limit.ITER, however, is predicted to operate at very low collisionality in the pedestal region, ν * ped ∼ 0.2 [9].So it is likely that the ITER pedestal stability is limited mainly by the kink-peeling branch of the peeling-ballooning instability.
A pellet mass threshold for ELM triggering is observed experimentally in DIII-D in plasmas with low-collisionality pedestals [10].Cryogenically frozen deuterium pellets were injected from the low-field side outboard midplane.The pellet injection frequency was intentionally kept below the frequency of natural ELMs in order to study individual triggering events, and to maintain low-collisionality plasma conditions.Over the range of 10 discharges with a total of more than 180 injected pellets it was found that pellets reliably triggered ELMs only above a certain mass threshold.
Pellet size thresholds for ELM triggering in a plasma with peeling limited pedestal were simulated for the first time [11] using the M3D-C1 code [12].A linear workflow was conceived and applied to DIII-D, showing how thresholds vary with plasma conditions and injection locations.In this workflow linear M3D-C1 simulations provide growth rates and determine the most unstable mode.Three-dimensional nonlinear M3D-C1 simulations confirmed that a pellet larger than the predicted linear threshold successfully triggers an ELM by first destabilizing intermediate-n peeling modes, which then in turn destabilize high-n ballooning modes shortly after.
Using linear and nonlinear M3D-C1 simulations, we predict the critical pellet mass threshold for ELM triggering in low-collisionality, peeling limited discharges for DIII-D and ITER.The latter is expected to operate in such a regime, as would other future devices like SPARC or a tokamak fusion pilot plant.Therefore, it is critical to assess the performance of pellet ELM triggering.In the next section we compare experimentally observed thresholds with simulated ones to provide confidence in our models.In section 3 we introduce a new nonlinear workflow that includes important pellet ablation physics in contrast to previously used linear simulations.In section 4 thresholds for peeling and ballooning limited pedestals are determined.Linear and nonlinear workflows are then applied to ITER in section 5 to find pellet mass thresholds for various poloidal injection locations.Final conclusions are summarized in section 6.

DIII-D experimental threshold and model validation
In DIII-D plasmas with low collisionality, peeling-limited pedestals, cryogenically frozen deuterium pellets were injected at the outboard midplane directed towards the magnetic axis [10].Two sizes of cylindrical pellets were used.A typical pacing-sized pellet is followed 20 ms later by a typical fueling-sized pellet, which has about 3.8 times the mass of the pacing pellet.The true mass of each pellet can vary, due to differences in pellet quality produced by cutting pellets out of an extrusion as well as differential ablation along the guide tube; the latter is only a very weak effect.Note that only pellets with a clear mass measurement by the microwave cavity detector are included in the database.In order to maintain the low-collisionality pedestal conditions pellets are injected every 200 ms, which is larger than the natural ELM cycle.Hence, the experiment is designed to study individual trigger events instead of ELM pacing.Since the injection frequency is much lower than the ELM frequency, pellets arrive at the plasma at various times relative to the previous natural ELM, which impacts the possibility of the pellet to trigger a new ELM.However, a clear separation of successful trigger events from failed ones is found.
As shown in figure 1, using data from [10], of 43 pacingsized pellets (with a relative mass less than 1) in this database that arrived at the plasma within the initial 50% of the natural ELM cycle, only one triggered an ELM.On the other hand, pellets above the observed mass threshold all reliably trigger ELMs.The blue bars in the figure count the number of pellets that fail to trigger an ELM, while the red bars count successful trigger events for each mass.The black dashed line shows the experimentally observed mass threshold.Counting also pellets that arrive later than 50% into the natural ELM cycle following the previous ELM (not shown here) adds more successful trigger events below and above the mass threshold, but failed ones only below the threshold.While small pellets (mass ⩽ 1) become more successful at later arrival times than earlier, failures to trigger are found at all arrival times.The threshold still marks the critical mass needed to reliably trigger an ELM.Keep in mind that ELM pacing aims to trigger ELMs within much less than 50% of the natural ELM cycle.This was demonstrated in [1], where ELMs were paced at up to 60 Hz compared to a natural ELM frequency of 5 Hz, hence ELM triggering occurred within less than 10% of the cycle.Also the pellet velocity was measured during the experiment and it was found that the probability of a pellet triggering an ELM is determined more by pellet mass than by velocity, as all pellets above the mass threshold were successful at triggering ELMs, regardless of velocity, which ranges from 40 m s −1 up to 400 m s −1 .At small pellet mass (⩽1) successful and failed trigger events occur at any velocity, which ranges between 100 m s −1 and 300 m s −1 .More details can be found in [10].
Comparing with the DIII-D experiment, pellet mass thresholds were determined using the linear M3D-C1 workflow for multiple times during discharge 178555 and compared against ELM triggering events at the respective time in the discharge.In M3D-C1 the pellet is introduced as a localized, Gaussian density perturbation added as a source term to the continuity equation in the code.The pellet location, distribution width and height are taken from a pellet ablation model and used to calibrate the input parameters in M3D-C1 so that the density perturbation in the simulation matches the pellet density from the ablation model.Note that in the linear workflow the perturbation is fixed and profiles do not evolve.More details are found in [11] and appendix.In all cases, thresholds simulated by M3D-C1 successfully separate (within error bars) events where a pellet fails to trigger an ELM in the experiment from events where another pellet successfully triggers an ELM. Figure 2 shows multiple pellets injected in shot 178555 as a function of relative pellet size on the x-axis and injection time on the y-axis.The ones marked in blue did not trigger ELMs while the red ones successfully triggered.The pellet size is normalized to the simulated M3D-C1 threshold size at the respective time.Note that the simulated threshold is different for each trigger event and was modeled independently for each event, using five different equilibra, each taken right before an ELM even occurs in the discharge, thereby reflecting the plasma pedestal conditions close to the ELM crash.Yet by normalizing each pellet size to its respective threshold, all modeled thresholds line up at the black dashed line, showing the clear separation of events as seen in the experimental data above.This good agreement provides confidence that the model is capable of predicting thresholds for pellet ELM triggering.Further validation is certainly needed.New experiments at DIII-D are planned, in particular for pellet injection from other poloidal locations than the midplane.

Stability analysis using 2D nonlinear simulations
We use the extended magneto-hydrodynamics code M3D-C1 [12], which solves the single-fluid MHD equations, to simulate ablating pellets in the edge of tokamak ELMy Hmode plasmas with low-collisionality, peeling limited equilibrium pedestal conditions.The code uses a refined mesh that increases significantly in resolution from the plasma core towards the separatrix, then drops in resolution in the scrapeoff layer and extends beyond a resistive vessel wall into an outside region with very low resolution.This outside domain is treated as a vacuum region.The plasma is initialized with an axisymmetric EFIT equilibrium, electron density and temperature profiles and an optional plasma rotation profile.A pellet is introduced into the simulation as an additional spatially localized source in the continuity equation for the density.Linear simulations hold the equilibrium plasma and pellet density fixed and solve for the time evolution of the least stable eigenmode for a given toroidal mode number n. From this the growth rate (positive for unstable modes, negative for stable) and mode structure can be obtained.Such linear simulations are computationally inexpensive and return results relatively quickly.However, the ablation of the pellet and the subsequent evolution of plasma profiles are neglected.Therefore, here we have developed a new procedure that utilizes 2D nonlinear simulations to determine growth rates and pellet size thresholds.This new procedure is then applied to the same DIII-D case as the previous linear workflow, discharge 178555 at 3055 ms.
In a nonlinear M3D-C1 simulation the plasma and pellet are initialized in the same way as in the linear ones, but then evolve over time.Typically, we continue the simulations for up to 100 Alfvén times.During the simulation the pellet moves radially inwards, ablates over time at a relatively constant rate and its density cloud spreads poloidally, while moving with the pellet and drifting due to ∇B drift.Note that in 2D simulations the pellet density is only poloidally localized while kept constant in the toroidal direction.Furthermore, the background plasma evolves with the ablating pellet.The temperature drops along the flux surface at the position of the ablating pellet due to dilution cooling and equilibrates at a value lower than without pellet ablation.The high density but relatively constant temperature at the location of the pellet makes the pressure peak, due to the comparatively fast time scale of electron thermal conduction relative to the ion sound speed of the expanding pellet mass.Further details are published in section 5 of [11].The main difference to a fully 3D nonlinear simulation is that the pellet is not toroidally localized, but such 2D nonlinear runs are still computationally suitable for systematic scans within reasonable CPU-time and cost.

Finding a threshold for pellet ELM triggering
Because the nonlinear simulation evolves the plasma and pellet over time it is insufficient to analyze the stability of the system at any single time.For each toroidal mode the growth rate also evolves over time, and the long-term stability can be estimated by a time average.If the growth rate is positive on average that means the mode grows more than it decays and is therefore unstable.This can be seen in figure 3, which shows the time evolution of the growth rate for an n = 7 mode for various pellet sizes.
The nonlinear simulation is analyzed for stability after each time step τ A = 6.48 • 10 −7 s (the normalized time in M3D-C1 which matches the Alfvén time for n e = 10 20 /m 3 and B t = 1 T), referred to as Alfvén time in the following.The growth rate strongly varies within the first 10-20 Alfvén times but then settles into a more long-term evolution.It is not suitable to determine the long-term stability from a single time slice.A typical ELM event occurs on the order of about 1-2 ms, which puts its initial growth up to its subsequent crash on a time scale of several tens of microseconds.The mode's growth rate needs to be estimated on this time scale.So, in the following we run the nonlinear simulations for at least 100 Alfvén times (about 65 µs), determine the stability for toroidal modes n = 2, 3, 5, 7, . . ., 23, 25 at Alfvén times τ A = 30, 40, 50, . . ., 90, 100 and average the respective growth rates over time.Figure 4 shows the result of repeating this procedure for various pellet sizes, given by their spherical radius.
For very small pellets with radii of 0.36 mm and 0.5 mm the time average growth rates are negative for all toroidal modes, which means the plasma edge remains stable.At 0.72 mm radius the plasma shows positive growth rates for n = 9 and n > 12.For larger pellets the plasma edge becomes more and more unstable.The most unstable mode can be found near n = 23.In the linear workflow, see [11], n = 9 was found to be the most unstable mode for this case, which can still be seen in the time average here as well by the increased growth rate around n = 9 compared to n = 7 and n = 11.However, during the time evolution, which is only captured by the nonlinear workflow, high-n ballooning modes grow significantly stronger, resulting in the large growth rates for n > 20.
Using the data points shown in figure 4 for the various pellet sizes at the most unstable mode we can determine the pellet size threshold for which the plasma edge would first become unstable, hence triggering an ELM.This is shown in figure 5.The top figure shows the interpolation of growth rates for the n = 9 mode with respect to pellet size and finds a critical threshold size, the size where the growth rate crosses into positive, of about 0.586 mm.However, contrary to the linear workflow n = 9 is not the most unstable mode.The bottom figure shows the same analysis for the most unstable mode n = 23, and finds a critical threshold size of 0.637 mm.The objective of this analysis is to determine the threshold for ELM triggering, which requires the most unstable mode to become unstable.Once a pellet is large enough to destabilize the most unstable mode, the instability will likely grow rapidly to an ELM crash.So in the nonlinear workflow the higher-n ballooning modes are much more unstable and will eventually become the dominant modes.Previous linear results determined the critical threshold size as 0.702 mm [11].So the new procedure using 2D nonlinear simulations confirms previous findings, which found that a linear analysis can provide a good first estimate despite its lack of pellet ablation and profile evolution physics.

Poloidal dependence of threshold in DIII-D
In previous work we also studied the dependence of the pellet size threshold on the poloidal injection location.Using linear simulations, it was found that the threshold varies significantly with the poloidal injection location [11].The outboard midplane has the lowest threshold, while the upper and lower x-point regions are the most unfavorable locations for pellet ELM triggering.
In figure 6 the previous linear results are compared with the results using the new procedure.Figure 6(a) shows the poloidal cross-section of the EFIT equilibrium of the respective DIII-D discharge.The red points mark the respective injection locations, which are the same as in the previous linear simulations.However, in the new procedure the pellets move, and the red arrows indicate the radial inwards motion of the pellets in the 2D nonlinear M3D-C1 simulation.Note that the size of the arrows does not reflect the absolute value of the pellet velocity, which can range from 40 m s −1 to 400 m s −1 in the experiment, just the respective direction.As discussed in section 2, the pellet velocity has little impact on the threshold in the DIII-D experiment.So, in the simulation a fixed pellet velocity of 100 m s −1 is used.Figure 6(b) then shows the previous linear results for the pellet size threshold in blue and the results of the new nonlinear procedure in green.On the low field side the nonlinear simulations return similar results as the previous linear ones.As mentioned in the previous section, the pellet size threshold is 0.637 mm at the midplane compared to 0.702 mm from the linear simulations.Near the x-points and even more so on the high field side the nonlinear simulations find significantly larger thresholds than the previous linear simulations.This suggests that the pellet ablation and subsequent drift of the pellet cloud material adds a small beneficial effect on the low field side, but significantly inhibits ELM triggering on the high field side compared to only a peak density perturbation.Hence, larger pellets are needed to trigger an ELM on the high field side.This is different from experiential observations in JET and JOREK simulations for higher pedestal collisionality [6,8].While the next section will show the significance of pedestal collisionality on thresholds for pellet ELM triggering, clearly, more studies are needed here to fully understand the dynamics.New experiments on DIII-D are planned, awaiting the installation of new injector hardware that would enable injection form multiple poloidal angles.

Thresholds of peeling and ballooning limited plasmas
The principle of pellet ELM triggering is that in order to trigger an ELM, the local pressure increase in the expanding pellet cloud pushes the equilibrium over the ballooning stability limit.Nonlinear MHD simulations confirm that the ablating pellet cloud raises the local pressure, despite the dilution cooling.Therefore, if the pedestal stability is naturally limited by ballooning, very small changes in pressure, hence small pellets, are already sufficient to trigger ELMs.This is the typical case in present day machines with medium to high collisionality.However, in future machines like ITER or SPARC the pedestal collisionality is estimated to be quite low, and the pedestal stability is more peeling limited than in present day discharges.The DIII-D discharge 178555 was designed to operate at lower collisionality to have a more peeling limited pedestal than other typical DIII-D discharges, which explains the pellet size thresholds found by M3D-C1 as well as observed during the experiment (see section 2).
To study the effect of peeling versus ballooning limited pedestals on the pellet size threshold we use VARYPED equilibria of discharge 178555 and determine the threshold for 3 different cases as shown in figure 7.All three equilibria are on similar iso-stability surfaces of the pedestal stability diagram calculated by the ELITE code; the inlet at the top left shows the stability diagram with the operating points of the  three cases.The red case is strongly peeling limited.The pellet size threshold for ELM triggering is significantly larger than for the other cases.In the red case a pellet with more than 4 times the mass is required to trigger an ELM than in the black case, which is the original equilibrium, as used in the previous sections.The blue case is ballooning limited and even a very small pellet is already sufficient to trigger an ELM.
This clearly shows that the distance of an equilibrium's operational point from the ballooning branch of the pedestal stability boundary strongly changes the threshold.For the same pellet size, ELMs may be triggered in equilibria closer to the ballooning branch, but are not triggered in equilibria further away.

Pellet mass thresholds in ITER
In this section we apply the linear, as well as the new 2D nonlinear workflow to ITER to find pellet mass thresholds in one of its 15 MA, D-T H-mode, Q = 10 scenarios.The equilibrium is taken from a set of standard scenarios provided by the ITER organization.Here we use the t = 100 s time slice, which is well into the flattop of this steady-state discharge.

Linear simulations
Since the linear workflow can provide a good estimate of thresholds on the low-field side, as has been seen in the DIII-D case, we will use its fast turnaround for a detailed scan of poloidal angles.However, it is likely that the linear results will underestimate thresholds on the high-field-side.The linear workflow uses the initial equilibrium, adds a poloidally and radially localized density source at a fixed position inside the separatrix and then determines the linear stability with respect to the density perturbation.Growth rates for toroidal modes between n = 1 and n = 25 are found.For the most unstable mode the critical pellet mass for which the growth rate first becomes positive marks the threshold for pellet ELM triggering.The results are shown in figure 8.
Figure 8(a) shows the poloidal cross section of the ITER equilibrium as well as the respective pellet injection locations used in the linear simulations.Note that in the linear workflow the pellet neither ablates nor moves; the background plasma profiles also do not evolve.The found thresholds are shown by the blue data points in figure 8(b).Note that for ITER we use pellet mass instead of pellet size to represent thresholds, which scales with the third power of the pellet radius.The outboard midplane is found to be the most effective location for ELM triggering.The orange and red circles represent the planned mass of pacing and fueling pellets, respectively.The thresholds found for the low-field-side are in between those.However, the simulated mass threshold is still more than twice the mass of planned pacing pellets.The current ITER design does not include a pellet launch tube at the outboard midplane though.At the planned injection locations in ITER the predicted mass threshold for ELM triggering largely exceeds the mass of pacing pellets (orange circle) and even the mass of fueling pellets (red circle).The locations are marked by the purple circles and arrows.Near the 220 • location, which is one of the planned injection ports in ITER, the mass threshold is 17 times larger than the pacing pellet mass (orange circle) and about 3 times larger than the maximum fueling pellet mass (red circle).

Nonlinear simulations
In the DIII-D case in section 2 we have seen that the 2D nonlinear workflow confirms the pellet thresholds on the low-field side but finds even larger thresholds on the high-field side.Since the nonlinear workflow requires more computational time than the linear one, we only apply it to selected poloidal angles while sufficiently cover points of interest like the midplane and the planned injection locations.
The green data points in figure 9 show the resulting pellet mass thresholds found by the 2D nonlinear workflow.For comparison the blue data points are found by the linear workflow and are the same as shown in figure 8.The nonlinear workflow again confirms the thresholds found on the low-field-side; some are even a bit lower than the linear ones.Also, again on the high-field-side the thresholds are significantly larger than the linear workflow indicated.The thresholds near the planned injection locations are confirmed to be exceedingly higher than planned pacing or even fueling pellets.On the other hand, fueling the plasma with pellets injected from around the high-field side would not be expected to unintentionally trigger ELMs in an otherwise ELM-free scenario.

Conclusions
In ITER, pellets are calculated to require more than 8 times the mass than currently planned to reliably trigger ELMs from the planned injection locations in its 15 MA, Q = 10, lowcollisionality discharges with peeling-limited pedestals.It is found that under such conditions mass thresholds are significantly larger than typical pacing pellets used in present day devices.The distance of an equilibrium's operational point from the ballooning branch of the pedestal stability boundary has a very strong impact on the threshold for ELM triggering.
The pellet mass threshold strongly depends on the poloidal injection location in DIII-D and ITER.The outboard midplane is the most effective location for triggering ELMs in all studied cases, while the simulated mass threshold is more than 40 times larger in the shown ITER case than in DIII-D for shot 178555 at 3055 ms.The upper and lower x-point regions are the most unfavorable locations in all studied cases.
A new nonlinear workflow is introduced, and its results are compared to linear simulations.The general trend observed in both devices is that linear simulations underestimate thresholds on the high-field side while mostly match thresholds on the low-field side.The nonlinear workflow includes pellet ablation, motion and spread of the pellet cloud, as well as evolution of plasma profiles.Even though simulated thresholds differ between the linear and nonlinear workflows, especially on the high-field side, the qualitative results remain unchanged.
The calculated pellet mass thresholds, comprised of M3D-C1 simulations for DIII-D and ITER as well as DIII-D experimental observation, show that pellet ELM-pacing will be challenging for operation in reactor like conditions.The pellet mass throughput would be incredibly demanding and an engineering challenge, both for injection and pumping, to use for high-frequency ELM pacing.On the other hand, injecting so much mass could change the pedestal, by raising density and lowering temperature, to a more ballooning limited regime, which in turn requires less pellet mass to pace ELMs.This raises the question, if the mass threshold depends on the pacing frequency.So far only individual triggering events have been studied in the regime with peeling limited pedestals.This work does not preclude pellet pacing as an actuator in ITER, but investigates constraints on conditions where it would be a useful tool.
Finally, if a new injection line at the outboard midplane is added to the design, ITER would still be able to trigger ELMs at a still reasonable pellet mass using this line, while using the existing lines, especially near 240 • , to fuel the plasma without triggering ELMs in an otherwise ELM-free scenario.The latter is certainly a very desirable feature.The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes.DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

Acknowledgment
Disclaimer: This report was prepared as an account of work sponsored by an agency of the United States Government.Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights.Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof.The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Appendix. Pellet perturbation in M3D-C1
This section provides additional information on how we set up the density perturbation for injected pellets in the M3D-C1 code.Since only 2D simulations are discussed here, we focus on this setup only.In the simulation the pellet is added as an additional main-ion density source S to the continuity equation, positioned at the injection location (R p , φ p , Z p ), and localized by a 2D Gaussian of width v p and amplitude A p in the R-Z plane referred to as pellet density or pellet cloud density throughout the paper.In this 2D setup the pellet density S is not toroidally localized, which implicitly turns the pellet into an axisymmetric density ring instead of an actually localized pellet.The two workflows discussed throughout the paper, linear and nonlinear, differ on how the pellet density S is initialized and evolved during the M3D-C1 simulation.In the linear case S remains static, so the only code inputs are the location (R p , φ p , Z p ),  width v p and amplitude A p , also called the 'pellet-rate' in literature, of the Gaussian distribution.Using the ITER case, discussed in section 5, as an example, the following figure demonstrates how we calibrate these inputs for spherical deuterium pellets of radius r p .Using a pellet ablation model [13] and its implementation in the XPELLET code [14] a spherical deuterium pellet of radius r p is injected into a fixed background plasma along a predefined trajectory with a fixed velocity.Here we use the ITER equilibrium as shown in figure 8(a), as well as n e and T e profiles as shown in figure A2.The XPELLET code ablates the pellet and returns the electron density deposited in the plasma.The solid lines in figure A1 show this result for three different pellet radii, using a pellet velocity of 100 m s −1 and a horizontal trajectory along the midplane starting outside the separatrix at R = 8.2 m.For each given pellet radius we fit the solid curves in the figure with a Gaussian using equation (A.1) to obtain the four M3D-C1 inputs.The dashed lines in figure A1 show the resulting density perturbation ∆n e = n e − n e0 , with n e the total density and n e0 the equilibrium density profile (see figure A2(top)), in M3D-C1 for each pellet radius.
In the nonlinear case S is initialized in a similar way, with some differences.A nonlinear M3D-C1 simulation evolves the plasma, including the density perturbation over time.The pellet density spreads poloidally, the pellet radius shrinks and the pellet position moves with a given pellet velocity.The M3D-C1 code internally applies the same ablation model as the XPELLET code.Inputs to M3D-C1 are now the pellet position (R p , φ p , Z p ), pellet cloud width v p and the pellet radius r p at t = 0, as well as the velocity vector and the pellet material.The pellet radius and pellet cloud width then determine the amplitude A p self-consistently.Note that most of these parameters now evolve with time; only the velocity is kept fixed at 100 m s −1 .To find the initial settings, we again match the solid lines in figure A1 with the M3D-C1 density perturbation ∆n e , but now within a few Alfven times after the simulation starts, typically around t = 5-10τ A .
Throughout the paper we use pellet radius r p and pellet mass m p as threshold measures.These are equivalent and can be converted by with n D = 0.1967 g cm −3 for frozen deuterium.Furthermore, to calculate the number of atoms in a pellet, divide its mass by the molar mass for deuterium m D = 2 g mol −1 and multiply by Avogadro's constant N A = 6.022 • 10 23 mol −1 .So, a pellet of radius 0.72 mm has a mass of about 0.3 mg and contains about 9.3 • 10 19 deuterium atoms.

Figure 1 .
Figure 1.Number of pellets arriving at the plasma within the first 50% of the natural ELM cycle sorted by their relative mass.Pellets marked in blue do not trigger ELMs in the experiment; the ones marked in red successfully trigger ELMs.The black dashed line shows the experimentally observed mass threshold.

Figure 2 .
Figure 2. Measured spherical radius of injected pellets, normalized by the M3D-C1 simulated threshold, for multiple times during discharge 178555.The black line shows the simulated normalized threshold.Blue markers represent pellets that do not trigger ELMs in the experiment; red markers represent successful trigger events.

Figure 3 .
Figure 3. Growth rate over time for a toroidal mode with mode number n = 7 in DIII-D discharge 178555.The colors represent various pellet sizes given by their respective spherical radius.

Figure 4 .
Figure 4. Time averaged growth rates for multiple toroidal modes and various pellet sizes.Note that the colors do not refer to previous figures.

Figure 5 .
Figure 5. Interpolation of growth rates for given pellet sizes to determine the zero crossing.Top: n = 9.Bottom: n = 23.

Figure 6 .
Figure 6.(a) Poloidal cross-section of DIII-D discharge 178555 at 3055 ms.The points mark the pellet injection locations and the arrows their respective injection direction.(b) Comparison of linear (blue) and nonlinear (green) pellet size thresholds, given by spherical radius in mm, with respect to poloidal angles as shown in (a).

Figure 7 .
Figure 7. Pellet size thresholds for three different VARYPED equilibria of discharge 178555.The inlet at the top left shows the respective operating points of the equilibria in the ELITE pedestal stability diagram.

Figure 8 .
Figure 8.(a) Poloidal cross section of the ITER equilibrium.The red markers show the respective pellet injection locations used in the linear workflow.(b) Pellet mass thresholds for multiple poloidal angles shown in blue.The orange and red circles show the upper limits of planned pacing and fueling pellet masses, respectively.The purple markers and arrows indicate the planned pellet injection locations in ITER.

Figure 9 .
Figure 9. Pellet mass thresholds for multiple poloidal angles.The green data points are obtained by the 2D nonlinear workflow, while the blue ones are from the linear workflow and are the same as shown in figure 8.
This material is based upon work supported by the US Department of Energy, Office of Science, Office of Fusion Energy Sciences, using the DIII-D National Fusion Facility, a DOE Office of Science user facility, under Award(s) DE-AC05-00OR22725, DE-AC02-09CH11466, DE-FC02-04ER54698, and Center for Tokamak Transient Simulation (CTTS) SciDAC (DE-SC0018109).This research used resources of the National Energy Research Scientific Computing Center (NERSC), a US Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231.This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC05-00OR22725 with the US Department of Energy (DOE).

Figure A1 .
Figure A1.Pellet cloud density S for various pellet radii, injected radially inwards from the outboard midplane in ITER.Solid lines are directly from the pellet ablation model, dashed lines are the M3D-C1 density perturbation.

Figure A2 .
Figure A2.Electron density (top) and temperature (bottom) profiles for the ITER case used in section 5.