A general stellarator version of the systems code PROCESS

We present modifications of the fusion reactor systems code Process that allow for a description of a general class of stellarator power plants, based on a stellarator coil-set and the respective MHD plasma equilibrium. For this, we modify Process such that each stellarator configuration enters the systems code via a set of effective parameters which can be calculated in advance before using them in new scaling models in Process. Further, we show two applications of the new Process version: firstly, we apply the code to three reactor-size stellarator devices with different aspect ratios, and secondly, to three coil-sets optimized for the same equilibrium with varying coil numbers.


Introduction
Stellarators are attractive candidates for a fusion power plant: they operate in steady-state and can be optimized for minimal plasma current thus avoiding current driven instabilities. Further, they do not necessarily rely on large poloidal field coils or a central solenoid. Stellarators also benefit from large connection lengths in island divertor configurations, easing power exhaust. Finally, the highly dimensional design space can be utilised to optimise the configuration according to relevant physics and engineering requirements at the cost of geometrical complexity.
The recent start of operation of the prototype advanced stellarator Wendelstein 7-X (W7-X) has shown that such config- * 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. urations can be realized with sufficient engineering accuracy [1,2], providing further incentive to study fusion power plant designs based on stellarators.
As of now, attractive stellarator configurations are developed through the process of stellarator optimisation [3][4][5][6][7], where a computational framework optimises a threedimensional MHD equilibrium and a corresponding coil set to fulfill a set of mostly physics-related figures of merit. To our knowledge, there exists so far no systematic framework that checks a configuration achieved by stellarator optimisation for a broader range of engineering constraints specific to fusion reactor design such as superconductor or neutronic limitations. Also, there is currently no framework available that is capable of quickly exploring a larger design space around a reference design point, while simultaneously judging the technological and economical feasibility. Systems codes can fill this gap between the conceptional magnetic configuration and the reactor technology, as visualized in figure 1.
Systems codes are coherent, holistic computational frameworks that aspire to model the crucial features of an engineered system. They typically consist of a set of simplified models that depict the governing design parameters and constraints. In the context of fusion power plants, the use of systems codes has several advantages: (a) They can check the feasibility of a given fusion reactor design point in a holistic way by taking physics, technological and economical constraints within one framework into account. (b) They can be used to find technologically or economically more suited design and operation points of a fusion reactor. (c) They can easily adapt technological advances due to their modular structure and thus allow a fast re-iteration of fusion reactor design points, when new technology becomes available. and which serves as an interface between stellarator optimisation and Process, using the MHD equilibrium and coil filaments to prepare a set of effective parameters for the systems code models. Secondly, Process itself is modified in a way to include this set of parameters in newly implemented or modified models, allowing to perform calculations involving different stellarator configurations. The two steps are explained in more detail in the next section. The outline of this paper is as follows: in section 2 we describe the structural changes in Process that were necessary to include more general stellarators. In section 3 we describe the newly developed models and their implementation. Finally, in section 4 we employ the Process framework with the implemented stellarator-specific changes for three example studies: first, the new magnet system model is benchmarked against a tokamak reference case. Secondly, we model three different stellarator configurations with distinct aspect ratios, using a 3, 4 and 5 periodic Helias configuration from [20]. Thirdly, we vary the number of coils for a specific W7-X equilibrium, scale the machine to reactor size with Process and study the impact of different coil numbers on the coil properties.

New workflow for stellarator-PROCESS
Stellarators, by their 3D geometry, impose non-trivial physics and engineering constraints on a fusion power plant design. For example, in contrast to tokamaks, the magnetic field strength on the inboard side of the coils can be different for every coil, the divertor area depends on the location of the magnetic islands, or the neutron wall load has large variations not only in poloidal, but also in toroidal direction. Further, stellarators can have vastly different coil and plasma boundary shapes. Thus, an accurate representation of systems codes relevant features at low computational cost is quite challenging for general stellarators. To mitigate this issue, we introduce an additional, automatized, calculation step between the output stemming from stellarator optimisation and the inputs that go into the systems code, as schematically shown in figure 2. In practice, the work-flow then is as follows. 'Stellarator optimisation' provides a 3D MHD equilibrium and a set of corresponding, as fixed considered, coil filaments at a reference point in major radius and aspect ratio. This reference point (equilibrium and coils) we denote with the symbol C from here on, which serves as input for the detailed calculations. The newly introduced intermediate calculation step (essentially the first part of the systems code models), involves accurate, but comparatively slow computations at this reference point. The result of these computations are a set of configuration-dependent effective parameters a i (C), which serve as input for newly implemented exact, fitted, or empirical scaling equations in the systems code.
The general idea behind this approach is to separate computationally heavy operations from the systems code. This means that every stellarator-specific systems code model consists of essentially two parts. The first part entails the detailed modelling of a sub-system outside the systems code. The second part, in turn, involves an associated (fast) scaling equation within the systems code that makes use of the results from the The workflow of the pre-calculation step: a configuration C (coil filaments and flux surfaces) is assumed as input from stellarator optimization. A set of Process relevant parameters a i is calculated based on a reference point C, which Process uses to calculate and optimize an iteration vector x for a reactor design point, according to an objective function f and according to the applied constraints. The found design point can be used again as feedback for stellarator optimization. detailed calculations. An example here would be the computation of the maximal coil force density f max (C) as effective parameter from 3D calculations for a reference coil set. In this example the scaling of f max within the systems code then is a linear scaling law in B max and the current density j, both parameters that the systems code optimizes for.
We implement the systems code models in a way that they reflect extrapolations of the reference point C in the following macroscopic design parameters: the major overall size of the machine (coil and plasma size), the minor plasma radius a at constant coil radius, and the total magnetic field strength on axis B t . For the plasma design, the implemented scaling parameters are the plasma density, temperature, and the ISS04 'renormalization' factor (a measure for the configuration-dependent quality of energy confinement [21]). The stellarator-Process version is capable of optimizing for devices by scaling these parameters as a part of the optimization vector now. In addition to the above listed set of iteration parameters, Process also optimizes in the engineering parameter design space, with the 'usual' parameters such as winding pack size, coil quench times, critical current density safety margins in the superconductor, copper fractions in the winding pack, net electricity output, etc, also see [8,9].
Note that by this prescription the coil number and the coil shapes are considered fixed by stellarator-Process and only the overall size of the coils is scaled. A broader device scan in different stellarator configurations or different coil-sets can be done by sampling different configurations C using stellarator optimization codes.

