NBI optimization on SMART and implications for scenario development

The SMall Aspect Ratio Tokamak (SMART) under commissioning at the University of Seville, Spain, aims to explore confinement properties and possible advantages in confinement for compact/spherical tokamaks operating at negative vs. positive triangularity. This work explores the benefits of auxiliary heating through Neutral Beam Injection (NBI) for SMART scenarios beyond the initial Ohmic phase of operations, in support of the device’s mission. Expected values of electron and ion temperature achievable with NBI heating are first predicted for the current flat-top phase, including modeling to optimize the NBI injection geometry to maximize NBI absorption and minimize losses for a given equilibrium. Simulations are then extended for a selected case to cover the current ramp-up phase. Differences with results obtained for the flat-top phase indicate the importance of determining the plasma evolution over time, as well as self-consistently determining the edge plasma parameters for reliable time-dependent simulations. Initial simulation results indicate the advantage of auxiliary NBI heating to achieve nearly double values of pressure and stored energy compared to Ohmic discharges, thus significantly increasing the device’s performance. The scenarios developed in this work will also contribute to diagnostic development and optimization for SMART, as well as providing test cases for initial predictions of macro- and micro-instabilities.


Introduction
Future fusion power plants based on the tokamak concept will arguably operate within a limited range of plasma and 5 Present address: Ecole Polytechnique Fédérale de Lausanne, Swiss Plasma Center, CH-1015 Lausanne, Switzerland. 6Present address: ITER Organization, 13067 St Paul Lez Durance Cedex, France.* Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.magnetic configuration parameters.Those parameters, optimized for safe steady-state operations, need to be reliably predicted during the design phase of the reactor, since they pose engineering constraints for designing the device.Optimization of future reactors is still undergoing design efforts to assess the merits and prospects for performance (see e.g.[1] and references therein).Two considerations defining the optimal configuration are arising as key indicators for the optimization process: the choice of the so-called aspect ratio (ratio of major to minor plasma radius), and the target scenario (or configuration) for steady-state operations.The latter includes negative triangularity regimes [2][3][4], for which plasma shaping is reversed from what is commonly used in most experiments.
Experiments in present tokamaks are expected to provide the basis for identifying the optimum regime(s) in which future reactors will operate.Given that present experimental devices are typically smaller than future reactors, and operate at considerably lower fusion performance, extensive validation of models is required for reliable, quantitative predictions.For example, negative triangularity scenarios have been studied in conventional aspect ratio devices, but the potential advantages for compact Spherical Tokamak (ST) devices operating at aspect ratio ⩽2 has not been extensively explored yet.Transport models, typically developed based on positive triangularity plasmas, also require more extensive validationand possibly development-when transferred to negative triangularity scenarios [5].To that end, the SMall Aspect Ratio Tokamak (SMART [6]) device is being commissioned at the University of Seville (Spain).The main goal of SMART is to test experimentally the potential advantages of negative triangularity scenarios in a ST configuration [6][7][8][9][10], and contribute to the overall understanding of the physics and transport scaling of ST-based fusion devices in terms of plasma confinement [11].Experimental data will also provide valuable information for further development and validation of predictive transport models.
The SMART Program includes three phases.After the initial commissioning of the device with Ohmic plasmas (Phase I), additional heating by Neutral Beam Injection (NBI) and enhanced capabilities for controlling the magnetic configuration are planned for Phase II.Phase III further expands the device capabilities in terms of achievable plasma current, magnetic field, auxiliary power and pulse duration.This work explores the optimization of NBI injection in SMART through the TRANSP and Orbit codes [12][13][14], in view of its commissioning in Phase II.The main goal is to estimate the range of plasma parameters that would be achievable in Phase II, and test a simulation workflow for further scenario development.These are necessary steps in the early phase of SMART, for example to provide estimates of plasma parameters for the design and optimization of diagnostics [15], as well as for an initial assessment of macro-and micro-instability properties for Phase II [16,17].
Results are compared to those from previous optimization studies [10,18] conducted with different codes than those used in this work.This initial exercise allows to verify the consistency of results, especially on NBI optimization, across codes that employ different approaches to solve for NB ion dynamics, as briefly discussed in section 3.
The target parameters for most simulations are conservative choices that do not require achieving the maximum values planned for Phase II, to demonstrate the device capabilities at the beginning of Phase II.Also, scenario predictions focus on the current ramp and flat-top phases.Little effort has been spent to optimize the ramp-down phase, which is typically much more difficult to predict since its evolution heavily depends on the previous plasma trajectory.
After an introduction to the design parameters of SMART in section 2, the main tools used for simulations are briefly summarized in section 3. Optimization of the NBI parameters is discussed in section 4. Section 5 contains the main results of this work.Section 5.1 presents the results for NBI optimization in the current flat-top, based on two target configurations.Section 5.2 expands the analysis to an entire discharge, including the current ramp-up phase.The main results from this work are finally summarized in section 6, which concludes the paper.

