Bifurcation-like transition of divertor conditions induced by X-point radiation in KSTAR L-mode plasmas

Density ramps with ion grad B drift directed into lower single null KSTAR L-mode plasmas are associated with a simultaneous and abrupt reduction of the divertor particle flux on both low- and high-field-side targets when the mid-plane line averaged electron density reaches a given level. Target embedded Langmuir probe signals show a clear ‘cliff edge’ behavior similar to that observed in the divertor target electron temperature in DIII-D H-mode plasmas (Eldon et al 2017 Nucl. Fusion 57 066039; McLean et al 2015 J. Nucl. Mater. 463 533–6). The collapse of the particle flux is observed along the whole divertor target area (from private flux region to the far scrape-off layer (SOL)). The critical upstream density of this target flux cliff is invariant under fuel gas throughput modulation. The transition along the cliff occurs in tens of milliseconds. With the cliff, carbon impurities and deuterium neutrals transported through the X-point to the core produce a strong radiation spot near the X-point, seen on bolometric signals, and increase the upstream density. The experimental observations are consistent with time-dependent SOLPS-ITER simulations, which also demonstrate an abrupt transition of the target flux and upstream density with the increase in X-point radiation. The timescale of the cliff predicted by SOLPS-ITER is consistent with the experiment, although, it is influenced by gas throughput or time-dependent numerical methods. In the L-mode phase space of separatrix electron density and temperature, branches are divided based on target temperature, because the latter is strongly coupled to the radiation front and ionization front due to the monotonic characteristic of the parallel electron temperature distribution. Since the H-mode condition operates at a much higher upstream density and electron temperature in phase space, dissipation from sputtered carbon alone leads to the density limit before reaching the X-point radiation condition. This is therefore consistent with the fact that cliffs have never been observed in H-mode KSTAR experiments.


Introduction
Detached divertor operation will be mandatory in ITER and all future reactor devices based on solid metal, actively cooled plasma-facing components (PFCs) [1]. Fundamentally, the detached condition is obtained by decreasing the electron temperature, T e along scrape-off layer (SOL) flux tubes in the divertor region so as increase momentum transfer and eventually achieve volumetric recombination in critical areas where the parallel heat flux is highest (in the separatrix vicinity) [2,3]. Divertor conditions in tokamaks are often characterized by the target plasma parameters such as target particle fluxes Γ t or T et , usually measured using embedded Langmuir probes (LP). Whilst the target T e measured in this way can often be prone to error in detached, or even high recycling conditions (see e.g. [4,5]), the probe ion saturation current, a direct measure of Γ t , is an unambiguous sign of detachment. It will be the primary experimental measure used throughout this paper.
Abrupt, or bifurcation-like transitions of the divertor plasma state or target parameters under conditions close to detachment have been observed in several devices. The divertor oscillations observed in the JET Mark I divertor configuration (carbon walls) were attributed to the coupling/decoupling of sputtered carbon [6]. The electron temperature 'cliff' in DIII-D [7,8], observed in single null, forward field Hmodes, again with carbon PFCs, was explained by the role of the E × B drift [9]. Such transitions have also been observed in machines with metal PFCs. An abrupt drop of the target flux and increase of the time derivative of the line averaged electron density was observed with the onset of X-point radiation (XPR) in C-Mod [10]. Similar behavior has been found in density ramp experiments with a horizontal outer target (OT) in JET with the ITER-Like Wall [11]. Generally, however, the transition to detachment is often a smoother process, without abrupt change of divertor parameters (see e.g. JET [12], ASDEX Upgrade [13] and TCV [14]). The strong XPRs seen in ASDEX Upgrade [1] are not associated with an abrupt transition of the divertor conditions, but fluctuations are observed when the outer target is detached [13]. The current convective instability [15], which is enhanced when there is asymmetry between the divertor target plasma temperatures, can explain these radiative fluctuations.
As described in this paper, bifurcation-like transitions of the divertor condition have been observed in KSTAR L-mode density ramp experiments with forward toroidal magnetic field, B t . Once the line averaged electron density,n e reaches a critical level, LP saturation currents across the entire (from private flux region (PFR) to the SOL) inner and outer targets are observed to decrease abruptly with the onset of the XPR.
In this paper, we primarily focus on understanding the cause of abrupt particle detachment in L-mode experiments using numerical plasma boundary simulations with the SOLPS-ITER code [16,17] in a time-dependent manner. The code reproduces the trend from the L-mode experiment and suggests a possible cause for the bifurcation-like divertor condition transition. Both experiment and simulation indicate that the major contributing factor appears to be the movement of the detachment front. Although our main focus is on the Lmode experiments and their modeling, we also provide some speculations on what might happen if the work could be extended to H-mode in future experiments. However, it is important to note that the H-mode discussion is not based on experimental results, but rather on theoretical considerations and possible extensions of the current study.