Models
Below, we introduce the newly developed stellarator-specific systems code models that aim to describe a general class of stellarators with a modular coil set, irrespective of their shape. The stellarator modifications to Process are comprehensive in the sense that they allow an equivalent modeling stellarators compared to the tokamak treatment [8,9].
For each model we describe both the external procedure of calculating the effective parameters as well as the systems code internal scaling equations. The effective parameters that are calculated in the external step are distinguished into two categories. The first type are configuration-specific quantities, that are used directly in follow up calculations and these are denoted by a i (C). The second type of parameters are those that are calculated as a reference point for the scaling equations and these are denoted as hatted values,â i (C), where C represents the configuration stemming from stellarator optimisation (3D MHD equilibrium and associated coil filaments).

Plasma volume and surface
The plasma volume V and the plasma surface area S are basic properties in Process. For example, subsequent calculations of the fusion power, fuelling rates, or material loads depend on the plasma volume. Similarly, the surface area is an important quantity to approximate the first wall area and to scale the heat flux densities.
The spatial location of stellarator-symmetric flux surfaces can be parameterized by a set of Fourier coefficients R c m,n and Z s m,n , where m and n are the poloidal and toroidal mode numbers respectively. The cylindrical coordinates for each flux surface can be obtained by Here, u describes a poloidal coordinate, v the polar toroidal coordinate, and s is a flux surface coordinate [22]. Equations (1) and (2) hold for stellarator symmetric configurations with a field period symmetry of N f . The volume enclosed by the last closed flux surface can be calculated for a reference size (R,â) according tô The surface area of a flux surface can be calculated bŷ √ g is the Jacobian determinant and |.| is the Euclidean norm.
The valuesV(C) andŜ(C) are calculated in the pre-processing step for a reference point in major radiusR and minor radiuŝ a. Within Process, the plasma volume and surface area is then simply obtained by the following scaling equations,

0D-transport
The 0D-transport model in Process imposes a power balance as an equality constraint, The left-hand side includes contributions from confinement loss P conf Loss , from bremsstrahlung P br , line radiation P line and synchrotron radiation P sync . The right-hand side includes heating from fusion alphas P α , a term of charged non-alpha particle heating P ¬α (e.g. in D-D fusion) and a term for auxiliary heating P aux . Writing these expressions explicitly, equation (6) becomes Here, f α is the fraction of the alpha particle energy that is deposited in the plasma, which is an input parameter in Process and depends on the configuration. Similarly f ¬α accounts for the particle confinement fraction of non-alpha particles. PROCESS' model for radiation losses (P br , P line , P sync ) is described in [23,24]. For P conf Loss , Process uses the effective energy confinement time τ E to determine the effective power transfer P conf where W is the total plasma energy. The energy confinement time τ E is obtained via empirical scaling laws. The used scaling law for stellarators in Process is the so-called ISS04 scaling [21], where a is the minor radius, R 0 is the major radius, n is the line averaged electron density, B t the toroidal magnetic field, ι 2/3 ≡ ι 2/3 (C) is the rotational transform (at s = 2/3), P is the combined effective plasma heating, and f ren is a proportionality factor that measures the magnetic configuration dependent deviation from the ISS04 scaling law. In principle, f ren is determined by C directly, although a reliable a priori method of calculating this factor is not available up to date. Instead, Process can iterate f ren within user set boundaries and return a needed configuration factor for the optimized power plant design point. The stored energy W in equation (8) is obtained from the imposed profiles for particle species averaged density n and temperature T: In stellarators, ρ is usually chosen as the effective radius, which fulfills g(ρ) ∼ ρ.
The temperature and density profile shapes for the electrons are input parameters in Process and can be specified using the parametric form Process implements the ion profiles as (user defined) multiples of the electron profiles. These profiles are taken to compute the radiation terms in the left-hand side of equation (7), see [23].
It should be noted that the imposed profile shapes are not per se consistent with the implied heating schemes or transport properties. However, in practice, the profile shapes can be determined by transport simulations independent of the systems code. Results from such simulations can then be used as input for Process, e.g. in profile shapes or heating source. Equation (7) serves as equality constraint in Process.

0.5D neoclassical transport model for stellarators
As Process lets the user choose T 0 , n 0 , α n and α T in equation (13) freely, we introduce a 'sanity check' of the confinement time here against a neoclassical model. The energy balance equation in steady state is Here, q is the flux surface averaged energy flux and p stands for the flux surface energy density sources and sinks. If one assumes constant energy flux on a flux surface, integrating equation (14) over a volume up to a radius ρ x yields where S(ρ x ) is the surface area at a radius ρ x . P rad is the radiation power and P heat is the heating power as specified in equation (7), both integrated values in the of S(ρ x ) enclosed volume. In Process, we choose ρ x = ρ core , where ρ core is an input parameter in Process, which determines the radius of a binary 'core' treatment [8]. ρ core is usually chosen in the order of ∼ 0.6 (ρ = 1 matches with the last closed flux surface).
The new model in Process now calculates a maximal allowable q max with the calculated heating and radiation power as Here, p α,rad V denotes the power density averaged over V(ρ core ). The volume over surface ratio at ρ core can be obtained approximately by scaling of equation (5).
Equation (16) can be compared against heat fluxes q neo from neoclassical theory, e.g. [25,26]. In Process we compare equation (16) against a neoclassical electron flux [26] where we take the profile shapes as given by Process and further assume the electrons to be in the 1/ν collisional regime and neglect the effect of the radial electrical field. The collisionality ν(n, T) can be calculated from classical statistical theory [26]. eff ≡ eff (C) is the averaged effective helical ripple and is an input parameter, which is calculated for every configuration C. q e,neo serves as an order of magnitude check for q max , as a design point with 2q e,neo ∼ q max indicates that profile gradients at the found design point cause similar purely neoclassical transport fluxes to q max and would not allow for an unknown turbulent heat flux q turb .
Using this model we try to circumvent the consistency issues and restrict profile gradients in the 0D transport model of Process for stellarators.