The SMART device
The main parameters of the SMART device are summarized in table 1 for the three planned phases of operations.Transition from one phase to the next will require achievement of engineering and physics milestones, as well as upgrades of some of the subsystems such as power supplies.
In terms of physics, Phase I aims at commissioning the machine with a limited set of diagnostics and modest plasma performance.Installation and commissioning of a NBI system for auxiliary heating and developing a feedback control system are the main physics upgrades foreseen for Phase II, which will extend the range of flexibility of achievable scenarios in terms of plasma parameters.Further upgrades of power supplies and of the NBI system are foreseen during Phase II, in parallel with physics studies.Phase III is expected to achieve the full machine capability and scenario flexibility with extended pulse duration and the addition of electron-cyclotron injection as an additional auxiliary heating system.Engineering upgrades mostly include power supplies, plasma-facing components (from bare stainless steel to C tiles), and diagnostics.
A central solenoid is used to drive inductive current [19].The main (toroidal) magnetic field is generated through 12 rectangular toroidal field (TF) coils [9].Initial operation (Phase I) will be limited to on-axis field B 0 ⩽ 0.1T.Upgrades of the TF power supplies in the subsequent phases II and III will enable operations up to B 0 = 1T.
Three pairs of in-vessel poloidal field (PF) coils complement two out-vessel coils to provide the required flexibility in the achievable magnetic configurations [7,8,20].According to initial modeling, a triangularity in the range δ = −0.6 → +0.6 can be achieved.The same PF coils are also used for vertical stability and divertor control.The vertical stability coils are the only ones planned to be controlled by a real-time feedback system, whereas all other coils will be controlled in feed-forward.
The cylindrical vacuum vessel has a radius of 0.8 m and height of 1.6 m (figure 1) [7,19].The vessel and coils geometry are expected to accommodate plasmas with major radius R 0 = 0.45 m and minor radius at the midplane a = 0.25 m.Several ports on the top/bottom of the vessel, as well as on the outboard wall, have been designed to host diagnostics, feedthroughs for the PF coils and access for the NBI system.

Simulation methodology
Previous simulations in preparation for SMART operations used a combination of codes for magnetic configuration  (FIESTA [21]) and transport analysis (ASTRA [22]).The first studies focused on Ohmic L-mode discharges [20].Next, NBIheated scenarios were considered [10], including scenarios for Phase II for assessment of transition from L-mode to H-mode.Separate works also explored start-up scenarios to assess the requirements on the PF coils' currents and power supplies [23], and the optimization of the NBI system starting from Phase II [10,18].This work builds upon previous work by the SMART Team to expand the analysis to include NBI for timedependent scenario development, and to predict the impact of NBI on the achievable plasma parameters within a common framework provided by the TRANSP code [12].As the simulations described here use different numerical tools and boundary conditions with respect to previous work, differences are expected in the results, especially in the case of predictive simulations where the magnetic equilibrium is evolved selfconsistently with transport models and during transient phases, like the plasma current ramp-up and ramp-down.These differences are indeed very valuable to identify uncertainties in the models and which boundary conditions have larger effect on the achieved global plasma parameters.

TRANSP code for scenario development
TRANSP [12] is an equilibrium and transport code for tokamak experiments, currently used for interpretive and predictive transport simulations.In predictive mode, turbulence transport models are used to derive density, temperature and rotation profiles for given diffusion coefficients or fluxes.TRANSP includes a dedicated module, NUBEAM [13], for fast ion physics based on neoclassical physics.This work is limited to fixed-boundary simulations using LCFS information from previous work and the TEQ equilibrium solver [24] in TRANSP to re-compute the full equilibrium.
For the predictive runs discussed in this work, the Multi-Mode Model (MMM [25]) for thermal transport is used as it has been shown to reproduce profiles in Spherical Tokamaks [11,26,27].Comparing the results from multiple thermal transport models will be the subject of future work.MMM contains physics modules to describe thermal transport induced by several classes of instabilities including the effects of sheared flows [25].As detailed in [25], MMM includes physics models for instabilities driven by electron/ion temperature gradients, trapped electrons, kinetic ballooning and peeling effects, collisionless and collision dominated magnetohydrodynamics drive, as well a model for drift resistive inertial ballooning modes.(The investigation of the mechanisms that dominate thermal transport by micro-instabilities in predicted SMART plasmas is beyond the scope of this work, and will be discussed in future publications once initial parameter ranges are assessed.) Of note, results from previous predictive analysis, e.g. reported in [10], were based on a-priori assumptions on transport coefficients inferred from scaling laws and/or other similar devices.This aspect is important to understand differences between the results previously reported and what is discussed in the following sections.Simulations discussed in section 4 use analytic profiles for density, temperature and rotation, based on previous simulations or inferred from similar experiments.This is done to benchmark the NBI calculations.Next, when optimizing the NBI configuration, the density profile is prescribed and the electron and ion temperature profiles are predicted inside ρ = 0.9 (see section 5), where ρ is the normalized toroidal flux used as radial variable.
Finally, simulations that predict both the ramp-up and the flattop phase (section 5.2) evolve the values of density and temperature at the separatrix, based on a two-point Scrape-Off-Layer (SOL) model [28].The two-point SOL model extends previous work [29,30] to calculate the electron/ion temperature and plasma density in the SOL, given the particle and power fluxes from the core through the LCFS calculated by TRANSP and given recycling coefficients.This approach is used for time-dependent predictions in section 5.2.