Experimental setup
The KSTAR L-mode density ramp experiment discussed here was conducted during the 2018 campaign. The deuterium molecule (D 2 ) gas puffing rate Γ D2 , is directly controlled either by density feedback (based onn e measured by millimeter-wave interferometer (mmWI) [18], beam path shown in figure 1) or feedforward during the experiments. Two key discharges (pulse numbers #20265, #20267) will be the focus of the analysis here and are characterized by total heating power 0.93 MW (mostly from neutral beam injection), flat-top plasma current I p = 0.6 MA, B T = 2.5 T, both in the forward direction and thus with the ion B × ∇B direction downwards into the lower divertor. The magnetic equilibrium shape is lower single null, with the inner and outer strike points positioned on the inboard (vertical) and central (inclined) divertor target, respectively.
In case of #20265,n e increased linearly from 2.0 − 3.0 × 10 19 m −3 with feedback control of Γ D2 during the I p flattop. For #20267,n e was ramped with feedback control of Γ D2 during the interval 0-5 s, followed by a fixed gas puff rate at Γ D2 = 1.0 × 10 21 s −1 without feedback control for 5-10 s. Embedded LP arrays [19] (see figure 1 for the location of LP in poloidal cross-section) are mainly used to determine divertor condition based on the measured saturation current converted to target parallel ion flux Γ ∥t . In the 2018 campaign, LP #41-44 were available on the inboard divertor, and LP #51, 53, 54 and LP #106-111, 113, 114 were available on the central divertor (outer or low-field-side target). The LP #51-54 and LP #106-114 are located in different toroidal sectors of the central divertor. A two dimensional radiation profile was obtained from the tangentially viewing infrared imaging video bolometer (IRVB) [20].  (figure 2),n e was feedback controlled to yield a linear density ramp. Whenn e reaches 2.8 − 2.9 × 10 19 m −3 (corresponding to 0.37-0.38 n GW where n GW is the Greenwald density limit [21]), there is an abrupt increase of dn e /dt andn e deviates from the linear level previously input to the feedback density control system. With the increase of dn e /dt, the inner target (LP #42-44) and outer target (LP #51-54) probe ion saturation current signals decrease simultaneously within 10-30 ms by 40%-70% compared to the peak value. Probe LP #41, located in the inner target (IT) far-SOL, does not show such a large decrease. The feedback density control system operates immediately after the increase of dn e /dt and reduces Γ D2 to restoren e to the desired linear phase (figure 2(b)). As a result, n e decreases again to the demand value and the target ion saturation current abruptly rises. This process repeats itself four times (shown as dashed vertical lines in figure 2) during the latter part of the density ramp, after which the LP signal remains at a reduced level.

Global description of the experimental results
Although the plasma control system is programmed to keep the strike point positions constant during the discharges, as the plasma density increases, there was an unintended, steady inward movement of 2-3 cm of the radial location, R OSP of the outer strike point throughout the density ramp (figure 2(f )). This is of the order of the spacing between the target LPs. There is, however, no significant change of R OSP during then e fluctuations, so that the abrupt change of the divertor LP signals is not caused by bulk motion of the magnetic equilibrium. The implication is that the rapid change of the target ion flux is the result of a transition of the divertor condition, i.e. particle detachment and re-attachment, caused by the density ramp.
Although not shown here, these ion saturation current fluctuations were also observed in similar, but less optimized L-mode density ramp experiments with density feedback control in the 2017 campaign (discharges #19086, #19088, and #19089, with slightly different discharge conditions total heating power 0.7 MW, I p = 0.5 MA and B T = 1.8 T). Shot sequence of the density ramp experiment with feedback control of density (#20265). (a) Line-averaged electron densityne, (b) deuterium gas puff rate Γ D2 , (c)-(e) parallel target particle flux measured by LPs (see [19] showing the location of the specified channel number as LP ##) and (f ) radial position of the outer strike point. The LP channels LP #41-44 are positioned on the inboard divertor target and LP #51-54 are on the central divertor target. The dashed line indicates the time at which the probe signal falls rapidly.
To isolate the effect of the gas modulation on the observed abrupt detachment, then e feedback control was deactivated in #20267 beforen e reached the value of 2.8 − 2.9 × 10 19 m −3 . The shot sequence of #20267 is shown in figure 3. The increase in dn e /dt occurred atn e = 2.8 − 2.9 × 10 19 m −3 (figure 3(a)), as in the discharge with density feedback control. In this case, the cliff transition of the divertor condition was observed without fueling modulation, but the saturation current oscillations do not occur because the gas puff strength is fixed. Similarly to discharge #20265, all probe signals drop simultaneously with increase of dn e /dt, except for at the inner target far-SOL (LP #41). The target ion flux reduction is largest in the near-SOL (LP #42 and LP #53) for both targets.
When studying the dynamic link between target profiles and upstream density, there is always a trade-off in these types of detachment experiments. Unlike the L-mode density scan experiment conducted in the 2017 KSTAR campaign [22], there is no gap in the density level during a density ramp experiment, making it easier to identify dynamics or bifurcations that are highly sensitive to upstream density. Sincen e was fixed during each of the discharges in these earlier density scan experiments, only fourn e levels were available, and the cliff target flux feature could not be captured. As shown in figure 4, the rollover of the target particle flux, characterized by the start of a decrease in signal after the peak, begins before the cliff (or bifurcation-like) drop occurs (see section 2.3 for a discussion on the equivalence of labeling this behavior as cliffor bifurcation-like). Therefore, it is necessary to distinguish between the 'onset of particle flux rollover' and the 'cliff drop (further detachment)'. The former, which serves as the criterion for the beginning of detachment, has a slower timescale than the latter. Although there is an asymmetry in the rollover of the particle flux [22], the bifurcation-like transition occurs symmetrically at both targets.
Finally, figure 5 shows the normalized beta and stored magnetic energy, which serve as measures of core confinement. As the target flux rollover progresses, both variables gradually decrease. However, after the 7 s, when the cliff is encountered, only minor fluctuations or slight drops are observed, and no significant trend difference arises. Thus, the impact of the target flux cliff on core confinement can be considered negligible.

Bifurcation-like transition of divertor condition
The bifurcation-like transition of the divertor condition is not a consequence of the fueling modulation since it occurs both with (#20265) and without (#20267) density feedback. This implies the existence of a critical upstream density which significantly changes the divertor condition. Figure 6 focuses on then e threshold for the 3-4 main Γ ∥t 'detachment/reattachment (cliff)' events in #20265 and the rapid Γ ∥t reduction in #20267. The thresholdn e values are taken at the point of maximum rate of change of dΓ ∥t /dt. The blue and orange bars indicate detachment start and endn e threshold values when Γ ∥t is on the cliff (end of the attached branch) and below the cliff (beginning of the detached branch), respectively. The red bar shows the re-attachmentn e threshold when Γ ∥t is under the cliff (end of the detached branch). There is no significant difference (∆n e,cliff start /n e,cliffstart < 3%) on the detachment threshold density between the five abrupt detachment/reattachment events in figure 6. Therefore, ∆n e,cliff start appears invariant under modulation of the fuel gas throughput (waveform of Γ D2 (t)).
The re-attachment thresholdn e values are higher than the detachment thresholdn e values by about 4-7 %, which is significantly higher than the measurement error level of mmWI, i.e. less than 0.25 % (n e ± 0.5 × 10 17 m −3 equivalently). This demonstrates that there is some slight hysteresis of the 'critical upstream density' or lagging ofn e compared to the change happening in the divertor region. Although there may be a better local metric (e.g. neutral pressure/density or target electron temperature) as a detachment or re-attachment cliff indicator, were unavailable on KSTAR the time of the experiments described here and so the study of these abrupt transitions relies almost exclusively on the upstream line averaged density and the response of the local divertor target ion flux.
Both the cliff detachment and re-attachment occur within 10-30 ms. The implication is that the intermediate state in the middle of the target ion flux cliff is unstable. Therefore, the abrupt transition between the detached and re-attached divertor condition can be referred to as a bifurcation-like transition in which the target ion flux selects one of the detached or attached branch (these will be referred to as XPR and non-XPR branch later in this paper). The timescale of the abrupt detachment is similar to, or longer than the parallel SOL transport timescales, but is faster than the diffusive processes which occur in the edge of the core plasma [12]. The evolution of Γ ∥t during the abrupt transition in KSTAR is almost linear, so the process does not appear to be enhanced by a feedback mechanism. There were no significant differences in the cliff fall and rise time. Figure 4 focuses on the behaviour of the Γ ∥t profiles at the inner and outer targets in pulse #20267 (the discharge with  feedforward density control) before and after the cliff transition. For the central divertor target, additional target probes (LP #106-114) in different toroidal locations have also been included to improve spatial resolution ( figure 4(b)). The strike point at each time was found from the magnetic equilibrium reconstruction. The 'cliff' of target particle flux at t ∼ 7 s (profile No. 11-12 in figure 4) is observed across the whole profile and therefore indicates complete detachment at both inner and outer targets. The target flux drop during 'rollover' (slow timescale behavior from peak to just before cliff transition) is larger for the outer target than the inner target. After the cliff transition, the target flux of outer target is 30% of the peak, while the inner target flux is 60% of the peak value. This behavior is translated to the asymmetric detachment observed in Lmode density scan experiments in which the outer target flux rollover occurs at a lower upstream density than at the inner target [22]. Thresholdne for the four detachment and three re-attachment events occurred in #20265 and a detachment event occurred in #20267. Line averaged electron density value was taken when the slope of the probe signal curve changed the most. Figure 7 shows the time series of #20267 including the degree of detachment (DOD, [12]) calculated using the near SOL LP and the radiation power measured by the IRVB [20,23]. An abrupt increase of DOD in near SOL of the both targets was observed around 7.0 s (figures 7(c) and (d)). The decrease of the outer target DOD from 3.5 to 5.0 s is due to the inward  figure 3(f ). Therefore, the outer target DOD was normalized at 6.0 s, during which time the strike point deviation from the probe measurement area was negligible. Since the strike point was not significantly changed along the inner target (vertical direction), the inner target DOD was normalized at 3.5 s.

Coupling of X-point radiation and abrupt transition of divertor conditions
The total radiated power (figure 7(e)) increases by 20% during the 200 millisecond interval before the abrupt increase of the DOD. A similar 10%-20% increase of the total radiated power was also observed in #20265 for all four successive bifurcations. The distribution of the local radiation at the 9 selected channels of IRVB changes significantly before bifurcation happens (figure 7(f )). The position of the selected channels are shown in the figure 8(a) as cross markers (the channels from top to the bottom correspond to core 1-4, X-pt. and SOL 1-4, respectively). Before the cliff transition, the channels just above the X-point (core 3-4) and on the X-point increase 50-100% and the channels below the X-point decrease over 200 millisecond (figure 7(f ) 6.8-7.0 s). After the cliff transition, the radiation at the innermost channel (core 1-2) increases and the channels in vicinity of the X-point (core 4, X-pt.) decrease, which demonstrate further movement of radiation zone from the X-point to the core (figure 7(f ) after 7.0 s).
The two-dimensional distribution of the radiation power density obtained by tomographic inversion of the IRVB signals [20] is shown in figure 8. The radiation distribution does not change significantly from 5.0-6.7 s. To see the change of the radiation zone clearly, the averaged radiation distribution shown in figure 8(a) was subtracted from given radiation distributions in each time instant shown in figures 8(b)-(h). The abrupt drop of the target flux occurred between t = 6.97-6.99 s (figures 8(d) and (e)). An additional radiation zone is formed at the X-point, and the radiation power increases (figures 8(b)-(d)). After the cliff transition, the X-point radiation slightly reduces and the center of the radiation zones moves further into to the core region (figures 8(e)-(h)). Since the carbon radiates dominantly in the temperature range of 6-11 eV [24], the increased radiation near the X-point implies cooling in the Xpoint vicinity close to 6-11 eV. Due to the limited divertor spectroscopy, it is not possible to determine which lines are dominant among the total radiation from IRVB, but based on the experiments in other carbon devices such as JT-60U [25] and TCV [26,27], it can be safely assumed that most of the radiation is coming from the carbon impurity. The increase of the radiation power density just above the X-point remains stable, so the results indicate an XPR regime as observed in [28] and not an X-point MARFE, [29][30][31] which is MHD unstable. Therefore, in this paper, the branch before the cliff transition will be referred to as the non-XPR branch, and the branch after the cliff transition will be referred to as the XPR branch. This avoids confusion with onset of rollover, which is often considered one of the criteria for detachment.

Setup of SOLPS-ITER simulation
The SOLPS-ITER 3.0.8 master branch has been used for all simulation results presented in this paper. The simulation grid is the same as that used in [22], but the pump and gas puff locations are modified as shown in figure 9, corresponding to the experimental conditions. In figure 1, 10 contours intersect the lower outer target SOL, and the last 2 of these touch the upper outer target. The number of contours can be varied by the level settings, but based on the EFIT contours and accurately converting this to the outer midplane (OMP), the upper target begins to receive flux from a flux tube that is approximately 8 mm further away in terms of distance at the OMP, compared to the 1st flux tube in the lower target. It can therefore be safely assumed that this particular case can be simulated as a pure single null configuration.
The effective pumping speed of the vacuum vessel (VV) pump and two cryo pumps are 9.7 m 3 s −1 and 25 m 3 s −1 , respectively [32]. Using the approach in [33], this yields pump absorption coefficients for the SOLPS-ITER simulations of 0.0018 and 0.017 respectively. The changes in the pump configuration and absorption coefficients do not significantly affect the divertor performance, but do affect the local plasma conditions near the gas injection location and pumps. The power crossing the core boundary is set to 0.8 MW, equally distributed between the electrons and ions. The plasma species consist of deuterium and sputtered carbon, with a chemical sputtering coefficient of 1% for both the target and first walls. The core particle source is given by a fixed flux boundary condition of 8 × 10 19 s −1 considering the source from the neutral beam. It should be noted that fluid drifts are not included in these simulations. Steady-state solutions are found with fueling throughput scans in the range of 5 × 10 20 − 3 × 10 21 atom s −1 . Time-dependent simulations are carried out for a fixed fueling throughput ({3, 4, 5, 8} × 10 21 atom s −1 ) starting from an initial condition obtained from the steady-state solution with the highest density (n e,sep = 3 × 10 19 m −3 ).
In this study, the quasi-steady state (QSS) EIRENE approach [34] is applied to time-dependent SOLPS-ITER simulations since the neutral flight time does not significantly affect the X-point radiation-induced cliff. This means that while the B2.5 time step is set to be 1 × 10 −4 s, the EIRENE time step is set to a larger value to obtain a fully relaxed neutral distribution and corresponding source terms at every EIRENE call. A B2.5 time step of 1 × 10 −4 s is large enough to ensure a correct solution and remains practical for performing timedependent simulations over a period of several seconds of Figure 9. SOLPS-ITER simulation grid. Magnetic equilibrium is the same as in [22], which is a typical KSTAR single null L-mode discharge with strike point on the central divertor target. Pump and gas injection locations are shown. The colored selected areas marked with numbers indicate x core region of the OMP y core region just above X-point z SOL region just above X-point { 1st SOL ring from OT to OMP. plasma time. In the simulations used in this paper, the time dependent mode in EIRENE is enabled (NFILEJ = 1 in the code). Consequently, trajectories exceeding the EIRENE time step are cut off (based on flight time), and information regarding these long-lived neutrals is recorded in a census array. However, to run the code in a QSS fashion, the maximum number of time-dependent neutrals is suppressed to less than 1% of the total number of the Monte Carlo trajectories. The EIRENE time step is set to 1 s for some time-dependent gas throughput ramp-up and down simulations in this paper, which is sufficient for the neutrals to be fully relaxed. However, other timedependent simulations in this study use an EIRENE time step of 1 × 10 −3 s, which is the default for the code. The effect of using different time steps is discussed later in this paper. Throughout this paper, the aforementioned method is simply referred to as the QSS EIRENE scheme. If the time-dependent  figure 11. The lowest density case in steady-state solutions is marked with an 'x', and the peaked particle flux cases are indicated with triangle markers for both IT and OT. mode in EIRENE is turned off (NFILEJ = 0 in the code), then the EIRENE time step is irrelevant and the EIRENE trajectories are only terminated by an event (e.g. ionization) or the maximum allocated CPU time. Although this method is perhaps best referred to as the QSS method, the solutions are not significantly different between the two methods because the EIRENE time step of 1 s is much larger than both the B2.5 time step and the mean ionization timescale.
In order for the simulation time evolution to have a physical meaning, numerical acceleration schemes and other settings that affect the timescales must not be used.
Here the under-relaxation parameter (rxf), internal iterations (nstg2), and numerical acceleration factors (defined in b2.numerics.parameters) are all set to 1. If the dynamic mode of B2.5 is used (changing the value of the 'b2mndt_style' switch from the default value of 1 to 2), the equations are solved in a different order from the default setting. In dynamic mode, the density equation is first solved, followed by the parallel momentum, electric potential, and electron and ion internal energy equations, so a more accurate discretization can be applied to the momentum and energy equations, and more accurate time evolution can be expected. Both the default time advance scheme of SOLPS-ITER and the dynamic mode were applied here, and the effect will be described later.
Overall divertor performance can be evaluated by changing the divertor condition through time-dependent simulations which linearly ramp the fueling throughput. In order to ensure that the time-dependent solution at each time is sufficiently quasi-steady state (QSS), 1 s was used for each ramp up and down. There are various methods to guarantee that each moment is sufficiently QSS (or converged), for example, by adding fixed gas throughput sections in between ramp up stages, forming a stair-like progression, and checking if the solution change is small enough in the fixed sections, or by comparing the time-dependent trace with a few steady-state solutions plotted in the phase space of n e,sep and T et [35]. In Lmode cases, simulations were performed with EIRENE time steps of either 1 × 10 −3 s or 1 s. For the H-mode-like condition, only the heating power was changed to 6 MW, and the rest of the conditions were kept the same, without including a pedestal or other H-mode specific features. The case of D only (EIRENE time step of 1 × 10 −3 s) and the case of D+C (EIRENE time step of 1 s) were also simulated.

Time evolution of the plasma parameters during XPR
In SOLPS-ITER simulations, the global particle balance is achieved by a fixed core boundary source representing the neutral beam source, fueling throughput, and pumping. The actuator that mainly changes the divertor condition is fueling (in the absence of the seeded impurity or pellet injection, etc), and most of the upstream and downstream parameters have a monotonic dependence with respect to fueling throughput.
Target particle flux is one of the exceptions. In a low recycling condition, it increases for gas throughput, but then turns to a decreasing trend at some point. This rollover is often judged as the starting point of particle detachment. The change in particle flux with respect to throughput is continuous even after rollover occurs, but when the throughput exceeds a certain value, the steady-state solution gives a drastically reduced target flux. That is, it shows a similar behavior as the target flux cliff observed in the experiment. Steady-state solutions before and after the cliff share the same gas throughput, but depend on the initial condition over a certain range. Therefore, there is hysteresis or an unstable intermediate state, resulting in bifurcation-like aspect in which one of the two branches is selected. The intermediate state along the cliff is only accessible via time-dependent simulation.
To clarify once again, the following discussions will involve a mixture of steady-state solutions (denoted as SS or only density) and time-dependent solutions (denoted as only time or density and time together), as described in section 3.1. In other words, we first obtain steady-state solutions in the range of n e,sep ∼ 5 × 10 18 − 3 × 10 19 m −3 , and then show the results of time-dependent fixed gas throughput simulations using the solution with the largest n e,sep of 3 × 10 19 m −3 as the initial condition. Figure 10 shows the evolution of the simulated inner and outer target fluxes mapped to the OMP. Rollover occurs with increased n e,sep with increased fueling throughput. The asymmetry in rollover between the two targets is consistent with observations from a previous KSTAR L-mode experiment where the outer target begins rollover at a lower n e,sep than the inner target [22]. Once n e,sep exceeds 3.01 × 10 19 m −3 , a rapid decrease in target flux occurs. Therefore the time-dependent snapshots are shown along this cliff in the absence of the steady-state solution. The case of fueling throughput of 3 × 10 21 s −1 with the default time-advance setting ('b2mndt_style = 1') are shown here. The drop of the target flux occurs within 100 ms before it reaches the steady-state solution. The value of the particle flux peak is consistent with the experimental value within a factor of two and qualitative trends are matched: the peak after the cliff is larger for the inner target than for the outer target. The target flux value after the cliff is significantly lower than the measured value, which is due to the absence of the core volume in the simulations, and excessive radiation due to the effect of boundary condition selection. This will be covered in detail later.
Density and radiated power over time are shown in figure 11. Density in all regions increases with gas throughput in steady-state. The total radiated power increases monotonically with gas throughput, but the main radiation zone moves from the divertor to the upper SOL above the X-point. As the cliff transition occurs from t = 0, the core density and total number of electrons in the core monotonically increase with time, while n e,sep decreases andn e remains almost constant. The increase in the total number of electrons in the core is consistent with the experimentally measured dn e /dt increase, but the discrepancy with the simulatedn e comes from the absence of the core region in the simulation. Since the OMP density profile is monotonic, the core region makes a dominant contribution to the experimentaln e . This is not captured in the SOLPS-ITER simulation, which models only the edge, SOL, and divertor regions. The decrease in n e,sep can be attributed to the following three factors: (a) a rapid decrease in the dominant neutral source (net contribution, other than neutral to neutral reflection) due to the target flux cliff, leading to neutral source depletion for SOL ionization; (b) a decrease in the quantity of ionization due to SOL power depletion caused by an increase in radiation power (mainly in the core region); (c) a decrease in the ionization rate due to excessive cooling in the SOL. These three factors compete with the opposing effect of an increase in the core ionization source and the associated radially transported source. However, n e,sep ultimately decreases. The radiated power continues to increase during the cliff transition and saturates when it reaches steady-state at t = 0.1 s. The main radiation zone shifts from the upper SOL to the core, and the radiated fraction there becomes 75% of the P SOL in the steadystate (t = 0.1 s). Figure 12 shows calculated two-dimensional distributions of radiated power density, T e and ionization source just before the cliff transition (highest density SS solution) and in the middle of the cliff transition with XPR (at t = 50 ms). The radiating zone shift from the upstream SOL to the core is confirmed to occur mainly through the X-point (figures 12(a) and (b)), as observed experimentally. The T e distribution is also consistent with the radiated power distribution since the carbon radiation is very sensitive to the background plasma temperature (figures 12(c) and (d)). In particular, the 5 eV front moves from the X-point (before XPR) to the core, accompanying the ionization front and corresponding ionization sources (figures 12(e) and (f )). This provides additional ionization sources to the core from the X-point. Thus the impact on the Figure 11. Simulated density and radiated power at various locations (see figure 9). 'SS' corresponds to the steady-state solutions from the density scans, plotted in increasing order of ne,sep or gas throughput. The time-dependent scan begins just before the X-point radiation at t = 0, starting from the highest density steady-state solution (ne,sep = 3 × 10 19 m −3 ). ne,core and ne,sep are measured at the OMP. Ne,core represents the total electron particle number in the core region. The inner/outer SOL and upper SOL boundaries are divided by their respective inner/outer poloidal locations at the X-point. core density and total number of electrons in the core is immediate, while it takes time to increasen e which is provided from further upstream ( figure 11).
The time evolution of the various parameters depending on the gas throughput level and the default/dynamic mode of SOLPS-ITER is shown in figure 13. As for figure 12, the initial condition for each was set to the highest density case of the SS solution. The first observation is that the larger the gas throughput, the shorter the cliff transition timescale. The simulation gas throughput can be similar to the experimental value only by the order of magnitude [36][37][38]. The experimental conditions of a toroidally localized gas injection location and pump are simulated as toroidally symmetric systems in SOLPS-ITER.
As shown in figure 11, n e,sep is decreasing in time but the trend changes at some point (t = 80 ms for the case of lowest throughput) when the core ionization source is provided through the X-point. This process is accelerated with the higher gas throughput cases. Both the inner and outer target temperatures decrease steadily with time, and are saturated to 0.2 eV and 0.1 eV respectively. The same is also true of the target flux with the timescale becoming shorter when the throughput is large. In the dynamic mode of SOLPS-ITER ('b2mndt_style = 2') the cliff drop timescale is significantly reduced. This is presumably because the more accurate discretization better accounts for the effect of the density change with each iteration in the momentum and energy equations. The choice of default or dynamic mode affects the timescale, but does not change the simulation trend. However, the dynamic mode requires a smaller B2.5 time step (less than 1 × 10 −5 s) for numerical stability, while the default mode remains stable with a time step of 1 × 10 −4 s. Considering that the time taken for the target flux cliff transition in the experiment is around 20-30 ms, it agrees well with the timescale obtained from the dynamic mode simulation.

Coupling of fronts with target electron temperature
As figure 12 indicates, the radiation front, T e = 5 eV front, and the ionization front appear to be coupled. Since KSTAR is a carbon device, sputtered carbon serves as the main radiating impurity species for cases without seeded impurities. Since the impurity emissivity strongly depends on T e , the radiation front can be characterized by the dominant radiating species and the dominant radiating temperature. The average of T e weighted by the carbon radiation power and its standard deviation is defined as follows ⟨T e ⟩ P rad,C ≡ ∑ T e P rad,C ∆V ∑ P rad,C ∆V (1) δ⟨T e ⟩ P rad,C ≡ √ ∑ (T e − ⟨T e ⟩ P rad,C ) 2 P rad,C ∆V Here, the sum is performed over the entire grid, and M is the number of non-zero weights, that is, the total number of cells (96 × 36 in case of the simulations described here). ∆V represents the cell volume. Higher charge state ions (C 4+ , C 5+ , C 6+ ) have a much higher ⟨T e ⟩ P rad,C than the other lower charge state ions ( figure 14). However, the radiated power contribution of the higher charge state ions is negligible. For steady-state solutions, as n e,sep increases the total radiated power also increases. During the cliff transition with XPR, the radiated power steadily increases in time, but decreases at the end of the transition (0.1 s). This is because C 2+ , C 3+ cannot effectively radiate due to excessive core cooling. Focusing on the lower charge state ions as the main radiator, the main emission temperature ranges can be represented by ⟨T e ⟩ P rad,C , which are 1 − 3, 5 − 8, 10 − 11 eV for C 1+ , C 2+ , C 3+ , respectively. The main radiation temperature of each lower charge state is quite insensitive to changes in the divertor condition (including the cliff) because it only changes by a few eV. The only exception is for the lowest steady-state density case (n e,sep = , ⟨T e ⟩ P rad,C due to the absence of the low temperature (1-20 eV) region in this attached condition. A similar deviation comes from the end of the cliff (t = 0.1 sec) with lack of the high temperature zone due to the limitation of the extent of the simulation grid and the core boundary condition. The ionization zone (or front) is also coupled to the T e , in a similar manner as the carbon radiation zone. The ionization potential of the main ion is 13.6 eV, but since the ionization front is determined by several factors (the electron density and temperature-dependent cross-section, atom density, and electron density), it is mainly located in the several eV region. The ionization front and radiation front can be defined in terms of the average location in the normalized parallel coordinate s ∥ , from the target to the mid-plane, weighted with the ionization source and emissivity respectively, as follows [22] Figure 15. 5 eV front (1st row), EI front (2nd row) and radiation front (3rd row) on several individual flux tubes are shown as a function of Tet at the inner (left column) and outer target (right column). s1-s5 in the legend indicate 1st to 5th SOL rings and the yellow shaded region corresponds to X-point location (appears to be a band due to the variation of the X-point location on s ∥ depending on the radial location: s1-s5).
Here, integration is performed from the target (t) to the upstream (u, selected as OMP for the outer SOL and IMP for the inner SOL), S iz is the ionization particle source from EIRENE, and P rad is the radiated power summed over species. The electron temperature and its parallel gradient are characterized by T et , considering the monotonic distribution of T e along the parallel coordinate. Therefore, the 5 eV front is strongly coupled with T et for both targets ( figure 15). The ionization front and radiation front are also well matched with 5 eV front location, and are hence correlated with T et for both targets. Since T et is strongly correlated with various divertor parameters [2,3,39], it also well characterizes the front location in parallel coordinates. Therefore, T et is a key parameter for characterizing the cliff induced by the XPR. The relationship between the 5 eV front location and T et is affected by whether the inner or outer target is considered, geometry, and the volumetric dissipation mechanisms. Details are beyond the scope of this paper. Electron density, radiated power, and ionization sources peak near the target for the low recycling condition (figure 16). As the divertor detaches, the electron temperature gradually decreases but still maintains a monotonic distribution along the parallel coordinate. The 5 eV front continues to move upstream and then crosses the X-point when n e,sep reaches 2.06 − 2.71 × 10 19 m −3 . This observation is consistent with the outer SOL 5 eV front in figure 15 (1st row, 2nd column), where the front can be seen crossing the X-point at these n e,sep values. The peaks of radiated power and the ionization source also move upstream and cross the X-point with the formation of the electron density peak slightly below the ionization source peak. In the XPR regime, the electron density, temperature, radiated power and ionization source of the 1st SOL ring gradually decrease with time.
Decrease of the electron density, temperature and radiated power peak at the 1st SOL ring during the XPR is due to the shift of these peak towards core by the X-point cooling and perpendicular transport ( figure 17). The radial electron temperature profile at the X-point is almost monotonically sustained from the core to the SOL by conduction. As the SOL continues to cool, the separatrix temperature also decreases, and when it reaches ∼5 eV, the ionization zone and radiation zone also pass through the X-point (separatrix) from the SOL, and the ionization source and radiated power peak at the X-point. The peaking of the radiated power (energy sink) causes the radial electron temperature profile to become locally non-monotonic at the X-point. From this point, core cooling proceeds rapidly both by the ionization source and the radiation that penetrates into the separatrix. In the XPR regime, the area right above the X-point has a temperature near 5 eV and becomes a zone where ionization and radiation are concentrated. To summarize, as the detachment progresses, the 5 eV front first moves from the target to the vicinity of the X-point, then penetrates into the core by perpendicular transport, accompanied by an ionization and radiation front.

Phase space behavior of KSTAR L-mode divertor SOL plasmas
In the absence of a seeded impurity and at constant input power, fueling is the only actuator which can change the divertor condition. By scanning the fueling throughput in a QSS manner, the operation space of the divertor SOL system can therefore be obtained assuming that the system is isolated from the core transport. The ramp up and down of fueling through- requires an order of magnitude higher fueling throughput compared to other cases for scanning from attached to detached conditions. Since the overall density level is low in L-mode, the neutral mean free path and corresponding ionization timescales are relatively larger than those in H-mode. This increases the probability of the existence of long-lived neutrals beyond the default EIRENE time step of 1 × 10 −3 s (see section 3.1). In the QSS EIRENE scheme used in this paper, most of the long-lived neutral information is cut-off, so the source from these neutrals (namely the ionization source) is underestimated. Therefore, the fueling requirement to obtain a similar ionization source is much larger than in the case in which the EIRENE time step is set to 1 s. In the H-mode cases for this simulation, since the density is generally high, the neutral path is shortened, and thus the effect of the EIRENE time step is relatively small. The change in the phase space due to the difference in EIRENE time step is not very large. However, if the time step is too small, the local ionization source by the gas puff stratum near the fueling location can be overestimated because the fueling demand may be large enough to be comparable to the recycled flux. Therefore, the distribution of ionization sources (and other sources originating from EIRENE neutral) can be distorted.
In these KSTAR SOLPS-ITER simulations, the temperature of both targets exhibits a similar trend in response to fueling throughput. It is therefore sufficient to examine only one phase space for either the inner or outer targets. In H-mode, n e,sep steadily increases as the fueling throughput increases. In L-mode, this trend continues, and when the T e,OT reaches about 0.7 eV, it leads to a rapid decrease in n e,sep . This coincides with the point at which the 5 eV front passes through the X-point shown in figure 15, so it is an XPR induced transition. In H-mode, the phase space curve is monotonic and hysteresis is insignificant. H-mode essentially has a phase space curve located in a higher n e,sep and higher T et region than L-mode. Similarly, in order to drive an XPR induced transition in Hmode, the X-point temperature must drop to several eV, and accordingly, the target temperature must drop to at least less than 1 eV. However, excessive fueling for this is limited by the density limit before T et drops below 1 eV. In fact, even in the SOLPS-ITER calculations, n e,sep steadily increases beyond 1 × 10 20 m −3 , which is a numerical solution. This is presumed to be the reason why the target flux cliff has not observed in the KSTAR H-mode experiment (see section 4.3 for further discussion). In D+C H-mode, it can be seen that the phase space curve shifts downward compared to D only H-mode, which is the effect of cooling by sputtered carbon radiation. It is expected that the same effect can be obtained through impurity seeding, and that a regime that shows a cliff could be found under H-mode conditions through appropriate level of impurity seeding.
The phase space of L-mode showing the XPR induced cliff structure is analyzed in more detail in figure 19. First, in the ramp up phase, n e,sep rises and, when n e,sep reaches ∼3.2 × 10 19 m −3 , T e,OT = 0.7 eV and the XPR proceeds.
Here, the first cliff transition begins, and n e,sep drops sharply to 1.3 × 10 19 m −3 while T e,OT changes by only ∼0.1 eV. This cliff transition is unstable, so no steady-state solution exists within the region indicated by the unstable zone in figure 19. This means that the QSS solution converges by selecting one of the non-XPR branch and the XPR branch. After the first cliff transition, the phase space point leaves the unstable zone and enters the XPR branch. Here, as fueling increases further, T e,OT steadily decreases, but n e,sep shows rollover. This is because, as shown in figures 11-13, the X-point ionization source gradually affects the mid-plane. At T e,OT = 0.2-0.3 eV, as the fueling starts to decrease again, the phase space point reaches the boundary of the unstable zone with a small hysteresis, and moves to the non-XPR branch as the second cliff transition occurs. The following two states, 'the state just before entering the 1st cliff transition' and 'the state immediately after experiencing the second cliff transition', lie on the same non-XPR branch and have the same T e,OT . However, the L-mode phase space has a hysteresis characteristic in that n e,sep can be very different for a given T e,OT (or T e,IT ) value. When fueling gradually decreases from the state after the second cliff transition, the hysteresis gradually disappears and returns to the initial state, while tracing the same phase space trajectory as the ramp up phase. The cause of hysteresis is the very different plasma parameters of the non-XPR and the XPR branches, so the timescale for pressure balance to QSS is different.
Finally, the directional property of the phase space was investigated. As can be seen from the H-mode phase space in figure 18, hysteresis is insignificant during ramp up and ramp down in H-mode. Therefore, any phase space point on the curve has bidirectional property and can be moved towards the desired direction on the curve freely through the control of the actuator (here, fueling). However, in Lmode, due to hysteresis, it is difficult to determine the directional property with only this information. Moreover, since the L-mode phase space curve in figure 19 is a trajectory due to a single fueling throughput ramp up and down, it is impossible to determine which trajectory will be followed if the sign of dΦgas dt is changed at an arbitrary point on the curve.
In order to elucidate the directional property of the L-mode phase space curve, the two most extreme points in non-XPR branch, 'the state just before entering the 1st cliff transition' and 'the state immediately after experiencing the second cliff transition' were selected and the gas level ramped down and up from there (figure 20). For convenience, the before and after cliff transition states will be referred to as states 1 and 2, respectively (denoted by red and yellow stars respectively in figure 20). First, in state 1, the new trajectory (red curve in figure 20) shows a slight hysteresis, but follows almost the same trajectory as when ramping up. This is similar to the small hysteresis seen in the XPR branch. On the other hand, since T e,OT in state 2 is the temperature at which a cliff transition can occur immediately, it can be considered that if a fueling ramp up is performed in this state, it will immediately enter the unstable zone and move to the XPR branch. However, in reality, the state 2 point returns to the curve at the first ramp up phase along the green curve in figure 20, and then passes through state 1 again, causing a cliff transition. Therefore, in state 2, whether fueling is reduced or increased, it still has the characteristic of returning to the QSS line shown in the ramp up on the non-XPR branch, but cannot pass through the unstable zone immediately from the state 2. What can be seen here is that the phase space point in the same branch always has the characteristic of returning to the QSS line regardless of the presence of hysteresis. In addition, since entry and exit to other branches through the unstable zone is directional, transition is possible only when a specific state at the boundary of the unstable zone is reached. Also, in the non-XPR branch, the two densities coexist at the boundary of the unstable zone. In the XPR branch, only a single state exists at this boundary.

Quantitative inconsistency of SOLPS-ITER solutions with experiment
The XPR branch solution can be considered as a kind of collapsed state solution shown in the case of ITER SOLPS modeling, in which a code bifurcation is caused by the core penetration of recycled neutrals with a radiation zone that yields no steady-state solution in between two branched solutions [40]. This is consistent with the KSTAR experimental observation that the intermediate state on the cliff is unstable and should go through a transition to a non-XPR or XPR branch.
In terms of timescale and qualitative representation of the phenomenon, the SOLPS-ITER simulation results for KSTAR demonstrate good consistency with experiments. However, there are quantitative inconsistencies, such as excessive core boundary cooling and an excessive target flux drop. Since the core region is not included in the SOLPS-ITER simulation grid, the code only solves a portion of the core just inside the separatrix. As a result, the radiation zone continues to cool with the formation of the XPR and penetrates into the core, eventually reaching the core boundary. At this point, most of the input power provided as the core boundary condition is emitted as core radiation, resulting in a collapsed solution ( figure 11). To prevent excessive core cooling, a deeper core grid can be used, or artificial heating power can be added to the core boundary (e.g. by applying a fixed temperature boundary condition) to maintain the temperature gradient and prevent the carbon radiation zone from approaching the core boundary. Additionally, since fluid cross-field drifts generally result in lower electron density and higher electron temperature on the high field side [41], including drifts will compensate for the bias of the radiation spot toward the low field side core boundary during collapse.

Two-point model formatting analysis of the target ion flux cliff drop
Since there is a decrease in P SOL due to the core radiation, there can be two possible scenarios which explain particle detachment according to the two-point formatting (2PMF) proposed in [3,42]: (i) a decrease in the basic two point model (2PM) term, which is proportional to p 2 tot,u /q ∥u , and (ii) a decrease in the volumetric loss term f Vol−loss Γ e∥t of the target flux in 2PMF, primarily resulting from momentum loss along the SOL due to plasma-neutral friction. Here, p tot,u is the upstream total pressure and the q ∥u is upstream parallel heat flux. These two effects are coupled because cooling of the SOL enhances the momentum and power losses due to the plasma-neutral interactions [22]. In the KSTAR experiment, the lack of measurements in the upstream and downstream SOL, mean that it is difficult to experimentally determine which one has the most significant effect on particle detachment. Therefore, 2PMF analysis was performed for the time-dependent SOLPS-ITER solution as shown in figure 21. The volumetric loss factors f mom−loss , f cooling , f Vol−loss Γ e∥t and the basic 2PM term of target flux Γ basic e∥t are defined in [3].
Since there is a decrease in P SOL due to the core radiation, there can be two possible scenarios to explain particle detachment according to the two-point model formatting (2PMF) proposed in [3,42]: (i) a decrease in the basic two-point model (2PM) term, which is proportional to p 2 tot,u /q ∥u , and (ii) a decrease in the volumetric loss term f Vol−loss Γe ∥ t of the target flux in 2PMF, primarily resulting from momentum loss along the SOL due to plasma-neutral friction.
2PMF analysis should be performed for a specific flux tube in the SOL, not for the integrated flux along the target. The outer SOL flux tube on which the peak target ion flux is located at t = 0 was therefore selected as representative. In addition, the OMP was selected as the upstream location. If the upstream is designated as the X-point, the difference in the resulting trend is not significant. In figure 21, the basic 2PM term decreases with time and then saturates to ∼0.6. The absolute magnitude of the volumetric loss term starts at less than 0.05 at t = 0 because of the momentum dissipation that occurs in the SOL during the continuous rollover of target flux. It then decreases monotonically with time, and since it drops by more than 80%, the decrease is larger than that of the basic 2PM term.
The 2PMF volumetric loss term of the target flux has a nonlinear relationship with the momentum loss factor and cooling factor.
(1 − f cooling ) is nearly uniform over time. This is because q ∥u decreases due to an increase in core radiation, but the SOL radiation also decreases because the radiation zone moves from SOL to core ( figure 11). On the other hand, (1 − f mom−loss ) continues to decrease with time due to further cooling of the SOL flux tube (mainly from molecule-plasma interactions), and its square contributes to f Vol−loss Γ e∥t . Therefore, the target flux cliff is the result of the simultaneous contribution of (i) the basic 2PM term corresponding to the upstream change and (ii) the momentum loss that occurs in SOL. The contribution of the momentum loss more significant. On the other hand, the continuous target flux drop before bifurcation (corresponding to the rollover in experiments) in the density scan simulations in [22] is primarily caused by the volumetric loss term associated with momentum loss. A notable difference between the simulation cases in [22] and those in this paper is that in the cited paper, Γ basic e∥t consistently increases with n e,sep without any decrease.

Exploring XPR-induced cliff behavior in KSTAR H-mode condition
The target flux cliff was observed in all density ramp experiments for KSTAR L-mode. However, in H-mode, with a heating power of 3 MW, the detachment pattern of the target flux in the D 2 puff experiment consistently appears smooth without cliff behavior. The line-averaged mid-plane density was raised up to 6 × 10 19 m −3 , corresponding to 0.7 of the Greenwald fraction, indicating that a high-density background was achieved through high fueling, as in the condition of the XPR regime proposed in [31]. However, neither X-point radiation nor a cliff were not observed. Since penetration of the radiation zone inside the X-point was not observed in the Hmode D 2 ramp experiment, it seems that the 5 eV front did not reach sufficiently close to the X-point. In H-mode, achieving further detachment with only D 2 and sputtered carbon is difficult due to the much higher heating power, background plasma density, and temperature compared to L-mode. To explore this, a detachment experiment with seeded impurity (N 2 , Ne) was performed, but the target flux cliff was not observed and only a smooth rollover was found [43]. These experimental observations are consistent with the H-mode phase space presented by the SOLPS-ITER simulations. In the case of including seeded impurity, considering the shift of the phase space curve to the lower T et direction by sputtered carbon in H-mode ( figure 18) and the effect of decreasing n e,sep for the increase in seeded impurity concentration [2], the phase space curve is expected to shift to the lower left. Therefore, it is expected that T et will fall below 1 eV before the Greenwald fraction rises excessively, allowing conditions for cliff and XPR to be reached. However, the parallel gradient of the electron temperature in the SOL becomes less stiff in the presence of strong SOL radiation, so the T et threshold for the cliff transition can be shifted in the phase space. In addition, since the main radiation zone of seeded impurities such as Ne and N 2 is formed at a higher temperature range than the carbon radiation zone and the ionization front of the main ion, a stable XPR may occur but a cliff may not be observed due to this mismatch.

Comparison with other devices and cliff event types
The divertor plasma parameter cliffs observed in other devices are summarized here and compared with those observed in this paper.
The bifurcation-like electron temperature drop observed in DIII-D [9] has been attributed to the role of E × B drift since it has been only observed at the low-field side target SOL in case of H-mode plasmas with a forward field configuration, where the ion B × ∇B is directed toward the divertor. The drift enhances cross-field transport from the near-SOL to PFR in the outer SOL. The positive feedback by the relation of E × B drift and local density accelerates the transition to the detached regime. Since the target flux cliff shown in KSTAR plasmas occurred simultaneously for both of the targets and regardless of PFR, near-and far-SOL (only excepting inner far-SOL), the cross-field drift would appear not be the main driver of the cliff observed in KSTAR.
The divertor self-oscillation observed in JET Mark I divertor [6,12] was attributed to the coupling/decoupling of the sputtered carbon sources with the radiated power. However, the oscillation of the n e and LP signals observed in feedback-controlled density ramp in KSTAR were due to the modulation of fueling throughput. Since the oscillations are not observed in the discharge without fueling modulation (#20267), the change of the sputtered carbon source cannot be a key driver of target flux cliff described in this paper.
There is an experimental observation in ASDEX Upgrade (with all tungsten PFCs) in which the radial decay length of the upstream density is increased as the divertor plasma detachs [44]. This implies that the effective perpendicular diffusivity is increased when detachment proceeds. If the effective perpendicular diffusivity, D ⊥,eff increased during the cliff, the sharp decrease in target flux in KSTAR can be understood [15]. But this should be accompanied by the depletion of the core density and increased upstream SOL density, so it is not consistent with the increased dn e /dt observed in KSTAR target flux cliff. Furthermore, the EDGE2D simulation shown in [45] indicates that an abrupt transition to the detached state was obtained when there is no correlation between D ⊥,eff and the edge collisionality. Therefore, the role of effective diffusivity on the bifurcation should be clarified with the support of the experimental upstream profile as a future work.
The X-point radiation or MARFE does not always accompany the abrupt change of the divertor condition, as shown in JT60-U [46] and ASDEX Upgrade [1]. Therefore, species or geometry dependence cannot be ruled out.

Leveraging time-dependent gas throughput scans and phase space diagrams for system identification
The QSS phase space diagram can be an important tool to determine whether or not a cliff might be expected to appear in a given system. Since the divertor-SOL system is very nonlinear, predicting system behavior through extrapolation is quite challenging. Therefore, the method of scanning the operation space through high fidelity simulation such as SOLPS-ITER has been mainly used. However, it is difficult to determine the dynamic behavior with a series of steady-state density scans that provide discrete points in phase space. Instead, timedependent gas throughput scans for a long enough time compared to the relaxation timescale gives continuous curve rather than the discrete points, providing more insight into existence of the discontinuity or abrupt changes in the operation space. In addition, time-dependent gas throughput scans can be utilized for system identification and for training reduced models [35,47].

Conclusions
In the KSTAR L-mode density ramp experiment, whenn e reached a specific critical density after target flux rollover, a cliff transition in which the target flux abruptly decreased or increased during a timescale of 20-30 ms was observed. The cliff transition occurred simultaneously across the entire inner and outer targets, regardless of gas input modulation. Since the radiation zone measured by the imaging bolometer system moved toward the core beyond the X-point while the target flux cliff was in progress, it can be seen that the cliff has a strong correlation with the SOL power depletion by the Xpoint radiation. In addition, the ionization front movement can be estimated through the rapid changes in dn e /dt.
Time-dependent SOLPS-ITER simulations qualitatively reproduces the KSTAR L-mode target flux cliff transition. The timescale was in good agreement with the experiment when the dynamic mode was used. As in the experiment, the state in the middle of the cliff was unstable, so a steady-state solution does not exist, and the transition can occur immediately above or below the cliff. Due to the characteristics of a carbon device, the main radiators are C 1+ , C 2+ , C 3+ , and their main radiation temperature coincides with 5 eV, where ionization of the main ion principally occurs. Therefore, the 5 eV front, ionization front, and radiation zone are almost identical, and as the detachment progresses, the 5 eV front rises to the upstream SOL and undergoes perpendicular transport to the core through the X-point. Therefore, as the radiation power steadily increases in the core region just above the X-point and the ionization source is supplied, the SOL is further cooled and the upstream density increases. The target ion flux drop is the result of a combination of upstream pressure drop and volumetric momentum dissipation. Due to the strong correlation between the T et and 5 eV front, T et can be used as a cliff discriminant parameter, and, in H-mode, it is difficult to meet the target flux cliff condition because the density limit occurs before T et becomes sufficiently low.
The practical and operational importance of this work lies in understanding how to avoid or achieve smooth detachment in the operation space by avoiding the cliff region. Since the target ion flux cliff occurs when the radiation front and ionization front simultaneously move beyond the X-point, it is speculated that the cliff can be avoided if impurity species other than carbon become the dominant radiator, and the detachment is driven by radiation rather than by plasma-neutral interactions. Shifting the detachment driver from fueling to impurity seeding could potentially circumvent the target flux cliff. This hypothesis will be tested in future work by creating an XPR condition that has not yet been achieved in KSTAR Hmode and investigating whether cliff transitions are accompanied in this case. An impurity seeding experiment under high background density conditions will also be conducted. Further understanding of the seed impurity dependence of the H-mode phase space curves in the time-dependent SOLPS-ITER simulations is required to investigate whether the cliff transition condition can be achieved in H-mode. Additionally, since the current KSTAR has an open divertor geometry, which can cause simultaneous cliff transitions in the entire inner and outer target, it is necessary to analyze the effect of the geometry on the divertor parameter cliff transition by comparing it with the operation space of the corner-shaped divertor geometry of KSTAR after upgrading the divertor.