Bifurcations of the magnetic axis and the alternating-hyperbolic sawtooth
BBifurcations of the magnetic axis and the alternating-hyperbolic sawtooth
C. B. Smiet,
1, 2
G. J. Kramer, and S. R. Hudson Princeton Plasma Physics Laboratory, Princeton University, Princeton, New Jersey, USA Huygens-Kamerlingh Onnes Laboratory, Leiden University,P.O. Box 9504, 2300 RA Leiden, The Netherlands
We present a sawtooth model that explains observations where the central safety factor, q , stayswell below one, which is irreconcilable with current models that predict a reset to q = 1 after thecrash. We identify the structure of the field around the magnetic axis with elements of the Lie groupSL(2 , R ) and find a transition to an alternating-hyperbolic geometry when q = 2 /
3. This transitionis driven by an ideal MHD instability and leads to a chaotic magnetic field near the axis.
The sawtooth oscillation in tokamaks consists of a slowrise in core temperature lasting several to hundreds ofms, followed by a rapid crash lasting 50 − µ s. It wasfirst observed in 1974 [1] and has been observed in almostevery tokamak since [2]. Understanding and mitigatingthis temperature-limiting instability is crucial for the suc-cess of future tokamak reactors such as ITER. Existingsawtooth models predict that the crash occurs when thecentral safety factor q ∼ q = 1 surface is removed. However, a significantly lowervalue of q is sometimes observed, and some experimentsindicate the q = 1 surface is not removed. In this paperwe present a sawtooth model that predicts a crash to oc-cur by a change in topology of the magnetic field in thecore when q = 2 /
3, and which is consistent with theseoutlier observations.The field lines in an axisymmetric tokamak lie on afoliation of nested flux surfaces that confine the plasma.The single field line at the center of this foliation is calledthe magnetic axis. The winding of field lines is quantifiedby the safety factor, q , (inverse of the rotational trans-form, ı ), which quantifies the ratio of toroidal to poloidalwinding of a field line on a surface.Sawtooth models predict a fast-growing ideal instabil-ity to occur at some value of q ( q at the magnetic axis)that rapidly mixes the plasma in the core, and q sub-sequently increases. Due to the temperature-dependent(Spitzer) resistivity, current re-accumulates near the axison a (slow) resistive timescale, which decreases q untilthe crash is triggered again. This leads to a characteristicsawtooth-pattern in diagnostics sensitive to temperaturevariations that measure near the magnetic axis.The current leading models, which are the Kadomtsevand the Wesson model, both predict the crash to occurclose to q = 1 due to a fast growing mode with 1/1mode numbers. In the Kadomtsev model [3], the crash istriggered when q is slightly below 1 and a q = 1 surfaceexists in the plasma. This configuration is unstable to ei-ther an internal kink mode, a resistive tearing instabilityof the q = 1 surface, or both [4]. The hot plasma withinthe q = 1 surface mixes with cold plasma outside andis deposited in a growing 1 / q = 1 surface and q = 1 after the crash. The q -profile is flat up to the inversion radius, which lies outsideof the pre-crash q = 1 surface [3].The Wesson model [5] states that the crash is triggeredby a quasi-interchange mode that occurs when q ∼ q = 1 after thecrash. Recent numerical simulations have shown otherphenomena in the q ∼ q just above 1 [6]. Thishas led to the formulation of a new explanation of thesawtooth [7], in which the saturated 1/1 poloidal flowmaintains a flat safety factor profile with q (cid:38) q ∼ q stays below 1, and where the q = 1surface is not removed. On TEXT Lithium fluorescenceimaging was used [11] to measure q = 0 . ± .
1, and onTokapole II the differential flux was directly measuredthrough probe insertion [12] giving q = 0 . − .
8. MTXalso reported q = 0 . − . ± . . ± .
1, whilston on TFTR motional stark effect was used [17, 18] toalso measure q = 0 . ± .
1. Both of these techniques usedsimultaneously on JET gave the same result of 0 . ± . q = 1surface which is in contradiction to Kadomtsev. A clearpattern emerges (for a review see [21]): the sawtoothcrashes that occur when the safety factor is below one allcluster around the value 0.7.On Jet, long-lived ( m, n ) = (1 ,
1) density perturba-tions, ’snakes’, were observed where a deuterium pelletcrosses the q = m/n = 1 / q = 1 surface. These snakes can sur- a r X i v : . [ phy s i c s . p l a s m - ph ] J u l vive sawtooth crashes [8, 22], indicating that the q = 1surface is not removed during the crash. A model fortheir persistence has been formulated postulating a split-up and re-formation of the snake during post-crash sec-ondary reconnection [23], but a simpler explanation isthat the q = 1 surface is simply not removed. Addi-tionally, Mirnov coil measurements of ellipticity inducedAlfv´en eigenmodes (EAEs) on the JT60U tokamak showmodes localized on the q = 1 surface that are presentboth before and directly after the crash [24], also indi-cating a process that leaves the q = 1 surface intact.Many of the crashes analyzed in Ref. [24] show EAE ac-tivity consistent with removal of the q = 1 surface, butsome clearly do not.The above issues are often referred to as the problem of‘incomplete reconnection’, suggesting that the Kadomt-sev process does not run to completion. Diamagnetic ef-fects can, in theory, stabilize reconnection when the dia-magnetic drift speed exceeds a threshold [25], and sucha process would result in a poloidally asymmetric post-crash profile, as can sometimes be observed [26, 27]. Adifferent model that predicts q < q < q = 0 . not consistent with the Kadomtsev andWesson models, namely when q is not set equal to 1, arecaused by an entirely different mechanism.We use mathematical group theory to identify a hith-erto overlooked bifurcation of the magnetic axis, use idealmagnetohydrodynamcs (MHD) stability calculations toidentify the associated unstable mode, and demonstratethat this mode produces magnetic stochasticity. Fromthis emerges a new model of the sawtooth crash thatdoes not remove the q = 1 surface, that predicts the crashto occur at q = 2 /
3, and that explains the observationsthat are irreconcilable with the models of Kadomtsev andWesson.We consider a typical tokamak-like magnetic field incylindrical coordinates, (
R, φ, Z ). We take the toroidalcomponent to be always positive, so that B ·∇ φ >
0, andwe assume that there is at least one closed flux surface(a surface where B · ˆ n = 0 with ˆ n the surface normal) somewhere in the domain of interest. We assume that themagnetic field is continuous, differentiable, and changescontinuously in time.The field line mapping, also called the Poincar´e map-ping, is constructed by integrating once around the torusalong the magnetic field. Integration is started from apoint on a prescribed surface, called the Poincar´e sur-face, that is transverse to the magnetic field. We choosethe Poincar´e section to be the plane φ = 0, and writethe mapping as f ( R , Z ) = ( R , Z ). The region of thePoincar´e section that is bounded by a closed flux surface(a region that is topologically a disk), is mapped to itselfunder the field line mapping. The field line mapping iscontinuous and differentiable. Because the magnetic fieldin a tokamak constitutes a Hamiltonian dynamical sys-tem [29], the field line map is the Poincar´e map or firstreturn map of this dynamical system.Brouwer’s fixed point theorem states that any contin-uous map from a disk to itself has at least one fixedpoint [30]. The magnetic axis is an example of a fixedpoint, as are the points on an intact q = 1 surface.At a fixed point we construct the Jacobian matrix ofpartial derivatives, M = ∂R ∂R , ∂R ∂Z ∂Z ∂R , ∂Z ∂R . (1)The matrix M describes to first order the behavior ofnearby field lines [31], f ( x + δ x ) ≈ f ( x ) + M · δ x , (2)where x = ( R, Z ) T . This matrix can be used to accel-erate convergence of iterative methods for finding fixedpoints and to compute the Lyapunov exponent of themagnetic field lines.At a fixed point the divergence free condition on themagnetic field, ∇· B = 0, guarantees that this map is areapreserving, so det( M ) = 1. Therefore M ∈ SL(2 , R ), thegroup of 2 × , R ) is a connected simple Lie group, whose elementscan be classified into elliptic, hyperbolic and parabolicsubsets by their action on the Euclidean plane.The eigenvalues of M are solutions to the characteris-tic polynomial λ ± = (Tr( M ) ± (cid:112) Tr( M ) − /
2. Eigen-values are complex when | Tr( M ) | < | Tr( M ) | ≥
2. When Tr( M ) > M is a hyperbolic element of SL(2 , R ). Thefixed point is surrounded by hyperbola that are left in-variant under this mapping and the (real) eigenvectors of M determine the directions in which points, upon consec-utive iterations of the field line map, converge to (formingthe slow manifold) respectively diverge away from (form-ing the fast manifold) the fixed point. When | Tr( M ) | < FIG. 1. Classification of the elements of SL(2 , R ). (a), (b) and(c): illustration of a typical alternating-hyperbolic, ellipticand hyperbolic element element respectively. Eigenvectors of M are shown in black, an example surface invariant underthe mapping in blue, and the mapping of select points areillustrated by the colored vectors. (d): Properties of elementsof SL(2 , R ) as a function of Tr( M ). For the elliptic fixed pointsthe rotational transform ı (left axes) and safety factor q (rightaxes) is shown in dashed red. Greene’s residue (left axes) isshown in green. the configuration of the field is that of an O-point, asillustrated in Fig. 1(b), and M is an elliptic element ofSL(2 , R ). The fixed point is surrounded by invariant el-lipses, which constitute the magnetic surfaces surround-ing the fixed point. When Tr( M ) < − M are negative. This mapping sends points on invarianthyperbolic surfaces to the branch on the opposite side ofthe fixed point, as illustrated in Fig. 1(a). M is also ahyperbolic element of SL(2 , R ), but because its map in-cludes a reflection, we call the fixed point an alternating-hyperbolic fixed point. When | Tr( M ) | = 2 the configura-tions include the identity mapping and shear mappingssuch as occur at an intact q = 1 surface. In this case M is called a parabolic element of SL(2 , R ).The field near an O-point is mapped from an elliptic(magnetic) surface to itself with a certain average rota-tion. This determines the local safety factor q (up to themultiplicity of the solution) of this fixed point throughcos(2 π/q ) = Tr( M ) . (3)This relation was given by Greene as sin (2 π/ q ) = R [31, 32], where Greene’s residue R = − Tr( M ).There must always be one more fixed point with pos-itive residue than there are fixed points with negativeresidue[31, 33]. The structure of SL(2 , R ), as well asthe residue and average angle of rotation are shown in Fig. 1(d).Different configurations of the magnetic field corre-spond to different elements of this group, and the struc-ture of this group indicate configurations when the topol-ogy of the magnetic field can change.As the magnetic field changes continuously in time, theposition of fixed points and the elements of the matrix M change continuously [33]. Exceptions to this occurwhen fixed points appear or disappear (always in pairs)or when a continuous line of fixed points breaks into anisland chain [33]. M traces a path through SL(2 , R ) andcontinuous functions on the Lie group, e.g. Tr( M ) and q ,also change continuously. There are two values of Tr( M )where M can transition between the different subsets ofSL(2 , R ) and the topology of the field around the fixedpoint changes. For example, the Kadomtsev process in-volves fixed points on the q = 1 surface, with Tr( M ) = +2and with M constituting a shear mapping, which changeinto a hyperbolic fixed point, an X-point with Tr( M ) (cid:38) ,
1) islandwith Tr( M ) (cid:46)
2, which is created with a local value of q = 1. In the Kadomtsev model this point goes on tobecome the new axis.When Tr( M ) = − q = 2 / (1 + 2 n ). In thiscase, the field lines make an integer number and a half rotations around the axis. If this occurs at the magneticaxis, the axis becomes an alternating-hyperbolic X-point.This can happen when q = 2 /
3, which is very close to q = 0 . ± .
1, the value at which experiments have shownsome crashes to be triggered [16, 17].In order to explain the rapidity of the crash we re-quire a fast instability that drives this transition. Wewill now study this situation using ideal MHD theory,and we will present a calculation that illustrates thatalternating-hyperbolic fixed points are prone to the for-mation of chaos.We construct two families of MHD equilibria each witha circular-cross section, a minor radius equal to 1m, a ma-jor radius equal to 3m, a plasma β of 3% on the magneticaxis, and a toroidal field of 1T. The pressure profile isquadratic in minor radius with a maximum on axis andzero on edge, and uniquely defined by the constraint that β = 3% on the axis. The first family of equilibria hasa q -profile given by q = q + 2 . ψ p , where ψ p is the nor-malized poloidal flux function. The safety factor on axisis thus q , and since ( r/a ) ≈ (cid:112) ψ p , it increases quadrat-ically with a value on the edge of q a = q + 2 . ≈ . R = 2 . r = 0 . β = 1 −
8% anda safety factor on edge q a = 3 −
6. The magnetic field islower than on TFTR which was at maximum 6T.We construct a second family of equilibria with a shal-lower safety factor profile, quartic in minor radius, given
FIG. 2. Ideal growth rates as a function of q on the magneticaxis calculated by the NOVA-K code. The first family with q = q + 2 . ψ p is given in red solid lines, and the secondfamily with q = q +0 . ψ p in black dashed lines. Insets: radialcomponent of the displacement vector of the ideal modes. Thearrows indicate the direction of the displacement. by q = q + 0 . ψ p . This profile has a low value of q a that would be challenging to achieve experimentally, butthe shallow profile highlights modes that occur aroundthe axis. We investigate the ideal stability of these MHDequilibria with q ∈ [0 . , .
1] using the well-benchmarkedNOVA-K code [34].We perform NOVA-K calculations in the ideal mode,not including kinetic effects. The code finds ideal MHDmodes in equilibria by solving the normal mode formula-tion of the linearized ideal MHD stability equations. Thenormal mode equation is given by − ω ρ ξ = F [ ξ ] , (4)where ω is the mode frequency squared, ρ is the plasmadensity, F is the linearized MHD force operator, and ξ is the displacement vector. Solutions with ω < / q <
1, and a 2 / q ∼ / q <
1, as is expected from literature. Itis this mode that drives the Kadomtsev and Wesson pro-cesses, displacing the hot core within the q = 1 surface.In the second family of equilibria (black curves) the q -profile is much flatter and shallower, and the 1 / q = 1 surface. In both families, exactly when q reaches FIG. 3. Effect of the 2 / A = 4 × − the magnetic topology aroundthe axis changes into the alternating-hyperbolic configuration.(b): At a higher amplitude A = 0 .
01 a finite region aroundthe core stochastisizes. / q = 2 /
3, fieldlines make one and a half rotations around the axis, andthus end up on the exact opposite side of the axis fromwhence they start. A perturbation that causes this tran-sition would be one for which a set of field lines (i.e. theslow manifold) is mapped closer to the fixed point (endingup opposite the axis but closer), and another set (i.e. thefast manifold) is mapped further away. This is what theideal displacement of a 2/3 mode accomplishes. The dis-placement is directed towards the axis in two directions,and away from the axis in the two others, and rotatingas a function of toroidal angle with the same rate as thefield lines, such that a field line is resonantly displacedtowards or away from the axis.To demonstrate this we will show the effect of an an-alytic 2/3 mode on the magnetic topology by perturb-ing a tokamak equilibrium field. We construct the equi-librium with a quadratic safety factor profile given by q = 2 / . ∗ ψ p , and all other parameters iden-tical to the families of equilibria described above. Theanalytic displacement is calculated from:Ψ = A exp (cid:18) − ψ p σ (cid:19) cos(2 θ − φ ) (5)where θ is the poloidal angle, A is an overall scaling fac-tor, and the radial envelope has a characteristic width σ = 0 .
15. The displacement vector ξ is calculatedthrough ξ R = − (1 /R ) ∂ Z Ψ and ξ Z = (1 /R ) ∂ R Ψ. In thenormal mode formulation of the ideal MHD equations,the perturbed field that corresponds with a displacementis given by δ B = ∇ × ( ξ × B ) (6)where B is the equilibrium magnetic field.We analyze the structure of the magnetic field by cal-culating the trajectories of low-energy (1keV) electronsusing the SPIRAL code [35]. The drift orbit of theselow-energy electrons is so small that they are effectivelyfield-line following. We calculate the Poincar´e map 1000times per trajectory for 127 starting points equidistantlyspaced between the equilibrium axis and the plasmaboundary.Poincar´e plots of B + δ B are shown in Fig. 3. Thescaling factor A = 4 × − results in a perturbed field | δ B | / | B | ∼ × − . At this amplitude the nested fluxsurfaces near the axis are broken up through resonancewith the applied perturbation as shown in Fig. 3(a). Notethat the center of the two ‘islands’ are themselves notfixed points of the field line map: the one maps to theother. The axis is an alternating-hyperbolic fixed point.At a value of A = 0 .
01, which corresponds to | δ B | / | B | ∼ . × − , the field around the magneticaxis becomes chaotic, shown in Fig. 3(b). A 3/3 islandchain is seen at the q = 1 surface. The 1 / q is measured to be 0.7 still often show a 1/1 temperaturedistribution during the crash phase [18, 36]. Stochastisa-tion happens for a pure 2 / q decreases on a slow, resistive timescale.When the 1 / q to 1, the safety factor decreasesto near q = 2 /
3. The ideal 2 / q = 1 surface (consistent withan inversion radius smaller than the radius of the q = 1surface [16, 20]), leaving that surface intact during thecrash. The poloidal and toroidal fluxes are redistributed,resulting in an increased safety factor. This shifts thefield out of resonance with the mode and the core re-gion heals with 2 / < q <
1. After the crash, currentdiffusion again slowly decreases the safety factor until q = 2 / q sawteeth occur clustered around the value of q = 0 . ± . q = 2 / q = 1surface is unaffected by this crash. The same holds forEAEs [24] that live on the q = 1 surface and are presentimmediately after the crash.One may wonder how the tokamak reaches a state with q this significantly below 1, but direct measurements tellus it can [11, 12, 14–18, 20] There are several reportedmechanisms that stabilize the internal kink mode [39], in-cluding toroidal rotation [40], fast particles [41, 42], anddiamagnetic effects [43]. Stabilization is experimentallyconfirmed by the decrease of sawtooth repetition ratethrough fast particle injection [44, 45]. (In such stabi-lized regimes, sawteeth with very low repetition rate canalso be observed, so-called giant sawteeth, [45, 46], whichdo exhibit a re-set to q = 1. These can be triggered bythe onset of energetic particle modes that remove thestabilizing fast particles [46].)In this paper we have focused on the sawtooth eventsthat occur when q = 2 /
3, but as figure 1 shows, a tran-sition of the axis to a(n) (alternating-)hyperbolic pointcan occur when q = 1 /n , (respectively q = 1 / (1 / n )).The new explanation of the sawtooth proposed by Jardinet al. [7] posits that crashes occur in a discharge where q is clamped just above 1 through a nonlinearly satu-rated quasi-interchange flow, and are triggered by cross-ing the threshold of the 2 / q = 1 would drive the transitionto a (regular) hyperbolic geometry. The presented nu-merical simulation (see fig. 6 (b)) exhibits a stochasticcore region surrounding a 2/2 structure similar to fig-ure3(b) which causes the crash. Sawtooth-like relaxationevents can also be observed when q = 2 [47] and aretraditionally attributed to a 2/1 double tearing instabil-ity. The present work suggests that a transition to a(n)(alternating-)hyperbolic transition could form a unifiedtheory that explains all three phenomena.The alternating-hyperbolic sawtooth model can possi-bly be measured on tokamaks with sufficiently fast di-agnostics. SRX tomography (f.ex. [36]) shows that thecrash phase is often (but not always) dominated by a 1/1signature. This would suggest that the 2/3 mode can actas a very fast trigger, which in turn sets off the 1/1 mode.The SRX data in [36] does show a few chords with a fasteroscillation frequency just before temperature equilibra-tion is initiated, but this instant is not tomographicallyreconstructed. On a modern tokamak electron cyclotronemission imaging (ECEI) would be fast enough to discernthe structure of the mode that initiates the crash.Most measurements of q ∼ . q ∼ q = 0 . q = 2 / [1] S. Von Goeler, W. Stodiek, and N. Sauthoff, PhysicalReview Letters , 1201 (1974).[2] I. T. Chapman, Plasma Physics and Controlled Fusion , 013001 (2010).[3] B. B. Kadomtsev, Soviet Journal of Plasma Physics ,710 (1975).[4] B. Coppi, R. Galvao, R. Pellat, M. Rosenbluth, and P. Rutherford, Fizika Plazmy , 961 (1976).[5] J. A. Wesson, Plasma physics and controlled fusion ,243 (1986).[6] S. C. Jardin, N. Ferraro, and I. Krebs, Physical reviewletters , 215001 (2015).[7] S. Jardin, I. Krebs, and N. Ferraro, Physics of Plasmas , 032509 (2020).[8] A. Weller, A. D. Cheetham, A. W. Edwards, R. D. Gill,A. Gondhalekar, R. S. Granetz, J. Snipes, and J. A.Wesson, Physical review letters , 2303 (1987).[9] D. Wroblewski and L. Lao, Physics of Fluids B: PlasmaPhysics , 2877 (1991).[10] Y. B. Nam, J. S. Ko, G. H. Choe, Y. Bae, M. J. Choi,W. Lee, G. S. Yun, S. Jardin, and H. Park, NuclearFusion , 066009 (2018).[11] W. West, D. Thomas, J. DeGrassie, and S. Zheng, Phys-ical review letters , 2758 (1987).[12] T. Osborne, R. Dexter, and S. C. Prager, Physical Re-view Letters , 734 (1982).[13] B. Rice and E. Hooper, Nuclear fusion , 1 (1994).[14] H. Soltwisch, Review of Scientific Instruments , 1939(1986).[15] H. Soltwisch, Review of Scientific Instruments , 1599(1988).[16] H. Soltwisch and H. Koslowski, Plasma physics and con-trolled fusion , 667 (1995).[17] F. M. Levinton, S. H. Batha, M. Yamada, and M. Zarn-storff, Physics of Fluids B: Plasma Physics , 2554(1993).[18] M. Yamada, F. Levinton, N. Pomphrey, R. Budny,J. Manickam, and Y. Nagayama, Physics of plasmas ,3269 (1994).[19] R. Wolf, J. O’Rourke, A. Edwards, and M. Von Heller-mann, Nuclear fusion , 663 (1993).[20] B. Rice, Review of scientific instruments , 5002 (1992).[21] H. K. Park, Advances in Physics: X , 1633956 (2019).[22] R. D. Gill, A. W. Edwards, D. Pasini, and A. Weller,Nuclear Fusion , 723 (1992).[23] D. Biskamp and J. Drake, Physical review letters , 971(1994).[24] G. J. Kramer, C. Z. Cheng, Y. Kusama, R. Nazikian,S. Takeji, and K. Tobita, Nuclear fusion , 1135 (2001).[25] M. T. Beidler and P. A. Cassak, Physical review letters , 255002 (2011).[26] I. Furno, C. Angioni, F. Porcelli, H. Weisen, R. Behn,T. Goodman, M. Henderson, Z. Pietrzyk, A. Pochelon,H. Reimerdes, et al. , Nuclear Fusion , 403 (2001).[27] Z. Pietrzyk, A. Pochelon, T. Goodman, M. Henderson,J.-P. Hogge, H. Reimerdes, M. Tran, R. Behn, I. Furno,J.-M. Moret, et al. , Nuclear Fusion , 587 (1999).[28] Y. I. Kolesnichenko, Y. V. Yakovenko, D. Anderson,M. Lisak, and F. Wising, Physical review letters ,3881 (1992).[29] A. H. Boozer, Magnetic field line Hamiltonian , Tech.Rep. (Princeton Univ., 1985).[30] L. E. J. Brouwer, Mathematische Annalen , 97 (1911).[31] J. M. Greene, Journal of Mathematical Physics , 760(1968).[32] J. M. Greene, Journal of Mathematical Physics , 1183(1979).[33] C. B. Smiet, G. J. Kramer, and S. R. Hudson, PlasmaPhysics and Controlled Fusion , under review (2019).[34] C. Cheng, Physics reports , 1 (1992).[35] G. J. Kramer, R. V. Budny, A. Bortolon, E. D. Fredrick- son, G. Y. Fu, W. W. Heidbrink, R. Nazikian, E. Valeo,and M. Van Zeeland, Plasma Physics and Controlled Fu-sion , 025013 (2013).[36] Y. Nagayama, M. Yamada, W. Park, E. Fredrickson,A. Janos, K. McGuire, and G. Taylor, Physics of Plas-mas , 1647 (1996).[37] A. H. Boozer, Physics of Plasmas , 072907 (2014).[38] Y.-M. Huang, A. Bhattacharjee, and A. H. Boozer, TheAstrophysical Journal , 106 (2014).[39] I. T. Chapman, S. D. Pinches, J. P. Graves, R. J. Ak-ers, L. C. Appel, R. V. Budny, S. Coda, N. J. Conway,M. De Bock, L. G. Eriksson, et al. , Plasma Physics andControlled Fusion , B385 (2007).[40] I. T. Chapman, T. C. Hender, S. Saarelma, S. E. Shara-pov, R. J. Akers, N. J. Conway, M. Team, et al. , Nuclearfusion , 1009 (2006).[41] M. Cole, A. Mishchenko, A. K¨onies, R. Kleiber, andM. Borchardt, Physics of Plasmas , 072123 (2014). [42] F. Porcelli, Plasma Physics and Controlled Fusion ,1601 (1991).[43] J. W. Connor, R. J. Hastie, and A. Zocco, PlasmaPhysics and Controlled Fusion , 035003 (2012).[44] C. Angioni, A. Pochelon, N. N. Gorelenkov, K. G. Mc-Clements, O. Sauter, R. V. Budny, P. C. De Vries, D. F.Howell, M. Mantsinen, M. F. F. Nave, et al. , Plasmaphysics and controlled fusion , 205 (2002).[45] D. J. Campbell, D. F. H. Start, J. A. Wesson, D. V.Bartlett, V. P. Bhatnagar, M. Bures, J. G. Cordey, G. A.Cottrell, P. A. Dupperex, A. W. Edwards, et al. , Physicalreview letters , 2148 (1988).[46] S. Bernabei, R. V. Budny, E. D. Fredrickson, N. N. Gore-lenkov, J. C. Hosea, C. K. Phillips, R. B. White, J. R.Wilson, C. C. Petty, R. I. Pinsker, et al. , Nuclear Fusion , 513 (2001).[47] Z. Chang, W. Park, E. D. Fredrickson, S. H. Batha, M. G.Bell, R. Bell, R. V. Budny, C. E. Bush, A. Janos, F. M.Levinton, et al. , Physical review letters77