Orbit code for NBI optimization studies
Orbit [14] is a guiding-center, particle-following code mostly used to simulate the dynamics of high-energy particles (EPs) for a given magnetic field configuration.The code evolves trajectories of a selected number of markers based on the Hamiltonian formulation for particles' orbits in axi-symmetric configurations [31].The study of the effects of perturbation fields, e.g.induced by plasma instabilities, is not considered in this work.
The Orbit code is well-suited for a quick assessment of the EP confinement in a given configuration, which enables scans over a broad range of plasma conditions starting from a single scenario identified e.g. through TRANSP analysis.Markers for Orbit simulations can be extracted from a fast ion distribution obtained through TRANSP/NUBEAM.The distribution can be representative of either confined fast ions resulting from slowing-down and scattering, or of the deposited NB ions including possibly lost particles.
Recent code developments [32] include the possibility to follow particles in the vacuum region between the last-closed flux surface (LCFS) and a limiter.This feature allows to estimate the poloidal distribution of lost energetic particles, which is required to assess the needs for protective structures preventing localized heat loads on the wall or diagnostic windows.(Incidentally, TRANSP/NUBEAM flags a marker as 'lost' if its position corrected for finite Larmor radius (FLR) effects travels farther than one Larmor radius from the LCFS.This implies that losses may be overestimated if the LCFS is located far from the wall, see additional discussion in section 5).
Of note, a guiding-center approach has been often questioned when dealing with fast ion simulations for compact, low-field devices such as STs (see [33] and references therein).Under those conditions, the NB ion gyro-radius can be comparable to-or even greater than-typical scale length of the background fields, thus breaking one of the main assumptions behind the guiding-center formulation.As shown in the following Sections, a comparison between Orbit/NUBEAM and the full-orbit code ASCOT [34,35] gives confidence that sufficiently accurate results in terms of overall trends and global quantities can be achieved for the cases discussed in this work.More extensive discussion on this topic is left for future publications.

Simulations for NBI optimization
Neutral Beam injection, or NBI, is widely used as an additional heating and current drive system on present tokamaks.Its optimization in terms of geometrical and beam parameters is required to minimize losses of NB ions, in order to maximize absorption and resulting plasma performance, see [36][37][38] among many others.
In this section, scenarios from previous work are first used to cross-validate the NBI optimization results reported in [10,18] with results from TRANSP/NUBEAM and Orbit analysis.As discussed later in this section, good agreement is found among different codes about expected NBI absorption and losses.The analysis is then extended to explore more NBI optimization options based on the scenarios discussed in section 5, especially for time-dependent simulations.
The NBI optimization study is based on parameters from a NBI system whose parameters are comparable to those of the system in use on the COMPASS device [39,40].The system will be available on SMART for Phase II.For the initial operations in Phase II of SMART, the expected injection energies are 15-40 keV and currents in the range 8-12 A, resulting in an injected power up to ≈400 kW once transmission and neutralizer losses are accounted for.Depending on the timing of further upgrades of the NBI system, the available power may increase to 1 MW, although a conservative upper value of 0.4 MW is assumed for most of the following analysis.The working gas is hydrogen.Note that in [18] the NBI injection voltage for most cases was limited to V inj = 15 kV.Experimentally, however, there can be interest in maximizing the amount of NB power delivered to the plasma to favor access to high confinement operations (H-mode) and to improve NB current drive.Therefore, an increased V inj = 25 kV has been typically used for this initial study, although it results in increased NB ion losses, as discussed later.A higher injection voltage is also consistent with the more recent study discussed in [10], although in that work it was assumed that the NB system would have been capable to deliver up to 1 MW of power at V inj = 25 kV.(No specific NBI system was targeted in [10], and 1 MW is likely to be achievable only by installing a second NBI line).
The initial target configurations and plasma parameters used in this work from ASTRA simulations [10] are shown in figure 2 for a pair of matched positive vs negative triangularity discharges.Density and electron temperature are comparable for the two cases, although the negative δ case is predicted to achieve slightly higher electron temperature with larger values of the outer gap, where the the outer gap is the distance between the last-closed flux surface (LCFS) on the outboard mid-plane and the vessel wall.
Figure 3 illustrates the geometry for scans of NB injection parameters, namely the tangency radius inside the plasma and the elevation of the NBI source outside the vessel to enable off-axis NBI.TRANSP simulations are first used to identify the optimum range of tangency radii that minimizes first-orbit losses and shine-through, thus maximizing the NB power effectively coupled to the plasma.Figure 4 shows TRANSP/NUBEAM results of NBI losses for a tangency radius scan covering the entire radial range from high-to low-field side, R tan = 19 cm−64 cm.The NBI is injected horizontally at the midplane.The injection voltage is set to The power lost through first-orbit losses remains negligible for R tan ⩾ 40 cm, which roughly corresponds to the magnetic axis (figure 4(a)).For more perpendicular injection, 19 cm ⩽ R tan ⩽ 40 cm, first-orbit losses increase rapidly and up to 20% of the injected power is lost.Shine-through losses have a different dependence on R tan .Except for the smaller R tan = 19 cm, the lost power remains ≈5%-10% up to R tan ≈ 50 cm, then rapidly increases to >20% as the path of NB injection is reduced and plasma density and temperature decrease towards the edge.For the configurations shown in figure 2, shine-through losses are generally higher for the negative    triangularity case, most likely due to the smaller crosssection, which results in a slightly reduced NB path inside the plasma.
Prompt losses calculated by NUBEAM are compared to those from Orbit in figure 4(b).The comparison is shown in terms of fraction of markers that are lost, rather than power as in figure 4(a), as that information is more readily available from Orbit.
The two codes use different conditions to classify markers as 'lost'.NUBEAM considers a marker lost when it either gets closer than one Larmor radius to a wall/limiter structure, or when its guiding center position is farther out than one Larmor radius from the LCFS.Therefore, some particles that exit the main plasma due to drifts but would otherwise reenter the plasma may contribute to enhance first orbit losses.By default, Orbit marks particles as lost once their guiding center crosses the LCFS.Similar to NUBEAM, this could lead to over-estimated losses.An option to follow particles in the vacuum region between the LCFS and the wall has been recently implemented [32].In this case, particles are lost once their FLR-corrected position hits the wall.Results of the NUBEAM/Orbit comparison shown in figure 4(b) confirm the expectations.First orbit losses to the LCFS from Orbit exceed those from NUBEAM by 10%-15%.When particles are allowed to re-enter the plasma in Orbit, losses decrease by ≈10% with respect to results from NUBEAM.In general, however, the two codes agree on the main trend vs. R tan .(Note that losses due to charge-exchange with edge neutrals, included in NUBEAM, are not implemented in Orbit.According to TRANSP, their contribution is ⩽3% for the two cases analyzed here).
Another possible mechanism for NB ions loss are ripple perturbations to the static magnetic field caused by the finite number of TF coils.Based on the geometry and number of the coils, the TF ripple field amplitude is up to 1% at the plasma edge on the low-field side, and rapidly decreases to ≪1% inside the plasma, cf figures 7 and 8 in [9].Orbit simulations are performed to assess the NB ion losses associated with TF ripple for a worst-case scenario with negative triangularity, in which the plasma region on the low-field side is maximized.As shown in figure 5, losses remain limited to ≈0.5% for the nominal value of the TF ripple field.Only by assuming a ripple field 4-10 times larger than the nominal value we find that the losses become comparable to those caused by other mechanisms, such as charge-exchange and shine-through.Rippleinduced losses are therefore neglected in the following analysis.
The results shown in figure 4 are consistent with previous work, which indicated an optimum tangency radius R tan = 44 cm [18].The consistency of the results across several codes (ASCOT, Orbit and TRANSP/NUBEAM) gives confidence in their validity, as well as in the validity of using TRANSP to extend the study to time-dependent simulations.
As discussed previously, in [18] the NBI injection voltage for most cases was limited to V inj = 15 kV as opposed to the V inj = 25 kV used in figure 4. Figure 6 shows the variation of NB ion losses for R tan = 44 cm and positive triangularity as the injection voltage is varied between 20 kV and 40 kV, overlapping with the V inj values explored in [10].Shine-through losses steadily increase from 3% up to 10% for the maximum injection voltage, while first-orbit losses to the wall remain negligible.For the maximum available NB current of 12 A, the injected power increases from 204 kW to 408 kW for V inj = 20 → 40 kV.Once losses are taken into account, that translates into an increase in coupled NBI power from ≈200 kW to ≈370 kW, which provides significant flexibility in terms of additional heating capability.
Additional flexibility in the profiles achievable by NBI heating in Phases II and III on SMART can be achieved by considering a more elaborate implementation of the NBI system, in which the NBI injector can be tilted vertically.Such flexibility clearly involves complexity in terms of engineering, but it also enables opportunities in terms of physics research and scenario control as demonstrated on other facilities such as DIII-D [36].
Examples of deposited NB ions and their profile after a 1 ms orbiting time are shown in figure 7. It is assumed that the NBI source (about 2 m away from the vacuum vessel) can  be displaced vertically in the range Z tilt = ±50 cm keeping the same cross-over point at the vacuum vessel.The nominal R tan = 44 cm for midplane injection is assumed for all cases.A sample of NB deposited ions is first obtained through NUBEAM, then Orbit is used to evolve particles' trajectories.First-orbit losses to the wall are recorded as discussed above.
From the distribution of markers shown in figure 7, it is clear that the possibility of tilting the NBI source would provide a great advantage in terms of achievable NB ion profiles.That is further demonstrated in figure 8, which shows the NB ion density and resulting NB driven current profiles for Z tilt = −50 → +50 cm.In fact, the possibility of varying the radial gradient of the NB ion profile would provide an excellent tool for controlled studies of the role of NB ions in driving deleterious instabilities such as Alfvénic modes.
In addition to modifying the radial NB ion density profile, a tiltable NBI system would also provide good control over the NB ion distribution in velocity space.Figure 9 shows examples of the NB ion distribution vs energy and pitch (ratio of parallel to total velocity) for Z tilt = −50, 0, +50 cm and V inj = 40 kV.As the tilt is increased from more negative to more positive, thus leading to a better alignment of the NBI injection path with the local magnetic field, the fraction of NB ions shifts from more trapped particles (pitch near zero) to more co-passing particles (pitch near 1).This fine control of the EP distribution resulting from NBI would also provide great flexibility for fast ion studies on a device such as SMART, enabling studies that are not easily achievable/controllable on larger fusion devices.

