Isochronous bifurcations of magnetic islands in tokamaks

In the phase space of twist Hamiltonian systems, the emergence of different isochronous island chains on the same rational surface is observed when multiple resonant perturbations with different wave numbers and the same winding numbers are applied.^1–3 All the concerned islands in the chains, called isochronous islands,² have the same winding number, which means they have the same average speed around an invariant circle. The transformation from one of these chains to another occurs through bifurcations called isochronous. Such isochronous islands and bifurcations have been found in several systems—see Refs. 1–3 and references cited therein.

In tokamaks, the magnetohydrodynamics equations for the magnetic field lines can be, in general, described as a twist Hamiltonian system after defining a suitable field line Hamiltonian.⁴ Thus, the emergence of isochronous bifurcations can be expected for multiple resonant perturbations with the same winding number, as considered in Refs. 5 and 6. In a Poincaré section, the magnetic island chains can be obtained by integrating the field line equation.

Recently, studies with visco-resistive magnetohydrodynamic (MHD) simulations have investigated the role of resonant magnetic perturbations (RMPs) on the onset of magnetic islands in the isochronous bifurcations at the NSTX-U tokamak.^5,6 In Ref. 5, a series of Poincaré sections were calculated by following magnetic field lines for two different approaches: (i) superposing the equilibrium axisymmetric field with the externally applied magnetic perturbation, but without considering the plasma response to this perturbation (the so-called vacuum approach); and (ii) superposing the equilibrium axisymmetric field with both the externally applied magnetic perturbation and the field due to the response of a visco-resistive single-fluid plasma. Perturbations create the reported islands due to the application of resonant external coils on the plasma equilibrium perturbed by an internal mode described by a linear resistive MHD approximation.^5,6 These calculations show the existence of island bifurcations. Evans et al.,⁶ however, have shown that island elongation is followed by a sequence of isochronous bifurcations when the perturbation current is increased.

Experimental observation of isochronous island chains during DIII-D plasma discharges has also been reported, reinforcing that further studies on the impact of RMPs are required to control these heteroclinic structures.⁷ These studies described a bifurcation from heteroclinic to homoclinic topology in DIII-D plasmas.

To study these effects on magnetic topology, we use a model that can generate isochronous bifurcations similar to the ones above when two RMPs are superposed for a cylindrical plasma with first-order correction on toroidicity. The first perturbation is created by a pair of resonant helical windings (RHWs) located on the external wall. It is usually applied to reduce wall erosion by generating a chaotic layer in the plasma edge.^8–13 The second perturbation is created by a helical current sheet (HCS) located on a rational flux surface that acts as a plasma current perturbation.^14,15 Since a first-order correction on toroidicity is applied to this cylindrical plasma model, an RMP field will result in the formation of magnetic islands at the resonant rational surface (the primary mode associated with the RMP) and also on the neighboring rational surfaces (the secondary modes associated with the RMP due to poloidal mode coupling induced by toroidicity). Similarly to the study reported in Ref. 5, the resulting magnetic topology around a rational surface in this work is a combination of primary and secondary modes caused by RHWs and HCS. Here, we estimate the current amplitude required for the onset of isochronous islands around rational surfaces.

Additionally, the analytical expressions for the magnetic field are obtained for a cylindrical plasma in equilibrium with first-order corrections on toroidicity included and with the superposition of perturbations due to a pair of RHWs¹⁶ and an HCS.^17–20 Solving the field line equations allows us to calculate the Poincaré maps for a plasma cross section, considering different combinations of parameters for specified perturbations. We have identified isochronous bifurcations located around rational flux surfaces, all of them modifying the distribution of field lines similarly to the bifurcations observed in visco-resistive MHD simulations^5,6 and in experiments.⁷ 

This paper is structured as follows: Sec. II presents the MHD equilibrium for a cylindrical plasma with first-order correction due to toroidicity. The magnetic field generated by a pair of RHWs and an HCS is calculated in Secs. III and IV, respectively. Section V shows the Poincaré sections of field lines, while conclusions are presented in Sec. VI. The  Appendix details the numerical procedure used to calculate the Poincaré sections and to identify the isochronous bifurcations of Sec. V.