Density limit
The density in stellarators devices is, at least empirically, bound by the Sudo limit [27], which accounts for excessive impurity radiation at high edge densities. This limit is proposed in the parametric form as Stellarator-Process can enforce this limit or multiples thereof. However, equation (20) was exceeded in W7-X and LHD experiments [28,29] and is likely dependent on edge impurity concentrations which are not governed by equation (20). There is however another density constraint, which is imposed by operational boundaries of an electron cyclotron resonance heating (ECRH) scheme in ECRH heated stellarator devices [30]. For reactor scenarios, ECRH heating using the O1 mode appears to be most suitable as it heats the lowest resonance of the electron gyro-frequency and thus requires lower gyrotron frequencies than higher resonant heating schemes. O1 heating implies the operational constraint where ω pe is the plasma frequency, ω gyro the gyrofrequency and ω max the maximum available gyrotron frequency. ω max depends on the available gyrotron technology and can be set by the user as an input. The critical density is reached when the plasma frequency matches the electron cyclotron frequency. Thus, the central electron density n e is limited to: n e < n ECRH e,crit = m e 0 e 2 ω 2 gyro , subject to: ω gyro < ω max . (22) Figure 3. Central density limits due to different ECRH heating schemes: the blue region indicates where O1 heating can be applied, the green region where X2 is feasible. Each shape area indicates a minimum required gyrotron frequency. For context, dashed lines indicate ignition according to Lawson criterion with different volume V and volume averaged ion temperatureT i (assuming n peak /n = 3 and ISS04 scaling from W7X parameters). Figure 3 visualizes the heatable densities with O1-heating, and in comparison an X2-heating scheme, with different maximal available gyrotron frequencies at varying magnetic field strengths. Equation (22) is implemented as a constraint in Process and ensures that the found design point is ECRH heatable in O1 mode. Note that there are heating schemes, such as electron Bernstein waves [31] or an X1 heating scheme, which could be used to heat a plasma beyond equation (22), but their relevance as a heating scheme in a stellarator reactor are still up for discussion and are not taken into account by Process yet.

Island divertor
There are three studied divertor concepts available for stellarator reactors: an ergodic divertor concept, also called helical divertor, for high shear configurations [32], a resilient nonresonant divertor concept [33] and a resonant, island divertor concept [34][35][36][37]. For now, we include only a description for an island divertor concept in Process, closely following the (previously implemented) model as proposed in [16].
In a stellarator with an island divertor concept, the magnetic field is designed such that the rotational transform ι res at the edge coincides with a low order rational number N p k/n, where m is the number of poloidal resonances (islands), k is the resonance order and N p is the field period of the machine. k is determined by radial B-field harmonics on or shortly behind the last closed flux surface, and, if the respective resonant harmonics are not actively suppressed, is typically equal to 1. The underlying concept of the island divertor is to use the magnetic islands for diverting the heat load coming from the plasma core and then intersect the islands with discontinuous divertor target plates. While the full physics description of the stellarator scrape-off-layer (SOL) is still a challenging and contemporary topic, fundamental geometrical considerations can be used to estimate the heat load on the divertor target plates. It is the goal of the proposed model here to provide an estimation of the peak heat load, as this is the constraining engineering limit, due to material limitations. The heat load on the divertor target plates q div is the ratio of the power arriving at the divertor P div and the area over which this power is effectively spread, A eff . One of the major strategies to reduce the heat load arriving at the divertor is to introduce low-Z impurities that are effective at radiating substantial power in the SOL. Consequently, the power arriving at the divertor is the power coming from the plasma core P core less the radiation from the impurities: where f rad is the radiation fraction, which needs to be given as an external input parameter.
The wetted area A eff on the divertor plates usually has the form of a strike-line with a total length L tot across all divertors and a width λ int . The heat load is then where P core is provided by the Process' plasma core model. Assuming that the heat load is distributed in equal shares across all divertor plates, then the total length L tot is simply the sum over all divertor targets L i , Here n = kN p , as defined previously. The strike-line length L strike on a single divertor plate can be estimated from the field line geometry. To this end, one needs to introduce the pitchangle Θ = dr/dl, which describes the radial displacement of a field line in the SOL along its arc-length and depends on the specific magnetic configuration C, but it is typically in the range of 10 −3 -10 −4 for stellarators. The strike-line is limited by the field line that just passes the divertor plate at the front and then after one toroidal turn (Δl ≈ 2πR) hits the target plate on the far side. Using the definition of the pitch-angle, the radial projection of the strike-line is Δr = 2πRΘ. The length of the strike-line on the divertor plate itself is then determined by the angle α lim = Δr/L strike under which the field line hits the target plate. The strike-line length on the divertor is then simply where F x is an additional broadening of the flux channel caused by diffusive cross-field transport. A model for this factor is given below in equation (30). A small intersection angle α lim helps to increase the strike-line length and reduce the heat load density. However, α lim is limited by the engineering accuracy under which target elements can be arranged, typically around ∼ 2 • . Generally, stellarators with an island divertor feature much longer connection lengths than tokamaks [38]. Consequently, the energy and particles have a longer dwell time in the SOL leading to a substantial cross-field broadening of the transport channel compared with tokamaks. We assume here that the cross-field transport is mostly of diffusive nature, allowing us to describe the strike-line width (also referred to as power decay width) by [39], Here, χ ⊥ is the perpendicular diffusion coefficient, which is an user-defined input, but usually taken in the order of ∼ 1 m2 s −1 [40]. τ is the characteristic dwell time of the particles in the SOL before reaching the target. As the particles follow the field lines, the dwell time τ depends on the connection length L c of the field line and the average speed of the particle, namely the ion sound speed c s = 2T/m (m here being the ion mass), and thus τ = L c /c s . The ion temperature (in the SOL) T is again a user-defined input, however since mostly detached scenarios are considered for a reactor design point for divertor protection, T must be on the order of 5-10 eV [40].
The connection length L c can be geometrically estimated by using again the definition of the pitch-angle Θ. If we define Δ as the radial distance from the LCFS to the target plate, then the connection length is simply The typical radial scale length Δ of the system is for the island divertor the radial extent of the magnetic islands w i . However, as the island is intersected by the divertor plates, only a fraction f of the island width is effectively used Δ = f · w i . Usually, the divertor plates are placed at the half radius of the islands, thus f is normally in the order of f ∼ 0.5. The full width of the island can be estimated from analytic theory [41], where ι = dι/dr is the magnetic shear at the edge, which is given by the magnetic configuration. Generally, stellarators with an island divertor need a comparably low magnetic shear in order to form sufficiently large magnetic islands. Finally, the previously mentioned flux channel broadening F x can be derived following the same diffusive ansatz, but for only one toroidal turn, which then becomes In conclusion, we have provided equations for all introduced parameters. Consequently, all the here derived relations can be consolidated in order to arrive at a heuristic scaling for the divertor heat load. By replacing the poloidal mode number m in terms of ι, k and N p , using equation (23) one obtains Here, ι (C), ι (C), N p (C), k (C) and Θ (C) are specific to the considered magnetic configuration and are easily obtained in the pre-processing step. χ ⊥ , α lim , f and T depend on the specific physics regime or the engineering design and must be provided by the user, but usually take values as indicated in the text above.
It is planned to validate this model against experimental results from W7-X in the future. Due to the analytic nature of the model, it will be possible to quickly adapt and test new findings and advances.
It should also be noted, since the heat load is usually limited by material constraints, the divertor model is also useful in reversing the parameters. For example, for a fixed design point and heat load limit, one can estimate the required radiation fraction that would be needed to make the design point feasible.