NBI-heated scenarios for Phase II of SMART
This section discusses initial results on scenario development for NBI-heated plasmas in Phase II of SMART operations.Based on the results on NBI optimization from section 4, simulations are first limited to predicting expected values of electron and ion temperature for a prescribed plasma rotation in the current flat-top phase.To reduce uncertainties in the predictive runs, electron density and plasma rotation are prescribed, as their predictions are presently unreliable when multiple transport channels are simultaneously activated in TRANSP predictions [41].Density and rotation values are based on available results from other compact (or spherical) tokamaks with parameters similar to SMART, such as LTX [42,43] and Globus-M [44].
Once reliable profiles are obtained for the current flat-top phase, predictions of the discharge are then extended backwards in time for an initial assessment of the plasma evolution during the current ramp, cf section 5.2.

Profile prediction in the current flat-top
As a starting point, the same magnetic configuration shown in figure 2 for the positive triangularity case is selected.The plasma is assumed to remain in low-confinement mode (L-mode) for the entire simulation, thus the achieved performance can be interpreted as a lower limit.The initial guesses of plasma profiles for the transport solver in TRANSP are shown in figure 10.As mentioned above, density and rotation profiles are prescribed.Based on ASTRA [22] simulations, cases  with central density n e = 2 → 4 × 10 19 m −3 at the beginning of the predictions are tested, with density increasing in time to mimic additional fueling by NBI (not included in ASTRA runs) and outgassing from the wall.Rotation is varied from 0 up to v t,0 ≈ 70 km s −1 on axis, corresponding to rotation frequencies up to ≈20 kHz commonly observed in other devices.Boundary conditions are imposed for toroidal flux ρ ⩾ 0.9 and profiles are predicted inside ρ = 0.9.Different boundary conditions for n e , T e and T i are tested to verify the robustness of the predictions.Simulations for NBI-heated discharges show no significant changes in the core values as edge T e , T i are varied within 30-100 eV.Ohmic simulations without NBI are also included as reference to assess the impact of NBI on the achievable range of electron and ion temperature.
The duration of the predictive simulations is 150 ms.That corresponds to several NB ion slowing down times and energy confinement times (≈10 ms), resulting in relaxed fast ion distribution and fast ion profiles.It is also sufficiently long to achieve a nearly stationary safety factor profile (q-profile), as discussed below.
The main results for predictions of T e and T i are shown in figures 11 and 12. Simulations of Ohmic plasmas result in peak electron and ion temperatures of 500 eV and 400 eV, respectively.Higher temperatures are predicted for lower densities.The flat region in temperature observed for 0.8 ⩽ ρ ⩽ 0.9 in figure 11 suggests that edge parameters should actually be lower in the Ohmic case.Examples of more self-consistent predictions through a reduced scrape-off layer model [28] are discussed in section 5.2.The addition of NBI and finite plasma rotation brings T e up to 700 eV on axis.The increase of T i is even more significant, with T i = 750 − 850 eV on axis compared to the 400 eV in Ohmic discharges.
Varying the electron density has a limited effect on the predicted profiles of T e and T i , as shown in figure 11.The main effect of increasing n e is to favor NBI deposition slightly more towards the low-field side, but the effective broadening of the temperature profiles remains limited.
The results about ion temperature show a dependence on toroidal rotation, cf figure 12(c), with higher temperatures achieved at the maximum rotation tested in these simulations.This is arguably caused by a reduction in transport associated with increased rotation shear [11], although more detailed analysis is required to identify specific instabilities responsible for thermal transport in these scenarios.
One important note on the above results is on the assumptions made to achieve nearly stationary conditions to infer 'converged' temperature profiles.For the cases shown above, the plasma boundary was kept fixed, although the safety factor profile is allowed to evolve in time.In reality, the discharge will evolve from a small plasma at start-up to the target configuration in the current flat-top, cf section 5.2.When the safety factor profile is evolved in TRANSP, central values q 0 < 1 are achieved after 40-50 ms.That implies that global n = 1 modes such as kinks or sawteeth would likely be destabilized, which would limit the achievable values of pressure and temperature.(n: toroidal mode number).Those effects are neglected in the analysis discussed above and will be the subject of future studies.An initial analysis, though, indicates that n = 1 kinks and n > 1 infernal modes are likely to be (linearly) unstable over a broad range of plasma parameters and may clamp the achievable values of temperature below those found here.Further considerations on the effects of large MHD events such as sawteeth are briefly discussed in section 5.2.3.