The parameters used in Secs. III–VI are related to the operation of the TCABR tokamak:²²  I_p = 100 kA, I_e = 4 MA, a = 0.2 m, R₀ = 0.6 m, and γ = 3. Figure 1 shows the safety factor radial profile considering Eq. (6) and the TCABR parameters. This figure also shows that radial positions (r/a) 0.610, 0.800, and 0.925 are approximately associated with safety factor values of 2, 3, and 4, respectively. In the following sections, trajectories of magnetic field lines around these radial positions are analyzed.

A pair of RHWs generates an external perturbation that is superposed to the equilibrium magnetic field. It consists of two helical wires located on the surface r = b separated by 180° and with current I_h in opposing directions, as indicated in Fig. 2.^16,23,24

Although the current distributions in this work are located on a surface—in Eq. (9), for instance, the surface is r = b—considering a current layer of finite thickness would result in more realistic values for the resonance amplitudes and the corresponding Poincaré maps. However, these alterations should not change the resonance sequence and the observed bifurcations.

In the low beta plasma considered in this article, the solution of Eq. (10) determines approximately the perturbing magnetic field in the vacuum approximation, obtained by neglecting the effect of viscous and resistive plasma reaction.

This perturbation excites the main mode (m_h,n_h) on the magnetic rational surface associated with the safety factor q = m_h/n_h and other sideband modes due to the toroidal correction.

In this work, a perturbation to the plasma current is represented as an RMP generated by a HCS located at the equilibrium rational surface r = r_s, with r_s < b and q(r_s) = m_s/n_s (m_s and n_s here are the mode numbers associated with this perturbation).¹⁵  Figure 3 shows a schematic diagram of the surfaces r = b and r = r_s with a representation of the helical direction in which the current density is applied.

In Eq. (16), the parameter $α s= B 0 θ / B 0 z$ defines the helical direction of the equilibrium magnetic field on the radial position r = r_s.

This perturbation excites the main mode (m_s,n_s) on the magnetic rational surface associated with the safety factor q = m_s/n_s and other sideband modes due to the toroidal correction.

The pair of RHWs is resonant to mode (m_h,n_h) = (2,1) and has a perturbation parameter I_h/I_p = 0.4%, while the plasma current perturbation is resonant to mode (m_s,n_s) = (8,2) and has amplitude j. The chosen wave numbers are among those observed in usual tokamak discharges. The toroidal correction of our cylindrical approach introduces secondary sideband resonant modes interacting with the dominant modes.

In this article, the “vacuum approximation” has been used to calculate the magnetic field structure, although plasma response effects may result in a screening (or even an amplification) of the external perturbations. However, for low beta plasma and non-marginally stable modes, as considered here, this effect is not so noticeable.^27,28

A set of initial values of radial and angular coordinates is chosen at the plane z = 0, and Eq. (25) is integrated numerically to obtain the Poincaré sections showing the intersection of field lines with this plane. The chosen wave numbers are among those observed in usual tokamak discharges. The  Appendix provides a detailed explanation of the numerical procedure used to solve Eq. (25), generate the Poincaré maps, and identify the bifurcations.

A complete toroidal description would make the identification of bifurcations more difficult than our cylindrical approach with toroidal corrections, which introduces the sideband modes but underestimates their amplitudes. In fact, for large-amplitude perturbations, the chaotic area increases while the island size is reduced, making the identification of bifurcations more difficult.

In the Poincaré sections, we consider the isochronous island chains (islands with the same winding number) created by the primary and secondary modes. Thus, we analyze the island chains on the q = m_h/n_h magnetic surface due to the primary mode (m_h,n_h) generated by the external RHWs and the secondary mode (m_s,n_s) generated by the internal plasma current perturbation. The competition between these primary and secondary modes gives rise to a sequence of isochronous bifurcations, on the considered rational surfaces, as the mode amplitude changes. These bifurcations modify the number of island chains formed by isochronous islands. As identified in Secs. V B–V E, the changes in the island chains occur due to local bifurcations²⁹ that change fixed point stabilities and create new fixed points,³ changing the number of island chains.

Initially, we analyze the distribution of field lines across the z = 0 plane when the perturbation is generated only by RHWs (in this first case, the perturbation due to the HCS is not considered) and then by an HCS (in this second case, the perturbation due to the RHWs is not considered).

Figure 4(a) shows the Poincaré section associated with the first case, considering I_h/I_p = 0.4% and I_s = 0. Two primary magnetic islands associated with the RHWs primary mode are formed at the q = 2 surface (r/a = 0.610). The toroidal correction applied to the axial component of the equilibrium magnetic field causes smaller secondary islands to appear on the neighboring rational surfaces: four magnetic islands are located at the q = 4 surface (r/a = 0.925), while three islands are observed at the q = 3 surface (r/a = 0.800).