Breeding blanket
To model the lithium blanket in a fusion reactor, PROCESS contains an helium-cooled pebble bed (HCPB) model developed at Culham Centre for Fusion Energy (CCFE) [9] and an HCPB model developed by Karlsruhe Institute of Technology (KIT) [42]. For the CCFE HCPB model, the energy deposited in the armour and first wall, blanket and shield are calculated using parametric fits to an MCNP neutron and photon transport model for a sector of a tokamak. The blanket contains lithium orthosilicate (Li 2 SiO 4 ), titanium beryllide (TiBe 12 ), helium and Eurofer steel. The energy multiplication by nuclear reactions in the blanket is given as 1.269. The KIT HCPB model allows for the energy multiplication factor, shielding requirements and TBR to be calculated self-consistently with the blanket and shielding materials and sub-assembly thicknesses. It also allows constraints to be set to meet engineering requirements. The blanket is split into subassemblies: the breeding zone, box manifold and back plate. Three breeder materials can be selected from: lithium orthosilicate (Li 4 SiO 4 ), lithium metatitanate (Li 2 TiO 3 ) and lithium zirconate (Li 2 ZrO 3 ). Together, the three sub-assemblies make up the total blanket thickness. Constrains can be set on the TBR, maximum allowed toroidal field (TF) coil fluence, maximum allowed heating of the TF coils and/or the maximum allowed helium concentration in the vacuum vessel. Through these constraint, the code can determine the thicknesses of the sub-assemblies and the overall blanket thickness.
For now, we assume these models to hold to first approximation also for stellarator devices. However, in contrast to tokamaks, stellarators can have significant variation of the neutron wall load in toroidal direction. This can be accounted for when adding a neutron peaking factor f peak to the models, which measures the inhomogeneity of the neutron load along the blanket area. For a given configuration, this factor can be calculated by Here, q max is the maximum and q avg the average neutron load in the blanket. When one constructs an intermediate, first wall like, hyper-surface between plasma and coils, one can approximately calculate q on this surface via Here, θ and φ are poloidal and toroidal coordinates on the surface, x S and x W are the position vectors of the source and the wall respectively, V S stands for the volume of the source and n is the normal vector of the wall. E n is the energy carried by a neutron in a D + T reaction (14.1 MeV). f S is the neutron fluence at the source point x S , which can be obtained using the Bosch-Hale fit [43] for a reference density and temperature profile, θ and ξ are fit functions and C 1 is a fit parameter, see [43] for their explicit form. An example calculation of the neutron wall load using equation (33) for a wall in a Helias 5 device is shown in figure 4. Equation (33) simplifies the geometry vessel by neglecting 'shadowed' regions in the vacuum vessel and it further does not account for neutron scattering, but it is a method to compute the peaking factor f peak computationally fast. More sophisticated values for f peak can be obtained [44] by dedicated 3D Monte-Carlo codes such as MCNP [45], which can include neutron scattering and further are able to resolve in detail vessel and blanket geometries at the cost of computational time. Equation (33) can be substituted with results from an MCNP run in Process, if more accuracy is needed.
The effect of the neutron inhomogeneity was implemented in the HCPB models in PROCESS now, using a calculation of f peak in the pre-processing step.
Future improvements of this model should replace the used tokamak-specific blanket models by stellarator specific models based on stellarator reference calculations, as conducted e.g. in [46].

Stellarator coils
For a given, averaged, toroidal magnetic field strength B t along the magnetic axis, Process should calculate the required coil current in the pre-defined coil filaments. This is achieved by using a simple linear scaling from a pre-calculated value for the averaged norm of the TF along the magnetic axis, B t axis , which can be obtained by integrating along the magnetic axis, where is the length of the magnetic axis, B t the magnetic field on the axis and s is a coordinate parameterizing the axis (not to be confused with the flux surface coordinate). Once determined for a reference point, the scaling of the coil current with respect to B t and R is of course linear, The needed coil current I 0 (C) for the respective B t at the reference design point can be calculated using the Biot-Savart equation, which is done numerically in the pre-processing step. The (vacuum) axis can be obtained by a field line tracer, e.g. [47], or as output by the equilibrium code VMEC [48].
Another important parameter for the coil design in a systems code is the maximum magnetic field on the coil surface B max , which is crucial for the superconductor material constraints. B max depends on the coil cross-section area and for its calculation at the reference design point withR,B, and the winding pack thicknessÂ WP we proceed as follows.
For stellarators in Process, and for the calculation of B max only, we approximate the winding pack to be of rectangular shape and to be homogeneously filled with a current carrying material. With these assumptions, Biot-Savarts volume integral can be in good approximation reduced to a Riemann sum of analytically solvable integrals of the magnetic field due to homogeneously filled straight cuboid beams [49,50]. For reasonable accuracies, each coil is discretised into O(100) straight beams, each producing a magnetic field B Beam i at position x. The total contribution of a coil to the magnetic field at a position x can then be approximated by The derivation and an explicit formula for B Beam i is given in appendix A. B max then becomes This descriptions allows, for our purposes, an accurate calculation of the magnetic field at the surface of the coils and in the current carrying material. The latter will be important for the force calculations that will be described further below.
B max depends on the winding pack cross-section. To reflect this scaling in the systems code, we calculate equation (38) for varying winding pack sizes in the pre-processing step and parameterized B max in Process via a fit function, which we choose here in the form of The first summand approximates the ideal part (due to an ideal toroid), the second summand includes the fitted scaling with changing winding pack size. a coil is the average minor coil radius, N the number of coils, and A wp the cross-sectional area of the winding pack. a 0 and a 1 are fit parameters that are obtained in the pre-processing step by varying A wp . The electromagnetic forces that act on the coils are important output and constraint parameters, as the integrity of the structural material is limited by the stress, which again scales with the force magnitude. This fact is especially limiting for compact devices at higher magnetic field, as those typically imply high operating current densities resulting in high force-magnitudes. The force density, as other effective parameters before, is calculated for a reference coil size and then scaled within Process. For this purpose, the magnetic field B is calculated inside the winding pack, using the finite winding pack Biot-Savart approximation introduced in equation (37). The Lorentz force density at a point x in the winding pack is then simply if the magnitude of j, the current density, is assumed to be constant and homogenous across the coil cross section and points along the tangential direction t of the coil. Figure 5 shows an example calculation of the force density distribution in a stellarator coil: in every poloidal cross section of the coil we discretize the winding pack cross-section into N × N volume elements dV for which we calculate a force density f using equation (40). f can be integrated over A wp to obtain a force densityf in N m −1 , or over the whole coil volume V coil to obtain a force F in N. The maximal value of f needs to be supported by the structural material in the winding pack,f result in coil jacket and coil insulation stresses and F is relevant for the outer coil support structure.
We calculate the effective parameter as the maximum of each of these forces according to f max (C) ≡ max θ,i |f| (θ a poloidal coil coordinate, i indicates the coil number) for every configuration and scale it in Process according to Here, j is the current density, I the coil current, the length of the respective coil (in Process). Hatted values again denote the values at the reference point where f max (C),f max (C) and F max (C) are calculated.
It should be noted that one needs to make an assumption about the orientation of the winding pack in order to calculate the force density. To this end, we choose the normal vectors of the winding pack to point into the cylindrical toroidal and radial direction respectively. In a realistic winding pack, which is optimized with respect to torsion and stresses, this normal vector might deviate from this assumption, however, f max will most likely not be affected significantly by this choice.
As stellarators can have significant lateral forces, Process also returns lateral and radial projections of equation (40) which are scaled analogously to equation (42). Figure 6 shows the order of magnitude of lateral projection off in a Helias 5 coil set.
To estimate the stress on the ground insulation of a coil set we use a simple model and only consider normal uniaxial stresses which depend on the poloidal coil coordinate θ, namely We assume that the forces F(θ) point orthogonal towards the outer boundary of the coil and thus create a pressure on the radially outward area of the coil A, which depends on the winding pack size. Assuming a fixed outer coil boundary condition, the maximal stress on this area, induced by the winding pack forces, then is where d WP is the radial thickness of the winding pack as calculated by Process from equation (48). This stress is subject to the elastic limit of the material under pressure. If a coil design as in [51] is assumed, this stress exerts on the ground insulation and its upper limit will be in the order of ∼ 100 MPa. In our implementation, the maximum allowable stress is a user defined parameter and if set, Process will optimize the design to fulfill this constraint.
It should be noted that we ignore stresses in the coil structural material for now, as accurate values for the peak stresses would require a detailed design of the coil support structure. Possibly, some simplifications of the support structure could be made, like a thin massive inter-coil shell, which could provide an idea about stresses in the coil support structure with the help of finite element calculations, but this is beyond the scope of this paper.
Another stellarator-specific output parameter of the coil module is the maximal curvature in the coils. This parameter is especially relevant for stellarators as the non-planar coils can have small bending radii that might not be in line with limitations imposed by the superconductor material. Again, in Process, the maximal curvature is implemented by a scaling equation, using a reference value κ max has been obtained in the pre-processing step, Here, d WP is the radial thickness of the winding pack. The term 1 estimates the curvature increasing effect of a radially extended winding pack. The reference value for the maximal curvature κ max (C) is calculated in the pre-processing step according to where γ i : I ⊂ R → R 3 parameterizes the ith coil in the set and θ is a local coil coordinate. κ max can be used to model the bending strain in the superconductor, which has direct implication on the critical current density of the superconductor. A bending strain model based on equation (47) is not yet implemented in Process for stellarators, but can easily be added, once respective models are available.