Examples of discharge evolution from ramp-up to Ip flat-top
Information on achievable parameters for phase II from the 'time-slice' approach used in section 5.1 can be used as starting point for further scenario exploration.In this section, examples of whole discharge simulations are discussed.A main difference for these simulations is the inclusion of a selfconsistent model to compute boundary conditions for T e and T i predictions [28], which leads to adjustments in the edge temperature based on the particle and power flux through the separatrix.As shown in the following Sections, resolving the thermal plasma and fast ion dynamic is key for further discharge optimization exercises.

Scenario evolution.
As shown in figure 13, it is assumed that the plasma configuration shifts from inner-wall limited shortly after the breakdown phase (not included here) to a diverted configuration during and after the current rampup.The magnetic configurations for the different phases, based on previous work from [20], are used as input for predictive TRANSP simulations.
The plasma is evolved from a inner-wall limited configuration, at I p = 40 kA, and the plasma current is ramped up to its flat-top value of I p = 200 kA within 50 ms.The plasma boundary evolves from an inner-wall limited configuration to a diverted configuration just before the I p flat-top phase.During that time, the plasma density on-axis is set to increase to n e = 3 × 10 19 m −3 , with an additional increase during the I p flat-top phase to account for NBI fueling and outgassing from the vacuum chamber walls.Consistent with simulations in section 5.1, plasma density and rotation profiles are prescribed, with zero rotation assumed for the Ohmic case in the absence of external torque from NBI.However, contrary to those simulations, the value of the density at the separatrix is calculated from a two-point SOL model and used as a boundary condition for the calculation of the analytic density profile.The nominal NBI parameters R tan = 44 cm, V inj = 25 kV and P inj = 250 kW, expected for the beginning of phase II, are used unless mentioned otherwise.Simulations at the higher NBI power P nb = 1 MW achievable after upgrades of the NBI system are not addressed here.
Figure 14 shows the calculated profiles of T e and T i for the Ohmic case and for two scenarios with NB injection.The difference between the two NBI cases is the time at which the beams are turned on, namely at the beginning of the current flat-top (t = 50 ms) or earlier during the current ramp (t = 30 ms).The Ohmic scenario achieves a peak electron temperature of ≈ 450 eV by t = 80 ms, then it saturates for the remaining part of the pulse.The ion temperature remains lower, with peak values of ≈250 eV.In general, the predicted profiles are more peaked and with lower edge values than those initially imported from ASTRA, cf figure 2. Although the profiles predicted by TRANSP appear closer to the expectation of L-modelike profiles, the discrepancy suggests that additional validation against real data, perhaps including other STs, is still required.
The addition of NB injection results in an increase in both ion and electron temperatures, as also shown in figure 15 for the on-axis and edge values.As expected, T e and T i rapidly increase on a time scale comparable to the NB ion slowing down time.Peak values with NBI are T e ≈ 650 eV and T i ≈ 1 keV, representing a nearly five-fold increase for T i compared to the Ohmic case by the end of the NBI pulse (t = 130 ms) and a ≈40% increase for T e .Contrary to T e , the core ion temperature does not saturate during the flat-top.This is caused by a decrease in the predicted ion thermal conductivity from the MMM model as the q-profile steadily drops, as confirmed by simulations in which the q-profile is maintained constant after q 0 drops to 1.
The simulations also show an increase in temperature after the beginning of the current ramp-down, t > 150 ms.This is likely due to unrealistic assumptions in how fast the current is decreasing compared to the rapid decrease in the plasma volume, which is simply assumed to shrink to the same innerwall limited boundary shown in figure 13.As the plasma shrinks, the still relatively high current results in large values of the Ohmic current density near the axis, hence the resulting Ohmic heating of ions and electrons.(Virtually no NB ions are present after t ≈ 150 − 160 ms).Further optimization of the ramp-down will require free-boundary TRANSP simulations, which are currently being implemented in the TRANSP equilibrium solvers based on realistic coils' geometry and constraints on the power supplies.
Figure 16 shows some performance indicators for the cases shown in figure 12.The plasma β (volume-averaged ratio of plasma to magnetic pressure) increases from β th = 1.8% to β th = 3.3% − 3.6% for the thermal components after NB injection.NB ions add an additional β fast ≈ 0.9%, bringing the total β close to 4.4%.Similarly, the plasma stored energy doubles from W MHD = 2 kJ to W MHD = 4 kJ, with an additional 1 kJ stored in the NB ion population.
An important result from the simulations with selfconsistent boundary conditions is the difference from the 'time-slice' simulations in which edge values were prescribed based on results from other devices, cf panels (c) and (d) in figure 15.The computed values of edge temperature through the two-point model [28] are considerably lower than the ≈100 eV initially assumed, e.g. in figure 12. Differences are larger for T i than for T e .The different boundary conditions result in predicted changes over the entire  radial profile, with self-consistent profiles being-generally speaking-more peaked and with lower values approaching the LCFS.
Similar considerations hold when comparing the results in figures 12-15 with previous results from predictive modeling reported in [10], in which thermal transport coefficients were prescribed based on a gyro-Bohm model [45] rather than computed self-consistently through anomalous thermal transport models [25], as in the approach adopted for this work.In general, values of electron and ion thermal conductivity computed through TRANSP/MMM are lower than those computed via the gyro-Bohm model, which explains the higher temperatures obtained in this work.The large differences between predicted Ohmic and NBIheated scenarios are promising in terms of flexibility for scenario control that can be achieved by installing an NBI system for phases II and III of SMART operations.A further increase in performance is expected for the planned 1 MW NBI system (cf table 1), which implies four times higher power than in the conservative predictions discussed in this work.The next section highlights additional features mostly relevant for NBI physics.

