Helical magnetic self-organization of plasmas in toroidal pinches with transport barriers formation

Nonlinear MHD modeling of toroidal pinch configurations for hot plasma magnetic confinement describes several features of the helical self-organization process, which is observed in both reversed-field pinches and tokamaks. It can also give a hint on why transport barriers are formed, by far one of the more interesting observations in experiments. The work tackles these two topics, helical self-organization and transport barriers formation - adding further information and examples to the results already presented in [Veranda, et al, Nucl.Fus. 60 016007 (2020)]. Regarding the topic of helical self-organization, a synthesis of the results obtained by a 3D nonlinear viscoresistive magnetohydrodynamics model will be presented. Modelling predicts a technique to “channel” reversed-field pinches into a chosen macroscopic helical shape and also predicts that the features of such helical self-organization, studied in the RFX-mod experiment in Padova, depend on two parameters only: plasma dissipation coefficients and edge radial magnetic field. They can be exploited to calm the natural tendency of reversed-field pinches to a “sawtoothing” dynamics, i.e. by decreasing visco-resistive dissipation and using helical edge fields not resonating with the plasma safety factor. Regarding the MHD description of the process of formation of transport barriers by magnetic chaos healing, we will describe the computation of Lagrangian structures, hidden in the weakly stochastic behaviour of magnetic field lines, acting as barriers to the transport. The radial position of such structures is observed to correspond to higher gradients of magnetic field lines connection length to the edge: this provides a further indication of their possible role in the formation of electron temperature barriers.


Introduction
This paper discusses self-organized states in currentcarrying pinches by focusing on the numerical solution of a 3D nonlinear visco-resistive magnetohydrodynamic (MHD) model, which considers the plasma as a quasineutral electrically conducting fluid. The paper is focused on the reversed-field pinch (RFP) configuration [1]. Numerical solution of such visco-resistive MHD model applied to the RFP predicted that the magnetic field can acquire a helical state in a self-organized manner [2][3][4][5][6], a feature later observed in high-current plasmas in various RFPs devices [7][8][9]. More recently, the MHD model has improved its descriptive capability of the typical experimental intermittency [10], predicted new ways for channeling helical self-organization towards a selected state highlighting the role of small edge helical magnetic field [11], a set of results then validated by experiments [12] in RFX-mod in Italy. This work focuses on the two key simulation parameters that are most involved in the emergence and control of 3D quasi-helical states in RFPs.A first parameter is the plasma visco-resistive dissipation, and defines a transition * e-mail: marco.veranda@igi.cnr.it between stationary helical equilibria and a cyclical dynamics. The other parameter is the amplitude of the helical magnetic boundary conditions, which control the level of intermittency and can force the self-organization towards a state chosen by an external feedback control system [13], with positive effects on plasma performances. Then, this work also discusses magnetic topology and confinement of magnetic field lines in RFP quasi-helical states.
The Lagrangian Coherent Structures (LCS) technique is used, borrowed from the study of dynamical systems (for a definition of LCSs see the review papers [14][15][16][17] and the applications in [18,19]). LCSs allow the discrimination of regions with different transport qualities: they are found in resilient bundles sorrounding the helical core during the whole dynamical evolution of the configuration, reinforcing the magnetic chaos healing effect provided by the formation of perfectly conserved helical magnetic surfaces in the core and safeguarding it, especially during relaxation events, whose periodicity is also analyzed in this work. Moreover we show a clear relationship between large gradients in connection length of magnetic field lines to the edge and the presence of LCSs -an important clue to their role in the formation of internal transport barriers.
As a final comment of this introduction we remind that, even though the paper is centered on the reversed-field pinch (RFP) configuration, similar modelling could be applied also to other current-carrying pinches like the tokamaks [20][21][22][23], the spheromaks [24], or the fieldreversed configurations [25] and that, on a much larger spatial scale, helical self-organization is an ubiquitous feature in astrophysical plasmas [26,27].
This work is structured as follows: in Sec. 2 we describe the employed MHD model. In Sec. 3 we focus on the pervasive presence of helical states in the solution of the nonlinear magnetohydrodynamics model. In particular, in Sec. 3.1 and 3.2 we describe the dynamical transition to quasi-helical states in the RFP, influenced by viscoresistive dissipation and by helical boundary conditions for the magnetic field. In Sec. 3.3 we describe the features of the periodicity of the relaxation events which tend to affect QSH regimes in the RFP. In Sec. 4, we present the studies about magnetic field line transport in the RFP. There we show how a resilient bundle of Lagrangian Coherent Structures (LCS) can give indications of local reduced transport of magnetic field lines, together with an enhanced magnetic chaos healing effect associated with the formation of quasi-helical states, in particular the ones based on non-resonant MHD modes. In Sec. 5, we draw some final remarks. We remark that the results described in this work are strongly related with the ones presented in [28], adding new material that reinforces the results obtained therein regarding the two topics, helical states dynamics and magnetic chaos healing, described in the following sections.

