Development and Application of SONIC Divertor Simulation Code to Power Exhaust Design of Japanese DEMO Divertor

: An integrated divertor simulation code, SONIC, has been developed in order to predict a self-consistent transport solution of the plasma, neutral and impurities in the scrape-off layer (SOL) and divertor. SONIC code has contributed to determining the divertor design and power handling scenarios for the Japanese (JA) fusion demonstration (DEMO) reactor. Radiative cooling scenario of Ar impurity seeding and the divertor performance have been demonstrated to evaluate the power exhaust scenarios with P sep = 230–290 MW. The simulation identiﬁed the decay length of the total parallel heat ﬂux proﬁle as being broader than the electron one, because of the ion convective transport from the outer divertor to the upstream SOL, produced by the plasma ﬂow reversal. The ﬂow reversal also reduced the impurity retention in the outer divertor, which may produce the partial detachment. Divertor operation margin of key power exhaust parameters to satisfy the peak q target ≤ 10 MWm − 2 was determined in the low n esep of 2 − 3 × 10 19 m − 3 under severe conditions such as reducing radiation loss fraction, i.e., f* raddiv = ( P radsol + P raddiv )/ P sep and diffusion coefﬁcients ( χ and D) . The divertor geometry and reference parameters ( f* raddiv ~ 0.8, χ = 1 m 2 s − 1 , D = 0.3 m 2 s − 1 ) were consistent with the low n esep operation of the JA DEMO concepts. For either severe assumption of f* raddiv ~0.7 or χ and D to their half values, higher n esep operation was required. In addition, recent investigations of physics models (temperature-gradient force on impurity, photon transport, neutral–neutral collision) under the DEMO relevant SOL and divertor condition are presented.


Introduction: SONIC Code and Application to Japanese DEMO Designs
An integrated divertor simulation code, SONIC (SOLDOR, NEUT2D and IMPMC codes), has been developed over almost two decades by QST (Fusion energy directorate was united from JAEA: Japan Atomic Energy Agency) [1,2], in order to predict a selfconsistent transport solution of the plasma, neutral and impurities in the scrape-off layer (SOL) and divertor, where magnetic field lines are attached to the divertor target plates at both ends. SONIC code consists of a two-dimensional multi-fluid model for the SOL and divertor plasma (SOLDOR), and neutral and impurity transport codes using Monte-Carlo (MC) techniques (NEUT2D and IMPMC, respectively). Three codes are calculated selfconsistently to achieve a steady-state solution. Fundamental equations and modeling of the SONIC code are shown and explained in Reference [2]. Transport equations of plasma particles, momentum, and ion and electron energies in SOLDOR are basically identical to B2-code [3]. Transport of neutral and molecular particles is traced inside the whole vessel by NEUT2D, and distributions of their densities are provided to calculate source and

Divertor Simulation and Geometry for JA DEMO
SONIC calculation meshes of all plasma areas and enlarged divertor areas are shown in Figure 2. Plasma transport in the SOL and divertor is simulated in the flux surfaces of the outer midplane SOL radius (r mid )  3.2 cm, which covers the most SOL area, connecting between the inner and outer targets. Here, the plasma transport is not calculated in the most outer mesh area, connecting to the vessel wall, except in the divertor target. The SOL area is divided into 24 radial layers along the magnetic flux surfaces, and the radial width at the outer midplane (r mid ) is increased from 0.7 mm (near the separatrix) to 1.3 mm (at the most outer flux surface). The inner and outer divertors cover all divertor plasma volume below the X-point, similar to ITER divertor, and the leg length (Ldiv) of 1.6 m is longer than that of ITER (1.0 m). The compression of neutral particles and efficient formation of the plasma detachment will be expected particularly near the strike-point. The poloidal angle between the separatrix and target surface at the strike point ( div ) is designed as 30 and 25 for the inner and outer targets, respectively, where the total (magnetic and geometrical) flux expansion at the inner and outer targets, i.e., fexp div /sin div where fexp div = (Bp/Bt) mid /(Bp/Bt) div , is similar (~12). The reflector is installed both at the inner and outer divertors, with the poloidal angle of 60° for the divertor design and operation window studies [31,38]. It is noted that the inner and outer reflector angles were recently increased to 90° and 80°, respectively, i.e., the V-shape corner geometry became more open. At the same time, they were extended behind the dome roof in order to protect the end-part of the cooling pipes against fast neuron irradiation from the main plasma. Influence of the "open reflector" geometry on the plasma detachment was small since the strike-point was relatively far (25 cm) from the V-shape corner.

Divertor Simulation and Geometry for JA DEMO
SONIC calculation meshes of all plasma areas and enlarged divertor areas are shown in Figure 2. Plasma transport in the SOL and divertor is simulated in the flux surfaces of the outer midplane SOL radius (r mid ) ≤ 3.2 cm, which covers the most SOL area, connecting between the inner and outer targets. Here, the plasma transport is not calculated in the most outer mesh area, connecting to the vessel wall, except in the divertor target. The SOL area is divided into 24 radial layers along the magnetic flux surfaces, and the radial width at the outer midplane (∆r mid ) is increased from 0.7 mm (near the separatrix) to 1.3 mm (at the most outer flux surface). The inner and outer divertors cover all divertor plasma volume below the X-point, similar to ITER divertor, and the leg length (L div ) of 1.6 m is longer than that of ITER (1.0 m). The compression of neutral particles and efficient formation of the plasma detachment will be expected particularly near the strike-point. The poloidal angle between the separatrix and target surface at the strike point (θ div ) is designed as 30 • and 25 • for the inner and outer targets, respectively, where the total (magnetic and geometrical) flux expansion at the inner and outer targets, i.e., f exp div /sinθ div where f exp div = (B p /B t ) mid /(B p /B t ) div , is similar (~12). The reflector is installed both at the inner and outer divertors, with the poloidal angle of 60 • for the divertor design and operation window studies [31,38]. It is noted that the inner and outer reflector angles were recently increased to 90 • and 80 • , respectively, i.e., the V-shape corner geometry became more open. At the same time, they were extended behind the dome roof in order to protect the end-part of the cooling pipes against fast neuron irradiation from the main plasma. Influence of the "open reflector" geometry on the plasma detachment was small since the strike-point was relatively far (25 cm) from the V-shape corner.  (b) Enlarged divertor region. Exhausted power (Pout) and particle (out) are given at the core-edge boundary (r mid /ap = 0.95). Ar is injected at the outer divertor above the SOL plasma region. [38] (c) Ar exhaust fluxes at the inner and outer divertor slots (ex in and ex out ) and Ar backflow fluxes (bk in and bk out ) are calculated once in every 10 IMPMC runs. In the following 9 IMPMC runs,bk in and bk out are provided proportional to puff Ar in "backflow model".
Exhaust plasma power (Pout) and particles (out) from the core region are given at the core-edge boundary (the magnetic surface corresponding to r mid /ap = 0.95). Pout is equally distributed to ion and electron energy sources (Pout i = Pout e = Pout/2). Psep is slightly smaller than Pout since the radiation loss in the main plasma edge (Prad edge ) is small in the JA DEMO concepts, i.e., 4-7% of Pout as shown in Section 3. Cross field diffusion coefficients of ion and electron energies and plasma particle are given over all plasma simulation area as χi = χe = 1 m 2 s −1 and D = 0.3 m 2 s −1 , respectively, which are the same as "standard" values in the previous ITER simulation [39]. Diffusion coefficient for the seeding impurity is given as the same value, i.e., Dimp = 0.3 m 2 s −1 .
For the parallel energy transport terms in the SOL and the divertor, electron conduction dominates when compared to ion conduction and their convections. In the multi-fluid plasma transport code (SOLDOR), heat flux limiters for the parallel conduction (a factor limiting the conduction energy transport to a fraction of the free streaming energy) are 0.2, 10 and 0.5 for electron, ion and viscosity, respectively. In the MC transport codes, NEUT2D and IMPMC, deuterium atoms (D0) and gas (D2) and impurity atoms (Ar) are traced in all mesh areas, as shown in Figure 2, including the most outer mesh in the vessel wall and above the dome, and in the sub-divertor. D2 is injected near the device outer midplane and Ar impurity seeding is at the upper part of the outer target. They are exhausted from the sub-divertor to the outboard exhaust opening, instead of the bottom in the early design [31]. This is because of reducing the fast neutron irradiation to protect the bottom part of the vacuum vessel. Distribution of the detachment plasma was not influenced, while the pumping system should be increased due to a reduction in the conductance of the exhaust route.