NBI performance.
Given the differences in plasma profiles as predicted by a 'time-slice' approach and those  Time evolution of the on-axis q-profile value, q 0 for Ohmic (black) and NBI heated scenarios with NBI starting at t = 30 ms (blue) and t = 50 ms (red).
obtained by evolving edge parameters self-consistently, cf section 5.2.1, this section revisits some of the previous results from section 4 related to NBI physics and performance.
The first question of interest is whether NB ions are expected to be sufficiently well confined for the predicted thermal plasma parameters, including during the current ramp-up phase, and whether NB ion losses might endanger some of the in-vessel components.Figure 17(a) shows the predicted NB ion losses as a function of time for the scenarios discussed in figures 14 and 15.Consistent with the 'time-slice' analysis, and in spite of the different thermal plasma profiles, losses remain below ⩽8% of the injected power during the I p flattop phase.Shine-through is the main loss channel, whereas first-orbit losses are negligible.Losses increase if the NBI is injected earlier during the current ramp, as expected for a smaller cross-section plasma at reduced current.Overall, the lost power remains ⩽ 15%, or ⩽40 kW during the short ramp-up transient phase.The fact that shine-through accounts for almost all the losses suggests that a properly cooled NB armor intercepting the escaping neutrals would be enough to protect the device from any significant damage.Scenarios for phase III, which are expected to achieve higher density and temperature than in phase II, have not been explored in this work.They are likely to be less prone to negative effects of shine-through losses, given the higher NB deposition rate, except perhaps if the NBI is fired very early during the current ramp.
The coupled NBI power is transferred to both electrons and ions through slowing down and thermalization.As shown in figure 17(b), for V inj = 25 kV about 60% of the power goes to electrons and 40% to ions.The relative fraction can be varied by varying the NB injection voltage, resulting in the ability to heat preferentially electrons at higher V inj or ions at lower V inj .
Another issue of interest for the NBI vs. Ohmic scenarios, as well as for potential upgrades of the facility, is how much flexibility NBI can provide for scenario control.Results shown in figure 18 indicate that-at least for the parameters predicted for phase II-NBI will have little impact on the overall evolution of the current profile.The total current driven by NBI, I NB ⩽ 17 kA, is at most 7.5% of the total current for I p = 200 kA scenarios, see figure 19.Given the NB injection geometry with R tan = 44 cm, the current density profile is highly peaked on-axis.As a comparison, the bootstrap current profile is much broader and accounts for at least twice the NB-driven current, sustaining more off-axis current than NB current drive.Nevertheless, NBI is still projected to be a major tool to provide flexible profile control for thermal plasma parameters in phases II and III.