Winding pack design
For tokamaks, Process is capable of optimizing the winding pack constituents (copper and superconductor fractions) with respect to the figure of merit. In [17] this degree of freedom was not implemented for stellarators, which we now enable using the following prescription. For the stellarator version of Process, we model the winding pack with N squared turns, surrounded by a coil jacket and some user defined ground insulation thickness on this coil jacket. Each of the N turns has a composition as shown in figure 7. The inner part of the conduit contains an approximate squared conductor area. The structure and helium fraction as well as the insulation thickness in the conduit cross section are user defined parameters, whose values are subject to external specifications. Especially the fraction for the structural material needs to match the inner winding pack stress constraints, which are non-trivial in 3D coils and require a sophisticated treatment. The copper-and superconductor fractions, in contrast, are subject to quench protection and can be calculated by Process, as will be addressed later in this section. The overall dimension of the turn area is a user defined parameter.
For stellarator coils, Process now optimizes the copper and the superconductor fractions according to the consistency equation Here, f j 1, is an iteration parameter and is bounded by user defined values. j crit is a parametric form for the critical current density of the superconductor which depends on T, the temperature in the superconductor, B max (A wp ) as given from equation (39) and the maximal strain in the superconductor. Currently, the implemented superconductor material parameterizations in Process cover Nb 3 Sn, NbTi, Bi-2212 and a REBCO-material [9].
The superconductor fraction f scu in the winding pack is a resulting parameter from the winding pack material area fractions, where f case is the case and insulation fraction of the whole turn area, f He is the helium fraction in the conduit area and f Cu and f oth are copper and other material fractions in the conductor area. Process finds the appropriate winding pack dimensions then by solving equation (48) for A wp , which is a simple root finding problem and is solved by Newton's method within Process. In equation (49), f Cu is an iteration parameter in Process and is bounded by quench protection arguments, which we will address below.
In the case of a coil quench, the internal TF coil current needs to be dumped into external resistors. The exponential decay time of the coil current during the quench is parameterized in Process by τ Q . This value is an iteration parameter, subject to the constraints: The first constraint restricts τ Q by the maximal allowable voltage across a coil and during a quench which is, for large resistances, approximately given by [9] E stoTF is the approximative average stored energy per coil, L the inductance of the coil set, N TF the number of coils, and I is the average coil current. The inductance of a stellarator coil set is calculated in the pre-processing step (e.g. by assuming a filamentary 3D curve approximation of the coils [52,53]) for a reference point and can be scaled in Process according to This equation is based on an ideal toroid, where a coil is the minor average coil radius. The restriction for τ Q is then The second constraint for τ Q due to the temperature rise during a quench can be quantified using an energy conservation argument leading to In appendix B we provide a short derivation of this equation.
Finally, the third constraint considers the fact that the changing current in the coils during a quench induces a stress in the vacuum vessel. The maximum allowable force density in the vacuum vessel during a quench f VV puts another lower bound on τ Q . We use a scaling equation to calculate the maximum force density based on a reference value according to where d VV is the vacuum vessel thickness, R VV the (approximate) major radius of the vacuum vessel and B the average toroidal magnetic field on axis. For now, we choose a sophisticated Ansys simulation from W7-X as a reference value as illustrated in figure 8, where 2.54 MN m −3 is the maximum value of the force density. Note that this step is not done in every pre-processing step, but instead is only provided once for the W7-X vacuum vessel. Due to lack of available models for generic 3D vacuum vessel, we assume for now that, in first approximation, this value also reflects the general inhomogeneity for any type of stellarator vacuum vessel. However, the reference value can be easily adapted for designs where more detailed simulation results exist. With values from W7-X, equation (54) becomes Also note that this constraint could in principle be overcome by a poloidal electric break, e.g. [54]. In Process, f VV is then bound to a user defined parameter and serves as an inequality constraint.

