Magnetic reconnection in 3D fusion devices: non-linear reduced equations and linear current-driven instabilities

Magnetic reconnection in 3D fusion devices is investigated. With the use of Boozer co-ordinates, we reduce the non-linear resistive magnetohydrodynamic equations in the limit of large aspect ratio and finite pressure fluctuations, to obtain a set of non-linear equations suitable for magnetic reconnection studies in stellarators. Magnetic flux unfreezing due to a finite electron mass is also considered. Equations that govern the linear regime and some of their general properties are given. We emphasise the role of magnetic geometry and identify how some aspects of stellarator optimisation could have an impact on reconnecting instabilities, in particular by exacerbating those enabled by electron inertia. The effect of 3D coupling on the linear reconnection rates and the mode structure is quantitatively addressed in the case in which the equilibrium rotational transform has one specific resonant location for which one mode can reconnect while coupled to an arbitrary number of non-resonant harmonics. The full problem is rigorously reduced to an equivalent cylindrical one, by introducing some geometrically modified plasma inertial and dissipative scales. The 3D scalings for the growth rates of reconnection instabilities and their destabilisation criteria are given.


Introduction
The high level performance in the recent operational phase of Wendelstein 7-X (W7-X) [1][2][3][4] has been a key step in the development of the international stellarator research programme [5][6][7][8][9][10]. While, on the one hand, W7-X has set new standards in stellarators and magnetic confinement fusion devices, on the other it has provided a wealth of stimulating physics results. One of the most intriguing is the presence of current driven instabilities, externally driven by electron cyclotron resonance heating, which results in sawtooth-like plasma oscillations, and even electron temperature collapse [11]. Since the physics of sawtooth oscillations in tokamaks 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. [12] is understood on the basis of magnetic reconnection theory [13][14][15], the explanation of such phenomena in W7-X requires a new look at the theory of magnetic reconnection in stellarators. Some fundamental aspects are, of course, expected to remain unaffected by the absence of axisymmetry. One of these is the boundary layer character of reconnecting instabilities, which occur as a local response acting to regularise current singularities otherwise driven by some special types of perturbed ideal magnetohydrodynamics (MHD) plasma displacements. Some others are, as we shall see, different. In the literature, we can already find studies of magnetic reconnection in stellarators [16,17]. In most of them, one sees that the rotational transform (poloidal/toroidal winding number) plays an important role. The reason lies in the fact that the singular currents regularised by the instability are radially located where the rotational transform is rational. Some authors (see Matsuoka et al [18], for instance) have produced stability diagrams for reconnecting instabilities, highlighting the role of the rotational transform of the helical field, ι δ ,ι δ , and that generated by the plasma current ι σ . It is established that the value of the derivative of the rotational transform at rational surfaces is also a critical parameter for the destabilisation of sawteeth in tokamaks [13,15], as well as for neoclassical tearing modes in quasihelically symmetric stellarators [19].
The effect of resistivity in hydromagnetic linear instabilities in general geometry was first studied by Johnson et al [20]. Here the authors used the so-called stellarator expansion [21], i.e. their analysis is limited to axial equilibria. A more general approach was presented by Glasser et al [22], who used Hamada [23] co-ordinates but, for the stellarator case, only considered high mode-numbers, thus not fully addressing general reconnecting instabilities in stellarators (which are mostly low mode-number instabilities). In these works, an ordering, specific to the resistive interchange instability, is implemented for the linearised equations. In our work, we propose a more general ordering in large aspect ratio and small β = µ 0 p 0 /B 2 , (the ratio of kinetic to magnetic pressure) which is independent of the linear instability we consider and generates a non-linear reduced model for magnetic reconnection in stellarators.
The pervasiveness of magnetic reconnection in many phenomena occurring in fusion plasmas makes its understanding in 3D devices indispensable, not only for stellarators, but also for externally perturbed tokamaks and conceptual devices that combine features of stellarators and tokamaks.
A pictorial understanding of 3D magnetic reconnection is a notoriously hard exercise, as fascinating as arid of quantitative results. For this work, we find it useful to exploit the concepts of 'small' and 'large' solutions introduced by Newcomb [24]. That is, one first solves the equation for the reconnecting magnetic flux, χ, in two distinct, asymptotically separated regions: one at the resonant surface (ψ ∼ ψ res , inner region), where radial derivatives are relatively large, and one far from it (outer region), where reconnection is negligible. Here ψ is a radial co-ordinate. Then, the study of the asymptotic behaviour of inner and outer solutions, in the overlapping domain, generates an indicial equation that gives the leading order terms of a power series in x = ψ − ψ res . For small enough pressure gradients, from both regions, two solutions with different asymptotic behaviours are generated: , the 'small' solution. In the external region, the complete solution is a linear combination of the small and large ones, and for x → 0 is χ ext = A L + B s x. The ratio 2B S /A L is related to the so-called tearing mode instability parameter, ∆ ′ . In the simplest case, 2B S /A L = ∆ ′ . For ∆ ′ of order unity, the relative perturbation is nearly constant, while for large ∆ ′ it is not. The first type of solutions are called 'tearing modes', the second are dissipative kink modes. We generally refer to 'reconnecting modes' encompassing both types. For the external region solutions, first in axisymmetry, Glasser Green and Johnson discovered a stabilising effect due to toroidal geometry [25], which sets a critical ∆ ′ for tearing mode destabilisation. Other authors have stressed the relation between ∆ ′ coefficients at different resonant locations in axisymmetric devices, for low [26,27] and moderate high [28] mode numbers. One of the main consequences of the work of Connor et al [26] is that the usual distinction between constant-χ and non-constant-χ fails, since the final eigenfunction is the result of a mixing of these two 'flavours'. In general geometry, resonant modes contribute to magnetic reconnection at their respective rational surfaces, but they do interact. In this article we address some aspects of this problem. We show that mode coupling is important not only in the external regions, but also right at the locations where magnetic perturbations can be resonant and reconnection occurs. In general geometry, for low β, the scaling of reconnection rates with key geometric factors is derived. We find that collisionless instabilities scale more unfavourably with the metric element that regulates how close magnetic surfaces are. In a subsidiary low-β limit, for a rotational transform that allows for one resonant reconnecting mode and an arbitrary number of non-resonant ones, it is shown that non-resonant modes have an order O(1) effect on the growth rate of the resonant one by re-defining inertial and dissipative scales. The full external region solution is also given in a special case, and the singular behaviour of the non-resonant mode, induced by 3D coupling to the resonant one, is quantified.
This article is organised as follows. In section 2, the nonlinear reduced magnetohydrodynamics equations for large aspect ratio and finite-β 3D devices are derived. In section 3, the linear equations governing reconnecting modes in 3D are presented, and some of their general properties are discussed. In a subsidiary low-β expansion, explicit growth rates for magnetic reconnecting instabilities in 3D are derived. Conclusions are in section 4.

