Review of magnetic islands from the divertor perspective and a simplified heat transport model for the island divertor

Magnetic islands in toroidal confinement devices are reviewed from the viewpoint of their divertor potential. Divertor-relevant geometric parameters are derived analytically, and the relationships among them are revealed. We explain how the island geometry limits the target length and demonstrate the importance of an appropriate numerical tool to minimize the risk of thermal overload of plasma-facing components in the divertor design. The currently available three-dimensional (3D) models are briefly discussed, and their strengths and weaknesses are evaluated. The highlight will be the introduction of a new energy transport model recently developed within the framework of the EMC3 code (Feng et al 2004 Contrib. Plasma Phys. 44 57)—the so-called EMC3-Lite version—primarily for the design and optimization of 3D divertors involving thermal overload concerns. While still undergoing experimental validation with the current graphite divertor of W7-X, it is already being used to develop a subsequent tungsten divertor for W7-X.


Introduction
The formation of magnetic islands in toroidal confinement devices is well described in the literature (see e.g. [1]). However, a systematic discussion of magnetic islands in terms of their divertor potential seems to be lacking. The viability of utilizing magnetic islands as a divertor solution for plasma exhaust-the so-called island divertor concept-has * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. been demonstrated on W7-AS [2,3] and W7-X [4,5] as well as on large helical device (LHD) [6,7]. While W7-AS and W7-X are low-shear stellarators and utilize their intrinsic magnetic islands, the LHD has a large magnetic shear and explored the n/m = 1/1 island created by external magnetic fields-the so-called local island divertor. The experimental results obtained so far with the island divertors are promising. Particularly noteworthy are the HDH mode found in W7-AS [8] and the stable complete detachment achieved in W7-X [5, 9, 10, and 11]. It has been shown that the geometry of the magnetic island can affect the performance of the island divertor [10, 12, and 13]. Therefore, an accurate description of the island topology and a quantitative understanding of the correlation among the key geometric parameters are essential for improving the island divertor concept.
As with all divertor concepts, the basic design principle of the island divertor is to minimize the peak heat load on the heat targets while maximizing the pressure of the recycling neutral gas at the downstream side, which involves the placement and arrangement of a series of heat targets and neutral gas baffles. Unlike the axisymmetric divertor of tokamaks, the island divertor is fully 3D, not only due to the helical structure of the magnetic island but also due to the discrete divertor plates and the resulting local plasma-surface interactions. The resulting scrape-off layer (SOL) consists of regions with different connection lengths [11] and the distribution of the incident angles of the field lines on the targets is complex [14], both of which change with the magnetic configuration. Therefore, the shaping and positioning of the divertor plates (targets and baffles) in an island divertor require a comprehensive optimization procedure for many free parameters in the magnetic configuration space. For example, the coil system of W7-X comprises ten identical modules, each consisting of five nonplanar coils, two planar coils, and one island control coil. They can be polarized and powered independently, creating a magnetic configuration space with a high degree of freedom, which is further complicated by plasma currents [15].
Given the geometric complexity and configuration flexibility, a proper numerical tool is desired in divertor design to avoid thermal overloads of plasma-facing components not intended to receive heat. The development of the tungsten divertor for W7-X, which is already underway, prompted a re-examination of the existing 3D models for their suitability. The current W7-X divertor in operation was designed with the help of a so-called field-line diffusion (FLD) method [16][17][18][19]. It successfully predicted the heat load observed at the targets, but not the heat flux on the adjacent baffles, where permanent hotspots were detected in some magnetic configurations of interest with a maximum heat load an order of magnitude higher than the design value [20]. On the other hand, although capable of calculating the heat load [20], the 3D edge plasma transport code EMC3-Eirene [21,22] is too slow and complex to be used for this purpose.
Given these circumstances, we have extracted and simplified the heat transport module of the EMC3 code and created an 'EMC3-Lite' version utilizing the existing EMC3 facilities, in particular the reversible field line mapping (RFLM) technique [23]. While being validated experimentally with the W7-X divertor, this light version is already being used in the tungsten divertor design. This paper provides a brief description of EMC3-Lite, including the heat transport model and some application examples. A detailed report of the experimental validation results will be presented elsewhere.
The rest of the paper is organized as follows. In the next section, we begin with a mathematical description of the magnetic islands, focusing on divertor-relevant geometric parameters and their correlations. We then use the W7-X island divertor to show the plasma-surface interactions in a typical 3D device. Section 3 reviews and evaluates the FLD model and the EMC3-Eirene code, with the aim of explaining the motivation for developing EMC3-Lite. There, we present the main physical aspects of EMC3-Lite, while leaving the technical description to the appendix. Section 4 discusses the lessons learned and remaining uncertainties in further developing the island divertor.