Structure mass
As shown in the previous section, large lateral forces can act on the non-planar stellarator coils. However, the details of the force distribution depend very much on the coil shapes and winding pack. This puts not only great demands on the support structure, but also makes it difficult to design an appropriate structure. Consequently, such designs for large stellarators are scarce. There exist only a few design concepts for a stellarator reactor, such as a bolted or welded plates [19] or support elements with 'stiffeners' [55]. Instead of implementing a specific design in Process, we choose to model only the total structure mass, while not being sensitive to the details of support structure. The total mass is a good proxy, both for the cost and the support structure complexity. As introduced already in [16], we stick here to an empirical scaling law from existing machines, as described in [56] to calculate the structure mass in Process based on magnetic energy W mag in the coil-set, Although equation (56) sees good empirical agreement, it does not show whether the design point has local unsupportable forces. In reality, the optimisation of the support structure is a difficult task to ensure the integrity of the device while avoiding local overloads.

Build consistency and port sizes
Scaling in R and the winding pack requires that Process checks the inner coil-coil distances in toroidal direction to prevent that coils come too close in azimuthal direction. We incorporate this constraint via an effective parameter of the minimal distance between two central coil filaments d min (C), which is calculated in the pre-processing step. This distance scales linearly with the major radius and is subject to the constraint where w WP denotes the toroidal width of the winding pack as calculated by the routine described in section 3.8 and w case is the implied coil casing width in toroidal direction. Furthermore, the radial distance between the plasma and the coils is also subject to build constraints. For stellarators, the most critical location is the point where the coils come closest to the plasma. One value for this distance at a reference device size is calculated in the pre-processing step and defines an effective value as d pc (C). In Process we then implement the scaling Here, f geo = ∂d pc ∂a accounts for how much the plasma wall distance changes when decreasing the minor radius in the same configuration. A is the (scaled) aspect ratio andÂ the aspect ratio at the reference point.
In Process, d pc is then subject to the constraint where d coil is the radial thickness of the coil (winding pack plus coil jacket and insulation), d VV is the thickness of the vacuum vessel, d shield of the thermal shield, d blanket the thickness of the blanket, d fw the thickness of the first wall and d SOL describes the width of the scrape-off layer. g ap accounts for the left available space.
Note that by this prescription, PROCESS only ensures radial build consistency along one radial line in the stellarator geometry and in general the gap g ap is a function of a poloidal and toroidal angle, g ap = g ap (φ, θ). Equation (59) is implemented via a stellarator specific inequality constraint in Process.
Finally, we calculate a maximal rectangular vertical port size area A max Port (C) in the pre-processing step for a reference point. Each dimension is then scaled linearly with the major radius within Process. The maximum port size limits the maximum size of blanket segments and is thus an important information to judge the feasibility of remote maintenance.

Concluding remarks
We listed the implemented changes in which Process' prescriptions of a stellarator power plants now differs from the tokamak prescription. For this, we identified important reactor relevant stellarator-specific features and implemented them to sufficient accuracy in Process using an additional pre-calculation step. However, there are more stellarator specific constraints in a power plant which are not included yet. For example, alpha particle damage on the wall and inhomogeneous radiation loads are approximated by the (axi-symmetric) description of Process. Proper stress and strain calculations for stellarator devices are ignored for now in Process. Capturing these Figure 9. The used tokamak DEMO TF-coil set for the comparison (output of tokamak-Process). The winding pack cross section shape is simplified as rectangular in stellarator-Process. modifications require more detailed calculations and stellarator design studies and need to be added in future publications.

Application
In this section we apply the modified Process code to three different scenarios: first, we carry out a benchmark of the newly developed stellarator coil module against the (established) tokamak coil module. Secondly, we apply Process for the first time exemplarily to three distinct stellarator configurations with different aspect ratios. This is possible only due to the newly developed models. The scenario demonstrates a possible use-case of Process to stellarator plasma optimization, as it allows to find feasible reactor design points (in terms of outer radius, aspect ratio, density or temperature). Further, Process can help to identify limiting constraints for a given stellarator reactor design point, which we demonstrate in this study too.  . ∼3000 Process runs, scanning the major radius against the toroidal magnetic field strength for a 1 GW net electricity Helias 5 configuration when neglecting the ECRH constraint (equation (22)) to visualize three of the newly implemented constraints. Green points indicate a valid solution. Blue points were allowed when the coil quench constraint was relaxed, yellow points were accessible if the blanket could be build more compact or larger coil-plasma distances were possible, red points are accessible by improving the confinement.
Thirdly, we use Process to compare three different stellarator coil-sets for the same magnetic configuration to demonstrate the capability of Process to provide input for stellarator coil optimization, as it is possible to provide necessary coil-coil and coil-plasma distances, both parameters which depend on technology and are used as inputs in stellarator coil optimization.
The example studies shown below aim to demonstrate the new capabilities of the approach, highlighting in particular the impact of engineering constraints on the design space. A fully detailed reactor design study is subject to a future work.

Benchmark against tokamak PROCESS TF coil module
The stellarator models were designed to accommodate any type of stellarator. This flexibility allows to model also a tokamak-coilset within the stellarator-Process version.
In this section, we briefly benchmark the results of the new stellarator coil module in Process against the output of the independent tokamak-Process TF coil module. For comparing the models developed in section 3 to the tokamak models, the reader shall be referred to [9].
For the benchmark, we start from 16 D-shaped TF coils, as shown in figure 9. The coil shapes are produced by tokamak-Process and are then taken as input for the stellarator run. We now obtain effective parameters a i (C) for the coilset in the pre-processing step as described in the previous section and then run stellarator-Process in optimization mode, optimizing for capital costs, which is equivalent to minimizing material masses. We fix the magnetic field strength on axis to 5.72 T and the aspect ratio to 3.1 and let Process find a consistent design point while optimizing for engineering parameters, such as the copper fraction in the conductor area, the winding pack dimensions, and the exponential coil quench dump time.
The result of the benchmark is displayed in table 1. Stellarator-Process converges to a similar design point as the tokamak-Process version. The winding pack dimensions and the copper fraction are optimized to similar values. The maximal magnetic field on the coils deviates by 5%, which is within the model accuracy. Generally, we find very good agreement of our stellarator coil model with the tokamak case, providing confidence in the implementation of the developed coil model.