Effect of large MHD events such as sawteeth.
As mentioned at the end of section 5.1, the results discussed above neglect the possible effect of plasma instabilities on the predicted performance for Phase II.As such, they should be considered as upper limits e.g. for the maximum electron and ion temperature achievable for given Ohmic and auxiliary power.Although information on the dominant instabilities that may be present in SMART plasmas is not yet available, a qualitative assessment of the effects of macroscopic instabilities can be performed with TRANSP by including a model for sawteeth.The latter would cause a redistribution of thermal plasma and NB ions in the region inside the q = 1 surface whenever the minimum of the q-profile drops below q = 1.
Two scenarios are tested in TRANSP, as shown in figure 20.In the first case, it is assumed that sawtooth crashes are triggered every 10 ms, which is comparable to the fast ion slowing down time over which the NB ion profile (and associated pressure) is restored after the previous crash.The q-profile on axis drops to 0.7 before the crash before being restored to q = 1 inside the q = 1 surface around normalized toroidal flux of 0.4, see figures 20(a) and (b).The electron and ion temperature drops from 650-700 eV and 900 eV to ≈550 eV and ≈750 eV on axis (figures 20(c) and (d)), respectively.T e and T i profiles flatten inside q = 1, as expected from the sawtooth model implemented in TRANSP assuming full reconnection after a crash.
In the second case, sawteeth are more frequent and the minimum of the q-profile is maintained around 0.85-0.9before a new crash occurs.The spacing between sawteeth is such that the electron temperature profile recovers in between crashes, although the ion temperature remains lower by 100-150 eV compared to the cases of no or infrequent sawteeth, see figures 20(e) and (f).
As a qualitative conclusion, the presence of large MHD events such as sawteeth would result in an overall decrease of up to 20%-25% of the central thermal plasma parameters in Phase II, compared to the case with no instabilities.Future work based on scenarios discussed in section 5.2 will target a more quantitative assessment of the possible macro-and micro-instabilities that could develop in Phase II of SMART operations [16,17].Example of the effects of large MHD events such as sawteeth on the predicted Te, T i profiles.(a) Evolution of the on-axis q-profile withouth sawteeth (black), with infrequent sawteeth (blue) and frequent sawteeth (red).(b) Safety factor radial profile for the case without sawteeth (black) and for the two other cases.Solid (dashed) lines refer to profiles just before (after) a crash.(c) and (d) Evolution of the on-axis Te and T i values withouth sawteeth (black), with infrequent sawteeth (blue) and frequent sawteeth (red).(e) and (f) Radial profiles of Te and T i for the case without sawteeth (black) and for the two other cases.Solid (dashed) lines refer to profiles just before (after) a crash.Profiles are taken from the time range 100 ms ⩽ t ⩽ 120 ms, corresponding to the current flat-top.

Conclusions
The SMART tokamak under commissioning at the University of Seville (Spain) will explore confinement properties of negative triangularity configurations in a compact (spherical) tokamak.After an initial phase limited to Ohmic heating, installation of a NBI system for auxiliary heating is planned to extend the achievable plasma regimes and performance.
Numerical optimization of NBI-heated scenarios indicates an optimal value for the NBI tangency radius, R tan = 44 cm, that minimizes first-orbit losses and shine-through for the plasma parameters expected in Phase II.Simulations that assume a vertically steerable NB injection show the potential advantage of such flexible system in terms of scenario and fast ion distribution control, which would compensate the increased complexity (and cost) from the engineering point of view.Predictions based on the increased flexibility in NBI will inform possible upgrades of the NBI system from Phase II to Phase III of SMART operations.
Based on initial simulations for the current flat-top phase with prescribed density and rotation, electron and ion temperatures are predicted to increase by 35%-50% for the nominal NBI parameters of R tan = 44 cm and V inj ≈ 25 kV compared to purely Ohmic discharges.This confirms the added value of a NBI system capable to deliver 0.25-1 MW power to SMART plasmas to extend the range of achievable parameters for thermal and fast ion transport studies.
The extension of predictive simulations for T e and T i to the entire discharge, including the current ramp-up and a self-consistent estimate of pedestal parameters, provides more insight on the predicted parameters.It is clear that adding selfconsistent boundary conditions through a reduced scrape-off layer model [28] largely modifies the simulations' outcome, for instance by suggesting more peaked profiles with lower edge temperature than what initially predicted in [7].Future work will focus on this issue, as well as on the comparison across different thermal transport models to assess the reliability of present numerical tools for scenario predictions.

Figure 1 .
Figure 1.The SMART device under commissioning at University of Seville, Spain.Reproduced with permission from [19].CC BY-NC-ND 4.0.

Figure 2 .
Figure 2. Plasma scenarios for NBI tangency radius scan.(a) Plasma boundary in a poloidal cross-section for positive (blue) and negative (red) triangularity, in addition to the wall boundaries and in-vessel PF coils shown as thick black lines (PF coils outside the vessel on the low-field side are not shown).Electron density (b) and temperature (c) are shown for the target positive and negative triangularity scenarios from panel (a).Solid (dashed) lines in panel (c) refer to electron (ion) temperature.The abscissa in (b) and (c) is the toroidal flux normalized to its value at the LCFS.

Figure 3 .
Figure 3. (a) Schematic of SMART at the midplane as seen from above, with lines repesenting the vessel wall at R = 0.8 m (black), the magnetic axis at R = 0.45 m (dashed), and three radii corresponding to Rtan = 20, 44, 64 cm.Solid lines show the NBI path corresponding to the three tangency radii.(b) Vertical cross section showing the NBI path outside the vessel (solid lines) and the corresponding deposition locations at guiding center for three elevations of the NBI source, Z tilt = −50, 0, +50 cm.