Geometric aspects of magnetic islands from the divertor perspective
The formation of magnetic islands in toroidal confinement devices is well known in the literature (see, e.g. [1]). For completeness, we begin with a mathematical description of magnetic islands in a toroidal device with a major radius R in which smooth, nested magnetic flux surfaces exist, labelled by r. The path of field lines on each flux surface can be described by a rotational transform ι − (r) = dθ/dϕ, where ϕ is the toroidal angle and θ is the poloidal angle in Boozer coordinates. Now, add a radial field b r (normalized to the toroidal field) that is resonant to a rational flux surface with ι − = n/m, namely where b rm is the magnitude of the resonant field, which, in reality, can be a function of the minor radius. For simplicity, it is assumed here to be constant. This resonant field will tear apart the nested flux surfaces near the rational surface to form a magnetic island. Define ∆θ = mθ−nϕ. The rotational transform can be rewritten as Let r a be the radius of the rational surface, and expand the rotational transform at r = r a as ι − = n/m + ι − ′ ∆r where ∆r = r − r a and ι − ′ is the magnetic shear at r a . In comparison with equation (2), we see immediately that Combining equations (1) and (3), we obtain with Ignoring the dependences of R on ∆θ and ∆r, integrating equation (4) from ∆θ 0 at ∆r = 0 yields which describes a closed curve labelled by ∆θ 0 with 0 ⩽ ∆θ 0 ⩽ π, here called 'island surface'. The continuous black curves in figure 1 show ten such island surfaces with  (6). Ten island surfaces are generated by varying ∆θ 0 from 0 to 0.9π with a step of 0.1π. The O-point is at ∆r/r i = 0 and ∆θ = π, while two X-points are located at ∆θ = 0 and 2π, respectively. The color contour in the top half shows the distribution of the internal field-line pitch normalized to its maximum value, while the lower part shows the lengths (normalized to 4Wr/brm) of the field lines to complete a poloidal cycle around the O-point.
∆θ 0 = 0, 0.1π,…, 0.9π, respectively, which, according to equation (6), possess the up-down and left-right mirror symmetry. The outermost surface with ∆θ 0 = 0 is known as the separatrix and the singularities with b r = b p = 0 (see explanation of b p below) are called the O-point (at ∆r = 0, ∆θ = π) or X-points (at ∆r = 0, ∆θ = 0 and 2π, respectively). W r = 2r i is known as the radial width of the island. Dividing equation (5) by r a , we can rewrite equation (5) as meaning that the island size r i scales linearly with the machine size. In addition, a low magnetic shear and a small poloidal mode number of the perturbation field favor the formation of a large magnetic island. For divertor purposes, the magnetic island must be thick enough to retain most of the recycled or sputtered neutral particles from plasma-facing components and to allow low-Z impurities to remove most of the heating power in the island under detached conditions. It is difficult to quantify how large an island must be for it to be used as a divertor, as the neutral penetration lengths depend on the plasma state and the plasma transport within the island, which in turn depend on many geometrical parameters of the island, not just the radial width of the island. Identifying the key geometric parameters and elucidating their relationships are the main objectives of this paper, while understanding their impact on the island divertor performance is a long-term subject for divertor programs along the low-shear-stellarator line. Regarding the effect of the island size on the detachment performance, dedicated experiments were conducted on W7-AS [12] supported by EMC3-Eirene modeling [13]. It was found that a stable (partial) detachment requires a radial island size of about 3-5 cm, depending on the b r strength. It was argued that the island must be wider than the radial width of the radiation layer [24]. Like W7-AS, W7-X is a low-shear stellarator with a five-fold toroidal periodicity (n = 5). The lowest-order resonances falling within the designed operating window of the rotational transform for W7-X have the poloidal mode numbers of m = 4, 5, 6. Radial field components resonant with these modes are inherently present in the field spectrum of the 3D coils and can produce magnetic islands at the plasma edge that can reach a W r of ∼10 cm. These 'natural' magnetic islands are being utilized at W7-X as a potential divertor solution for plasma exhaust, following the low-shear-stellarator line to a reactor.
Besides b r and W r described above, there are other geometric parameters, which are also important for the island divertor. We start with the island-divertor-relevant poloidal field component b p . The total poloidal field component b θ (normalized to the toroidal field) around the rational surface can be split into two parts, namely b θmn and b p , as where b θmn represents the resonance nature of the magnetic islands and is irrelevant for divertor purpose. According to equation (3), the second part b p can be expressed as which describes the poloidal movement of a field line relative to the X-points. In fact, the interplay of b r and b p causes a rotation of the field lines around the O-point, on which the island divertor is based. With the help of equations (3) and (5), we can rewrite b p as being the maximum magnitude of b p . It is interesting to see that the ratio of b pm /b rm linearly scales with the ratio of r a /r i . Defining W p = 2π·r a /m as the poloidal island width on average, the ratio of b pm /b rm can also be expressed as The pre-factor 4/π arises from the fact that b r is a sine function of ∆θ, while b p linearly scales with ∆r. Using their respective average values ofb r = b rm · 2/π andb p = b pm /2 instead, we can rewrite equation (12) where ε i describes the poloidal elongation of the magnetic island. Under the condition of r i /r a ≪ 1, equation (11) shows that b pm is much larger than b rm , unless the poloidal mode number m is significantly larger than four. Now, we define the internal pitch Θ i as Θ i = b 2 r + b 2 p , which determines the projection of the parallel heat or particle flux on the plane of [∆r, ∆θ] shown in figure 1. Combining equations (1), (10), and (11), we have The upper half of figure 1 show the contours of Θ i /b pm with r i /r a = 0.1 and m = 5. Since b pm is much larger than b rm , the resulting Θ i roughly follows the ∆r-scaling of b p .
The field lines inside the island rotate around the O-point due to Θ i . The length of the field-lines to complete a poloidal turn (with respect to the O-point), L c , can be calculated by combining equations (3), (5) and (6) as where the symmetry condition of the island has been used. The resulting L c -distribution is illustrated in the lower part of figure 1, where L c is normalized to 4W r /b rm .
For clarity, figure 2 shows its profile, i.e. L c as a function of ∆r from the inner separatrix at ∆r = −1 to the O-point at ∆r = 0 and ∆θ = π. Note that the L c shown in figure 2 is normalized to W r /b rm . The separatrix passing through the X-point has an infinite L c . Moving from the separatrix into the island, the L c value quickly falls in a narrow ∆r-range beyond which the further change of L c is no more than a factor of two. As will be discussed later, this narrow ∆r-range is usually much smaller than the typical radial decay lengths of the plasma in the island. L c has a minimum at the O-point. Setting ∆θ 0 → π in equation (14), we obtain the L c of the O-point as where x = ∆θ/ (π-∆θ 0 ), and the cosine function is expanded to the second order in a Taylor series. Equation (15) shows that the ratio of W r /b r is a good measure of L c for the island divertor. For given R, m and ι − ′ , equation (5) shows that W r ∝ √ b rm and therefore L c ∝ 1/ √ W r . This means that the connection length will decrease if external perturbation fields are introduced to increase the island size.