Numerical tools
The simulations presented here are performed with the SpeCyl code [29], recently verified numerically against the PIXIE3D code [30] with excellent results [31]. SpeCyl deals with the simple visco-resistive approximation which combines Faraday law, momentum equation (viscous Navier-Stokes) and the single fluid resistive Ohm's law (in Eqs. 1a,1b,1c) to evolve the magnetic B and velocity u fields in time t. Plasma density is assumed to remain uniform and constant, ρ = 1. Pressure is neglected. This yields a model with two parameters: the dimensionless resistivity η and viscosity ν. The model has been extensively used for compressible laboratory plasmas, capturing major physical effects observed in experiments (see for example [7,[32][33][34][35][36][37]). The equations written in dimensionless units are: Numerical simulations are performed in cylindrical geometry with aspect ratio R 0 /a = 4, where a = 1 is the cylinder radius. Simulations reported in this paper employ resistivity increasing towards the edge, η(r/a) = η 0 (1 + 20(r/a) 10 ), and uniform viscosity ν 0 , with η 0 ∈ [10 −7 , 10 −4 ] and ν 0 ∈ [10 −5 , 10 −2 ]). The relevant time scales of the system are the Alfvén time τ A = a(µ 0 ρ) 1/2 B −1 0 with B 0 on axis field, the resistive time τ R = µ 0 a 2 σ 0 , with σ 0 on-axis electrical conductivity and the viscous time scale τ ν = a 2 /ν. We can then define the Lundquist number S = τ R τ A , the viscous Lundquist number M = τ ν τ A , the Prandtl number P = S M and the Hartmann number H = S P −1/2 (the last two quantitites can be obtained by properly rescaling Eqs.1, as done in [6]). Magnetic boundary conditions are defined as follows. The edge radial magnetic field is either zero, corresponding to an ideal conducting wall, or helically modulated through MP [38]. The MP helical twist is defined by n MP with poloidal periodicity m = 1, and MP intensity is measured by the quantity MP%= b r (a)/B θ (a)%. The simulations start from an axisymmetric, non-reversed, unstable Ohmic equilibrium [39] with pinch parameter Θ = B θ (a) B z = 1.6. Fig. 1 presents the main features of ohmic axisymmetric equilibria for RFP. An uniform induction electric field E = E 0ẑ is imposed to sustain the plasma current, with E 0 /η 0 = 4.2. A slight perturbation nonlinearly induces the reversal and starts the MHD dynamics, typically driven by modes of resistive-kink/tearing nature [38].