When only the HCS is present, as shown in Fig. 4(b), eight primary magnetic islands are formed at the q = 4 surface (r/a = 0.925). Secondary islands also form on different rational surfaces due to toroidicity induced poloidal mode coupling. Note that two sets of isochronous islands appear on both the q = 4 surface (r/a = 0.925—green and magenta islands), associated with the HCS primary mode, and on the q = 3 surface (r/a = 0.800—red and blue islands), associated with an HCS secondary mode. In fact, as indicated in Fig. 4(b), field lines located on the green islands do not intercept the magenta ones at the q = 4 surface, while blue islands do not intercept the red ones at the q = 3 surface.

When the superposition of perturbations due to RHWs and HCS is considered, as shown in the Poincaré section in Fig. 5, two main magnetic islands are formed at the q = 2 surface (r/a = 0.610). Smaller islands are formed on neighboring rational surfaces, indicating that the trajectory of field lines is similar to the case without the HCS perturbation, as shown in Fig. 4(a). The magnetic islands located around r/a = 0.925 and r/a = 0.800, associated with safety factor values 4 and 3, respectively, are only slightly distorted compared to the Poincaré section in Fig. 4(a). However, as will be shown in the next figures, these islands bifurcate, generating new islands if the amplitude I_s is large enough. In these two rational surfaces, we observe isochronous islands when the perturbation is generated only by an HCS located at the q = 4 surface (r_s/a = 0.925), as shown in Fig. 4(b).

In this section, we consider the primary modes (m_h,n_h) = (2,1) and (m_s,n_s) = (8,2) and present examples of topological bifurcations that are generated when the current I_s is increased and I_h is kept constant. In all cases considered here, the perturbation is generated by RHWs and HCS, with the phase angle (⁠ $ϕ p$⁠) equal to zero. The examples correspond to bifurcations observed at two different rational surfaces with q = 3 and q = 4.

When the I_s/I_p is increased from 0.116% to 0.232%, while I_h/I_p is kept constant at 0.4%, a bifurcation is observed at the q = 3 surface (r/a = 0.800), as shown in the Poincaré section in Fig. 6(a). In this case, each of the three main magnetic islands (represented in red), which also appear in Figs. 4(a) and 5, now contains two additional smaller isochronous islands (green and blue), with the magnetic topology resulting in a combination of secondary modes associated with the RHWs (three red islands) and HCS [six isochronous islands in blue and green, also seen in Fig. 4(b)]. Field lines on the green region do not access or intercept the blue one. The observed elongation and bifurcation of magnetic islands is consistent with isochronous bifurcations observed experimentally⁷ and numerically investigated with resistive MHD simulations.^5,6 Figure 6(b) shows an enlarged view of the Poincaré section around r/a = 0.800. When I_s/I_p is increased to 4.65%, as shown in Fig. 6(c), the green and blue islands become larger and, instead of the main magnetic islands observed in Figs. 6(a) and 6(b), there is a thick chaotic layer around r/a = 0.800, highlighting the predominant role of the HCS secondary mode (six isochronous islands) in shaping the magnetic topology around this rational surface.

Figure 7 shows a series of Poincaré sections of field lines around the q = 4 surface (r/a = 0.925), considering that I_s/I_p is modified in the range from 0.291% to 0.407%, while the I_h/I_p is kept at 0.4% for all cases. In Fig. 7(a), four magnetic islands (represented in magenta) are located around r/a = 0.925, indicating that the RHW secondary mode (four islands) is predominant on topology. The distortion of field lines observed around r/a = 0.925 is similar to the modification of magnetic islands around rational surfaces observed in resistive MHD simulations (see Fig. 3 of Ref. 5).

When I_s/I_p is increased to 0.314%, as shown in Fig. 7(b), a bifurcation is observed in the q = 4 surface: two small islands (green and blue) are formed inside the main ones (magenta) around r/a = 0.925, which is also the region where two sets of isochronous islands are located when only the perturbation due to the HCS is considered, as shown in Fig. 4(b). In Fig. 7(c), it is shown that, when I_s/I_p is modified to 0.320%, another bifurcation is observed: the secondary (green and blue) islands are completely separated from the main islands (magenta), and they are delimited by a magnetic surface (red) that does not connect to the main islands. Thus, we have one chain with 4 islands and the other with 8, a total of 12 islands. However, when the I_s/I_p is increased to 0.407%, only two sets of isochronous islands are observed around r/a = 0.925, represented in magenta (4 islands) and red (4 islands) in Fig. 7(d). In this case, the influence of the HCS secondary mode, characterized by two sets of 4 isochronous islands, is predominant on the magnetic topology [as observed in Fig. 4(b) for I_h = 0], without the secondary islands created by the RHWs in Fig. 4(a) for I_s = 0.

