{'fraction_non_alphanumeric': 0.03841430775528197, 'fraction_numerical': 0.01640479863975566, 'mean_word_length': 4.439667722874026, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
KNOSOS: a fast orbit-averaging neoclassical code for arbitrary stellarator geometry

KNOSOS (KiNetic Orbit-averaging SOlver for Stellarators) is a freely available, open-source code that calculates neoclassical transport in low-collisionality plasmas of three-dimensional magnetic confinement devices by solving the radially local drift-kinetic and quasineutrality equations. The main feature of KNOSOS is that it relies on orbit-averaging, which removes the dependence on the coordinate along the magnetic field line, and allows to solve the drift-kinetic equation very fast. KNOSOS treats rigorously the effect of the component of the magnetic drift that is tangent to magnetic surfaces, and of the component of the electrostatic potential that varies on the flux-surface, {\varphi}_1. Furthermore, the equation solved is linear in {\varphi}_1, which permits an efficient solution of the quasineutrality equation. As long as the radially local approach is valid, KNOSOS can be applied to the calculation of neoclassical transport in stellarators (helias, heliotrons, heliacs, etc.) and tokamaks with broken axisymmetry. In this paper, we show several calculations for the stellarators W7-X, LHD, NCSX and TJ-II that provide benchmark with standard local codes and demonstrate the advantages of this approach.


Introduction
Stellarators are non-axisymmetric devices in which the magnetic field is created basically by external magnets, without the need of any mechanism to drive current within the plasma. This provides them with an inherent capability for steady state operation and makes them less prone to plasma magnetohydrodynamic instabilities, but it also generally produces larger energy losses: at low collisionalities, the combination of magnetic geometry and particle collisions leads to a variety of stellarator-specific neoclassical transport regimes, which usually give a large contribution to the radial energy and particle transport in the core of the device (Dinklage et al. 2013;Dinklage & the W7-X Team 2018). Of special relevance are the 1/ν, the √ ν and the superbanana-plateau regimes (Ho & Kulsrud 1987;Beidler et al. 2011;Calvo et al. 2017), in which the energy transport coefficients show a positive temperature dependence, much more unfavourable than the negative T −1/2 scaling of the banana regime of the axisymmetric tokamak. The fundamental reason for this behaviour is that in a generic stellarator, unlike in an axisymmetric tokamak, trapped particle orbits have non-zero secular radial drifts. The exception are omnigenous stellarators: in these magnetic configurations, the secular radial drifts vanish (Cary & Shasharina 1997;Parra et al. 2015), and the level of neoclassical transport is low, similar to that of the tokamak. Quasisymmetric stellarators are a particular family of omnigenous stellarators, see e.g. (Landreman & Catto 2012).
The two world's largest stellarators in operation, Wendelstein 7-X (W7-X) (Klinger et al. 2017;Wolf & et al. 2017) and the Large Helical Device (LHD) (Takeiri & the LHD team 2017), have relied on optimization of neoclassical transport for their design and operation. The magnetic configuration of W7-X has been designed to be close to omnigeneity with poloidally-closed contours of the magnetic field strength; one of the goals of the project has been to prove the constructability and reliability of such design (Sunn- Pedersen et al. 2016). In LHD, the plasma column can be shifted inwards so that the minimum values of the magnetic field along the field line have approximately the same value (see figure 2 of )), a well-known geometric property of some omnigenous magnetic fields (Mynick et al. 1982;Landreman & Catto 2012); discharges performed using this magnetic configuration consistently show better energy confinement (Yamada et al. 2005). Finally, a particular kind of quasisymmetry, quasiaxisymmetry, was the design criterion of the National Compact Stellarator Experiment (NCSX) (Zarnstorff et al. 2001). Power reactor designs exist for these three stellarator concepts (Sagara et al. 2010).
It is then clear that optimization of neoclassical transport is a crucial issue towards a stellarator reactor. One of the most common goals of stellarator optimization efforts is the minimization of the so-called effective ripple, a figure of merit that provides information of the level of transport in the 1/ν regime. While there is little doubt that minimization of this quantity should be a design criterion in any future stellarator, it has important limitations. On the one hand, empirical studies of the energy confinement time of several devices aimed at obtaining a unified International Stellarator Scaling law (ISS04) have not shown a very strong correlation between reduced effective ripple and improved energy confinement (Yamada et al. 2005;Fuchert et al. 2018); on the other hand, self-consistent neoclassical transport simulations performed in the configuration space of W7-X, complemented with simplified anomalous modelling (accounting for nonnegligible turbulent contributions to transport), have shown mild increases of the energy confinement time for configurations of significantly reduced effective ripple . This points towards one of the obvious limitations of the effective ripple: it is only an appropriate figure of merit for neoclassical transport if the plasma species are in the asymptotic 1/ν regime. However, bulk particles are distributed close to a Maxwellian that typically spans across several transport regimes. Even in cases in which the collisionality is low and the neoclassical predictions of the radial energy flux agree with the experiment, the parameter dependence of the experimental energy flux does not follow the scaling expected for any specific neoclassical transport regime, see e.g. , because the flux is caused by a combination of transport regimes.
The reason for choosing the effective ripple as a figure of merit is that the 1/ν regime is the low-collisionality regime of stellarators in which the effect of the magnetic geometry on transport can be encapsulated in a straightforward manner in a single quantity that is independent of density and temperature, and radial electric field. Furthermore, this quantity can be efficiently calculated by solving the bounce-averaged drift-kinetic equation, e.g. with the NEO code (Nemov et al. 1999). None of this has been possible so far for other low-collisionality regimes for arbitrary stellarator geometry.
Moreover, for other regimes such as the √ ν and the superbanana-plateau regimes, the effect of the electric field (radial and tangential to the flux-surface, the latter associated to the variation of the electrostatic potential on the flux-surface, ϕ 1 ) has to be considered (Calvo et al. 2017), and this quantity is determined by imposing ambipolarity of the neoclassical particle fluxes and quasineutrality, which in turn depend on the plasma profiles, and specifically on the gradients. In order to address this issue, selfconsistent neoclassical transport simulations have started to be performed in the last few years: the neoclassical fluxes are calculated with the DKES code (Hirshman et al. 1986) and then the ambipolar and energy transport equations are solved (the latter with a prescribed energy source) (Turkin et al. 2011;Geiger et al. 2015). Although we will see that DKES makes use of the so-called monoenergetic approximation, which allows to solve the drift-kinetic equation for any regime in a three-dimensional space (instead of a four or five-dimensional one, as the codes mentioned at the end of this paragraph), using DKES to self-consistently solve neoclassical energy transport is still computationally expensive at low collisionality. Moreover, it is inaccurate at sufficiently low collisionality due to some additional approximations that it makes: DKES uses an incompressible E × B drift (Beidler et al. 2007) and does not include the tangential magnetic drift or the radial E ×B drift caused by the variation of the electrostatic potential within the flux-surface (the latter makes the fluxes depend non-linearly on the plasma gradients (Calvo et al. 2018b)). Some or all of these approximations are absent in more recent codes such as SFINCS (Landreman et al. 2014), EUTERPE (García-Regaña et al. 2013García-Regaña et al. 2017) or FORTEC-3D (Satake et al. 2006), but at the expense of higher computational cost.
We have developed a new code, the KiNetic Orbit-averaging Solver for Optimizing Stellarators, KNOSOS, based on the analytical techniques developed in a series of papers (Calvo et al. 2013(Calvo et al. , 2014(Calvo et al. , 2017(Calvo et al. , 2018b. It solves local drift-kinetic equations that will be summarized in the next section and that accurately describe neoclassical transport in the 1/ν, √ ν and superbanana-plateau regimes. The equations include the effect of the magnetic drift tangential to flux-surfaces and the radial E × B drift due to the variation of the electrostratic potential within the flux-surface; the radial electric field E r and ϕ 1 are obtained by imposing ambipolarity and quasineutrality, respectively. Local drift kinetic equations are valid for large-aspect-ratio stellarators or configurations close to omnigeneity (see the discussion before equation (2.23) in §2). Unlike preliminary versions of KNOSOS (Velasco et al. 2018;Calvo et al. 2018b), this version does not require a explicit split of the magnetic field magnitude into omnigeneous and non-omnigeneous pieces. The goal of this code is to be, at the same time, accurate and fast, so that it allows to perform comprehensive parameter scans and to provide input to other codes or suites of codes. Generally speaking, the goal is to improve our confidence in neoclassical predictions, in light of recent theory developments, and to be able to fully exploit these predictive capabilities. To facilitate this objective, the code is freely-available and opensource.
The rest of this paper is organised as follows. §2 presents the drift-kinetic and quasineutrality equations solved by KNOSOS. Then, §3 summarises how the equations are solved: the drift-kinetic equation is written in terms of a few integrals along the magnetic field lines in §3.1, and these integrals are discussed in §3.2; the parameter space and discretization of the drift-kinetic equation is discussed in §3.3 and §3.4, and the consistent solution of the drift-kinetic and quasineutrality equations is presented in §3.5. §4 shows several calculations for real magnetic confinement devices and comparisons with widely benchmarked neoclassical codes: the monoenergetic transport coefficients are compared with DKES in §4.1; the effect of the tangential magnetic drift on the energy flux is discussed in §4.2; the variation of the electrostatic potential along the flux-surface is compared with EUTERPE in §4.3. Finally, §5 summarizes the conclusions. Additionally, there are three appendixes: appendix A discusses the collision operator, and appendices B and C describe algorithms employed to accelerate the calculation of the bounce integrals.

Equations
In this section, we briefly present the equations solved by KNOSOS. Their derivation and further details can be found in previous work by (Calvo et al. 2017(Calvo et al. , 2018b. We first define the coordinate system that we will use. The flux-surfaces are labelled by the radial coordinate where 2πΨ t is the toroidal magnetic flux. The magnetic field lines on the surface are labelled by an angular coordinate where θ and ζ are poloidal and toroidal Boozer angles, respectively, and ι is the rotational transform. Finally, l is the arc-length along the magnetic field line. In these coordinates, the magnetic field B can be written as where primes stand for derivatives with respect to ψ, and Ψ t = ±1 depending on whether the magnetic field is parallel or antiparallel to the direction of the Boozer toroidal angle (i.e. depending on the sign of B · ∇ζ).
As velocity coordinates, we choose the particle velocity v = |v| , (2.4) the pitch-angle coordinate (2.5) and the sign of the parallel velocity For each species b (i will denote bulk ions and e electrons) we need to calculate the deviation of the distribution function from a Maxwellian for trapped particles, that we denote by g b (ψ, α, l, v, λ, σ). The Maxwellian distribution function reads where n b is the density, T b the temperature and m b the mass. Trapped particles are those for which v = 0 at some point along their trajectories. For them, 1/B max λ 1/B min , where B max and B min are the maximum and minimum values of the magnetic field strength on the flux-surface, respectively. The equation for g b (ψ, α, v, λ) is complemented with the boundary condition (2.10) and the condition 2π 0 g b dα = 0 . (2.11) The coefficients of equation (2.9) are integrals over the arc-length between the bounce points l b1 and l b2 , i.e., between the points where the parallel velocity of the particle is zero (see a sketch in figure 1). On the right-hand side of equation (2.9) C lin b [g b ] is the linearized pitch-angle-scattering collision operator: For the ions, since m e /m i 1, this single-species collision operator is correct, but electron-ion collisions need to be retained in the electron drift-kinetic equation. For both species, we follow the common practice (see e.g. ) of using equation (2.12) with an effective collision frequency accounting for inter-species collisions. This is discussed in more detail in appendix A. On the left-hand-side of equation (2.9), is a combination of thermodynamical forces (Z b is the charge number and the elementary charge is denoted by e) and the drift velocity, is the sum of the magnetic and E × B drifts: Here, ϕ is the electrostatic potential, that can be split as ϕ(ψ, α, l) = ϕ 0 (ψ) + ϕ 1 (ψ, α, l) , which means that ϕ 0 and ϕ 1 will be the dominant contribution to the radial and tangential components of the electric field, respectively (and in turn to the tangential and radial components of the E × B drift, respectively). They can be determined by solving two additional equations. The component of the electrostatic potential that varies on the flux surface, ϕ 1 , is obtained from the quasineutrality equation, which for a pure plasma (i.e., composed of electrons and one ion species) reads The sum is done over kinetic species, and λ = 1/B max at the boundary between passing and trapped. Here, we have used that, in terms of our coordinates, velocity space integrals are of the form that g b is even in σ and that g b = 0 for λ < B −1 max . We note that, since ϕ 1 and g b appear in equations (2.9) and (2.18), both equations need to be solved consistently.
The radial electric field is given by the radial derivative of the piece of the electrostatic potential that is constant on the flux surface, where r = a ψ/ψ LCF S , ψ LCF S being the flux-label at the last closed flux-surface and a the minor radius of the device. The radial electric field is set by the ambipolarity of the neoclassical radial particle fluxes, In our variables, where ... denotes flux-surface average. Finally, the radial energy flux is given by KNOSOS solves equations (2.9) and (2.18), together with equation (2.21). These equations have been rigorously derived in (Calvo et al. 2018b) under the hypotheses of low collisionality, large aspect ratio and closeness to omnigeneity (we note that large aspectratio is a common characteristic of real stellarators  while, as noted in the introduction, closeness to omnigeneity is a property sought in present and future devices). At low collisionalities, the motion of particles along the magnetic field is much faster than collisions, and the distribution function does not depend on the coordinate along the field line. Closeness to omnigeneity makes neoclassical transport describable by a radially-local equation for the deviation of the distribution function of trapped particles from a Maxwellian. In particular, it guarantees that the bounce-averaged radial drift is small enough so that even in situations of small E × B drift, and terms proportional to ∂ ψ g do not appear in equation (2.9). Finally, the large-aspect ratio approximation allows us to use the pitchangle collision operator, equation (2.12).
Let us finally discuss the neoclassical regimes that equations (2.9) and (2.18) can describe. The second term on the left-hand side of equation (2.9) includes the radial magnetic and E × B drifts caused by the inhomogeneity of the magnetic field strength and of the electrostatic potential on the flux-surface, respectively. This means that equation (2.9) can model the 1/ν regime and the transport caused by ϕ 1 . The first term of the left-hand side includes the precession tangential to the flux-surface caused by the radial variation of the electrostatic potential (i.e. the radial electric field E r ) and of the magnetic field strength. This implies that we can model the √ ν and superbanana-plateau regimes. As discussed previously, radially global effects are not accounted for.

Solution of the equations
In this section we provide an overview of how equations (2.9) and (2.18) are solved. We first give an explicit expression for equation (2.9) in §3.1 and we discuss how to calculate its bounce-averaged coefficients in §3.2. We then devote §3.3 to build the grid in which we will evaluate the distribution function, and §3.4 to discuss the discretization of the equation. Finally, the solution of quasineutrality, equation (2.18), is addressed in §3.5.

Final expression of the drift-kinetic equation
Using the expressions of the pitch-angle scattering collision operator described in equation (2.12) and of the magnetic and E × B drifts in right handed Boozer coordinates, equation (2.9) can be written in terms of a few bounce integrals: We note that only v d,b , ν λ,b , F M,b and Υ b depend on the species: the bounce-integrals are only determined by the magnetic configuration and the electrostatic potential. The magnetic shear appears explicitly in I v M,α . Equation (3.1) is a differential equation in two variables only, α and λ, which is the origin of the fast performance of KNOSOS that will be demonstrated in §4. The radial coordinate ψ is a parameter, since we are solving radially local equations; v is a parameter as well, since ϕ 1 ϕ 0 ; and finally l has disappeared since the coefficients are bounceaverages of certain quantities. The calculation of these coefficients is described in §3.2.

Calculation of the coefficients of the drift-kinetic equation
The integrals in l are done using a extended midpoint rule (see e.g. Press et al. 1996, subroutine midpnt). This open formula is appropriate for integrals that are improper in the sense that they have an integrable singularity at the integration limits. This is our case, since by definition λB(l b1 ) = λB(l b2 ) = 1. The number of points that we use is not pre-defined: starting from being one, it is tripled until the integral converges.
Let us now note that integrals such as those of equations (3.2) may be difficult to converge if the numerator does not go to zero in the integration limits, or it does, but slower than the denominator. This may happen, first, if λ is such that l b1 or l b2 are close to a point l T where B(l) has a local maximum B(l T ) for fixed α; second, if the interval (l b1 , l b2 ) contains a point l B where B(l) has a local maximum and λ is close to λ B ≡ 1/B(l B ). In such cases, I(λ) may become very large; if the inverse of λ is equal to the corresponding maximum of B, the integral diverges logarithmically. We can physically identify these situations in the example of figure 1: divergences happen at bifurcations, where orbits go from being trapped in a particular region in l to be trapped, for smaller λ, in a wider region (the boundary between passing and trapped particles, is a particular case of this).
One can ease the convergence, and thus make the calculation faster, by removing the divergence and solving it analytically as explained in Appendix C of (Calvo et al. 2017). This is described more in detail in our Appendix B. Additionally, in Appendix C we will discuss how the fact that field lines are straight in magnetic coordinates is used to accelerate the evaluation of the magnetic field strength at each point α, l without loss of accuracy.

Spatial and velocity grid
In §3.2 we have seen how the integrals of equations (3.2) are calculated. These integrals will be evaluated at the points (α, λ) in which we want to determine the distribution function g b . The selection of these points constitute the subject of this subsection.
Let us start with the spatial grid. We have seen that ψ is a parameter, and l does not appear in the bounce-averaged drift-kinetic equation, which leaves us with the field line label α. There are, however, two complications: first, at a given α and λ, several wells may exist (in other words, several pairs of l b1 and l b2 ), which means that we need to use an integer label w for them (as we will discuss more in detail in the following subsection). Second, even if g b does not depend on l, its integrals over velocities (needed e.g. to compute ϕ 1 , see equation (2.18)) do, so we must define a two-dimensional angular grid. As a general rule, when doing so, we try to minimize the number of points at which g b needs to be solved, in order to save computing resources. With this in mind, we make use of periodicity and align the grid points with the field lines. The grid points also aligned with ζ = 0.
More details are discussed along figure 2 (top), which shows precisely one example (for the W7-X case discussed in §4.1) of how the angular grid is built: by following several field lines until they have completed a full poloidal turn. The distance between two consecutive field lines ∆α is taken to be an integer fraction (one sixth, in the plot) of 2πι/N , where N is the number of toroidal periods. This is how the green points are located, with uniform spacing along the toroidal angle. Along the field lines, several maxima of the magnetic field are found, plotted with magenta circles. It is observed that we are dealing with a relatively optimized configuration, in the sense that most trapped particles are so in a major well that coincides with one field period (black continuous arrow). In other words, their bounce points l b1 and l b2 are two consecutive magenta points, separated toroidally by a characteristic angular distance ∼ 2π/N (smaller for large values of λ, close to the bottom of the magnetic well). In the example, several ripple wells are found (grey arrows). For small enough values of λ, trajectories trapped in more than one field-period exist: in this example, particles may move between ζ = 0 and ζ = 6 2π N (black dashed arrow); trajectories with smaller λ (that is, trapped in more than 6 toroidal periods) are ignored in this case; this procedure effectively sets the boundary between passing and trapped particles.
Periodicity allows us to project all these grid points onto the first period. The result is a bidimensional grid in α and l, with N α and N l points in each direction. N α is the integer quantity such that N α < 2π ∆α N α + 1. The N l points along the field line are distributed uniformly in the toroidal angle along a toroidal period, and N l is the largest power of 2 that is smaller than or equal to N α . This will be useful for a fast computation of the Fourier transform, needed when solving quasineutrality. Toroidal periodicity is also enforced at the corners of the grid: for instance, in figure 2 bottom, point α = α Nα−4 , ζ = 0, is not contained in the wells marked in magenta (we have not looked for maxima of B at smaller l). Using periodicity, the value of the distribution function at this point will be taken to be equal to the value at α = α 1 and ζ = 2π/N . The number of points where this has to be done can be minimized by putting one of the corners of the grid close to the global maximum of B on the flux-surface. For each of the nodes of this grid (and for each of the possible values of λ) the points along the trajectory and the bounce points of particles trapped in one or several field-periods are now clearly identified, and the integrals of equation (3.2) can be evaluated.
Let us turn our attention to the velocity grid, where we are using λ and v as coordinates. Since we have seen in §2 that only trapped particles need to be calculated, an obvious choice for the former is a uniform grid †, with N λ + 1 values between λ 1 ≡ 1/B max † When the particles are in the 1/ν regime, special attention should be paid to bifurcations, where g b has discontinuous first λ-derivatives (Nemov et al. 1999;Calvo et al. 2014), and a non-uniform grid, adapted to the structure of maxima and minima at fixed α, is a more efficient choice (Kernbichler et al. 2016). The same applies to very low collisionalities, when the and λ N λ +1 ≡ 1/B min . The distribution function will not be evaluated at λ N λ +1 , which will be ghost points employed for imposing the contour conditions at the bottom. When integrating in λ, we will use the extended trapezoidal rule (Press et al. 1996). Finally, v is a parameter in our calculations: equation (3.1) will be solved for several values v i of the velocity and the solution will be numerically integrated in v. Since the integrand of equations (2.19) contains an exponential coming from the Maxwellian distribution, we will use Gauss-Laguerre of order 64 (Press et al. 1996): being v th,b the thermal velocity of species b, and w i a set of tabulated real numbers. This procedure requires solving the monoenergetic drift-kinetic equation for n = 64 values of v/v th,b , typically from ∼ 10 −2 to ∼ 10 2 . However, the contribution of the largest v i to the integral can tracked and usually neglected, and this allows for an important reduction of computing time. Let us finally note that this is a standard and well-tested choice in neoclassics and gyrokinetics (e.g. Velasco et al. 2011;Barnes et al. 2019), although other velocity-space discretization methods have been proposed in recent years (Landreman & Ernst 2013) that could be easily implemented in KNOSOS.

Discretization of the drift-kinetic equation
In §3.3 we have built a grid in variables α and λ. Three integers can be used to label any point (α i , λ j , w) of this grid: i runs from 1 to N α , j from 1 to N λ and w = I, II... is an integer that labels wells for a given α and λ. At a given point, we define and so on (in order to ease the notation, g i,j,w does not contain a species index). The final step in the discretization of the drift-kinetic equation is the determination of the way in which we approximate the derivatives of g b of equation (3.1) at each point of this grid.
Let us start with the collision operator, which divided by (3.4) and can be expanded into two terms We represent the λ grid at fixed α in figure 3. Here, λ 1 is the boundary between passing and trapped particles. In this example, only one complete well is plotted at λ 2 , labelled I. If one moves to larger λ, a bifurcation appears in the vicinity of λ j0 , with two wells labelled I and II. At a larger value of λ, there are the bottoms, where the wells have their minimum magnetic field (different in I than in II) and beyond which no orbits are allowed. At a generic point, we make use of equation (3.5) and then employ central finite differences with second-order accuracy with ∆λ = λ j+1 − λ j . Derivation is done at fixed α and well-label w. Nevertheless, contribution to the flux is concentrated on very thin λ layers. For the wide parameter range that will be studied with KNOSOS, the uniform grid is considered appropriate. equation (3.6) relies on ∂ λ g being continuous, which is not fulfilled at bifurcations if the effect of the tangential drift is small. There, we use finite differences with second-order accuracy directly over equation (3.4) and summing over wells, This discretization is best suited to obtain the expected relation between different values of ∂ λ g at the bifurcation for the 1/ν regime (Nemov et al. 1999;Calvo et al. 2014). Finally, we have two kinds of boundary conditions: one at the boundary between passing and trapped particles, corresponding to equation (2.10), and one at the bottom, corresponding to regularity (Calvo et al. 2013), Here we have employed a ghost point λ N λ +1 at exactly the bottom of the well, where ∂ λ g b | i,N λ +1,w = 0. One precision must be made: while in omnigenous magnetic fields the values of the maxima and minima of B are the same when moving in α, and equation (3.9) can be used as such for all α, this ceases to be true in a generic stellarator. For instance, the distance from λ N to the local bottom will be exactly ∆λ for one field line and smaller elsewhere (it may even happen that the contour condition must not be applied . This requires introducing straightforward corrections to equations (3.6), (3.7), (3.8) and (3.9). Let us now turn our attention to the terms with the first derivative in α in equation (3.1), which we multiply by (3.10) We represent the α grid at fixed λ in figure 4 (top). In this example, there is only one well at α 1 , labelled I. If one moves from smaller to larger α, a bifurcation appears in the vicinity of α i0 , with two wells labelled I and II. At a larger value of α, the wells merge into a single region labelled again I. The last point of the grid, α Nα , is close to α 1 + 2π. Non-centered finite differences with second-order accuracy are used. For a given fluxsurface, for each solution of the drift-kinetic equation, the sign of the coefficient in front of ∂ α g b (i.e. the direction of the flow in the α direction) indicates whether forward or backward differences should be used, with ∆α = α i+1 − α i . We thus find it useful (in terms of computing time devoted to filling the matrices of the linear problem) to discretize separately the terms I v E ,α ∂ α g and I v M ,α ∂ α g. In the first case, the direction is given by the sign of Ψ t ∂ ψ ϕ 0 , i.e., by the sign of the radial electric field; in the second case, by the particle charge and, in general, by the value of λ and w. With this in mind, four different matrices are created at the beginning of the simulation for a given flux-surface: one corresponding to forward differences being used everywhere; one corresponding to backward differences everywhere; and two λ (and w)-dependent discretizations, in which forward (backward) differences are used according to the sign of I v M,α . Any matrix appropriate for describing equation (3.10) will be a linear combination of two of these four pre-calculated matrices, and a simulation including ions and electrons and/or different values of the radial electric field will generally make use of the four of them.
Periodicity in α is easily imposed by replacing equation (3.11) at i N α − 1 with (3.14) respectively, and equation (3.12) at i 2 with respectively. We note that, since ι is generally irrational, 2π + α 1 − α Nα will e.g. be slightly smaller than ∆α We also note that bifurcations do not pose a problem for α-derivatives, due to g b being continuous in α. For example, in the vicinity of α i0 , see figure 4 (top), the forward derivative is discretized where continuity of g b has been used in the first two expressions of equation (3.17). Equivalent expressions can be obtained for the backward derivative.
One final caveat has to be made. In an omnigenous magnetic field, the contours of minimum B on a flux-surface must encircle the plasma (toroidally, poloidally, or helically). This is not true for a generic stellarator, in which local minima of B exist on the flux-surface. Close to these minima, moving in α at constant large λ is not always possible, as these trajectories may not exist. This situation is illustrated in figure 4 (bottom), at α i0 . There, instead of equation (3.12), we use (3.18) and we have implemented two models: in one, λ j0 is the value of λ closest to λ j in which trajectories exist for all α; in the second model, λ j0 is the closest value of λ in which trajectories exist at α i0−2 . The relative differences between the two models are smaller than e.g. the error bars of DKESin figure 7. We note that (with different manifestations for other choices of velocity coordinates) an incorrect treatment of this kind of particles is common to all existing radially local codes. For each of the species b, we end up with an equation that is linear in g b and can be written as a linear problem in matrix form. The matrix that represents is square with approximately N λ × N α elements per row, and sparse, with ∼ 6 non-zero elements per row: between 3 and 5 for the α derivatives, and typically 2 additional points for the collision operator. Although their relative weight varies with ν λ,b , v d,b and ∂ ψ ϕ 0 , the non-zero elements are always at the same position for a given flux-surface, which can be used to save computing time, by using the four pre-computed matrices described above.
We solve the linear problem with a direct solver from the PETSc library (Balay et al. 1997(Balay et al. , 2019a based on LU factorization. The reason is that the matrix is not large enough to require iterative methods, and reusing the LU factorization greatly accelerates the solution of the quasineutrality equation, as discussed in §3.5.

Solution of the quasineutrality equation
We will solve the quasineutrality equation by means of a response matrix approach (similar methods are used in gyrokinetics for the calculation of the electrostatic potential fluctuations (Kotschenreuther et al. 1995)). Let us first rewrite equations (3.1) and (2.18) making explicit the dependence on ϕ 1 : It can be observed that equation (3.20) is linear in ϕ 1 , and therefore the response of the distribution function g b (and of its velocity integral) of species b to certain ϕ 1 can be calculated as a superposition of the responses to a complete set of harmonics that parametrize ϕ 1 (θ, ζ). We can perform this parametrization efficiently thanks to the Fast Fourier Transform, using N = 2(2N n + 1)(N m + 1) coefficients: (the grid defined in §3.3 is not uniform in θ, so an interpolation is done before the Fourier transform). We can now denote u k (θ, ζ) each of the N basis elements (e.g. cos(θ + 2N ζ)) and the combined system of drift-kinetic and quasineutrality equation can be symbolically written as where ϕ 1 ϕ 1 ϕ 1 is a vector whose N components are the coefficients of the expansion of ϕ 1 according to equation (3.22) and A is a generally dense N × N matrix. In this linear, system, the right-hand side ϕ 0 and A have been filled, the new linear system can easily be solved, e.g. using a new LU decomposition, to obtain ϕ 1 ϕ 1 ϕ 1 , i.e., the set of coefficients ϕ (c) mn and ϕ (s) mn that parametrize the solution to quasineutrality. Finally, since the response of g b to every basis element has already been computed, a simple linear combination yields the distribution function that is solution of the drift-kinetic and quasineutrality equations, without requiring an additional solve of the former.
In summary, the drift-kinetic equation is solved a total of N +1 times (for each species), but LU factorization is done once (for each value of v). The linearity of the system of equations due to the smallness of ϕ 1 , together with the method that we have chosen for solving the drift-kinetic equation, yields a large reduction of the computing time needed for solving the system of equations: the code is roughly N faster (with N ranging from 100 to 1000), with respect to an equivalent code that allowed the particle orbits be modified by ϕ 1 .

Results
In this section, we show calculations for a variety of three-dimensional magnetic configurations in order to compare KNOSOS with widely-benchmarked codes and to illustrate its performance. In §4.1, we will solve a simplified drift-kinetic equation, without the magnetic drift and electric field components tangent to the flux-surface, and we will compare our results with bidimensional databases of DKES monoenergetic transport coefficients. The effect of the tangential magnetic drift in the energy flux, calculated for realistic kinetic profiles, will be discussed in §4.2. Finally, solutions of the quasineutrality equation will be compared with EUTERPE calculations in §4.3.

DKES-like monoenergetic transport coefficients
In this subsection, we will show that KNOSOS can be used for creating a DKES-like database of monoenergetic transport coefficients at low collisionalities. We will compare our calculations with DKES, both in results and computing time. Let us first discuss the rationale behind the monoenergetic approach, which is not specific of DKES, and the particular simplifications involved in DKES. More details can be found in the overview paper .
Predictive transport simulations solve the energy transport equation for every species: 3 2 where P b is the net energy source to species b and the energy flux Q b contains a turbulent contribution, at least close to the edge, that is currently provided by simplified models (Turkin et al. 2011). Calculating the time evolution of the energy, as in (Pedersen et al. 2015), or finding the steady-sate solution as in , requires evaluating the neoclassical contribution to Q b a large number of times. The monoenergetic approach, together with some simplifications to the drift-kinetic equation, provides a way out of solving the drift-kinetic equation many times. Strictly speaking, monoenergetic transport coefficients can always be calculated if the velocity v is a parameter in the drift-kinetic equation that is being solved, as in the case of equation (3.1): one can rewrite as a convolution of monoenergetic transport coefficients where g b is the solution of equation (3.1). Up to this point, the reduction in computation time associated to the monoenergetic approach has not derived from this rewriting, but from the fact that v is a parameter in equation (3.1), which is then easier to solve than a drift-kinetic equation with energy diffusion in the collision operator. Additionally, some fundamental simplifications are done by DKES: (4.5) Here,ĝ b is the solution of a modified version of equation (3.1), simplified aŝ In other words, the effect of the tangential electric field and the tangential magnetic drift is ignored, and an incompressible E × B tangential drift is used (this last simplification is specific of DKES and is not used by other codes in ). While it is well known (Calvo et al. 2017) that these effects need to be kept in the drift-kinetic equation for an accurate computation of the radial fluxes, there is a range of situations in whicĥ Q b ≈ Q b (this will be discussed in detail in §4.2) and this inaccuracy allows for a very large reduction of the computing time. The reason is that, for a given flux-surface, when normalized by the plateau valueD * the transport coefficientsD * 11 only depend on two v-dependent dimensionless parameters, the collisionality (4.10) and the normalized radial electric field v E * = E r vB 0,0 . (4.11) Here, R 0 is the major radius, and the main Fourier mode of B (see appendix C) is B 0,0 ∼ 1 T in all the simulations presented in this paper. Since there is no species dependence, in the rest of the subsection we follow the common practice of dropping the species index when discussing monoenergetic calculations. A predictive transport simulation thus requires to precompute a so-called database of (DKES-like) monoenergetic coefficientsD * 11 (ν * , v E * ). Once this is done, the calculation ofQ b for given n b , T b and E r using equation (4.4) requires a few bidimensional interpolations and an integral in v. The problem then lies in the computation of the databaseD * 11 (ν * , v E * ) for every new magnetic configuration, which typically takes hours, due to the poor convergence of DKES (and most neoclassical codes ) at low collisionalities. We will show that the bounce-average technique greatly reduces the computing time by using in KNOSOS equation (4.6) and comparing the results with DKES. Calculations without the simplifications made by DKES are left for §4.2.
In order to illustrate the performance of KNOSOS in a variety of three-dimensional configurations, we choose four very different types of stellarators. Figure 5 shows the map of the magnetic field strength on the flux-surface ψ/ψ LCF S = 0.3 of the highmirror configuration of the helias W7-X (top left), the R ax = 3.75 m configuration of the heliotron LHD (top right), an equilibrium of NCSX close to quasiaxisymmetry (bottom left) and the standard configuration of the heliac TJ-II (bottom right) . Figure 6 shows the first comparisons between KNOSOS and DKES, in which the normalized monoenergetic transport coefficientsD * 11 are calculated for several values of the collisionality and the normalized radial electric field. Figure 6 (top left) contains data for the W7-X high-mirror configuration, which we discuss in more detail. The expected 1/ν dependence is observed at the highest collisionalities and, due only to the absence of tangential magnetic drift, for small values of v E * . There is √ ν characteristic behaviour elsewhere, with smaller levels of transport for larger |E r |. The comparison between KNOSOS and DKES is satisfactory, with agreement within the error bars of the DKES calculation, and only at the highest collisionalities, and for the largest values of E r , there are very small differences. The calculation for all the points of this case was made with N α = 32 and N λ = 64, and it took 2.0 seconds in a single standard CPU. Of this time, around 0.7 seconds were used for setting the grid and performing the bounceaverages, and then it took less than 0.04 seconds to calculate each point. This number may be reduced even further using smaller N λ for the cases of largest collisionality and smallest radial electric field: in the √ ν regime, transport is given by a small layer close to the boundary between passing and trapped particles, the size in λ of this layer is proportional to ν λ /E r (Calvo et al. 2017), and this determines the required number of grid points N λ in each case. Similar results can be seen for LHD in figure 6 (top right). N α = 32 and N λ = 64 grid points were used, and the total computation time was 2.1 seconds. For NCSX, figure 6 (bottom left) the agreement is good except for the higher collisionalities, where the 1/ν regime should connect with a banana regime (see figure 15 of ). This regime could be easily implemented in KNOSOS following (Landreman & Catto 2012). N α = 32 and N λ = 64 grid points were used, and the total computation time was 1.0 seconds. Finally, figure 6 (bottom right) contains the results for TJ-II, the hardest case due to its complicated magnetic geometry, see figure 5 (bottom right). N α = 32 and N λ = 128 were used, and the simulation took 157 seconds. Here, one can see more  clearly that the cases deeper in the √ ν regime are more difficult to converge, and indeed the points corresponding to v E * B 0,0 10 −3 T and ν * < 10 −4 are barely converged. The rest of the simulations agree with DKES and reach even lower collisionalities than those typically required for describing a TJ-II plasma, whose ion temperature never exceeds a few hundred eV. Figure 7 contains, for each of the four configurations, radial profiles of the transport coefficientD * 11 for two cases, v E * = 0 and v E * B 0,0 = 3 × 10 −4 T, for a given collisionality. They are meant to represent the level of transport in the 1/ν regime (there,D * 11 is by definition proportional to 3/2 ef f , being ef f the effective ripple) and the √ ν regime, respectively. We choose ν * = 2×10 −5 for W7-X (top left) and LHD (top right), ν * = 10 −5 for NCSX (bottom left) and ν * = 3 × 10 −5 for TJ-II (bottom right). It can be observed that the good agreement holds for all cases at all radial positions. The comparison of the different parts of figure 7 provides additional information that may be relevant when devising a stellarator optimization strategy: in general, configurations with lower 1/ν transport show lower √ ν transport as well. This is not surprising considering that both quantities are connected to the bounce-averaged radial component of the magnetic drift, which appears in the source of the drift-kinetic equation (2.9) in both regimes, and which is in turn proportional to the variation of the second adiabatic invariant on the flux-surface, ∂ α J. Inasmuch as the optimization procedure actually reduces the size of ∂ α J, both the 1/ν and √ ν (and superbanana-plateau) regimes will generally be optimized. Nevertheless, using directly the effective ripple as figure of merit of neoclassical transport does not automatically guarantee a reduction of ∂ α J, and the √ ν transport may remain unoptimized. Figure 7 (top) may represent an example of this situation: while this W7-X configuration is designed to have low level of 1/ν transport at an intermediate radial position (where the plasma volume is relatively large and neoclassical transport is expected to be at least comparable to anomalous transport), the √ ν transport is smallest exactly at the magnetic axis. We will argue the relevance of optimizing transport regimes of collisionality lower than the 1/ν, which has now become possible with KNOSOS, in the next section.

4.2.
Effect of the tangential magnetic drift on the radial transport of energy In §4.1, we have shown solutions of equation (4.6), a simplified drift-kinetic equation that is not accurate when the tangential components of the magnetic drift and of the electric field play a role, one needs to solve equation (3.1). In this section, we will demonstrate the importance of solving equation (3.1) instead of equation (4.6), i.e., of computing Q b and notQ b , when calculating the radial energy flux in real plasmas. It must be noted that the solution of equation (3.1) with KNOSOS is not computationally more expensive than that of equation (4.6): in the superbanana-plateau regime, that may arise in the presence of the tangential magnetic drift for certain values of E r , transport is dominated by a resonant layer whose size decreases with (ν λ /E r ) 1/3 , i.e., slower than the boundary layer that determines the √ ν transport (Calvo et al. 2017). Calculating Q b instead ofQ b does not require a larger value of N λ in general.
In this section, we focus on characterizing the effect of the tangential magnetic drift for the particular case of ϕ 1 = 0. We advance one of the salient results: this effect will be non-negligible even at not very low collisionalities. The reason is that the calculation of the energy flux for a given plasma, characterized by the kinetic profiles, requires the solution of the drift-kinetic equation for several values of the velocity, see equation (3.3), with the normalized particle energy (v/v th,b ) 2 spanning several orders of magnitude. This means that, even if the thermal particles are in 1/ν regime, there are particles with higher v that are in lower collisionality regimes. Figure 8 contains simulations for the high-mirror configuration of W7-X at ψ/ψ LCF S = 0.25, which corresponds to r/a = 0.5. We choose a pure hydrogen plasma, with n e = 8.0×10 19 m −3 , ∂ r n e /n e = −2.0 m −1 , T e = T i = 4.0 keV, ∂ r T e /T e = ∂ r T i /T i = −3.0 m −1 . These are values comparable to those measured in high-performance OP1.2 plasmas of W7-X (Klinger & the W7-X Team 2019) in the region of crossover between positive and negative radial electric field, corresponding to electron and ion root solutions of the ambipolarity equation (Pablant & et al 2019). In these plasmas, neoclassical transport calculated neglecting the tangential magnetic drift typically accounts for around half the total experimental transport. Figure 8 (top) contains an overview, in logarithmic scale, of the ion and electron radial energy flux as a function of the radial electric field. Empty and full blue boxes correspond toQ i and Q i respectively, both calculated with KNOSOS. We immediately see thatQ i overestimates the radial energy flux at small values of the radial electric field, specially at E r = 0 (strictly the only point of the figure whereQ i is proportional to ε 3/2 ef f ). The tangential drifts make the ion flux decrease, differently in the case ofQ i and Q i , as we will discuss below. Finally, empty and full red boxes correspond toQ e and Q e calculated with KNOSOS. In this plot, is difficult to notice any difference between the different electron calculations. Figure 8 (top) contains additional black lines that are the result of combining calculations with DKES and KNOSOS. We will leave the discussion of these results at the end of the section. Figure 8 (bottom) contains a blowup of the most relevant range of the previous data in linear scale. Here, the effect of the tangential magnetic drift on the energy flux can be observed more clearly: the size of the peak at small |E r | is reduced and displaced to positive (negative) values in the case of electrons (ions). The effect is larger for the ions due to their larger normalized Larmor radius ρ i * , which makes them leave the 1/ν regimes at relatively higher collisionalities. We have mentioned that these plasmas are close to the crossover between ion and electron root, and this figure can help us discuss some features of transport in both situations. In electron root, the radial electric field is expected to be positive and large, and the electrons are expected to give the largest contribution to energy transport. According to figure 8 (bottom), Q e provides a minor, although systematic, correction toQ e , below 10% for this plasma profiles and configuration. The situation is different in ion root, typically characterized by a negative radial electric field that is small in size, and dominant ion transport. Here, including the tangential magnetic drift can lead to large corrections, above 50% in some cases.
For the sake of completeness, figure 9 contains two more cases. In figure 9 (top) we repeat the calculation for a W7-X plasma of much higher collisionality, choosing n e = 1.6×10 20 m −3 , ∂ r n e /n e = −2.0 m −1 , T e = T i = 2.5 keV, ∂ r T e /T e = ∂ r T i /T i = −3.0 m −1 . We first note that the electrons are deep in the 1/ν regime, since ν e * = 3.4 × 10 −2 and ρ e * = 1.4 × 10 −5 . Nevertheless, Q e (E r ) does not show the linear dependence expected when the 1/ν dominates. This is an indication of what we advanced at the beginning of this section: even in plasmas nominally in the 1/ν regime, the contribution of the √ ν regime is not negligible, and should not be neglected in the optimization procedure. For the ions, even at these higher collisionalities and low temperatures, ν i * = 1.6 × 10 −2 is not much larger than ρ i * = 6.0 × 10 −4 divided by the inverse aspect ratio. This means that, for ions slightly more energetic than the thermal ions, the tangential magnetic drift is relevant at small values of |E r | (Calvo et al. 2018b). Figure 9 (top) shows indeed systematic differences betweenQ i and Q i . Finally, figure 9 (bottom) contains a calculation with the same kinetic profiles of figure 8 (top) for the inward-shifted configuration of LHD. It can be observed that the effects discussed in figure 9 (top) are even more pronounced, to the extent of changing qualitatively the Q b (E r ) dependence (and making it more similar to that reported in (Matsuoka et al. 2015): while practically any increase of |E r | causes a reduction of Q i in W7-X, this is not the case for LHD. For finite ion-root values of E r , Q i (E r ) has a peak whose height is determined by superbanana-plateau transport.
In light of these results, two comments related to stellarator optimization can be made. First, the fact that the monoenergetic transport coefficients respond to small tangential E × B drifts differently in the inward-shifted LHD, with respect to other configurations, was already discussed in , and it can be observed more clearly when calculating the energy flux including the tangential magnetic drift. We also note that part of the neoclassical optimization of W7-X comes from its large aspect-ratio, which tends to make the tangential magnetic drift smaller, when compared with the E × B drift. It is then clear than a systematic study of the different lowcollisionality regimes, and their different configuration dependence, should be addressed when devising an stellarator optimization strategy. Second, a comprehensive optimization strategy will involve, at least, solving energy transport consistently with ambipolarity and quasineutrality. Along this subsection, we have compared Q andQ at fixed E r , but a more systematic study applied to real discharges of W7-X, including the experimental validation of E r predictions, is ongoing .
Let us finally discuss the black lines of figures 8 and 9, which correspond to combining simulations of DKES and KNOSOS. As we have argued at the beginning of this section, calculating the radial energy flux requires solving the drift-kinetic equation for velocities (v/v th,b ) 2 spanning from ∼ 10 −2 to ∼ 10 2 , typically. Similarly to what we discussed for v v th,b , this means that particles with v v th,b could be in the plateau regime, and they would not be described by equation (3.1). In order to quantify this effect, and to show that it is negligible for the high-performance plasmas of W7-X, we perform calculations of Q i (E r ) and Q e (E r ) combining KNOSOS with DKES. This can be done by rewriting equation (2.22) as (4.12) where H is the Heaviside function, v 0 is a cut-off velocity,D * 11 (v) comes from DKES in this case and (4.13) from KNOSOS. The latter is calculating according to equation (4.3) solving the driftkinetic equation that is correct at low collisionalities with ϕ 1 set to zero. In other words, monoenergetic transport coefficientsD * 11 coming from DKES are used above certain collisionality when performing the velocity integral and monoenergetic transport coefficients D * 11 coming from KNOSOS are used below that collisionality. The cut-off velocity v 0 must correspond to particles in the 1/ν regime, which is correctly described by the two codes. This guarantees that both codes are employed in the parameter region where they are accurate (and fast).
In figures 8 and 9 (bottom), the black lines corresponding to using equation (4.12) barely separate from the solution of equation (3.1). This means that the contribution of the plateau regime to the energy flux is negligible. Only for ions in the presence of very negative values of the radial electric field, in the high-density W7-X calculation, starts the black line to separate from the blue signs. This is to be expected: due to the approximate (exact in the case ofQ i ) |E r | −3/2 -dependence of the transport coefficients in the √ ν, the contribution of low collisionalities to transport is reduced for very large values of |E r |, and therefore the contribution of the plateau becomes non-negligible.

Tangential electric field
Neoclassical physics gives rise to ϕ 1 , and the associated tangential electric field produces radial drifts in all species. This is the reason why we need to solve consistently the drift-kinetic equations of the bulk species and quasineutrality (Calvo et al. 2017), but the effect is more relevant for impurities, due to their larger charge number, changing even qualitatively transport (e.g. making it depend on the radial electric field in the so-called mixed collisionality regime (Calvo et al. 2018a;Buller et al. 2018)). With impurity transport in mind, simulations of ϕ 1 for the stellarators W7-X, LHD and TJ-II have been performed in the last years with three codes, EUTERPE, SFINCS and recently FORTEC-3D (García-Regaña et al. 2013;García-Regaña et al. 2017, 2018Mollén et al. 2018;Fujita et al. 2019). Nevertheless, the number of simulations remains small because they are computationally very demanding, specially at low collisionalities. A more comprehensive study, including dependence on the configuration, collisionality, and bulk plasma profiles thus remains to be done. In this section, we will show that KNOSOS can reproduce the results of EUTERPE (with adiabatic electrons and no tangential magnetic drift) and, by accounting for the effect of the tangential magnetic drift, describe stellarator regimes only simulated before for simplified geometries (Calvo et al. 2018b;Velasco et al. 2018). Since it can do so while keeping the computing time low, this opens the door to a number of new impurity transport studies.
We start by reproducing the results of (García-Regaña et al. 2017), specificaly of two low-collisionality plasmas of LHD and W7-X. These are expected to be the plasma conditions of largest eϕ 1 /T i so that, even in optimized magnetic configurations, the effect   on the radial transport of impurities may be large. It will be confirmed (as advanced in a previous work (Velasco et al. 2018) in a simplified calculation) that the inclusion of the tangential magnetic field leads to qualitative changes in ϕ 1 , making it larger. Figure 10 shows the variation of the electrostatic potential on several flux-surfaces of the inward-shifted configuration of LHD for a low-collisionality plasma (described in (García-Regaña et al. 2017), termed AIII, and characterized by a small negative E r ). Each row corresponds to a different flux-surface, and each column to a different calculation method. Let us start by comparing the left column, calculated with EUTERPE, with the center column, calculated with KNOSOS using equation (4.6) for describing ion transport. The two methods should give the same results, and it can be observed that, although there are some differences, reasonable agreement between the two codes is obtained. It should be emphasized that differences in calculated values of ϕ 1 similar but smaller to those reported here, have been shown to produce negligible differences in impurity transport (Velasco et al. 2018). If we now focus on the right column, we observe, as discussed in detail in (Velasco et al. 2018), that the inclusion of the tangential produces relevant differences (in particular, more important than those betweens columns left and center): the amplitude becomes larger, and the phase changes, with the angular dependence of ϕ 1 turning from being stellarator-symmetric (as expected for ions in the  √ ν regime), to not having definite symmetry (as corresponds to the superbanana-plateau regime (Calvo et al. 2018b)). Figure 11 contains a similar calculation performed for a low-collisionality plasma of W7-X (described in (García-Regaña et al. 2017), termed IV, and characterized by a larger negative E r ). Again, each row corresponds to a different flux-surface, and each column to a different calculation method. The comparison between EUTERPE and KNOSOS solving the same equation (left and center) is good close to the core, but it becomes slightly worse closer to the edge, where KNOSOS underestimates the amplitude of ϕ 1 . When the tangential magnetic drift is included (right), the results change very slightly in the core and do not change elsewhere. This feature is likely caused by the large radial electric field, which leaves the ions in the √ ν regime (instead of the superbanana-plateau).
The computing time for each of these KNOSOS simulations is of the order of a minute in a single processor. We note that including kinetic electrons (which may be necessary for high electron temperature) would roughly double the computing time. This is to be compared with the (m i /m i ) 1/2 ≈ 43 factor in Monte Carlo codes such as EUTERPE and FORTEC-3D.
Let us finally mention that the experimental validation of ϕ 1 predictions has drawn much attention in the last years: experimental measurements of ϕ 1 were first obtained at the edge of the TJ-II stellarator Pedrosa et al. (2015), and very recently in its core region . The validation of KNOSOS predictions, including finer scans in the magnetic configuration, is left for a forthcoming work.

Conclusions
KNOSOS is a freely-available open-source code that provides a fast computation of neoclassical transport at low collisionality in three-dimensional magnetic confinement devices, thanks to a rigorous application of the orbit-averaging technique to the the driftkinetic equation and an efficient solution of the quasineutrality equation. We have shown that, when solving equivalent equations, KNOSOS reproduces the calculations of DKES and EUTERPE in simulations that can be orders of magnitude faster. This makes it a tool that can be used for a variety of physics problems, that we summarize next. As a first obvious application, it can provide a fast calculation of the level of transport of a magnetic configuration for low-collisionality transport regimes not usually considered in stellarator optimization, such as the √ ν and superbanana-plateau regimes. Optimization programmes, however, are slowly starting to provide more a accurate characterization of transport by performing predictive simulations with prescribed sources and turbulent transport models. KNOSOS can also contribute to overcome two of the main limitations of this approach: the large computing time needed to create a database of mononoen-ergetic neoclassical transport coefficients and/or the lack of accuracy involved in the monoenergetic approach itself. But a fast neoclassical code can have uses beyond stellarator optimization. For instance, the transport of impurities caused by their interaction with the bulk ions (via ϕ 1 or through inter-species collisions) has drawn much attention in the last years; however, a systematic study of its dependence on the magnetic configuration, colllisionality, and bulk plasma profiles remains yet to be done, due to the large computing resources needed for the combined solution of the quasineutrality and drift-kinetic equations of the bulk species. This will be addressed in forthcoming papers, in combination with analytical formulas for the radial flux of impurities in a variety of neoclassical regimes (Calvo et al. 2018a. Finally even in situations in which turbulence is dominant, a fast neoclassical code may be required. Its output (the radial electric field, the tangential electric field or the complete distribution function of the bulk species) can be read by gyrokinetic codes when studying the effect of neoclassical transport on turbulence. This effect is expected to be largest in those low-collisionality regimes in which the specificities of KNOSOS, very small computing time and inclusion of the tangential magnetic drift, are most relevant.

Acknowledgments
This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. This research was supported in part by grant ENE2015-70142-P, Ministerio de Economía y Competitividad, Spain and by grant PGC2018-095307-B-I00, Ministerio de Ciencia, Innovación y Universidades, Spain. and erf(x) = (2/ √ π) x 0 exp(−t 2 ) dt (A 4) is the error function. and so on.