Quasi-Helical solutions of the visco-resistive MHD model
In this section, we describe the results of an extended set of simulations of the visco-resistive MHD model described in Sec. 2. The resulting global picture consolidates the importance of the Hartmann number and the key role played by edge magnetic perturbations in ruling the helical selforganization process in reversed-field pinches. This section is structured as follows: at the beginning we describe the dynamics of couples of RFP simulations, choosing extreme values of the Hartmann number and MP amplitude (Sec. 3.1). In this way the possible qualitative behaviour obtained when spanning the simulations parameter space are found. Then, quantitative results will be shown for a larger number of simulations, analyzing the energy associated with the dominant helical mode versus the energy of the other MHD modes, often used also in experimental analyses (Sec. 3.2). The results described in the following can be synthetized as follows. Several features of the helical regimes are ruled by two quantities only, the Hartmann number and magnetic boundary conditions. An increase of the Hartmann number (i.e. a reduction of visco-resistive dissipation) at fixed MP's amplitude allows a transition from a single-mode behaviour to a multiple helicity one, while a proper choice of MPs at fixed high Hartmann number can pace the sawtoothing dynamics and can stimulate quasi-helical states in between sawteeth, with its same helical twist. Figure 1. Radial profiles of the equilibrium magnetic, current density and velocity fields in the RFP cases. The B z component of the RFP magnetic field slightly reverses at the edge, giving the name to the configuration. Magnetic field is normalized to its on-axis initial value, velocity field is normalized to the corresponding Alfvèn velocity. Dotted lines represent the same quantities after the relaxation event following a slight perturbation of the initial equilibrium, resulting in the beginning of the 3D MHD dynamics. In this picture the modulus of all the quantities is plotted. In the bottom-right corner, corresponding to high values of H (H > H c ) and using MPs, we find a state with a high level of helical symmetry (with the same twist of the MPs, which in this case is m=1, n=-6, MP% = 2%).

RFP helical regimes
We start by considering a magnetic perturbation applied on the m=1, n=-6 mode, which corresponds to the helical twist of a mode non-resonating with the safety factor. Typically in RFX-mod operation at high current the m=1, n=-7 mode is observed as the spontaneously dominant one in quasi-helical states, even though it was shown in [12] that n=-6 states can be induced by the feedback control system [13,40]. We first show in Fig. 2 the four different regimes which can be obtained making extreme choices for the quantities H, MP%.
In Fig. 2 we plot the temporal evolution of the total energy associated to the most important m=1 helical modes. Consider the helical states found at low Hartmann number and zero MP (top left): at t ∼ 0.1τ R the stationary helical solution of the model is made up of a single MHD mode with m=1, n=-11 (SH Single Helicity states, with welldefined magnetic surfaces in the whole plasma volume). This mode is different from the most unstable ones, which at the very beginning of the simulation around t ∼ 0.02τ R create a sequence of modes with increasing periodicity number starting from n=-8, sequence that is typical of the RFP multiple helicity relaxations. However, at a certain time, the (1,-11) mode overcomes the others, which exponentially decay to vanishing amplitude. At high values of H and zero MP (top right) instead, a strong competition between MHD modes ends in the presence of a sawtoothing regime (MH Multiple Helicity states, with stochastic magnetic field lines). In particular, one notices that the energy associated to the modes is two orders of magnitude lower in the high H cases, where it reaches values compatible with experimental measurements in the RFX-mod device and that the m=1, n=-6 mode has negligible amplitude.
In the second row of Fig. 2, we turn on the helical boundary condition on the m=1, n=-6 mode with amplitude MP%= 2%. In the low H and MP-on quarter (bottom left) we observe that the system reacts to the new boundary condition by amplifying to large amplitude the one stimulated by MPs in the spectrum of helical modes. In the end the m=1, n=-11 mode remains the most energetic one, while the m=1, n=-6 mode remains stationary at a finite energy. Finally, in the high H and MP-on quarter (bottom right) we find a systematic repetition of quasi-helical states, with a clearly dominating MHD mode which has the same twist of the applied MPs. Thus, in the end, the applied MP is successful in overcoming the other modes thanks to the natural depression at higher H values of the competing modes. This last example is the most adherent to the experimental observations in RFX-mod, as was discussed in [10,12].