MC Calculation Technique in Divertor
In the SONIC simulation, computer resources are consumed primarily for impurity transport calculations (IMPMC). In order to save the impurity atom (Ar) transport (b) Enlarged divertor region. Exhausted power (P out ) and particle (Γ out ) are given at the core-edge boundary (r mid /a p = 0.95). Ar is injected at the outer divertor above the SOL plasma region. [38] (c) Ar exhaust fluxes at the inner and outer divertor slots (Γ ex in and Γ ex out ) and Ar backflow fluxes (Γ bk in and Γ bk out ) are calculated once in every 10 IMPMC runs. In the following 9 IMPMC runs, Γ bk in and Γ bk out are provided proportional to Γ puff Ar in "backflow model".
Exhaust plasma power (P out ) and particles (Γ out ) from the core region are given at the core-edge boundary (the magnetic surface corresponding to r mid /a p = 0.95). P out is equally distributed to ion and electron energy sources (P out i = P out e = P out /2). P sep is slightly smaller than P out since the radiation loss in the main plasma edge (P rad edge ) is small in the JA DEMO concepts, i.e., 4-7% of P out as shown in Section 3. Cross field diffusion coefficients of ion and electron energies and plasma particle are given over all plasma simulation area as χ i = χ e = 1 m 2 s −1 and D = 0.3 m 2 s −1 , respectively, which are the same as "standard" values in the previous ITER simulation [39]. Diffusion coefficient for the seeding impurity is given as the same value, i.e., D imp = 0.3 m 2 s −1 .
For the parallel energy transport terms in the SOL and the divertor, electron conduction dominates when compared to ion conduction and their convections. In the multi-fluid plasma transport code (SOLDOR), heat flux limiters for the parallel conduction (a factor limiting the conduction energy transport to a fraction of the free streaming energy) are 0.2, 10 and 0.5 for electron, ion and viscosity, respectively. In the MC transport codes, NEUT2D and IMPMC, deuterium atoms (D 0 ) and gas (D 2 ) and impurity atoms (Ar) are traced in all mesh areas, as shown in Figure 2, including the most outer mesh in the vessel wall and above the dome, and in the sub-divertor. D 2 is injected near the device outer midplane and Ar impurity seeding is at the upper part of the outer target. They are exhausted from the sub-divertor to the outboard exhaust opening, instead of the bottom in the early design [31]. This is because of reducing the fast neutron irradiation to protect the bottom part of the vacuum vessel. Distribution of the detachment plasma was not influenced, while the pumping system should be increased due to a reduction in the conductance of the exhaust route.

MC Calculation Technique in Divertor
In the SONIC simulation, computer resources are consumed primarily for impurity transport calculations (IMPMC). In order to save the impurity atom (Ar) transport calculation in the sub-divertor, a technique to assume the impurity backflow at the exhaust slots, i.e., "backflow model", was developed [18] and it has been incorporated in the IMPMC calculation. Ar atoms in the sub-divertor region are executed 1 in every 10 IMPMC runs, and the exhaust fluxes at the inner and outer divertor slots (Γ ex in and Γ ex out ) and backflow fluxes (Γ bk in and Γ bk out ) were calculated as shown in Figure 2c. In the following nine IMPMC runs during the SONIC conversion process, Γ bk in and Γ bk out are provided proportionally when the seeding rate (Γ puff Ar ) changes, where Γ ex in and Γ ex out are absorbed at the divertor slots. This "backflow model" technique can also reduce oscillation to a steady-state solution, which is mainly caused by the different time scales of impurity transport between the divertor chamber and the sub-divertor. As a result, the total IMPMC calculation time is reduced significantly, and a stable solution of the SONIC simulation can be obtained. to satisfy the Ar particle balance of the net exhaust flux Development of neutral-neutral (and neutral-molecular, molecular-molecular) elastic collision models (NNC) in NEUT2D, and influences on the plasma detachment and the divertor pressure, are described in Section 4.3. On the other hand, NNC was not activated for the divertor performance study.