Reduced magnetohydrodynamics equations for stellarators
Magnetic reconnection in plasmas occurs only if some physical mechanism allows for magnetic flux unfreezing. In this work we consider resistive MHD, and briefly discuss electron inertial effects. The full MHD equations, however, are not always ideal for numerical simulation, in particular when the aim is to simulate processes that are slower than magnetosonic waves propagating across the magnetic field. 'Reduced' MHD equations, where these waves are eliminated by a suitable ordering, are then more appropriate. Since these equations were first derived by Strauss [29], many derivations have been published, but their application to stellarator physics is very limited [30]. Most derivations are specific to tokamaks, and even those that are not, such as the one by Hazeltine and Meiss [31], tend to make assumptions that are not appropriate for stellarators. The aim of this section is to present a derivation that overcomes this problem.

Background and orderings
Reduced MHD relies fundamentally on a disparity between length scales parallel and perpendicular to the magnetic field B, Most quantities of interest are thus taken to vary on the length scale L ∥ in the direction of B and a factor ε −1 faster in the other two directions. An exception pertains to the equilibrium field B 0 , which varies on the length scale L ∥ in all directions. The magnetic field is written as a sum of this 'guide field' and a 'perturbation', and so are the gradients in the direction of B 0 and δB, thus allowing for a corresponding non-linearity in the equations.
For the case of a tokamak with inverse aspect ratio ε, Hazeltine and Meiss [31] choose B 0 to denote the equilibrium toroidal field and δB thus to represent both the equilibrium poloidal field and the perturbation due to instabilities. This choice is necessary for the study of magnetic reconnection because otherwise the gradient of the parallel equilibrium current would be smaller than that of the perturbed current, ∇J 0 ∇δJ, and the primary linear drive for the instability, dJ 0 /dr, would not appear in the appropriate order.
In a stellarator, it is not immediately obvious how best to choose B 0 . Indeed, several choices are possible, and since there is no unique way of identifying the toroidal component of the magnetic field, we shall identify B 0 with the full equilibrium field, which we assume satisfies a scalar-pressure equilibrium relation, J 0 × B 0 = ∇p 0 , with nested magnetic surfaces. It is then possible to write B 0 in Boozer coordinates [32,9], In the first of these expressions, the first term is assumed to be dominant, if the inverse aspect ratio is ϵ 1. If we thus take the major radius to be L ∥ and minor radius to be L ⊥ , we have In order to eliminate the fast magnetosonic wave but retain the shear Alfvén wave, we only consider perturbations that vary slowly with respect to the former wave and scale in the same way as the magnetic-field perturbation |δB|/B ∼ ϵ, i.e.
where v A = B/(µ 0 ϱ) 1/2 denotes the Alfvén speed, ϱ the plasma density, and V the plasma flow velocity. It follows that the plasma displacement is then comparable to the minor radius, so that the plasma pressure perturbation is of order unity, Our equations will therefore be appropriate for describing nonlinear instabilities. Finally, we need to choose an ordering for the normalised plasma pressure, β = 2µ 0 p 0 /B 2 . The natural choice follows from the equation of motion, which for the equilibrium (dV/dt = 0) suggests that, since the pressure varies on the length scale of the minor radius but B 0 only on the length scale L ∥ , we should choose . This ordering has consequences for the last term in equation (1), because of the equilibrium relation [9] In this equation, the only quantities that vary over a flux surface are K and B, so it follows that K must be of the order of where M = ∆B/B denotes the 'mirror ratio', i.e. the relative variation of B over the surface. Interestingly, it follows that the final term in equation (1) is in general not negligible, This ratio is not small unless M 1, since we have taken β = O(ε). We further note that the general expression for the parallel equilibrium current, obtained from the curl of equation (1), to lowest order in ϵ 1, because G ∼ ε −2 I, and