Statistical analysis of numerical simulations
We compute the time average of the magnetic energy associated to the various MHD modes involved in the dynamics, so that the considerations made in the previous sections can be generalized by considering a set of around 100 simulations with varying parameters. We will show the results of a subset of simulations with n MP = −6, and keep in mind that if we consider helical regimes with helicities close to this one (n ∈ [−5, −9]) the results do not change.
Thinking in the Fourier space makes clear that the presence of a QSH or of a SH state can be diagnosed when we detect a large positive difference between the magnetic energy of a dominant m=1 mode and the sum of the energy of other m=1 modes (this indicates the presence of a helical structure, with an extreme case corresponding to the absence of m=1 secondary modes (SH)). Another important indication of the presence of a QSH or of a SH state is a low magnetic energy of the m=0 MHD modes: in fact nonlinear interaction between m=1 modes in the  Fig. 2). Such states were deeply studied in the past, for example in Refs. [3,6]. In particular it was shown that the dominant helicity in highly dissipative regimes without MPs depends on the value of Θ [41]. The scaling of the secondary modes' amplitude in the first panel of Fig. 3 changes abruptly around H c ∼ 3 · 10 3 : this is the signature of a dynamical transition. In the region H < H c the scaling is W m=1 M,sec ∝ H 3.1 , calling for a reduction of the Hartmann number to access helical states at very high plasma dissipation (meaning low plasma current and/or high plasma density). In the region H > H c the scaling is W m=1 M,sec ∝ H −0.6 , calling for an increase of Hartmann number to reduce fluctuations (meaning high plasma current and/or low plasma density).  Fig. 2). For values of n MP ∈ [−9, −5] there are always two regions where long-lasting QSH regimes, i.e. obeying to the two requirements written before, can be found. Pure and robust and stationary helical regimes are present at H H c and are studied in [3,6,41], where the features of the intermediate turbulent regimes found at H ∼ H c are also described. The regimes detected at H H c display the typical intermittent behaviour observed in RFP experiments. They were studied in [10,11]. Though the low-H regimes would represent a stationary, stochasticity-free option for the operation of an RFP device, comparison with experimental results tends to rule out the high-dissipation solution. In fact a reasonable lower bound for the Hartmann number value in typical RFX-mod experiments is H > 10 5 (see Ref. [42] for an estimate of H in RFX-mod). In summary, the extended numerical study presented in this subsection provides a complete picture of the different helical regimes of the RFP and of the transition between them. The major conclusion is that there is just one region where long-lasting QSH state qualitatively similar to experimental observations can be found, i.e. only at high Hartmann number and only using MPs with the proper amplitude.

Scaling of the "sawtoothing" characteristic time
We name the characteristic time between a sawtoothing event in the RFP τ saw , and we compute it by looking for the time intervals between maxima of selected physical quantities, like the magnetic energy associated to the dominant mode. We define τ saw as the average of the time intervals in a single simulation (see the bottom right panel of Fig. 2 for a visual definition of τ saw ) and we compute the related error δ τ saw .  P): noting that the equation can be properly rescaled as described in [6], the previous relation can be analyzed, provided that the value of τ saw is rescaled as well (using τ saw →τ saw = P − 1 2 τ saw ). We can thus obtain a new scaling which is linked to the previous one by a simple relation obtained by substituting the definitions of (H, P) in terms of (S , M).This choice of parameters yields a scaling of the rescaled sawtoothing period asτ saw ∝ H 0.76±0.03 P −0.01±0.04 which depends mostly on H, see Fig. 4. The relatively strong dependence of τ saw on dissipation, already found in [29], may explain why RFP quasi-helical states are observed to become more persistent (i.e. with less sawtoothing events) when increasing current (i.e. decreasing resistivity): in Fig.5a of [43] a strong dependence of QSH-persistence is in fact observed with the Lundquist number, which in turn strongly depends on plasma current.