PROCESS for stellarators with different aspect ratios
Historically, stellarator reactor design studies were performed by individual calculations for single magnetic configurations. Such an approach is tedious for the large magnetic configuration space of stellarators. Our models allow for the first time to model different types of stellarator configurations within the same code framework within seconds of computational time. To demonstrate this capability, we showcase below exemplary studies for three different stellarator configurations, Helias 3, Helias 4 and Helias 5 as introduced in [19,20,57], with a field periodicity of 3, 4, and 5 respectively, as shown in figure 10.
For each coil-set we calculate a vacuum VMEC free boundary equilibrium and determine the effective parameters as described in the previous section. We then run Process in optimization mode where we optimize for minimal major radius at constant aspect ratio and a required net electricity output of 1 GW, which corresponds to approximately 3 GW fusion power. For this, we optimize the following parameters: the overall temperature and the overall density for fixed profile shapes, α T = 1.2, α n = 0.35, as defined in equation (13), following neoclassical transport simulations conducted in [16]. Further, we optimize the major plasma radius and the overall magnetic field strength. In addition, Process optimizes for coil current densities, winding pack dimensions and material fractions.
These optimization parameters are bound by several imposed constraints: for the radial build constraint, a fixed radial component thickness of 1.15 m is assumed, including vacuum vessel, breeding zone, blanket structure mass and neutron shielding, consistent with neutronics calculations conducted in [46]. The SOL width is taken as 15 cm. For every configuration we assume that 80% of the born fusion alpha particles heat the plasma. This value is expected to increase in future stellarator devices, as improved alpha particle confinement in stellarators is only recently addressed in stellarator optimization, but with promising results already [58,59]. We further impose an ECRH heated ignition point, using the prescription in subsection 3.4, and assume maximal available gyrotron frequencies of 200 GHz. For comparison, ITER operates with 170 GHz gyrotrons [60]. Requiring ECRH inhabitability constrains both, the density and the magnetic field from above. The ISS04 transport model as in section 3.2 is assumed and the 0D power balance equation including Bremsstrahlung and line-radiation terms is enforced as a constraint equation. The superconductor material is taken as Nb 3 Sn at 4.75 K operation temperature and the current density we limit to 80% of the critical superconducting density. Superconductor Table 2. A selection of Process' output parameters for the converged design point for each of the three Helias configurations. The design points were optimized with respect to minimal major radius and for a net electricity output of 1 GW, which approximately corresponds to 3 GW fusion power. strain is neglected for the critical current density. All devices assume an island divertor, which is described by the model in section 3.5. In this model, a radiation fraction of 85% in the SOL is assumed and radiating plasma impurities are neglected. The maximum allowable divertor heat load is set to 10 MW m −2 and a volume averaged upper beta limit of 5% is imposed. Coil current densities, winding pack dimensions and material fractions are subject to quench restrictions and ground insulation stress as described in subsection 3.8.
To visualize the restriction of the parameter-space by the newly implemented constraints, we conducted an R-B-scan of the Helias 5 device in figure 11: from this, one can identify the coil quench constraint, which forbids points at high Bfield and smaller major radius. Further, the new radial build consistency model through imposed blanket and shielding requirements dominantly rules out design points at smaller major radius and an increase in confinement properties (an enhanced ISS04 proportionality factor) would allow to access the region at lower B-field. Note that for this scan we neglected the ECRH constraint, which is very sensitive on the technological assumptions (gyrotron frequency and heating scheme).
The exact positions of the constraints also depend on other factors like imposed steel fraction in the conductor or insulation thicknesses. An eventual hoop or inner-winding pack stress limit would further limit the design space, but this was not developed yet for stellarators.
Using the above listed set of optimization parameters and constraints, important output parameters of the optimized reactor design points with respect to minimal major radius are shown in table 2. The study in [20] assumed NbTi superconductors, which we replaced by Nb 3 Sn superconductors here. This allows for higher field strengths as the found 7 T on axis for Helias 5. The possibility to switch superconductor material also demonstrates the advantage of technological flexibility of the systems code framework. The found densities in table 2 are in line with the ECRH heating constraint as described in subsection 3.4. Despite of the different aspect ratio and the different coil numbers, the total masses found by Process are comparable for all three devices and only increase slightly for the machine with higher aspect ratio. Comparing the coil masses of the Nb 3 Sn Helias 4 design against the previous NbTi study [57], we find good agreement (6.8 kt compared to 7.5 kt) for a similar design point.
The major radius of all three designs are limited by the plasma-coil distance, which needs to include space for blanket and shielding and the radial extension of the coils. This we found by relaxing the radial build constraint, which resulted in significantly smaller major radii for all three machines (also compare figure 11). This fact also explains the significantly larger plasma volume of 1560 m 3 of Helias 4 compared to the other two devices, as Helias 4 features a comparably small plasma-coil distance, namely 1.7 m at reference size compared to 1.9 m for Helias 5 at reference size [20]. We further observe a relevant design restriction by the imposed ECRH constraint, mainly given by the critical O1-mode density limitation. This provides further motivation to develop high power, high frequency gyrotrons or an X1 heating scheme for a stellarator reactor at high density.
While our models present substantial progress in the engineering feasibility of stellarator designs, there are still missing factors that are under development and not yet included in this study, such as high fast particle wall loads, superconductor strain or transport effects that are not covered by the ISS04 assumption, or stress limits in coils and in the structure mass. Modelling these effects for general stellarators is beyond the scope of this work, and can be added in future publications.