Magnetic-field perturbation
We now consider perturbations to this equilibrium. The equation of motion (2) simply becomes (to lowest order) ensuring that fast magnetosonic waves are absent. Furthermore, the relation 0 = ∇ · δB ∇ · δB ⊥ implies that we can write It follows that, for any function f (ψ, θ, φ) varying more quickly across the field than along it, In the last expression, the term proportional to I is negligible, being of order O(ε 2 ) compared with the first term. The Jacobian can similarly be approximated by 1/ with a relative error of order ε. The two first terms on the right represent the equilibrium current from equation (5), which in contrast to B 0 varies perpendicularly to the field on the length scale L ⊥ . The gradients of the equilibrium current and its perturbation are thus comparable.

The shear-Alfvén law and vorticity equation
We are now ready to consider the shear-Alfvén law [31] B · ∇ whereV = ∂V/∂t + V · ∇V, and the curvature vector, defined by κ = b · ∇b, is almost equal to its equilibrium value. The last term in the square bracket is relatively small since κ ∼ 1/L ∥ , and in the previous term we may make the approximation From Ohm's law, we conclude thatV × B d∇ϕ/dt, where ϕ denotes the electrostatic potential and the resistivity has assumed to be small. If we subtract the equilibrium version of the shear-Alfvén law, from equation (7), the remainder becomes where we have written B instead of B 0 since they are equal to the requisite accuracy. This is an appropriate form for the flute-reduced shear-Alfvén law in stellarator geometry. It is very similar to that of Hazeltine and Meiss [31] but has been derived under slightly different assumptions. In particular, the current associated with B 0 is allowed to vary rapidly with the minor radius and the pressure perturbation can be as large as the equilibrium pressure.

Ohm's law
It is not difficult to express Ohm's law, in this notation. The right-hand side expresses the effect of any non-inductive current drive, which is assumed to drive the equilibrium current. To proceed, we need only recall equation (6) and note that A −χB, so that

Equations for density and pressure
The density is governed by the continuity equation where the compressibility is relatively small, and can therefore be neglected. The continuity equation therefore reduces to The pressure obeys the entropy production law d dt where the right-hand side is usually neglected on the grounds that only phenomena occurring on time scales shorter than the energy confinement time are considered. The equation thus becomes ∂ ∂t where we have earlier written p = p 0 +δp. Equations (8)-(11) form a closed system of four non-linear equations for the four unknowns (χ, ϕ, ρ, δp). The density profile is however usually not of any great importance for the study of MHD instabilities. If it is assumed to be flat (∇ρ = 0), the pressure equation reduces to Note that all nonlinearities now appear in the form of Poisson brackets.