Topology of quasi-helical states: barriers to the transport of magnetic field lines
In this section we discuss Poincaré plots, Lagrangian Coherent Structures (LCS) computation and connection length L c to show that a suitable choice of the MP spectrum at the edge is an important ingredient for the attainment of higher magnetic order and that bundles of LCSs sorround the helical core of the configuration. The studies presented in this work are performed with the NEMATO field line tracing code (presented in [44] and numerically benchmarked with success in [45]).

Lagrangian Coherent Structures and the degree of magnetic order
Few preliminary words about the LCS tool, borrowed from the study of dynamical systems. LCSs can be used to distinguish regions with different transport qualities and their most important feature is that, for a meaningful finite time span which characterizes a LCS [18], magnetic field lines belonging to different regions cannot mix with each other. The LCS technique allows a more refined analysis than a simple inspection of a Poincaré map. As shown in [12,46], when a Poincaré map suggests only a stochastic behaviour, the LCS tool highlights some structures, which we color in blue in the panels a)-f) of Fig. 6.
If we look at Fig. 5 we can define a high magnetic order if there is a large area of the Poincaré plot occupied by conserved magnetic surfaces (colored areas) or, looking at Fig. 6, if we compute the presence of well-defined bundles of LCS (plotted in blue), or if we measure a large connection length of magnetic field lines to the edge (white areas of high connection length maps on the poloidal sections, plotted in Fig. 6) and if such positive features resist to the sawtoothing events.

Poincaré plots of quasi-helical states
The two simulation cases that we consider are characterized by the visco-resistive parameters S = shows that a core region of conserved magnetic surfaces is present when the QSH state is forming (red region in the first row, blue in the second). Surrounding the core is a stochastic region, whose radial extension depends on the relative intensity of the dominant helical mode and its perturbations (grey region). As a first consideration we notice that the area of conserved magnetic surfaces in the nonresonant case is higher than that in the resonant case, both during the formation of the helical state and its relaxation to a MH state (last column of Fig. 5): this was already found in [12], but here the greater resilience to magnetic chaos is confirmed during the whole cycle of evolution of a QSH state. Another difference between the first two rows is that in resonant cases the core-conserved region does not enclose the cylindrical axis at r/a = 0, while the non-resonant cases generally do (apart from the relaxation phase). This may be important for understanding the features of the electron temperature profiles measured in RFP experiments, which display clear internal transport barriers previously put in relation to magnetic chaos healing and helical q profile [7,47].

Lagrangian Coherent Structures in quasi-helical states
LCSs are computed in the two cycles of QSH states, and are colored in blue in Fig. 6. We considered just the LCS located in the core, i.e. at r/a < 0.7. LCS are computed using a finite time L z of 4 toroidal turns for each magnetic field line (for a justification of the finite time choice, a critical parameter of the whole method for computing LCSs, see [12,46]). From a topological point of view there are no relevant differences between the LCSs' shape in the two quasi-helical states: they tend to encompass the helical core as it widens during the evolution of the helical state. A certain degree of difference can be observed in the presence of more developed bundles of LCSs in the non-resonant case, panels a)-c) than in the resonant case. Interestingly, in the latter case the position of a bundle of LCSs can be correlated with the position of reversed shear of the helical safety factor profile (see Fig. 5 and Ref. [47]).
LCSs are observed to evolve on a time scale slower than the dynamical one, as their shape does not change significantly between the beginning and the middle of QSH states in the first two columns of Fig. 6), which are separated by τ dyn ∼ 10 4 τ A . Another important time scale is the one related to the LCSs' power of separating different spatial regions, which lies on a shorter timescale τ sep ∼ 10 3 τ A [12], after which LCSs become leaky.
An assumption at the basis of our use of the LCS technique is that the magnetic field line are used as a proxy for the particle trajectories along the lines. It is our plan to use tools beyond this ansatz to take into account the behaviour of charged particle in the magnetic field under scrutiny, like the one described in Ref. [48].