The magnetic island as a divertor
The basic idea of the island divertor concept is to utilize the internal field-line pitch Θ i to guide the edge plasma onto a set of pre-arranged solid plates to facilitate the handling of plasma exhaust. Here, one of the most important design targets is to distribute the heat flux as widely as possible over the solid surface to minimize the peak heat load.
The local heat flux density perpendicular to the surface of a target is given by where α is the incident angle of field lines on the target and q ∥ is the heat flux density along the field-lines (equation (16) may be invalid in the case of very shallow incidence angle, which is not considered in this paper). The parallel heat flux obeys the Bohm condition of where n t ≡ n et = n it , c ⩾ C s = ion sound speed, n et and n it , T et and T it are respectively the electron and ion density and temperature at the target, γ e and γ i are their energy sheath transmission factors, and E i represents an effective potential energy of the ions released during the surface recombination processes on the targets. In this paper, we assume pure hydrogen plasmas with where m i is the proton mass. While q ∥ is determined by plasma transport associated with the island geometry, equation (16) shows that q ⊥ can be reduced by decreasing the incident angle α. Unfortunately, the potential of target inclination cannot be utilized to an unlimited extent due to misalignment of the usually segmented target plates-an engineering issue. Assuming that there is a technical limit for α, α lim , (e.g. 2 • according to [25]), in the following, we explain why the targets in W7-X are discrete.
A continuous target would follow the ι − = n/m resonance. In this case, the maximum achievable incident angle is constrained by the maximum value of Θ i , i.e. b pm . The m = 5 island in W7-X has a r i /r a ratio of about 0.1, which, according to equation (11), results in a b pm = 8 × b rm . In fact, the internal field-line pitch shown in figure 1 sets an upper limit on the incident angle for any continuous target plates. The m = 5 island chain with b rm < 0.001 has a b pm that is smaller than 0.008, corresponding to an incident angle of about 0.46 • . This maximum incident angle is achievable at the outermost island tip, for example, by inserting a plate perpendicular to the island surfaces at ∆θ = π. Moving from the island tip into the island, Θ i reduces and will decrease further if the plate must be inclined with respect to the island surfaces for reasons of island stability due to plasma currents and of configuration flexibility. To increase the angle of incidence, the divertor plate must be inclined with respect to the helical field (the resonant poloidal field component + toroidal field). As a result, the plate has a limited toroidal length and thereby becomes discontinuous, as sketched in figure 3. Shown there is a portion of an island surface with m = n = 5. Under the assumption r a /R ≪ 1, we do not differentiate between the helical and toroidal lengths in this paper. The red dashed lines indicate a field line coming from upstream (top of figure 3) somewhere in the lower part of the island shown in figure 1. The red field line cuts the surface into a spiral band with a width of 2πR·sin(Θ i ) ≈ 2πR·Θ i , which, in fact, defines the poloidal (with respect to the O-point) width of the parallel transport channel of the plasma in the island. The blue dashed line represents a continuous plate that is resonant to the 5/5 island and would result in an incident angle of Θ i , which is much smaller than α lim . To meet the technical requirement, the continuous plate must be interrupted and inclined, as indicated by the solid blue triangle on the left side. Consequently, part of the plasma traveling along the field line to the local target on the left can leave the parallel transport channel through cross-field transport into the shaded region (target shadow in figure 3). The parallel channel has a projection length on the left target of L t = 2πR·sin(Θ i )/sin(α lim ) ∼ 2πRΘ i /α lim , which is by a factor of Θ i /α lim shorter than the blue dashed target in figure 3. In reality, the target is somewhat longer than L t , accounting for the plasma leakage into the target shadow region (TSR) mentioned above. This is also the reason for the emergence of the solid triangle on the right side of figure 3. Associated with resonant nature of the magnetic islands, the field lines in the TSR can reach a length of 2πR, e.g. ∼35 m in W7-X. As will be seen later, plasma transport in the TSR is the most important issue addressed in this paper.