Linear theory of low-β reconnecting instabilities
We now consider β ϵ, neglect pressure perturbations, and linearise equations (8) and (9), to obtain and This is our starting point for the investigation of linear reconnecting instabilities in 3D devices.

General properties of the external region
In the external region, as already discussed in the Introduction, the plasma is considered to be marginally stable and resistive effects are negligible. Thus, after taking ∂ t ≡ γ → 0, and η → 0 in equation (12), we obtain and Equation (14) is the fundamental external region equation for arbitrary 3D geometry in the limit of negligible equilibrium pressure gradients. As opposed to the cylindrical case, considerable complications come from the differential operator on the LHS, since, in arbitrary geometries, the Laplacian is where each element of the metric tensor, g ij , as well as the Jacobian 1/ √ g are functions of the Boozer angles, θ and φ.
This introduces a coupling of the Fourier harmonics of the magnetic perturbation, χ ext . Such coupling has been studied in axisymmetric geometry by several authors. Obviously, in this case, only poloidal harmonics are coupled. However, in non-axisymmetric geometry, toroidal coupling is also expected. Indeed, it is not difficult to show that the Fourier expansion of equation (14) is where and subscripts refer to Fourier indices, with perturbed quantities written as δQ(ψ, θ, φ, t) = exp[−i(nφ − mϑ) + γt]δQ(ψ), and we suppressed the subscript 'ext.' By using these results, equation (15) is then written in the following useful form where each coefficient can be inferred from the expressions above. With equation (18), we are aiming at describing a mode χ mn that is resonant at the radial location ψ mn so that ι(ψ mn ) = n/m. Then, if we introduce ι(ψ) = ι(ψ mn ) + ι ′ (ψ mn ) (ψ − ψ mn ) + · · · . , at the RHS of equation (18) we obtain the familiar singular term also found in cylindrical geometry. It is interesting to note that the Fourier coefficients of the metric element g ψψ play a key role in the classification of solutions of equations (18), since they multiply the highest order derivative. In fact, when the second order derivatives dominate, for ψ → ψ mn , equations (18) are a coupled system of equations first studied by Giovanni Plana [33,34] (see appendix A). A particular property of this system of equations is that, if for one given resonant mode χ mn = 0, the second order derivatives of any other, non-resonant, mode can be infinite! We will study this phenomenon in detail for the simplest case of one resonant and one non-resonant mode, but first we introduce the inner layer equations for an arbitrary number of modes.

General properties of the inner region
As usual, in the reconnection layer, the electrostatic potential cannot be neglected. The analysis is further complicated by the inclusion of a finite growth rate and resistivity. The only simplification allowed pertains to high order radial derivatives. We shall assume where m ≥ 0, to obtain two second order differential equations that are still tractable, and very useful to predict how growth rates are expected to scale with geometric factors. Thus, by using equation (19), we obtain from equations (12) and (13) and where we took care to combine the fields in such a way that the Fourier series of δJ ∥ /B, in θ and φ, can be introduced unambiguously. One remarkable feature of these equations is the fact that, due to the θ and φ dependence of B, g ψψ and 1/ √ g, perturbed quantities are coupled. This issue did not arise in previous works on coupled reconnecting modes in tokamaks.
The system of equations (20)-(21) is of fourth order. This can easily be seen by using equation (20), writing ϕ as a function of χ and χ ′ ′ , and inserting the result in equation (21). We neglect the coupling for simplicity, and find 2 , and x = ψ − ψ mn . Thus, we have two exponential solutions that, far from the reconnection region, We also find an exact solution where C is a constant. This solution has the so-called 'twisting parity', that is ϕ ′ (0) = 0, and does not provide magnetic reconnection. The correction to this solution, as we approach the ideal region, is evaluated from the full equation (12), which we re-write, again neglecting geometric couplings and taking constant coefficients for illustrative purposes, in the following form withγ = Gµ 0 ϱg ψψ γ/B 4 , v 2 = m 2 g θθ /g ψψ , µ = ι ′ mg ψψ and ν = −imI ′ constants. We remind the reader that, far from the reconnection region, the second order derivatives are no longer dominant, thus the m 2 terms must be retained. We then find where C ′ and C ′ ′ are constants of integration. An xlogx term of the external region solution is therefore driven by the layer twisting parity solution. The role of the xlogx term in the matching of the tearing parity solution is discussed in appendix A. We are then left with a tearing parity solution, which we investigate in greater detail.