Study of the connection length to the edge
LCSs' ability to confine magnetic field lines, a feature not to be argued from a simple inspection of the Poincaré plots of Fig. 5, can be analyzed by computing the connection length to the edge, i.e. the length L c that a magnetic field line starting at (r 0 , θ 0 ) covers to reach an edge radius taken as a reference for the computation r/a = 0.75. The results are shown in Fig. 6, where magnetic field lines have been integrated until a maximum normalized length of L c,max = 10 5 (corresponding to a time of τ conn ∼ 10 5 τ A , i.e. much greater than the duration of a QSH cycle). We observe first that the areas in white, characterized by L c = L c,max , are generally larger in non-resonant states. We also notice regions where the connection length drops sharply: these regions are in good correspondence with the presence of bundles of LCS, another evidence of their role in defining special regions for the transport of magnetic field lines. We also note that when the ordering role of the helical mode fades, in the third column of Fig. 6, a rich set of topological structures appear connecting the plasma core to the plasma edge with relatively low connection length. Consider for example the green spots at r ∼ 0.2, θ ∼ 0 in panel c), which are named "escape channels" in plasma literature [49,50]: there, they are shown  to possess a rich fractal structure which may explain some anomalous values of transport coefficients in the plasma edge. As a final point it is important to notice that a "barrier" structure persists even during the reconnection event (see panels c) and f)). The detailed impact of these properties on temperature profiles are under study, both from the numerical and the experimental point of view: the use of a numerical code for the solution of the anisotropic heat transport equation [51,52] will help determining whether LCSs are responsible for the high temperature gradients observed in correspondence of Internal Transport Barriers [53].

Summary and final remarks
Helical self-organization is an ubiquitous feature in current-carrying toroidal plasmas for magnetic confinement of fusion plasmas. In this work we focus on the reversed-field pinch configuration, where a global selforganization process systematically occurs at high plasma current [1], impressing a helical shape to the whole plasma column [54,55], with beneficial effects on confinement [7]. The first result of this paper is that there are two quantities ruling the nature of helical states and the possibility to interact constructively with self-organization, namely plasma dissipation (Hartmann number H, the inverse geometric mean of resistivity and viscosity) and helical magnetic boundary conditions. Synthetizing, there are two combinations of Hartmann number and helical boundary conditions leading to quasi-helical solutions for the RFP: the pulsating solutions at high H and the stationary helical solutions at low H. We have shown that it is possible to interact constructively with a helically self-organized RFP plasma by lowering dissipation (i.e. by increasing plasma current in the experiment) and favoring a selected MHD mode by modulating appropriately the edge magnetic field (i.e. by acting on the plasma-edge with advanced feedback control system). We also showed that the sawtoothing period in the simulations described in this paper is shown to depend on the Hartmann number only, i.e. on the simple product of resistivity and viscosity, a result to be validated experimentally.
The second result of this paper is given by the identification of bundles of Lagrangian Coherent Structures as responsible for the increased confinement of magnetic field lines close to the core of the helical plasma. LCSs surround the core of helical states and resist to the sawtoothing events in RFPs, even when magnetic topology gets less ordered. They are correlated to sharp gradients in the connection length of magnetic field lines to the edge, a fact that further supports the LCSs' role in the formation of internal transport barriers. In the future, an analysis of heat transport taking into account the strong anisotropy in the thermal conductivity of fusion plasmas and the presence of magnetic islands and stochastic magnetic fields is foreseen (using the tools described in Ref. [51,52]), aiming at clarifying the role of LCS in sustaining transport barriers evolution -following previous studies [56,57] and in the addition to already established mechanisms of microturbulence suppression caused by sheared flows [58,59]. Future numerical work will be dedicated to analyze the role of a finite pressure gradient and to the effect of a realistic boundary on the helical self-organization of the RFP, adding a thin resistive shell and a vacuum layer between the plasma and the ideal shell where MPs are applied.