Considering that a bifurcation is observed around the rational surface with q = 3 when I_s/I_p = 0.232% and that another set of isochronous islands is observed around the surface with q = 4 when I_s/I_p  $≥$ 0.314%, our results are consistent with those reported in Ref. 5, in which isochronous bifurcations are observed on different rational surfaces as the perturbation parameter was increased.

Figure 8 shows an enlarged view of the Poincaré sections in Fig. 7 in the region 0.923 ≤ r_s/a ≤ 0.929 and 0.25 ≤ θ/(2π) ≤ 0.35, considering the same perturbation parameters of Fig. 7. The described bifurcation process can be confirmed; in other words, Figs. 8(a)–8(d), clarifying the observed local bifurcations, are magnifications of Figs. 7(a)–7(d) around a fixed point.

In Sec. V B, we analyzed the effect of increasing I_s on the distribution of field lines, supposing that the phase angle in Eq. (12) was $ϕ p=0$⁠. We now study the effect on the magnetic topology of a perturbation created by an HCS, located at the surface r/a = 0.925 with $ϕ p=π$ and (m_s,n_s) = (8,2), superposed to a perturbation created by a pair of RHWs located at r = b with (m_h,n_h) = (2,1). While I_h/I_p is kept constant at 0.4%, we change the parameter I_s/I_p from 0.116% to 0.232% and analyze the Poincaré sections of field lines, as shown in Fig. 9.

Figure 9(a) shows that four secondary islands (associated with the RHWs) are formed, when I_s/I_p = 0.116%, at r/a = 0.925. When I_s/I_p is modified to 0.232%, a bifurcation is observed, as shown in Fig. 9(b): small islands (represented in green and blue) are formed inside the main magnetic islands (magenta). There are, however, two main differences between the bifurcations observed in Figs. 7(b) (null phase angle) and 9(b) (finite phase angle): first, the two magnetic islands represented as green and blue surfaces in Fig. 9(b), for finite phase angle, appeared when the current I_s was set at a value approximately equal to 74% of the one associated with Fig. 7(b); second, the small islands (green and blue) in Fig. 7(b) are located outside the main magnetic islands [represented in magenta in Fig. 7(b)], while in Fig. 9(b), for a finite phase angle, they are formed inside the main islands.

This example shows that the phase angle between the two perturbing currents, which can be adjusted by modifying the RHW position, alters the bifurcations in the plasma.

The current expressed in Eq. (26) has a reversed direction in comparison to the one shown in Eq. (9). Considering that I_h/I_p = 0.4%, and for RHWs with (m_h,n_h) = (2,1), Fig. 10 shows the Poincaré section of field lines for the z = 0 plane, initially without the effect of the HCS. In Fig. 10(a), two islands are formed around r/a = 0.610, and smaller islands are generated around the radial positions with rational safety factor values due to the inclusion of toroidicity correction on the axial component of the magnetic equilibrium field. These islands are similar to the results shown in Fig. 4(a), with the difference that the positioning of magnetic islands is shifted along the angular direction due to the modification of electric wires on the external wall. In Fig. 10(b), we show that four magnetic islands, represented in magenta, are formed around r_s/a = 0.925. In the following section, we analyze an isochronous bifurcation observed around this region when the HCS is considered.

The influence of two perturbations on the distribution of magnetic field lines is analyzed in this section: the first one is caused by a pair of RHWs with (m_h,n_h) = (2,1) whose current is described by Eq. (26), while the second one is generated by a plasma current perturbation modeled as an HCS with zero phase angle (⁠ $ϕ p=0$⁠) and (m_s,n_s) = (8,2).