Power Exhaust Scenario and Key Parameters
The JA DEMO concepts challenge plasma performance, being higher than ITER levels, with comparable f rad main (0.2-0.4), and a large divertor power handling factor of P sep /R p = 28-34 MWm −1 , corresponding to 1.8-2.1 times larger than that of ITER (16 MWm −1 ). Divertor volume is increased, i.e., R p and L div are 1.37 and 1.6 times of ITER (6.2 m and 1.0 m), respectively, which will contribute to increasing the radiation power and to reducing the plasma temperature.
A simple formula of the plasma heat load (q plasma ) is described by P sep , f* rad div , λ q// mid and geometry parameters, as following, q plasma = P sep ·(1-f* rad div )·(sinθ div /f exp div )· 1/(4πR p ·λ q// mid ·S det ), where f* rad div = (P rad sol +P rad div )/P sep , i.e., radiation power fraction in the SOL and divertor normalized by P sep , and power reduction in the plasma detachment is represented by an additional factor (S det ) in this simple formula. The total flux expansion is 12 as explained in Section 2.2. λ q// mid is e-folding length of the parallel heat flux (q // ) profile near the outer midplane separatrix, which is determined by the power balance between the plasma radial diffusion and parallel conduction. SONIC simulations evaluated profiles of the plasma density, temperature and heat load of the target as a self-consistent solution of the plasma, neutral and impurity transport processes. Particularly for the reactor design, the plasma detachment is required to significantly reduce both the plasma temperatures (T e div , T i div ) near the strike-points and the peak heat load (q target ) less than 10 MWm −2 level [40]. P sep , f* rad div and diffusion coefficients are important key parameters for the divertor performance, which have been systematically investigated by the SONIC code [38]. When compared to that of ITER (0.5-0.6) [24], larger f* rad div is required for the larger P sep /R p design.
High plasma density is preferable for the reactor design and power exhaust. On the other hand, restriction of n e at the main plasma separatrix (n e sep ) is a common issue for DEMO divertor designs [41,42], since Greenwald density (n GW = I p /πa p 2 [10 20 m −3 , MA, m]) corresponds to 6.7 x10 19 and 7.3 × 10 19 m −3 for the JA DEMO 2014 and higher-κ designs, respectively, both of which are lower than n GW of ITER (1.1 × 10 20 m −3 ). Thus, n e main is required to be 1.2 times higher than n GW , assuming a quadratic function-like peak profile, in order to achieve the expecting P fusion . At the same time, n e sep /n GW was reported to be~1/3 in the H-mode plasma experiments [43,44] and similar fraction in "standard" ITER simulations [39]. A recent experiment result was also lower than critical values of 0.4-0.5 [45], which was consistent with the edge ballooning stability model. As a result, the operation range of n e sep /n GW is expected to be 0.3-0.5, thus the divertor performance and the operation boundary for the JA DEMO were investigated in the low density range of n e sep = 2 -3 × 10 19 m −3 . SOL density (n e sep ) scan for the JA DEMO higher-κ case, corresponding to the lower P sep , was performed with increasing fuel gas puff rate from 25 to 150 Pam 3 s −1 for a given P out = 250 MW and Γ out = 1 × 10 22 Ds −1 at the core-edge boundary [31]. The total radiation power in the plasma edge (0.95 < r/a < 1), SOL and divertor (P rad edge + P rad sol + P rad div ) was adjusted to be 200 MW by feedback of the Ar seeding rate in the iterative SONIC calculation. The input radiation fraction corresponded to f rad in = (P rad edge + P rad sol + P rad div )/P out = 0.8, which was so far determined as the reference value, and this series of n e sep scan was referred to as Case-1 in Reference [38]. Here, P sep (= 235-240 MW) is slightly smaller than P out , thus f* rad div corresponds to a slightly smaller value (~0.78). Power exhaust parameters of the reference JA DEMO 2014 scenario are P out = 300 MW and P rad edge + P rad sol + P rad div = 240 MW (Case-2). In addition, lower radiation fraction cases of the JA DEMO higher-κ and JA DEMO 2014 (f rad inp = 0.7 for each P out ) were also investigated as Case-3 and Case-4, respectively. Lower limit of the divertor plasma operation is discussed in Section 3.4. Figure 3a,b show distributions of the radiation power density (W rad ) in the inner and outer divertors at n e sep = 2.0 × 10 19 m −3 for Case-1, corresponding to a near lower boundary of the n e sep range. Total radiation powers in the two divertors are comparable, i.e., 79 and 82 MW, respectively. The large radiation peaks near the separatrix are maintained at the upstream on both divertor legs. In the inner divertor, a large W rad is seen at 40-60 cm upstream of the target near the separatrix, and it is maintained far above the inner target. T e is decreased to~1 eV over most of the area of the target, which we describe as "full detachment", as shown in Figure 3c. Total heat load (q target ) is evaluated by including surface recombination of the ions (q t rec = n i div C s div E ion , where n i div , C s div and E ion are ion density, sound velocity at the divertor sheath and recombination energy, respectively), radiation power load (q t rad ) and neutral flux load (q t n ), in addition to the plasma heat flux (q t plasma ). The peak q target of 4.2 MWm −2 is seen near the separatrix, mostly attributed to q t rec as shown in Figure 3e. Here, locations of multi-peaks in the q target profile near the strike-point change in the converged solutions, which are attributed to peaks of n e div , n i div , T e div and T i div profiles. The time-averaged profile may be more appropriate to describe the peak q target value. operation boundary for the JA DEMO were investigated in the low density range of ne sep = 2 -3 × 10 19 m −3 . SOL density (ne sep ) scan for the JA DEMO higher- case, corresponding to the lower Psep, was performed with increasing fuel gas puff rate from 25 to 150 Pam 3 s −1 for a given Pout = 250 MW and out = 1 × 10 22 Ds −1 at the core-edge boundary [31]. The total radiation power in the plasma edge (0.95 < r/a < 1), SOL and divertor (Prad edge + Prad sol + Prad div ) was adjusted to be 200 MW by feedback of the Ar seeding rate in the iterative SONIC calculation. The input radiation fraction corresponded to frad in = (Prad edge + Prad sol + Prad div )/Pout = 0.8, which was so far determined as the reference value, and this series of ne sep scan was referred to as Case-1 in Reference [38]. Here, Psep (= 235-240 MW) is slightly smaller than Pout, thus f*rad div corresponds to a slightly smaller value (~0.78). Power exhaust parameters of the reference JA DEMO 2014 scenario are Pout = 300 MW and Prad edge + Prad sol + Prad div = 240 MW (Case-2). In addition, lower radiation fraction cases of the JA DEMO higher- and JA DEMO 2014 (frad inp = 0.7 for each Pout) were also investigated as Case-3 and Case-4, respectively. Lower limit of the divertor plasma operation is discussed in Section 3.4. Figure 3a,b show distributions of the radiation power density (Wrad) in the inner and outer divertors at ne sep = 2.0 × 10 19 m −3 for Case-1, corresponding to a near lower boundary of the ne sep range. Total radiation powers in the two divertors are comparable, i.e., 79 and 82 MW, respectively. The large radiation peaks near the separatrix are maintained at the upstream on both divertor legs. In the inner divertor, a large Wrad is seen at 40-60 cm upstream of the target near the separatrix, and it is maintained far above the inner target. Te is decreased to ~1 eV over most of the area of the target, which we describe as "full detachment", as shown in Figure 3c. Total heat load (qtarget) is evaluated by including surface recombination of the ions (qt rec = ni div Cs div Eion, where ni div , Cs div and Eion are ion density, sound velocity at the divertor sheath and recombination energy, respectively), radiation power load (qt rad ) and neutral flux load (qt n ), in addition to the plasma heat flux (qt plasma ). The peak qtarget of 4.2 MWm −2 is seen near the separatrix, mostly attributed to qt rec as shown in Figure 3e. Here, locations of multi-peaks in the qtarget profile near the strike-point change in the converged solutions, which are attributed to peaks of ne div , ni div , Te div and Ti div profiles. The time-averaged profile may be more appropriate to describe the peak qtarget value. Profiles of (c,d) ne div , Te div and Ti div , (e,f) integrating heat load components at the inner and outer divertor target, respectively [38]. Profiles of (c,d) n e div , T e div and T i div , (e,f) integrating heat load components at the inner and outer divertor target, respectively [38].