W7-X island divertor
Following the machine symmetry (five-field-period and stellarator-symmetry), the W7-X divertor is made up of ten identical discrete divertor modules. Each module consists of a horizontal and a vertical heat target, and a set of baffles for controlling the recycling neutrals, as shown in figures 4 and 5.  There, the standard and high mirror configurations are used as two typical examples. Both are bounded by the 5/5 island chain, generated by the non-planar coils alone, but with different combinations of these coil currents. In order to maintain the configuration flexibility of W7-X, the targets are poloidally (with respect to the plasma center) almost tangential to the basic flux surface form constrained by the 3D coils-see figure 4. As explained in the previous section, in this target arrangement, the plasma in the island is mainly guided by b r rather than b p from the periphery of the confinement region across the island towards the targets. Assuming Θ i ≈ b r = 0.001 and α lim = 0.035 (2 • ), the targets in W7-X with R = 5.5 m should have a minimum length of 2πRΘ i /α lim ∼1m, which indeed agrees with the toroidal length of the footprint of the long flux tubes on the targets, as shown in the lower part of figure 5.
The resonant radial field creating the magnetic islands in W7-X is not a constant, but decreases rapidly towards the plasma center. As a result, the outer half of the island (with X-point as a reference) is more extended in the radial direction than the inner half. The targets cut the outer half of the island. Thus, the target-to-target connection length, L ct , shown in figure 4 is not the same as those shown in figures 1 and 2. Nevertheless, the L c -estimation made in section 2.1 is largely applicable to the island SOL of W7-X (the region with L ct > ∼100 m in figure 4), including the L c -distributions shown in figures 1 and 2 and the absolute L co value given by equation (15). However, there are two other regions that do not appear in the model island depicted in figure 1, namely the private flux region (PFR) and the TSR as illustrated in figure 4. According to the connection length L ct , there are three SOL regions, namely the island SOL region with L ct > ∼100 m, the PFR with ∼100 m > L ct > ∼35 m (one toroidal turn of W7-X), and the TSR with shorter L ct . A detailed explanation of these regions can be found in [11]. Here, it should be emphasized that a proper estimation of the heat flux in TSR is beyond the capability of the FLD model applied so far to the island divertor. Figure 5 displays the distributions of the incident angle and L ct on the divertor plates of interest. The horizontal target plate is composed of a low-ι target, a so-called high-ι tail and a baffle-like plate in between (middle baffle). The incident angle changes its sign on a single plate or even within a section of a plate. The magnitude of the incident angle varies from zero to a few degrees. In the standard configuration, the footprint of long flux tubes in the island SOL is mainly located at the horizontal target, which moves to the vertical target in the high-mirror configuration. The sign of the incident angle defines from which toroidal direction the plasma comes to the target, while the connection distribution indicates from which SOL region the plasma comes. In both cases, the long flux tubes with L ct > 100 m are indeed terminated at the targets as designed, and the vertical and middle baffles for neutral compression are positioned in the TSR with short connection lengths, where no significant heat load was expected. However, the first island divertor experiments with the standard and high-mirror magnetic configurations showed local heat fluxes on the two baffles that are much higher than their design values-the main trigger for the work in this paper.

Heat load evaluation models
The complex island topology and the local plasma-surface intersection presented in the previous section demonstrate the need of an appropriate physical model and a corresponding numerical tool for the design of the island divertor. In this section, we introduce three models and assess their capabilities and limitations, namely the FLD model, the EMC3-Eirene code, and the EMC3-Lite version.