Figure 4 .
Figure 4. NBI losses vs. tangency radius for midplane NB injection.(a) Fraction of injected NBI power lost through first-orbit losses to the wall (solid) and through shine-through (dashed) calculated by TRANSP/NUBEAM.(b) Fraction of lost particles calculated by TRANSP/NUBEAM (black, green) and through Orbit for the positive triangularity case.As in panel (a), the dashed line correspond to shine-through losses from TRANSP.

Figure 5 .
Figure 5. (a) Contour of the nominal TF ripple field amplitude in a poloidal cross-section for a negative triangularity configuration.(b) Particle lost fraction as a function of the TF ripple amplitude, rescaled with respect to the nominal amplitude.

Figure 6 .
Figure 6.Fraction of injected NBI power lost through first-orbit losses to the wall (solid) and through shine-through (dashed) as a function of the NBI injection voltage.Losses to the wall are multiplied by 100 for clarity.First-orbit losses are negligible for V inj < 40 kV, hence the results are noisy due to poor statistics of the lost markers.

Figure 7 .
Figure 7. Examples of NBI deposited markers (green) and their relaxed profile after 1 ms orbiting (black).The NBI is tilted from Z tilt = −50 cm to Z tilt = +50 cm at fixed nominal tangency radius Rtan = 44 cm.The yellow markers represent deposited NB ions lost to the wall.

Figure 8 .
Figure 8.(a) NB ions density and (b) current density as the NBI is tilted from Z tilt = −50 cm to Z tilt = +50 cm for fixed nominal tangency radius Rtan = 44 cm.

Figure 10 .
Figure 10.(a) Density profiles at the beginning of the simulation, assuming the plasma current has already achieved its flat-top value Ip = 200 kA.(b) Initial guesses for Te (solid) and T i (dashed) profiles used in TRANSP predictions during the current flat-top phase.

Figure 11 .
Figure 11.(a) Prescribed density profiles near the end of the current flat-top.(b), (c) Results from Te and T i predictions assuming a rotation on axis ωtor = 120 krad s −1 .The lower, dashed curves refer to Ohmic simulations without NBI injection for central ne = 3 − 4.5 × 10 19 m −3 .Upper curves refer to runs with V inj = 25 kV and an injected NBI power of 250 kW.

Figure 12 .
Figure 12.(a) Prescribed toroidal angular velocity.(b), (c) Results from Te and T i predictions as the assumption on toroidal rotation is varied, cf panel (a).The lower, dashed curve refers to an Ohmic simulation without NBI injection.Upper curves refer to runs with an injected power of 250 kW.

Figure 13 .
Figure 13.(a) Poloidal cross-section of SMART with the wall and in-vessel PF/divertor coils.The two coils outside the vessel on the low-field side are not shown.The three contours show the LCFS at t = 0 (blue), t = 50 ms (cyan) and t = 150 ms (red).(b) Waveform of plasma current, Ip.(c) Waveforms of central (black) and edge (red) density.

Figure 14 .
Figure 14.Predicted profiles of Te (top row) and T i (bottom row).Panels (a), (d) refer to Ohmic simulations.Panels (b), (e) show results with NBI starting at t = 50 ms, at the beginning of the current flat-top.Panels (c), (f) show results for NBI starting at t = 30 ms, near the middle of the current ramp.Note the different color scale for T i between panels (d)-(f).

Figure 15 .
Figure 15.(a) Central (solid) and edge (dashed) electron temperature for Ohmic (black) and NBI scenarios.NBI is turned on at t = 50 ms (blue) and t = 30 ms (red).(b) Same as in panel (a) for the ion temperature.(c), (d) Radial profiles of Te and T i at t = 125 ms, near the end of the current flat-top.The green lines show the profiles obtained with fixed boundary conditions for Ohmic (dashed) and NBI plasmas (solid), cf figure 12.

Figure 16 .
Figure 16.Volume-averaged ratio of plasma to magnetic pressure, β, for (a) thermal plasma and (b) fast particles from NBI for the scenarios shown in figure 12. (c) Energy stored in the thermal plasma (solid) and in fast particles (dashed).

Figure 18 .
Figure18.Time evolution of the on-axis q-profile value, q 0 for Ohmic (black) and NBI heated scenarios with NBI starting at t = 30 ms (blue) and t = 50 ms (red).

Figure 19 .
Figure 19.Calculated total current from NBI (a) and bootstrap (b).The total current is Ip = 200 kA during the flat-top.(c) Radial profile of the NB-driven current density at t = 125 ms, showing a on-axis peaked profile for the nominal NB injection parameters.(d) Radial profile of bootstrap current density for the NBI-heated scenarios (blue, red) compared to the Ohmic case (black).

Figure 20 .
Figure 20.Example of the effects of large MHD events such as sawteeth on the predicted Te, T i profiles.(a) Evolution of the on-axis q-profile withouth sawteeth (black), with infrequent sawteeth (blue) and frequent sawteeth (red).(b) Safety factor radial profile for the case without sawteeth (black) and for the two other cases.Solid (dashed) lines refer to profiles just before (after) a crash.(c) and (d) Evolution of the on-axis Te and T i values withouth sawteeth (black), with infrequent sawteeth (blue) and frequent sawteeth (red).(e) and (f) Radial profiles of Te and T i for the case without sawteeth (black) and for the two other cases.Solid (dashed) lines refer to profiles just before (after) a crash.Profiles are taken from the time range 100 ms ⩽ t ⩽ 120 ms, corresponding to the current flat-top.

Table 1 .
SMART parameters for the three phases.