Radiation Loss and Plasma Detachment
In the outer divertor, a large W rad is also seen at the upstream (40-60 cm) near the separatrix, where local c Ar (= n Ar /n i ) is increased up to 2%. On the other hand, the W rad peak shifts toward the target at the outer flux surfaces as shown in Figure 3b, and it becomes smaller than 10 MWm −3 (lowest color bar) and located just above the target (a few centimeters). The plasma detachment is produced within~12 cm near the strikepoint as shown in Figure 3d, which we describe as "partial detachment". The peak q target of 5.5 MWm −2 is seen at the boundary of the attached region, where both T e div and T i div are increased from 1-2 to~15 eV, and n e div is decreased from~1.5 × 10 22 m −3 tõ 1.5 × 10 21 m −3 . At the same time, the W rad peak shifts toward the target, thus q t plasma and q t rad become dominant. In the partial detachment, the peak q target is sensitive to the plasma temperature and density profiles and location of the W rad peak. Figure 4 shows distributions of T e , W rad , c Ar , Mach number of the plasma flow (M // = V i// /C s , where V i// is parallel ion velocity and C s = [(T e +T i )/m i ] 0.5 ), the static plasma pressure (n e T e + n i T i ) and the total plasma pressure, i.e., static and dynamic pressures (n e T e + n i T i + m i n i V i// 2 ), along three representative magnetic field lines, corresponding to the second, sixth and ninth radial meshes outside the separatrix. The plasma detachment is produced at the former two target areas, and the peak q target is seen at the target in the third one (attached plasma). Along the magnetic field line near the separatrix (second radial mesh), the plasma detachment (T e and the total plasma pressure decrease to 1~2 eV and a half value of the upstream total pressure, respectively) is produced at a few 10 cm above the target as shown in Figure 4a,f. As described above, the peak q target is sensitive to the plasma profiles both across and along the magnetic field line. In order to simulate an abrupt change of the plasma parameters and impurity concentration near the divertor target, the mesh size is gradually reduced towards the target as shown in Figure 2b, and the minimum size (in a few cm above the target) is 1.5-3 mm in the poloidal direction (corresponding to 28-55 mm along the magnetic field line). Along the magnetic field line at the attached region (ninth radial mesh), T e (and T i ) decreased from 25 eV (15 mm above the target) to 15 eV (at the target) due to enhancement of W rad above the target. It is noted that the plasma flow changes toward the target in the narrow region, M // is enhanced tõ 0.82 at the target as shown in Figure 4d. While the static plasma pressure is decreased from 2.1 kPa to 1.5 kPa, the total plasma pressure is conserved along the magnetic field line as shown in Figure 4e,f. In the outer divertor, a large Wrad is also seen at the upstream (40-60 cm) near the separatrix, where local cAr (= nAr/ni) is increased up to 2%. On the other hand, the Wrad peak shifts toward the target at the outer flux surfaces as shown in Figure 3b, and it becomes smaller than 10 MWm −3 (lowest color bar) and located just above the target (a few centimeters). The plasma detachment is produced within ~12 cm near the strike-point as shown in Figure 3d, which we describe as "partial detachment". The peak qtarget of 5.5 MWm −2 is seen at the boundary of the attached region, where both Te div and Ti div are increased from 1-2 to ~15 eV, and ne div is decreased from ~1.5 × 10 22 m −3 to ~1.5 × 10 21 m −3 . At the same time, the Wrad peak shifts toward the target, thus qt plasma and qt rad become dominant. In the partial detachment, the peak qtarget is sensitive to the plasma temperature and density profiles and location of the Wrad peak. Figure 4 shows distributions of Te, Wrad, cAr, Mach number of the plasma flow (M// = Vi///Cs, where Vi// is parallel ion velocity and Cs = [(Te+Ti)/mi] 0.5 ), the static plasma pressure (neTe + niTi) and the total plasma pressure, i.e., static and dynamic pressures (neTe + niTi + miniVi// 2 ), along three representative magnetic field lines, corresponding to the second, sixth and ninth radial meshes outside the separatrix. The plasma detachment is produced at the former two target areas, and the peak qtarget is seen at the target in the third one (attached plasma). Along the magnetic field line near the separatrix (second radial mesh), the plasma detachment (Te and the total plasma pressure decrease to 1~2 eV and a half value of the upstream total pressure, respectively) is produced at a few 10 cm above the target as shown in Figure 4a,f. As described above, the peak qtarget is sensitive to the plasma profiles both across and along the magnetic field line. In order to simulate an abrupt change of the plasma parameters and impurity concentration near the divertor target, the mesh size is gradually reduced towards the target as shown in Figure 2b, and the minimum size (in a few cm above the target) is 1.5-3 mm in the poloidal direction (corresponding to 28-55 mm along the magnetic field line). Along the magnetic field line at the attached region (ninth radial mesh), Te (and Ti) decreased from 25 eV (15 mm above the target) to 15 eV (at the target) due to enhancement of Wrad above the target. It is noted that the plasma flow changes toward the target in the narrow region, M// is enhanced to ~0.82 at the target as shown in Figure 4d. While the static plasma pressure is decreased from 2.1 kPa to 1.5 kPa, the total plasma pressure is conserved along the magnetic field line as shown in Figure 4e,f.  Profiles of T e div and T i div in the attached region are relatively flat or rather increased toward the outer SOL, but they are still lower than those at the upstream SOL, as shown in Figure 5a, where the radial distances are mapped to the same magnetic surfaces at the outer midplane. The midplane SOL width of the plasma detachment and location of the peak q target correspond to ∆r mid~1 and 1.25 cm, respectively. Near the separatrix, T e and T i are significantly reduced from 380 and 812 eV at the X-point, respectively, to 1-2 eV at the target. On the other hand, Figure 5b shows that reduction in the total plasma pressure is only half in the simulation. Here, the dynamic pressure near the X-point is very small, while it is added to the static pressure, independently of the plasma flow direction. The plasma flow profile in the outer divertor and SOL is explained in Section 3.3. A significant reduction in n e div , i.e., "roll-over", was expected, at least, by further reduction of the plasma temperature to less than~0.5 eV, leading to enhancement of volume recombination process. Studies of atomic and molecular collision processes and/or enhancement of the plasma diffusion may be necessary to produce the further reductions in T e div and T i div , leading to further reduction in the particle flux and n e div as was explained in Section 1 [13,23].
Triangle corresponds to each field line location near the X-point. Positive Mach number represents the flow direction toward the outer target. The plasma detachment is produced near the target in the former two field lines, and the peak qtarget is seen at the target (attached plasma) in the third one.
Profiles of Te div and Ti div in the attached region are relatively flat or rather increased toward the outer SOL, but they are still lower than those at the upstream SOL, as shown in Figure 5a, where the radial distances are mapped to the same magnetic surfaces at the outer midplane. The midplane SOL width of the plasma detachment and location of the peak qtarget correspond to r mid ~1 and 1.25 cm, respectively. Near the separatrix, Te and Ti are significantly reduced from 380 and 812 eV at the X-point, respectively, to 1-2 eV at the target. On the other hand, Figure 5b shows that reduction in the total plasma pressure is only half in the simulation. Here, the dynamic pressure near the X-point is very small, while it is added to the static pressure, independently of the plasma flow direction. The plasma flow profile in the outer divertor and SOL is explained in Section 3.3. A significant reduction in ne div , i.e., "roll-over", was expected, at least, by further reduction of the plasma temperature to less than ~0.5 eV, leading to enhancement of volume recombination process. Studies of atomic and molecular collision processes and/or enhancement of the plasma diffusion may be necessary to produce the further reductions in Te div and Ti div , leading to further reduction in the particle flux and ne div as was explained in Section 1 [13,23]. Figure 5. Profiles of (a) Te and Ti at the outer divertor target and near the X-point, (b) plasma static pressure and total plasma pressure including the dynamic one at the outer divertor and near the Xpoint, (c) electron and total parallel heat fluxes near the X-point (q//e Xp , q//e Xp + q//i Xp ) and the heat load at the outer target. Distances from the separatrix at the outer target and near the X-point are mapped to the outer midplane SOL radius.
As a result, SONIC simulation showed that the peak qtarget for the reference power exhaust scenario (Case-1) can be reduced to lower than 10 MWm −2 for both inner and outer divertors even at the low ne sep , which is expected near the lower boundary of the JA DEMO operation window. At the same time, large Wrad and plasma detachment can be maintained in the long divertor leg without extending to the X-point, SOL and inside the main plasma.

Influence of Heat Flux Profile on Plasma Detachment
Figure 5c compares profiles of the heat load at the outer target (qt plasma , qt plasma + qt rec , total qtarget) and the electron and total parallel heat fluxes near the X-point (q//e Xp and q//e Xp + q//i Xp ). Both q//e Xp and q//e Xp +q//i Xp profiles are described approximately by a two-exponential function such as q//(r mid ) = q// near exp(−r mid /q// near ) + q// far exp(−r mid /q// far ). The decay lengths of "near-SOL" (q// near ) and "far-SOL" (q// far ) for the q//e Xp profile, mostly attributed by the electron conduction, correspond to 2.2 and 9.1 mm, respectively. Both Te and Ti near the As a result, SONIC simulation showed that the peak q target for the reference power exhaust scenario (Case-1) can be reduced to lower than 10 MWm −2 for both inner and outer divertors even at the low n e sep , which is expected near the lower boundary of the JA DEMO operation window. At the same time, large W rad and plasma detachment can be maintained in the long divertor leg without extending to the X-point, SOL and inside the main plasma.

Influence of Heat Flux Profile on Plasma Detachment
Figure 5c compares profiles of the heat load at the outer target (q t plasma , q t plasma + q t rec , total q target ) and the electron and total parallel heat fluxes near the X-point (q //e Xp and q //e Xp + q //i Xp ). Both q //e Xp and q //e Xp +q //i Xp profiles are described approximately by a two-exponential function such as q // (∆r mid ) = q // near exp(−∆r mid /λ q// near ) + q // far exp(−∆r mid /λ q// far ). The decay lengths of "near-SOL" (λ q// near ) and "far-SOL" (λ q// far ) for the q //e Xp profile, mostly attributed by the electron conduction, correspond to 2.2 and 9.1 mm, respectively. Both T e and T i near the separatrix become high in the DEMO plasma due to large exhaust power and low density when compared to those of ITER. Thus, the q // profile in the SOL generally becomes narrow when compared to the ITER simulation result (3.6 mm [39]), whereas χ i = χ e = 1 m 2 s −1 and D = 0.3 m 2 s −1 are the same values. The SONIC simulation result shows that large q // near the separatrix (shown as "near-SOL" in Figure 5c) is mostly reduced in the divertor due to large radiation loss and detachment.
Location of the peak q target (∆r mid = 1.25 cm) corresponds to a boundary between the "near-SOL" and the "far-SOL" in the q //e Xp + q //i Xp profile, and the peak q target is determined mostly by the local T e div and T i div , i.e., plasma heat load, as shown in Section 3.2.
λ q// near and λ q// far of the q //e Xp + q //i Xp profile are 2.9 and 13.2 mm, respectively, which are larger than those of the q //e Xp profile. Transition from the "near-SOL" to the "far-SOL" of the q //e Xp + q //i Xp profile is mostly attributed to the ion conduction and convection transports. Particularly, the convection component in q //i Xp is driven from the outer divertor to the upstream SOL in the wide SOL region (∆r mid ≤ 1.7 cm), as shown in Figure 6a, and it does not contribute to the ion heat flux towards the outer divertor. Therefore, a flat "shoulder" appears in the q //i Xp profile as shown in Figure 6b. The small flow reversal (M // = −0.1~−0.2) is produced even above the divertor target as shown in Figure 4d, and the plasma static pressure is gradually decreased toward the upstream SOL. These simulation results suggests that it is produced above the divertor target by locally increasing neutral ionization and plasma pressure, and it is extended to the low-field-side SOL as proposed in Reference [46]. As a result, the decay length of the q //e Xp + q //i Xp profile is increased when compared to that of the q //e Xp profile, and, at the same time, the transition from the "near-SOL" to the "far-SOL" is emphasized clearly. At the same time, as shown in Figure 4c, c Ar is reduced significantly above the outer target and it is gradually increased toward the upstream SOL, in particular on the magnetic field lines, where the relatively small flow reversal (M // = −0.2) is generated. This may sustain the attached plasma region at the outer target and produce the partial detachment in the simulation. separatrix become high in the DEMO plasma due to large exhaust power and low density when compared to those of ITER. Thus, the q// profile in the SOL generally becomes narrow when compared to the ITER simulation result (3.6 mm [39]), whereas χi = χe = 1 m 2 s −1 and D = 0.3 m 2 s −1 are the same values. The SONIC simulation result shows that large q// near the separatrix (shown as "near-SOL" in Figure 5c) is mostly reduced in the divertor due to large radiation loss and detachment. Location of the peak qtarget (r mid = 1.25 cm) corresponds to a boundary between the "near-SOL" and the "far-SOL" in the q//e Xp + q//i Xp profile, and the peak qtarget is determined mostly by the local Te div and Ti div , i.e., plasma heat load, as shown in Section 3.2. q// near and q// far of the q//e Xp + q//i Xp profile are 2.9 and 13.2 mm, respectively, which are larger than those of the q//e Xp profile. Transition from the "near-SOL" to the "far-SOL" of the q//e Xp + q//i Xp profile is mostly attributed to the ion conduction and convection transports. Particularly, the convection component in q//i Xp is driven from the outer divertor to the upstream SOL in the wide SOL region (r mid  1.7 cm), as shown in Figure 6a, and it does not contribute to the ion heat flux towards the outer divertor. Therefore, a flat "shoulder" appears in the q//i Xp profile as shown in Figure 6b. The small flow reversal (M// = −0.1 ~ −0.2) is produced even above the divertor target as shown in Figure 4d, and the plasma static pressure is gradually decreased toward the upstream SOL. These simulation results suggests that it is produced above the divertor target by locally increasing neutral ionization and plasma pressure, and it is extended to the low-field-side SOL as proposed in Reference [46]. As a result, the decay length of the q//e Xp + q//i Xp profile is increased when compared to that of the q//e Xp profile, and, at the same time, the transition from the "near-SOL" to the "far-SOL" is emphasized clearly. At the same time, as shown in Figure 4c, cAr is reduced significantly above the outer target and it is gradually increased toward the upstream SOL, in particular on the magnetic field lines, where the relatively small flow reversal (M// = −0.2) is generated. This may sustain the attached plasma region at the outer target and produce the partial detachment in the simulation.  Radial diffusion of the SOL plasma significantly affects profiles of the heat and particle fluxes. Effects of the smaller diffusion coefficients on the divertor operation were investigated. Simulations with reducing both χ and D to half values, i.e., χ e = χ i = 0.5 m 2 s −1 , D = 0.15 m 2 s −1 and T e , T i and n e profiles were compared in Figure 7a of Reference [38]. For the smaller χ and D case, gas puff rate was reduced to 50 Pam 3 s −1 in order to compare profiles at similar n e sep = 2.0 × 10 19 m −3 , since n e sep is increased by the reduction of D. Radial gradients of T e , T i and n e profiles were increased, and n e was decreased in the whole SOL region. T e and T i at the outer midplane separatrix (T e sep and T i sep ) are increased from 360 to 390 eV, and from 820 to 1190 eV, respectively, where increases in T i sep are significant. Therefore, λ q// near values of q //e Xp and q //i Xp profiles become smaller (from 2.3 to 1.6 mm and from 3.2 to 2.5 mm, respectively) due to reduction of the radial diffusion and enhancement of the parallel conductive transport. Figure 6c,d show the Mach number and the heat flux profiles, respectively, for the smaller χ and D case. Peak of the flow reversal shifts toward the separatrix and the "shoulder" in the q //i Xp profile becomes narrow (∆r mid ≤ 1.1 cm), thus the transition from the near-SOL to the far-SOL shifts toward the separatrix. As a result, λ q// near of the q //e Xp + q //i Xp profile is decreased from 2.9 to 2.2 mm. At the same time, transition from the near-SOL to the far-SOL becomes rather gradual. Figure 6e compares λ q// near values of q //e Xp and q //e Xp + q //i Xp profiles for the four cases (larger P sep and/or smaller f* rad div as explained in Section 3.1). All λ q// near values for the q //e Xp +q //i Xp profile are 1.2-1.6 times larger than those for the q //e Xp profiles, and they are decreased from 2.5-3.3 mm (for reference χ and D) to 1.7-2.2 mm.
It was reported that λ q// near was scaled by λ q// Eich = 0.7·B t −0.77 ·q 95 1.05 ·P sep 0.09 (mm, T, MW) [47], based on an experiment database of the heat load profile at the outer target under the attached divertor condition. The scaling predicts 1.3 mm for JA DEMO, respectively. These λ q// near values are still larger than the empirical scaling (λ q// Eich ). Further reductions of χ and D will be required to simulate the divertor performance for the empirical SOL heat and particle flux profiles.
It is noted that the near-SOL region in the q //e Xp + q //i Xp profile becomes narrow and that local T e Xp and T i Xp at the transition to the far-SOL are increased. Influence of the upstream heat flux profile was seen in the partial detachment profile at the outer target; the width was decreased from 12 to 7 cm for Case-1 as shown in Figure 7 of Reference [38], the peak q target was also shifted from r div = 15 to 10 cm, and the peak was increased from 5.8 to 9.5 MWm −2 due to large increases in the local T e div and T i div .

Divertor Operation in the Low Density SOL
Simulation results of the divertor plasma detachment and heat flux profile for the JA DEMO high-κ reference (Case-1: P sep~2 35 MW and f* rad div~0 .78), and effects of the diffusion coefficients are discussed in Section 3.2 and 3.3, respectively. Divertor operation in the low n e sep range and influences of the key parameters such as P sep , f* rad div and λ q// mid on the peak q target were systematically investigated [38]. A higher P sep~2 83 MW case with the same f* rad div (Case-2) is another reference for JA DEMO 2014. Divertor operations for extreme cases with reducing f* rad div to~0.7 (Case-3 and Case-4), and all four cases with reducing χ and D were also determined.
Operation boundary of the low n e sep for the standard χ and D is summarized in Figure 7c, where the peak q target are plotted as the function of n e sep and the outer q target is generally larger than the inner q target . The power exhaust in the divertor, i.e., P sep ·(1f* rad div ), for the four cases corresponds to 50, 60, 75 and 90 MW. Circles and squares show series of reference JA DEMO higher-κ and 2014 cases, respectively, where the peak q target in Figure 3f is marked by an open circle. The peak q target is generally increased with increasing P sep ·(1 -f* rad div ), and with decreasing n e sep . The two series of f* rad div~0 .8 are acceptable in predictably low n e sep (≥ 2 × 10 19 m −3 ) to reduce q target ≤ 10 MWm −2 . The JA DEMO higher-κ case has advantages to provide enough operation margin to the recrystallization temperature of the tungsten target (~1200 • C). Triangles and diamonds show the series of lower f* rad div (~0.7) cases. The peak q target is increased with reducing width of the plasma detachment, as shown in Figure 7a,b: profiles of T e div , T i div , n e div and heat load components at n e sep = 2.0 × 10 19 m −3 for the JA DEMO higher-κ case with the lower f* rad div (the peak q target is also marked by an open circle in Figure 7c). The detachment width is decreased to 7 cm at the outer target, and the peak q target is seen at~3 cm outside of each detach-attach boundary. The peak q target is significantly increased to 12.5 MWm −2 due to the increase in q t plasma , while q t rad is slightly reduced, where both the local T e div and T i div increase to 33 eV, while the local n e div is similar at~2 × 10 20 m −3 . In the outer flux surfaces (r div > 20 cm), n e div is decreased to a few 10 19 m −3 and T e div and T i div are increased to 20-30 and 70-100 eV, respectively. T i div becomes larger than T e div by a factor of three to four due to reduction in the energy exchange rate. These changes of the plasma profiles and the peak q target are consistent with those of the upstream SOL plasma as shown in Section 3.3, i.e., the flow reversal region is decreased in ∆r mid ≤ 1.1 cm, λ q// near of q //e Xp and q //e Xp + q //i Xp profiles are decreased to 1.9 and 2.7 cm, respectively, and the transition from the near-SOL to the far-SOL of the total q // Xp profile shifts toward the separatrix. As for the two lower f* rad div cases, higher n e sep operations, i.e., n e sep ≥ 2.3 × 10 19  by an open circle in Figures 7c). The detachment width is decreased to 7 cm at the outer target, and the peak qtarget is seen at ~3 cm outside of each detach-attach boundary. The peak qtarget is significantly increased to 12.5 MWm −2 due to the increase in qt plasma , while qt rad is slightly reduced, where both the local Te div and Ti div increase to 33 eV, while the local ne div is similar at ~2 × 10 20 m −3 . In the outer flux surfaces (r div > 20 cm), ne div is decreased to a few 10 19 m −3 and Te div and Ti div are increased to 20-30 and 70-100 eV, respectively. Ti div becomes larger than Te div by a factor of three to four due to reduction in the energy exchange rate. These changes of the plasma profiles and the peak qtarget are consistent with those of the upstream SOL plasma as shown in Section 3.3, i.e., the flow reversal region is decreased in r mid  1.1 cm, q// near of q//e Xp and q//e Xp + q//i Xp profiles are decreased to 1.9 and 2.7 cm, respectively, and the transition from the near-SOL to the far-SOL of the total q// Xp profile shifts toward the separatrix. As for the two lower f*rad div cases, higher ne sep operations, i.e., ne sep  2.3 × 10 19 and 2.7 × 10 19 m −3 , are required.  Figure 7c for comparison. Similarly, the peak qtarget is increased with increasing Psep·(1 -f*rad div ) and they are reduced with increasing ne sep for each case. The lower boundary of ne sep for qtarget  10 MWm −2 was determined to be 2.0 × 10 19 and 2.3 × 10 19 m −3 for Case-1 and Case-2, respectively. Therefore, the high f*rad div (~0.78) cases are acceptable in the low ne sep operation, while higher ne sep is required for Case-2. On the other hand, for the two lower f*rad div (~0.67) cases (Case-3 and Case-4), the peak qtarget and ne sep are significantly increased, thus the divertor operation is difficult in the low ne sep range (2 − 3 × 10 19 m −3 ). As a result, the reduction in  and D to half values significantly affected the divertor power exhaust, because of a reduction of the detachment width and increases in Te div and Ti div at the attached region, rather than due to increase of the peak q//, as shown in Section 3.3. In order to produce the wider detachment plasma and to reduce the acceptable peak qtarget ( 10 MWm −2 ), a high f*rad div of a level of 0.8 was necessary for the JA DEMO 2014 and JA DEMO higher- designs.

Development of SONIC Modelings for DEMO Design
Some proposals for improving physics models in the SONIC code have been investigated and influences on the divertor performance were evaluated under the DEMO SOL Results of the peak q target for the smaller χ and D cases are plotted by open symbols in Figure 7d. Four guidelines show results with the standard χ and D cases in Figure 7c for comparison. Similarly, the peak q target is increased with increasing P sep ·(1 -f* rad div ) and they are reduced with increasing n e sep for each case. The lower boundary of n e sep for q target ≤ 10 MWm −2 was determined to be 2.0 × 10 19 and 2.3 × 10 19 m −3 for Case-1 and Case-2, respectively. Therefore, the high f* rad div (~0.78) cases are acceptable in the low n e sep operation, while higher n e sep is required for Case-2. On the other hand, for the two lower f* rad div (~0.67) cases (Case-3 and Case-4), the peak q target and n e sep are significantly increased, thus the divertor operation is difficult in the low n e sep range (2 − 3 × 10 19 m −3 ). As a result, the reduction in χ and D to half values significantly affected the divertor power exhaust, because of a reduction of the detachment width and increases in T e div and T i div at the attached region, rather than due to increase of the peak q // , as shown in Section 3.3. In order to produce the wider detachment plasma and to reduce the acceptable peak q target (≤10 MWm −2 ), a high f* rad div of a level of 0.8 was necessary for the JA DEMO 2014 and JA DEMO higher-κ designs.

Development of SONIC Modelings for DEMO Design
Some proposals for improving physics models in the SONIC code have been investigated and influences on the divertor performance were evaluated under the DEMO SOL and divertor conditions. Further validations of simulations such as benchmarks with other divertor codes will be required to be incorporated as defaults. First, an improvement of the kinetic simulation model for the thermal force on impurity ions was proposed for the low collisionality plasma under the SOL condition in DEMO, such as T e sep of a few 100 eV, T i sep of several 100 eV and n e sep of a few 10 19 m −3 . Second, a photon transport model for hydrogen isotopes was investigated under the detached plasma condition of the DEMO divertor, where both neutral density and divertor size are increased when compared with existing devices and ITER. Third, a neutral-neutral collision (NNC) model was implicated since neutral-neutral elastic scattering plays an important role for distributions of neutrals, molecules and detachment plasma.

Improvement of Kinetic Model in Impurity Transport
SONIC simulation results showed that the midplane c Ar sol for the Case-1 series was decreased from 1.3 to 0.4-0.5% with increasing the gas puff rate (and n e sep ), and that the divertor shielding factor (c Ar div /c Ar sol ) was increased from 1.4 to 2.4-2.6 [31]. At the same time, c Ar sol was varied also in the poloidal direction, e.g., c Ar sol was decreased from~1% (near X-point) to 0.4% (at outer midplane). Parallel transport of the seeding impurity ion in SOL and divertor is determined mostly by modeling of the thermal (temperature-gradient) force (TF) and friction force (FF) from the bulk ions. Conventional formula of a kinetic model [48], which is based on the collisional plasma and proportional to ∇T i along the magnetic field line, was incorporated in the IMPMC code, and it has been used under the DEMO SOL plasma condition, i.e., low collisionality plasma. A kinetic TF model was proposed [49,50], which incorporated a collisionality parameter and ion heat flux limiter in the conventional TF formula as a collisionality-dependent term. The model reduces TF significantly in the low collisionality SOL. Self-consistent plasma simulations were carried out for the JA DEMO higher-κ reference case. Figure 8a shows distributions of c Ar along the separatrix magnetic field line from the outer target to the inner one. The conventional TF model (A) produces Ar accumulation peaks in some poloidal locations near the X-point of the low-field-side (LFS) SOL and along the highfield-side (HFS) SOL. Effects of the proposed TF model on the c Ar profile were evaluated for (B), assuming the same profiles of plasma temperature, density and flow velocity for the conventional model (A); c Ar is significantly decreased over all SOL regions (above the X-point). On the other hand, a self-consistent SONIC solution of the proposed model (C) shows that TF is reduced by 20-80% in the SOL when compared to the conventional TF model (A), and that the local accumulation peaks of c Ar are reduced. While the averaged c Ar level is reduced slightly in the overall SOL region, it is comparable near the outer midplane. A similar effect was also observed in the poloidal profile of W rad as shown in Figure 8b. Here, the total heat load on the target was comparable since the total P rad sol and P rad div were adjusted to be comparable (~190 MW and f* rad div~0 .78) for (A) and (C) cases by feedback of the Ar seeding rate. In the overall SOL region, P rad sol is reduced by 18%, which corresponds to the reduction in the c Ar sol level. At the same time, Ar influx into the main plasma is also reduced. In the divertor, c Ar div is decreased due to increase in n e rather than n Ar , thus P rad div is slightly increased from 175 to 180 MW. As a result, the proposed TF model affected to reduce the TF on impurity transport and local accumulation peaks of c Ar div under the low collisionality SOL condition. It also reduced P rad sol , but the value was relatively small. Processes 2022, 10, x FOR PEER REVIEW 14 of 19

Development of Photon Transport Model and Its Effect of Plasma Detachment
Under the detached condition of the DEMO divertor, i.e. high neutral density and large divertor size, photon absorption of Lyman and Balmer line radiations by neutral atoms will enhance the ionization process and change the distribution of the detached plasma. Plasma density becomes significantly high, such as a few 10 21 m −3 , as shown in Figure 3. Typical optical pass length of Ly- is decreased to ~1 mm, and ionization process is enhanced locally by the photon absorption. A local photon transport model is based on collisional radiative (CR) model calculation of neutral atom [51] and assumes photon absorption in the small optical pass (2 mm radius sphere from the photon source), where a database of effective ionization rate-coefficient is determined as a function of ne, Te, neutral density (n0), temperature (T0) and molecule density. With considering the photon absorption process, an effective ionization rate coefficient is enhanced in low Te (1-2 eV) and high (ne  10 21 m −3 ) plasma when compared to the recombination rate, as shown in Figure  9a. Therefore, profiles of ne div and Te div in the detached divertor will be modified. A selfconsistent iterative calculation of the photon transport model and the SONIC code was carried out for the SlimCS DEMO inner divertor [52], where the plasma density was significantly high up to several 10 21 m −3 , and volume and surface recombination processes contributed to the divertor heat load. Figure 9. (a) Effective ionization rate coefficient without and with the photon trapping for neutral temperature (T0 = 1 eV) and density (n0 = 10 20 m −3 ) and the three electron density cases [52]. Simulation for JA DEMO: total target heat load profiles at (b) inner and (c) outer divertors with and without photon transport model. Self-consistent photon transport simulation was also performed for the JA DEMO divertor as shown in Figure 3, where the maximum ne div is lower and Te div is comparable

Development of Photon Transport Model and Its Effect of Plasma Detachment
Under the detached condition of the DEMO divertor, i.e. high neutral density and large divertor size, photon absorption of Lyman and Balmer line radiations by neutral atoms will enhance the ionization process and change the distribution of the detached plasma. Plasma density becomes significantly high, such as a few 10 21 m −3 , as shown in Figure 3. Typical optical pass length of Ly-α is decreased to~1 mm, and ionization process is enhanced locally by the photon absorption. A local photon transport model is based on collisional radiative (CR) model calculation of neutral atom [51] and assumes photon absorption in the small optical pass (2 mm radius sphere from the photon source), where a database of effective ionization rate-coefficient is determined as a function of n e , T e , neutral density (n 0 ), temperature (T 0 ) and molecule density. With considering the photon absorption process, an effective ionization rate coefficient is enhanced in low T e (≤1-2 eV) and high (n e ≥ 10 21 m −3 ) plasma when compared to the recombination rate, as shown in Figure 9a. Therefore, profiles of n e div and T e div in the detached divertor will be modified. A self-consistent iterative calculation of the photon transport model and the SONIC code was carried out for the SlimCS DEMO inner divertor [52], where the plasma density was significantly high up to several 10 21 m −3 , and volume and surface recombination processes contributed to the divertor heat load.

Development of Photon Transport Model and Its Effect of Plasma Detachment
Under the detached condition of the DEMO divertor, i.e. high neutral density and large divertor size, photon absorption of Lyman and Balmer line radiations by neutral atoms will enhance the ionization process and change the distribution of the detached plasma. Plasma density becomes significantly high, such as a few 10 21 m −3 , as shown in Figure 3. Typical optical pass length of Ly- is decreased to ~1 mm, and ionization process is enhanced locally by the photon absorption. A local photon transport model is based on collisional radiative (CR) model calculation of neutral atom [51] and assumes photon absorption in the small optical pass (2 mm radius sphere from the photon source), where a database of effective ionization rate-coefficient is determined as a function of ne, Te, neutral density (n0), temperature (T0) and molecule density. With considering the photon absorption process, an effective ionization rate coefficient is enhanced in low Te (1-2 eV) and high (ne  10 21 m −3 ) plasma when compared to the recombination rate, as shown in Figure  9a. Therefore, profiles of ne div and Te div in the detached divertor will be modified. A selfconsistent iterative calculation of the photon transport model and the SONIC code was carried out for the SlimCS DEMO inner divertor [52], where the plasma density was significantly high up to several 10 21 m −3 , and volume and surface recombination processes contributed to the divertor heat load. Figure 9. (a) Effective ionization rate coefficient without and with the photon trapping for neutral temperature (T0 = 1 eV) and density (n0 = 10 20 m −3 ) and the three electron density cases [52]. Simulation for JA DEMO: total target heat load profiles at (b) inner and (c) outer divertors with and without photon transport model. Self-consistent photon transport simulation was also performed for the JA DEMO divertor as shown in Figure 3, where the maximum ne div is lower and Te div is comparable Figure 9. (a) Effective ionization rate coefficient without and with the photon trapping for neutral temperature (T 0 = 1 eV) and density (n 0 = 10 20 m −3 ) and the three electron density cases [52]. Simulation for JA DEMO: total target heat load profiles at (b) inner and (c) outer divertors with and without photon transport model. Self-consistent photon transport simulation was also performed for the JA DEMO divertor as shown in Figure 3, where the maximum n e div is lower and T e div is comparable to those for the SlimCS DEMO divertor. Larger effects were mainly seen in wide region of the inner divertor; the target heat load decreased by 20% near the separatrix and 30% at the peak q target location, as shown in Figure 9b. This was explained by decreasing T e div and increasing n e div , due to enhancement of the ionization by the photon excitation. At the same time, n 0 at the target (n 0 div ) and pumping particle flux in the divertor were reduced. Similar effects were also observed in the outer divertor as shown in Figure 9c; reduction of the target heat load was less than 10% near the peak q target location due to lower n e div of~1.5 × 10 20 m −3 . As a result, for the DEMO-level power and particle handling, energy transport accompanying with the photo-excitation will significantly affect the evaluation of the heat load profile in the fully detached divertor and very high n e div (≥ a few 10 21 m −3 ) as well as the plasma profile.

Neutral-Neutral Elastic Scattering Effect on Detachment and Neutral Pressure
In the NEUT2D code, tracking of a neutral particle (D 0 ) after either charge-exchange reaction or elastic scattering with an ion is calculated with selecting a new velocity, as explained in Reference [2]. A molecular particle (D 2 ) is generated with a small probability of D 0 particles incidental to the vessel wall and target, and it is also traced with the same technique with considering dissociation and ionization processes with electron and elastic scattering with ion.
In the high-recycling and detached divertor conditions, neutral-neutral elastic scattering (NNC) plays an important role for the neutral and molecular distributions, the particle exhaust rate and the resultant detachment formation. An NNC model was implicated in the NEUT2D code on the basis of the latest NNC cross-section database, where collision rates and momentum exchange rates were evaluated from differential cross-section database for D 0 -D 0 , D 0 -D 2 by Krstic [53] and D 2 -D 2 by Phelps [54]. The NNC model is a rather theoretical expression when compared to those used in EIRENE [55]. The preliminary result shows that neutral pressure in the sub-divertor region is increased by a factor of two to three, as shown in Figure 10. Here, exhaust port from the sub-divertor was located under the dome in this simulation case. The mean-free-path of neutrals is decreased to 0.01-0.1 m, where molecular transport can be treated as the transient regime between the molecular flow and the viscous flow. Neutral pressure in the sub-divertor region significantly affects the engineering design of the particle and He exhaust system. The NNC model was not used in the power exhaust study since the benchmarking with the other neutral transport code such as EIRENE and DEGAS was not done, and it requires large resources. Calculation speed should be increased to handle two species in the MC calculation. to those for the SlimCS DEMO divertor. Larger effects were mainly seen in wide region of the inner divertor; the target heat load decreased by 20% near the separatrix and 30% at the peak qtarget location, as shown in Figure 9b. This was explained by decreasing Te div and increasing ne div , due to enhancement of the ionization by the photon excitation. At the same time, n0 at the target (n0 div ) and pumping particle flux in the divertor were reduced. Similar effects were also observed in the outer divertor as shown in Figure 9c; reduction of the target heat load was less than 10% near the peak qtarget location due to lower ne div of ~1.5 × 10 20 m −3 . As a result, for the DEMO-level power and particle handling, energy transport accompanying with the photo-excitation will significantly affect the evaluation of the heat load profile in the fully detached divertor and very high ne div ( a few 10 21 m −3 ) as well as the plasma profile.

Neutral-Neutral Elastic Scattering Effect on Detachment and Neutral Pressure
In the NEUT2D code, tracking of a neutral particle (D0) after either charge-exchange reaction or elastic scattering with an ion is calculated with selecting a new velocity, as explained in Reference [2]. A molecular particle (D2) is generated with a small probability of D0 particles incidental to the vessel wall and target, and it is also traced with the same technique with considering dissociation and ionization processes with electron and elastic scattering with ion.
In the high-recycling and detached divertor conditions, neutral-neutral elastic scattering (NNC) plays an important role for the neutral and molecular distributions, the particle exhaust rate and the resultant detachment formation. An NNC model was implicated in the NEUT2D code on the basis of the latest NNC cross-section database, where collision rates and momentum exchange rates were evaluated from differential cross-section database for D0-D0, D0-D2 by Krstic [53] and D2-D2 by Phelps [54]. The NNC model is a rather theoretical expression when compared to those used in EIRENE [55]. The preliminary result shows that neutral pressure in the sub-divertor region is increased by a factor of two to three, as shown in Figure 10. Here, exhaust port from the sub-divertor was located under the dome in this simulation case. The mean-free-path of neutrals is decreased to 0.01-0.1 m, where molecular transport can be treated as the transient regime between the molecular flow and the viscous flow. Neutral pressure in the sub-divertor region significantly affects the engineering design of the particle and He exhaust system. The NNC model was not used in the power exhaust study since the benchmarking with the other neutral transport code such as EIRENE and DEGAS was not done, and it requires large resources. Calculation speed should be increased to handle two species in the MC calculation.

Conclusions: Progress Summary
An integrated divertor simulation code, SONIC, has been developed in order to predict self-consistent transport solution of the plasma, neutral and impurities in the SOL and

Conclusions: Progress Summary
An integrated divertor simulation code, SONIC, has been developed in order to predict self-consistent transport solution of the plasma, neutral and impurities in the SOL and divertor. The SONIC code has been applied for analysis of experiment results in JT-60U and design studies of JT-60SA. It has also applied to Japanese DEMO divertor designs for SlimCS and JA DEMO. While the roll-over of the ion flux and plasma pressure was reproduced at the onset of the plasma detachment, further developments of the plasma and neutral transport modeling are still necessary to quantitatively simulate the divertor detachment. Since the detached divertor operation is a crucial for the DEMO rector design, recent version of the SONIC were used for the divertor design for the DEMO operation scenario. This paper summarized applications of the SONIC code to the JA DEMO divertor design study and recent major achievements.
A new framework of MPMD approach and an MPI data exchange scheme was developed for SONIC code (version 4), which is suitable for multi-author code development. This SONIC revision contributed particularly to improve performance of multiple impurity calculations by IMPMC codes. In addition, an acceleration technique of "backflow model" was incorporated to save computer resource for the impurity atom transport in the sub-divertor.

•
Divertor design and reference operation scenarios for JA DEMO were determined JA DEMO divertor had a power handling of P sep = 230-290 MW and the power handling parameter of P sep /R p = 27-34 MWm −1 , corresponding to levels two times larger than that of ITER. While the JA DEMO divertor covers all divertor volumes below the X-point based on the ITER divertor, a longer leg length (1.6 m) was chosen to obtain compression of neutral particles and efficient formation of the plasma detachment under the low density and high temperature SOL plasma condition. The SONIC code has been used to evaluate the heat load and the divertor plasma for the two power exhaust scenarios, i.e., P sep~2 83 MW for JA DEMO 2014 design and P sep~2 35 MW for JA DEMO with higherκ plasma design. Ar seeding was used as a reference impurity, and series of n e sep scan solutions were obtained for the reference divertor radiation fraction of f* rad div~0 .8. The primary requirement of the divertor operation (peak q target ≤ 10 MWm −2 ) was satisfied in the low n e sep range of 2 -3 × 10 19 m −3 , while the lower P sep scenario had larger heat load margin.

•
Effect of SOL heat flux profile on the partial detachment in the divertor was clarified Large q // near the separatrix ("near-SOL") was significantly reduced in the divertor due to large radiation loss and detachment, and the peak q target at the outer divertor appeared at the boundary of attached plasma region in the partial detachment. The simulation showed ion convective transport from the outer divertor to SOL ("flow reversal") near the separatrix, which did not contribute to the ion heat flux towards the outer divertor and produced a "shoulder" in the q //i Xp profile. Therefore, the total q //e Xp + q //i Xp profile had 1.2-1.6 times larger e-folding length in "near-SOL". The flow reversal also reduced the impurity retention in the outer divertor, which may sustain the attached plasma region and produce the partial detachment.
Since the flow reversal from the outer divertor target was not found in the major tokamak experiments, i.e., it was measured to the target (forward) direction near the X-point, further studies of modeling and experiment are necessary even in existing devices.

•
Operation boundary of JA DEMO divertor was clarified for key power exhaust parameters Severe case studies of reducing f* rad div to~0.7 or/and diffusion coefficients to the half values showed that width of the partial detachment was decreased and the peak q target was increased with increasing the local T i div and T e div , accompanied by shifting the flow reversal peak toward the separatrix. Case studies of f* rad div~0 .7 showed that higher n e sep operations such as those larger than 2.6 × 10 19 and 2.3 × 10 19 m −3 were required for the JA DEMO 2014 and JA DEMO with higher-κ scenarios, respectively. For the smaller χ and D cases, the influence of reducing f* rad div~0 .7 was significantly seen both in the peak q target and n e sep , i.e., q target was increased to larger than 10 MWm −2 and n e sep also increased, and the divertor operation expected at n e sep higher than~3 × 10 19 m −3 .
Consequently, the proposed divertor geometry and reference power exhaust parameters (P sep = 230-290 MW, f* rad~0 .8, χ = 1 m 2 s −1 , D = 0.3 m 2 s −1 and n e sep = 2 -3 × 10 19 m −3 ) were so far consistent with the JA DEMO design concepts. Severe case studies showed margins from the reference parameters, and higher n e sep operation is preferable if a severe parameter is required.
Further improvements of the divertor geometry and operation options such as different seeding impurity will be necessary to extend the partial detachment width and to reduce local T e div and T i div at the attached plasma region. Feasible guidelines of the diffusion coefficient profile over the main plasma edge, near-SOL and far-SOL will be demanded to determine the divertor operation window for the DEMO design.

•
Physics models for divertor simulation were evaluated in DEMO SOL and divertor conditions Three physics models have been investigated in the SONIC code, and influences on the divertor performance were evaluated under the DEMO SOL and divertor conditions; Further validations of simulation models will be required to be incorporated as defaults.
(i) A kinetic model for the thermal force (TF) on impurity ions was evaluated for high temperature and low collisionality SOL plasma. The revised model affected to reduce both the TF on impurity transport and local accumulation peaks of c Ar div under the low collisionality SOL condition. It also reduced P rad sol , but the value was relatively small. (ii) Photon transport model showed that energy transport accompanying the photoabsorption/excitation significantly affected the evaluation of the heat load and the plasma profiles in the fully detached divertor under the DEMO-level high density condition (n e div ≥ a few 10 21 m −3 ). (iii) The neutral-neutral collision (NNC) model was implicated since neutral-neutral elastic scattering plays an important role for distributions of neutrals and molecules, and detachment plasma. The preliminary result shows that neutral pressure in the sub-divertor region is increased by a factor of two to three, which significantly affects engineering design of the particle and He exhaust system.