FLD
The FLD model is widely used in the stellarator community to estimate the strike positions of the plasma on plasmafacing components in three-dimensional environments. It is fast, easy to implement and can reasonably predict the locations where the main plasma-surface interaction will occur. It is a widely applied tool, but the physics involved seems to be less discussed.
The FLD method is actually a Monte-Carlo approach for the source-free continuity equation in a fluid context, including a parallel convective process and a perpendicular diffusion. Field lines, representing particles, are started from the last closed flux surface (LCFS) and assigned a positive or negative ion sound speed (C s ) of equal probability parallel to the magnetic field-lines. The parallel paths of the particles are actually the field lines, and parallel particle tracing is therefore equivalent to field line integration. After a certain parallel integration length ∆ℓ, which can be controlled by a time step τ as ∆ℓ = C s · τ , the field line (particle) is randomly displaced in the direction perpendicular to the field line with a step of where D is the diffusivity and ⃗ ξ ⊥ is a 2D random vector with mean 0 and variance 1. Both C s and D are assumed to be constant. The Monte-Carlo process described above can be found elsewhere [21,26], and it is mathematically equivalent to the following two deterministic equations: where b is the unit vector of the magnetic field, and n + and n − are the densities of the particles with positive and negative C s , respectively. Adding equations (19) and (20) yields where n = n + + n − and u ∥ = (n + − n − ) · C s /n. Equation (21) is known as the continuity equation. It should be mentioned that the boundary condition at the target, i.e.u ∥ = ±C s , is implied. Since C s and D are independent of n and we are only interested in the distribution of the particle flux density n t ·C s on the target, a computational mesh is not needed. The resulting particle deposition distribution is taken as representative for the heat load. Obviously, dividing equation (21) by a constant D would lead to the same heat load distribution. Thus, the ratio of C s /D is the only input parameter for the FLD model.
The FLD model captures an essential physical point of particle transport, namely that the time scale of parallel particle transport is determined by the ion sound speed as τ ∥ ∼ L c /C s . Assuming a simple SOL with a constant L c , the particles will spread out over a perpendicular length scale of which is well known in the literature [27]. The FLD model takes into account the realistic plasmafacing components and the magnetic configurations. The magnetic field is pre-calculated and gridded. Field-line tracing is usually performed using the Runge-Kutta method, where the local field components are obtained through interpolation. It is numerically easy to implement and was the workhorse in the development of the current W7-X divertor. Figure 6 shows the distributions of the power deposition resulting from the FLD model for the standard and high-mirror configurations with D/C s = 5 × 10 −6 m. We see that the horizontal and vertical targets receive most of the power. There is no remarkable heat load on the vertical and middle baffles.  The biggest concern here is not to what extent the particle flux can approximate the heat load, but rather the prescribed parallel flow velocity, which should actually be determined by the parallel momentum balance. Unphysical results will occur in the SOL domains where the plasma cannot enter with a prescribed u ∥ -field. To make this point clearer, we introduce a simple divertor geometry as sketched in figure 7. Particles are started uniformly at the interface between the closed region (red) and the main SOL region (green) and are assigned positive and negative C s with the same probability of 50%. Note that the initial velocity remains unchanged during the whole particle tracing procedure and therefore the destination (the left or right target) of a particle is actually defined already by the direction of its initial velocity. The resulting u ∥ varies from −C s to +C s in the SOL range in contact with the confinement region. Outside this range (divertor legs in figure 7), u ∥ remains at the ion sound speed. Consider a case where there is an in-vessel component, which emerges at the beginning of the divertor leg on the right, as indicated by the filled triangle in the far SOL region in figure 7. In this case, all particles diffusing into the far SOL region between the triangle and the target on the right will move all the way to the target without exception, independently of the parallel distance between the triangle and the target. No particles can reach the right shoulder of the triangle, which is obviously not true from either a kinetic or a fluid perspective.
The actual situation in W7-X is much more complicated than that illustrated in figure 7, but the principle is the same. It is therefore to be expected that the FLD model will underestimate the heat load on the blue triangle on the right in figure 3.

EMC3-Eirene
EMC3-Eirene is a fully 3D edge plasma transport code that includes reduced fluid models for electrons, ions and impurities, and a kinetic description of neutral gas. It self-consistently describes plasma, impurity, and neutral gas transport, their interactions, and PSI processes in realistic 3D magnetic configurations and divertor geometries. It can reproduce the heat load observed at the vertical and middle baffles, which is not predicted by the FLD model. In fact, its success in explaining the hotspot observed at the vertical baffle in the high mirror configuration [20] triggered the re-examination of the limitations of the FLD model presented in the previous section.
Despite its capability of simulating the heat load in 3D divertors, the EMC3-Eirene code is not ideally suitable for scanning a large number of magnetic configurations. One of the main difficulties is the construction of a computational mesh for each magnetic configuration, which is required to solve the nonlinear fluid equations coupled with neutral gas. There are too many constraints in the computational mesh to make an automatic mesh generator possible. Instead, a large amount of manual work is necessary.
The EMC3 code applies a Monte-Carlo method to solve the fluid equations, based on the same principle described briefly in the previous section for the FLD model. A more detailed formulation of the Monte-Carlo procedure developed for EMC3 can be found in [21,26]. Here, we would just like to point out one of the technical highlights of the EMC3 code in terms of handling magnetic field lines, namely the so-called RFLM technique [23]. It is a numerical method for interpolation of magnetic field lines in flux tube coordinates. Compared to the Runge-Kutta method used in the FLD model, the RFLM method is typically three orders of magnitude faster, but with comparable accuracy. This strength of the EMC3 code motivated us to simplify the EMC3 code into a version that is more suitable for developing 3D divertors-the topic of the next section.