Dispersion relation
The inner layer equations (20)- (21), for all types of modes (resonant and non-resonant) are and sets of (2M+1)(2N+1) equations. Among all the couples of integers (m, n), there is a special one, and multiples thereof (which we do not consider here) for which ι(ψ m0n0 ) = n 0 /m 0 . Thus, for (m, n) = (m 0 , n 0 ), equation (27) gives In (29) we see that, if χ m0n0 ∼ χ mn ∼ ϕ mn , then and, for ι ′ x 1, the LHS of equation (28) becomes since all the non-resonant harmonics of the electrostatic potential can be neglected. We now notice that, for (m, n) = (m 0 , n 0 ), equation (28) gives while, for (m, n) = (m 0 , n 0 ), we have (33) Since in the sum on the LHS of equation (33) there is at least one non-null term, then and, for x → 0, that is, for (m, n) = (m 0 , n 0 ), and the second order derivative of a non-resonant harmonic is a linear combination of second order derivatives of all other harmonics, including the resonant one. Let us recall that −N ≤ n ≤ N, −M ≤ m ≤ M, so equation (35) is in fact a set of (2M + 1)(2N + 1) − 1 equations. Each of them involves χ ′ ′ m0n0 once. We write all these equations in matrix form, emphasising such χ ′ ′ m0n0 terms . . .
where one procedure to obtain the matrix G is given in appendix B. Hence, we are finally able to obtain, by inverting the matrix G, all the non-resonant terms in the sum of the RHS of the resonant Ohm's law (29), which now reads where G −1 · g ψψ k,1 are the entries at the kth row of the column vector G −1 · g ψψ .
Basically, the role of non-resonant harmonics is to modify resistive and inertial scales by a factor of order unity. Then, the eigenvalue equation of the full 3D problem is obtained by introducing the generalised Fourier transform,χ(p) = dx exp(ipx)χ(x), and combining the transformed equations for the resonant mode. This gives d dp with and We notice that p should not be confused with the plasma pressure, since we are considering the ultra low-β limit β ϵ. Equation ( where and χ ext m0n0 is the (m 0 , n 0 ) Fourier harmonic of the solution of the external region equation equation (18).
We now have a dispersion relation that encapsulates, to leading order, the effect of 3D geometry, and yet yields some very familiar scalings, since the problem of reconnecting low mode-number instabilities in stellarators has been reduced to an equivalent cylindrical limit.
which is the Coppi scaling for nearly ideally marginal dissipative reconnecting modes, with a geometric O(1) correction. Equation (46) is obtained from equation (42) by taking the δ i,e /δ → 1 limit. Notice that the dispersion relation equation (42) is valid for both tearing and dissipative kink instability, therefore for constant and non-constant χ-perturbations across the resonant layer.
In the collisionless case, we simply replace resistivity by electron inertia η/µ 0 → γd 2 e [39,42] (where d e = c/ω pe is the electron skin depth, with ω pe = (ne 2 /m e ϵ 0 ) 1/2 the electron plasma frequency) to obtain the scaling [43][44][45][46] γ ∼ d eĝ ψψ 00 with a boundary layer width The collisionless equivalent of the Furth Killeen Rosenbluth [39,40] result is obtained in a similar way. The results are for the growth rate, and (49) for the boundary layer width. Thus, we can say that, depending on the regime of interest, several scalings with the modified metric elementĝ ψψ 00 (and the Jacobian) are possible. Collisionless modes feature growth rates with stronger scalings inĝ ψψ 00 than resistive ones. As a result, an increase inĝ ψψ 00 would exacerbate reconnection instabilities, just as it does for ion-temperature-gradient driven electrostatic turbulence [47,48], sinceĝ ψψ 00 ∝ g ψψ 00 , and g ψψ = ∇ψ · ∇ψ measures the strength of radial gradients. This result applies even in the absence of non-resonant modes coupling. Another result, which is stellarator-specific and new, is that all harmonics contribute to the actual growth of the instability! The quantity that measures the contribution of all harmonics to the growth rate is the ratio A quantitative assessment of the effect of non-resonant harmonics is a nontrivial task, but is possible through the evaluation of the matrix G. We note that the necessary condition for these modes to be simulated numerically is that the quantity δψ in , and the relevant physical scales that define it, must be resolved. We conclude this section by stressing that ∆ ′ m0n0 > 0 is the sufficient condition for instability in the presence of only one resonant location. The stabilising effect of the average magnetic curvature [25], which causes a threshold in ∆ ′ , has been neglected owing to the low β ϵ subsidiary limit used to derive equations (12)- (13). This would be just one effect, among a plethora of other kinetic effects, which are here neglected. For large Larmor radii compared with the layer width (which is the most relevant case is fusion devices) all such kinetic effects can be treated analytically [49][50][51]. They all generally cause some finite ∆ ′ threshold with which the potential geometric thresholds will compete. In any case, regardless of the layer model used in the analysis, the quantity ∆ ′ m0n0 must be evaluated from an ideal MHD calculation treating the exterior region. In the case in which two modes are considered, of which one is resonant and the other is not, an exact external solution has been found, and this is reported in appendix C.

Conclusions
In this work, we have presented several features of the theory of magnetic reconnection in 3D magnetic fusion confinement devices putting emphasis on the mathematical nature of the problems encountered.
A set of non-linear equations for reduced resistive megnetohydrodynamics in stellarator geometry has been introduced. The equations are derived in the limit of small inverse aspect ratios ϵ 1, and finite β ∼ ε, where β is the ratio of kinetic to magnetic plasma pressure. They describe: the plasma flow evolution, via a vorticity equation, equation (8); the evolution of the component of the vector potential parallel to the equilibrium magnetic field, equation (9); and the plasma pressure fluctuations through the entropy law, equation (11). All equations are derived by using Boozer co-ordinates and include the equilibrium magnetic field mirror term to leading order. They are suitable for general 3D devices: tokamaks, stellarators, as well as hybrid, as long as the aspect ratio is large.
We studied linear reconnecting instabilities supported by the 3D reduced MHD equations derived in the subsidiary limit of small β ϵ. Equations (12)-(13) are the governing equations for magnetic reconnection this regime. Far from the resonant locations where reconnection occurs, ψ ψ mn , the perturbed vector potential (external solution) is described by a set of coupled equations of the Bessel-Plana type. Here ψ is the radial flux co-ordinate, and ψ mn is the radial location where the rotational transform is rational, ι = n/m. A key role in the 3D coupling of modes is played by a matrix constructed with the coefficients of the Fourier series in Boozer angles of the metric element g ψψ = ∇ψ · ∇ψ. In the limit of ψ approaching the resonant location ψ nm , the Fourier components of the external solution show the same asymptotic radial behaviour as the solution of the cylindrical problem first studied by Newcomb. If one resonant mode is present, 3D coupling is responsible for logarithmic singularities in the radial derivatives of all other modes, even non-resonant ones. We identified two types of reconnecting modes that can be destabilised, and show how they relate to the standard tearing and dissipative kink mode. The influence of the metric elements and other important geometric factors have been established.
The effect of 3D mode-coupling has been investigated by rigorously reducing the 3D problem to an equivalent cylindrical one. This was done for an equilibrium rotational transform with one resonant location, thus allowing for one reconnecting mode and an arbitrary number of non-resonant harmonics. The effect of geometry is that of redefining plasma inertial and dissipative scales. Explicit formulae for the growth rate, in terms of the device geometric factors, and destablisation criteria were given. by defining z = ζ q /q, and u = vz −p , where q = 1/(2p + 1). The general integral of equation (A3) is known, and so that of equation (A2), which is v = ζ 1/2 C 1/(2q) (ciζ q /q) , where C ν is a linear combination of the Bessel functions. In our case, care must be taken in the choice of which combination, since we need to evaluate lim ψ→ψ ± mn χ ′ mn , and the two limits can result in spurious phase factors that must be excluded, since all coefficients of equation (A1) are real. Let us then introduce ψ − ψ mn ≡ s, and consider χ mn ≡ χ R , for s > 0. We then have and the general solution is χ R = √ sC 1 (2i √ κs). We now choose where A R and B R are arbitrary constants. We recognise the 'small', J 1 , and 'large', Y 1 , solutions in Newcomb language. Equation (A5) implies that with a logarithmic divergence for s → 0 + . For s < 0, we introduce χ mn ≡ χ L , ρ → −s, with ρ > 0, and equation (A2) becomes The general solution is now χ L = √ ρC 1 −2 √ κρ , with C a linear combination of Bessel functions. After choosing χ L = A L √ ρJ 1 −2 √ κρ + B L √ ρY 1 −2 √ κρ , we find [52,54] and The logarithmic singularities cancel when evaluating the difference χ ′ R − χ ′ L , but determine the constant on integration of the twisting parity solution derived in section 3.2.

Appendix B The geometric matrix G
One procedure to construct the matrix G is the following one.