PROCESS for stellarators with different coil sets
Stellarator coil optimization is, at least traditionally, carried out for a fixed magnetic configuration. For every configuration, however, there exists an infinite number of coil-sets producing (approximately) the same magnetic field [62] and choosing the right coil set is a trade-off between field accuracy and engineering constraints, such as the minimal curvature, port sizes required by remote maintenance, coil-coil distance, coil-plasma distance, engineering tolerances [63], or costs. In this section we demonstrate that stellarator-Process can help judging the reactor relevance of the coil-set by providing further details, as its material usage and forces, by including coil quench constraints and by considering other plant constraints at the same time. For this purpose, we generate three exemplary coil-sets targeting a W7-X like equilibrium, using the coil optimization code FOCUS [61]. The chosen coil-sets have 30, 50, and 60 coils respectively, and their corresponding Poincaré plots are shown in figure 12. Albeit similar flux surfaces and island positions compared to the W7-X equilibrium, further physics properties of the respective equilibriums were not checked here, as our purpose here is just an exemplary application of Process to different stellarator coil sets for the same equilibrium. The coil-set with 30 coils was not able to match the target iota at the boundary, which results in a lack of the desired island structure there, but for the sake of having a coil-set with significantly less coils, while still retaining a significant coil-coil gap, we neglect this fact and include this coil-set in the following study.
We use Process now to scale the overall size of the machine to reactor size (22 m major radius). The plasma size is then held fixed and the new Process implementation for stellarators is used to optimize for capital costs, which is equivalent to minimizing the coil masses in this case. Process then finds the required coil sizes, the copper and superconducting material fractions, assuming Nb 3 Sn superconductors, to match the build constraints in radial and toroidal direction, the coil quench protection constraints and the superconducting critical current density constraint. Relevant coil related Process output parameters are shown in table 3. As expected, local coil forces are significantly larger for a design with less coils. The stored magnetic energy scales with the coil minor radius which is, approximately 20% larger in the 30 coils device. The design with less coils allows for larger vertical ports, which could ease remote maintenance, however this seems to come at a cost, as the design with 30 coils is found at significantly higher total coil masses compared to the design with 50 or 60 coils, likely induced by significantly higher maximal B-field values at the coils which then again requires larger material masses to fulfill the critical superconductor current density and the quench constraints. Finally, the left over coil-coil gap vanishes for the device with 50 and 60 coils, which indicates that these two designs would benefit from larger coil-coil distances. This information can be used to re-iterate the imposed coil-coil distance in the coil optimization step to obtain coil-sets which feature larger coil-coil (and coil-plasma distances) to allow for a smaller overall machine size.

Summary and outlook
In this work we presented modifications of the fusion reactor systems code Process to model a general class of stellarators. For this, we modified Process in a way that it covers several stellarator specific features of a fusion power plant. Some stellarator specifics, such as an accurate description of the alpha particle wall load in stellarators, the inhomogeneous plasma radiation load on wall materials or a 3D stress model are left out for future publications.
As a result of the modifications, Process can now be used to obtain stellarator reactor design points that are, within the implemented model coverage, in line with current technology, taking as input solely the common output of stellarator optimization, namely a plasma boundary and the respective coil set. This new code modification allows for the first time to compare different stellarator configurations within the same wholistic systems code and thus can contribute to stellarator optimization, as it can help constraining the high dimensional design space, that stellarator plasmas and coilsets allow for. Process further allows high level optimization of design parameters with respect to economical figure of merits, such as component masses, which can help guide future stellarator reactor studies, as it allows a fast adjustment to new technological advances. The new framework further calculates coil thicknesses that are in line with superconductor and coil quench constraints, which again can be used to re-design coils.
The implementation of the presented models was done in two frameworks, in a pre-processing code and in Process itself, which allows to achieve rapid reactor design points with Process within seconds, once the effective parameters are computed. This timescale makes Process, or some of its implemented submodules, suitable to include in stellarator optimization routines in the future. We demonstrated applications of the new Process code for general stellarators in three use-cases: first, we benchmarked the results of the newly implemented coil model against the (independent) tokamak description of Process. We not only found a similar optimized design point but also sufficient agreement in the relevant parameters themselves. Secondly, we applied Process to three previously found configurations [20], and obtained an example reactor point with more detailed physics and engineering parameters which are in line with the newly implemented constraints. In the third application, we demonstrated that the technological constraints implemented in Process can be used to provide insights in important input parameters for stellarator coil optimization, such as the coil-coil distances or the coil-plasma distance in the coil-set, which are subject to non-trivial material constraints, as superconductor properties or coil quench considerations.
From our studies conducted here, we find that the major radius for nearly all examined machines is limited by the required blanket space and that coils situated further away from the plasma would likely be beneficial for reactor designs with smaller major radius for these configurations. In other words, the major radius of the used devices, and thus the major cost driver of a stellarator reactor device, was consistently found to be constrained from below by the plasma-coil distance, not by lack of confinement quality. To obtain valid compact designs, a major focus for coil and plasma optimization could lie on finding designs that allow for large coil-plasma gaps, while still retaining a tolerable field error and an acceptable coil-coil gap to allow for finite size coils. Overall, by the studies in section 4 we demonstrated that design window analyses of stellarator devices with different plasma shapes, TF periods, number of coils and coil shapes have become possible within Process.
The Process code is maintained at the CCFE in Culham, Oxfordshire, UK (A description of the code can be found here: https://ccfe.ukaea.uk/resources/process/.) The 'pre-processing' step uses the in section 3 presented calculations and was automatized and implemented as a python tool, which is maintained at the Max-Planck-Institute for Plasmaphysics in Greifswald, Germany.

Appendix A. Biot-Savart with finite conductor size
Here we derive the magnetic field B at a point p due to a current carrying rectangular cuboid (beam) as it is used in equation (37). The cuboid and used conventions in the following is shown in figure 13.
When a 3D stellarator coil is approximated by N such beams, this procedure allows a fast evaluation of the magnetic field near and, very useful for force calculations and superconductor constraints, within the conductor. This method was also used in [49].
Let b be the vector in longitudinal (y-) direction of the beam, while n points in normal (x-) direction. Define the functions: x p − x N(x − x p , y − y p , z) 3

(A.2)
N(x, y, z) = z 2 + y 2 + x 2 (A. 3) where x p are projections according to: x p = p · e x . 2b is the dimension of the beam in x and 2d in y direction. If the current density j in the winding pack is approximated as a continuous constant function across a rectangular cross section, pointing w.l.o.g. in Cartesian z direction, Biot-Savart's volume integral can be written as: The integral over F 1 and F 2 have an analytical form then, as it is shown below.
For convenience, define (A.6) Then (A.7) And analogously for F 2 it is (A.8) This simplifies equation (A.4) to a one dimensional integral along the z-direction, which can be solved numerically. However, using equation (A.6), the integral in z-direction can also be solved analytically, and the magnetic field B can then be written as

Appendix B. Quench protection
We shortly provide the derivation of the critical current density as limited by a simple coil quench protection argument as given in the final form in [9].
In thermal equilibrium and without losses the heat produced by the copper resistivity during a quench is equal to the heat needed to rise the temperature in the material by dT, dQ heat = dQ temp . Now, the quench restriction is to impose The integral on the left-hand side runs over the whole quench time while the integral on the right-hand side goes from the operation temperature T op to a maximal T max . The difference T max − T op is usually chosen in the order of 150 K. If one assumes an exponential decay of J after a quench detection time t d as: Using the definition of the relative winding pack material fractions f as in equation (49)  With this, one ends up with (identifying J 0 with the copper current J cu ) In terms of the total winding pack current density, equation (B.13) can be rewritten using 1 − f He = f cond and J WP = J Cu f Cu f cond (1 − f case ): Equation (B.14) constraints the winding pack current density by a temperature rise during a coil quench. This value is dependent on the chosen copper alloy, which enters in η and c i .