EMC3-Lite
According to equation (17), q ∥ is a function of n t , T et and T it . There are two ways to simplify the EMC3 model. The first is to calculate n t while fixing T et and T it . The second is to calculate T et and T it for a given n t . The first concept involves an iterative solution of the continuity and momentum equations for a selfconsistent calculation of the plasma density and the parallel velocity. Here, the neglect of ionization sources is concerning but necessary; otherwise, the Eirene code would be involved, making our simplification less meaningful. The second option is to extract the energy transport module from EMC3. The volume sources of energy due to plasma-neutral and plasmaimpurity interactions can be neglected at low plasma densities so that the Eirene code and the impurity transport module of the EMC3 code are not needed. In addition, the heat load distribution on the target is determined more by the energy transport than by the particle transport. For these considerations, we opt for the second strategy. Below we focus on the physical aspects of EMC3-Lite. Some technical details can be found in the appendix.
Ignoring the energy source terms due to the neutrals and impurities and assuming n ≡ n e = n i = constant, the energy transport equations for electrons and ions included in EMC3 reduce to where Γ ∥ = nu ∥ is the parallel particle flux density, κ e and κ i are the Spitzer heat conductivities for electrons and ions, χ e and χ i are the anomalous conductivities in the direction perpendicular to the magnetic field, and k is the energy exchange coefficient between electrons and ions. The parallel particle flux Γ ∥ = nu ∥ is closely related to the ionization source of the neutrals recycled from the targets. The associated convective heat flux is non-negligible downstream, but is usually much smaller than the other terms elsewhere. It will be considered in the boundary condition, but not in the transport equation. Note that the ions have a much smaller classical heat conductivity than the electrons. Neglecting the convective energy fluxes and the parallel conduction of ions, adding equations (23) and (24) yields where it has been assumed that T ≡ T e = T i and χ = χ e + χ i . At the target, the following Bohm boundary condition holds − κ e ∇ ∥ T| tartget = n · C s γT t (26) where γ = γ e + γ i . Compared to equation (17), (26) omits the potential energy E i , as it is much smaller than γT t in the cases where equation (25) is appropriate. We will assume that χ and n are both constant. Even so, equation (25) is nonlinear because C s scales with the square root of T and κ e depends on T as κ e = κ e0 T 5/2 , where κ e0 is a physical constant. This would mean that a computational mesh is required to represent the T-field, and equation (25) must be solved iteratively. Instead, we restrict equation (25) to low density plasmas where no significant temperature gradient exists between upstream and downstream, and assume that κ e = κ e0 T 5/2 0 where T 0 is an input constant. T 0 is also used to preset C s . Under these assumptions, equation (25) becomes linear, and the heat flux on the target can be simulated utilizing the Monte-Carlo method developed for the EMC3 code without the need of a computational mesh-an advantage of Monte-Carlo over the conventional finite difference or finite volume methods. Here, however, it is necessary to translate the boundary condition of equation (26) into Monte Carlo 'language'.
As the plasma density n is a constant, we divide equations (25) and (26) by n so that the resulting coefficient κ e /n has the unit of diffusivity. Since we are not interested in the spatial distribution of T, we ignore the parallel variation of the magnetic field strength, which would otherwise greatly limit the parallel tracing steps. Then, the parallel heat condition with constant κ e /n is known in probability theory as a Wiener process with the transition probability of: with δ = 2κ e τ /n (28) where δ is the space step and τ is the time step. In the Monte-Carlo algorithm used in EMC3, equation (26) is equivalent to a loss probability p loss , i.e. the ratio of the 'particle' flux defined by equation (26) to the total particle flux that hits the target according to equation (27). By mathematical manipulation, we obtain where is the parallel decay length of T at the target resulting from equation (26). Note that equation (29) is accurate only up to the first order of δ/λ T∥ . Thus, δ/λ T∥ must be kept much smaller than unity, which requires τ ≪ nλ 2 T∥ /2κ e . EMC3-Lite requires χ, n and T 0 as inputs. The total power into the SOL is only a scaling parameter. It is difficult to specify quantitatively the range of n and T 0 to which equation (25) applies. Here, we make an estimate of the parallel temperature gradient length for reference. The assumption that there is no significant T-gradient between upstream and downstream requires λ T∥ > L c /2 (from upstream to target). For hydrogen plasma with γ = 7, n = 10 19 m −3 and T 0 = 100 eV, equation (30) yields a λ T∥ of ∼130 m, which is comparable to, or even slightly smaller than, half of the L c values shown in figure 4. Note that λ T∥ ∝ T 2 0 /n. A more reasonable choice of these parameters will be determined by comparison with experimental results.
The EMC3-Lite version is currently undergoing extensive experimental validation at W7-X to assess its suitability for the tungsten divertor design of W7-X, which is already underway. Detailed model-experiment comparison results will be published elsewhere. Here, we use the standard and high mirror configurations as two typical examples to highlight the essential improvement over the FLD model in predicting the heat load distribution. The results are displayed in figure 8. Compared to the FLD results shown in figure 6, significant heat loads now occur on the middle baffle, especially in the standard configuration. In the high mirror configuration, a hotspot is created on the vertical baffle with a local heat flux similar in magnitude to the peak heat load of the targets. These results agree well with experimental observations. In particular, IR-cameras detected a hotspot at this location of the vertical baffle where the peak heat load was almost an order of magnitude higher than the design value, which severely limited the discharge time with the highly optimized high-mirror configuration in the past experimental campaigns.