Figure 11 shows the Poincaré plot of field lines for two cases: in the first one, shown in Fig. 11(a), I_s/I_p is equal to 0.116%, while in the second one, shown in Fig. 11(b), I_s/I_p is equal to 0.174%. In Fig. 11(a), the magnetic islands are distorted around r/a = 0.925 in comparison to the Poincaré section shown in Fig. 10(b), with the RHW secondary mode shaping the distribution of magnetic islands around this rational surface. When I_s/I_p is increased to 0.174% as shown in Fig. 11(b), the effect of the HCS secondary mode on magnetic topology is amplified, and a new set of four red magnetic islands around r/a = 0.925 is observed. This type of bifurcation, caused by the HCS, is similar to the ones reported in numerical and experimental studies.^5–7 

In plasmas confined in large aspect ratio tokamaks, represented as a cylindrical geometry with toroidicity included as a first-order correction, we considered magnetic islands created by a pair of resonant helical windings (RHWs), located on the external wall, and a helical current sheet (HCS), representing a plasma current perturbation near the plasma edge. Separately, these perturbations created resonant islands on the plasma's rational magnetic surfaces. When the two perturbations acted together, the resultant island's configuration changed through bifurcations as the current amplitudes were varied. In this work, we investigated isochronous bifurcations on a given rational magnetic surface for which the previous and final islands are isochronous, i.e., have the same winding number.

Poincaré sections of magnetic field lines were calculated by integrating the field line equation, and topological bifurcations were observed on sets of magnetic islands when one perturbing current was modified. Isochronous bifurcations were triggered when the current amplitude associated with the HCS was large enough while the current on the RHWs was kept constant, reorganizing field lines on magnetic islands near the radial position where the HCS is located. In these bifurcations, the number of island chains changes from the one determined by the RHW resonance to the other determined by the HCS resonance.

Isochronous bifurcations were also observed when the current in the HCS was increased and displaced 180° along the angular direction, resulting in pairs of small magnetic islands appearing inside the main magnetic islands where the HCS is located. These bifurcations also appeared when the current on the RHWs was reversed; in this case, sets of isochronous islands were observed around the surface where the HCS is located.

The isochronous bifurcations presented in this work are similar to the ones observed in visco-resistive MHD simulations of NSTX-U plasmas and observed in DIII-D plasmas. In further research, our results may yield additional insight into how plasma transport is affected by isochronous bifurcations in these and other machines.

The authors thank the São Paulo Research Foundation (FAPESP, Brazil) under Grant Nos. 2022/08699-2 and 2018/03211-6 for financial support. The authors acknowledge Dr. Marisa Roberto for her invaluable contributions to this work.

The authors have no conflicts to disclose.

Andre C. Fraile Jr.: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Methodology (equal); Software (lead); Writing – original draft (lead); Writing – review & editing (equal). Marisa Roberto: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Writing – original draft (supporting); Writing – review & editing (equal). Gustavo P. Canal: Formal analysis (equal); Writing – review & editing (equal). Ibere L. Caldas: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In this section, we present the numerical procedure used to calculate the Poincaré maps shown in Sec. V. The displacement along a magnetic field line was numerically calculated by solving Eq. (25) with the Dormand and Prince method³⁰ with order 5(4) and an adaptive step. In order to calculate the Poincaré map for the plasma cross section z = 0, we choose a set of N_IC initial conditions (r, θ, z) = (r_i,θ_i,0), with i = 1, 2, …, N_IC. Then, Eq. (25) is solved N_IC times from z = 0 to z = 2πR₀N_t, where N_t is the number of turns along the periodical length 2πR₀. Each time the magnetic field line crosses the z = 0 section, the corresponding position (r, θ) is computed, and the Poincaré map is plotted with the set of N_t⋅N_IC points.

As an example of this calculation, Fig. 12 shows the Poincaré map calculated with N_IC = 22 and N_t = 2000 considering the superposition of two perturbations: a pair of RHWs with (m_h,n_h) = (2,1) for I_h/I_p = 0.2% and an HCS with (m_s,n_s) = (8,2) for I_s/I_p = 0.058%.

As shown in Fig. 12, magnetic islands are formed around radial positions with rational safety factor values. After a Poincaré map such as that in Fig. 12 was obtained, we identified isochronous bifurcations (as those presented in Sec. V) by increasing the parameter I_s/I_p and observing the distribution of magnetic islands around rational surfaces while the parameter I_h/I_p is kept constant.

Figure 13 shows this stochastic parameter estimated for I_h/I_p = 0.4% and $ϕ p=0$⁠.

Since η < 1 in the range I_s/I_p ≤ 0.6%, neighboring separatrices do not intercept for the magnetic islands chains associated with q(r₄) = 4 and q(r_9/2) = 9/2.