Discussion
Magnetic islands are inherently present in non-axisymmetric toroidal confinement devices. If the islands are large enough to intercept most of the recycled neutral particles and sputtered impurities and to allow low-Z impurities to dissipate most of the heat power within them, they have divertor potential for plasma exhaust. The radial size of the islands scales with the square root of the ratio of the major radius to the magnetic shear. Therefore, large low-shear stellarators favor the formation of large magnetic islands. A divertor concept based on such magnetic islands-the so-called 'island divertor' is being explored on the low-shear stellarator W7-X as a potential solution of plasma exhaust under reactor conditions.
The geometric parameters relevant to the island divertor and their correlations have been derived analytically. Their effects on divertor performance have been partially investigated in the W7-AS and W7-X experiments, but there are still many aspects that remain unexplored. The most relevant parameters and their respective functions in the island divertor transport are listed below. (a) The strength of the radial resonant field determines the radial projection of parallel transport, thus controlling the relative weights between the perpendicular and parallel transport of the plasma across the magnetic islands. (b) The radial width of the islands is responsible for screening the neutral particles of both the working gas and impurities. (c) The poloidal width of the islands affects the viscous transport of the parallel counter-flows on the adjacent island legs. (d) The non-resonant part of the poloidal field component conducts heat toward the X-points of the islands, which influences the impurity radiation distribution under detached conditions [11]. (e) The connection length determines the parallel transport time scale of the plasma, and thus the perpendicular width of the plasma deposition profiles on the targets. (f) The internal pitch of the field lines in the islands constrains the practical realization of the divertor targets.
There are different reasons why the island divertor in W7-X is discrete. The most fundamental one is that the maximum achievable internal field-line pitch in the W7-X islands is too small to be fully utilized in practice to reduce the heat flux on the target. Note that the situation will not change as the machine size grows. The localized divertor targets, the high configuration flexibility, and the helical resonant structures of the magnetic islands present challenges for the island divertor design, for which suitable numerical tools are necessary.
Concerning the conventional field-line diffusion model, we have derived the transport equation, explained the main physics involved, and clarified why it failed to predict the baffle overloads observed in the first W7-X experiments. The field-line diffusion model does not include the momentum balance for a self-consistent calculation of the parallel velocity, producing unphysical results in the TSR. As for the EMC3-Eirene code, although capable of performing thermal load calculations, it is less practical to use to scan the vast space of magnetic configuration in W7-X. The large amount of manual work required in grid construction for each magnetic configuration makes the EMC3-Eirene unsuitable for this task. We have therefore simplified the EMC3 code into a new 'light' version, the so-called EMC3-Lite, for developing the tungsten divertor of W7-X, which is already in process. Its creation was triggered by, but its use is not limited to, the W7-X tungsten divertor. EMC3-Lite includes only parallel classical electron conduction and a perpendicular anomalous conductive process, which are the dominant heat transport processes at low plasma densities. It neglects particle transport but includes a parallel particle flux in the Bohm boundary conditions at the target. EMC3-Lite is undergoing extensive experimental validation using the current W7-X divertor as a test bed. It can reproduce the main features of the measured power depositions, including the hotspots on the neutral gas baffles not predicted by the field line diffusion model. Compared to EMC3-Eirene, EMC3-Lite is more automated, especially in grid construction. Thanks to the reversible field-line mapping technique of the EMC3 code, which speeds up the field-line tracing by three orders of magnitude compared to the conventional Runge-Kutta integration method, EMC3-Lite can complete a typical case study for W7-X within a few CPU minutes.
None of the 3D models addressed in this paper can handle the drifts, of which the E × B drift is the most important. Drift effects were indeed observed in the first divertor experiments of W7-X [28]. Nevertheless, no drift theory has been available so far to explain the observed baffle overload. By contrast, EMC3-Eirene and EMC3-Lite can both reproduce the observed baffle overload without drifts. While the extent to which the drifts can relocate the heat load remains an open question, the existing error fields pose another challenge. In these circumstances, the best we can do is to use the experimental data collected so far at W7-X to estimate the range of the uncertainties of the EMC3-Lite model, and then use them as tolerances when building the next tungsten divertor.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Acknowledgments
This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom research and Training Programme (Grant Agreement No. 101052200-EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
The corresponding author would like to thank Drs Dirk Naujoks, Joris Fellinger and Craig Beidler and Professor Per Helander for their valuable comments on the manuscript.

Appendix. Technical aspects of EMC3-Lite
One of the technical highlights of EMC3-Lite compared to the traditional field-line diffusion model is that the field lines needed for describing the parallel transport are interpolated rather than integrated. The numerical method for reconstructing continuous field lines from a number of pre-calculated field lines of finite length is known as RFLM [23]. Here, we briefly explain the basic idea.
A flux tube is defined by four field lines toroidally spanning ∆ϕ = ϕ 2 − ϕ 1 with ϕ 2 > ϕ 1 , and has a convex cross-section everywhere in this ϕ-range. In the flux tube, we introduce α and β as two local independent coordinates so that each field line in the flux tube can be uniquely identified using the two coordinates. Its trajectory in real space is given by where the coefficients on the right hand side are determined by the four given field lines. Note that the two coordinate systems share the common ϕ as the third coordinate. Obviously, the accuracy of this bilinear interpolation scheme depends on the cross section of the flux tube as well as its length [23]. As a consequence of the magnetic shear, the cross-section of a flux tube becomes increasingly deformed as its length increases, thus limiting the actual length of the flux tubes to satisfy the convex cross-section requirement. The resulting flux tubes may be too short to cover an entire toroidal period of the magnetic configuration. In this case, additional flux tubes must be constructed to extend the toroidal domain, e.g. from ϕ 2 to ϕ 3 , where ϕ 3 > ϕ 2 , and this procedure is repeated until the periodicity of the magnetic configuration is completed. To obtain continuous field lines, the following coordinate transformation is performed at the toroidal interface, e.g. at ϕ 2 : Where I ∈ [ϕ 1 , ϕ 2 ] and II ∈ [ϕ 2 , ϕ 3 ] indicate the two neighboring flux tubes through which a specified field line passes. The continuity of the field line in real space is ensured by maintaining its position at the interface during the coordinate transformation.
The actual length of the flux tubes depends on the magnetic shear. For most interesting configurations in W7-X, the flux tubes can be up to half the length of a field period. The edge domain of interest is filled with a large number of flux tubes without gaps and overlaps in real space. Using the standard magnetic configuration of W7-X as an example, figure A1 illustrates a portion of the 2D mesh in the bean-shaped crosssection, which is used to generate the flux tubes. The original 2D mesh is composed of 100 × 512 (radial × poloidal) flux tubes, corresponding to a spatial resolution of about 0.3 cm × 0.6 cm on average. The edge domain of EMC3-Lite starts from a closed magnetic flux surface and ends at the wall. A magnetic flux surface is chosen as the innermost domain boundary to facilitate the treatment of the boundary conditions there. It is not necessarily the LCFS and is automatically determined by EMC3-Lite for each individual magnetic configuration. Each grid point defines a field line terminating in the adjacent triangular plane (one half field period distant from the bean-shaped plane), obtained using the same Runge-Kutta method as the FLD model. Any field lines that extend beyond the wall will be removed, and the remaining field lines form the flux tubes, in which field lines are interpolated. The performance of the RFLM technique is evaluated in terms of accuracy and computational speed, compared to the Runge-Kutta method. Within a comparable accuracy (figure A2), RFLM is three orders of magnitude faster.
The computational mesh described above is only needed for the RFLM technique to reconstruct the field lines. The only constraint on the mesh is that it must be fine enough to meet the accuracy requirement of the interpolated field lines, which makes the automatic generation of the mesh possible. As in the EMC3 code, the Monte-Carlo procedure in EMC3-Lite is formulated in the flux tube coordinates and traces 'particles' without the need to know their positions in real space. However, plasma-facing components are usually specified in real space and are difficult to transform into the flux tubes. Therefore, calculating the intersection of the particle trajectory on the target requires frequent switching between the two coordinate systems, which can be quite computationally time consuming.
A technical solution to this problem is to discretize the targets and baffles of interest in the mesh. For this purpose, each flux tube is toroidally resolved by a certain amount of cells. All cells that intersect the targets or baffles are identified and marked as 'limiter' cells. The coordinate transformation is then carried out only when a particle comes into a limiter cell, which greatly speeds up the computation.
All plasma-facing components are represented in a socalled Kisslinger's format [17,18], in which all plates are cut into toroidal segments given in cylindrical coordinates, as briefly explained in figure A3. The actual tiles that make up the Figure A1. A portion of the 2D mesh in the bean-shaped plane for generating flux tubes. The innermost red curve is a flux surface, while the outermost represents the wall. The original 2D mesh has a resolution of ∼0.3 cm × 0.6 cm (radial × poloidal) on average. Figure A2. Poincare plots of field lines reconstructed using RFLM (green) and calculated by 4th order Runge-Kutta (red) with an integration step of ∼2 cm. The former is a factor 1000 faster than the latter. target panel in W7-X follow this segmentation to a greater or lesser extent. Compared to triangulation, which is commonly used for 3D visualization, Kisslinger's format is more suitable for resolving the heat deposition profiles of interest in this paper. Each tile is divided into a number of quadrangular elements in cylindrical coordinates, which are used to evaluate the heat load calculated by the Monte Carlo process.
The calculated distribution of power deposition suffers from Monte-Carlo noise. The statistical error is inversely proportional to the square root of the number of particles collected in each area element and therefore depends on the resolution and the total number of particles launched. The required computation time is proportional to the total number of the traced particles. For a typical resolution of 0.5 × 2.5 cm 2 (poloidal × toroidal), which is used to produce figure 8, Figure A3. A plate file in Kisslinger's format. The file starts with a note line for the user. The plate is defined by NT toroidal cut lines denoted by ϕ i in degrees with i∈ [1,NT], and each cut line consists of NP points given by [R i,j , Z i,j ] in cm with j∈ [1,NP]. Note that cut lines toroidally outside the EMC3-Lite domain are not recognized. The parameter NPER describes the toroidal periods of appearance of the plate. R_ and Z_shift in cm can be used to adjust the entire plate without changing its shape. where M is the total area elements and q ⊥ij is the heat flux density at the area element A j measured in the i-th realization. As one can see, a result with a statistical error of a few percent can be obtained within a CPU time range of a few minutes. Although primarily designed to evaluate the thermal load of plasma-facing components in 3D environments, EMC3-Lite also provides other capabilities of general interest for investigating the edge topology of magnetic configurations. For example, it includes a Biot-Savart calculator for the magnetic field from given current filaments; It can calculate the distributions of the connection length in real space and on the target and the distribution of incidence angles. Moreover, it also provides some functions to support the mesh generation for the EMC3-Eirene code, especially in determining the outermost boundary of the EMC3 domain-the biggest challenge in the mesh construction. EMC3-Lite is structured modularly and is therefore technically easy to expand.