Density limit studies in the tokamak and the reversed-field pinch

In the RFX RFP [14], data points in the plane (n_G, n₀), with n₀ the central electron density, seem to follow a Hugill–Greenwald scaling, with scarce data points for n₀ > n_G, as published several times in the past [15–21], but a real disruptive limit when approaching n_G has never been found [18, 20]. This can be highlighted by grouping data in classes, according to the value of the toroidal loop voltage, as shown in figure 1(a). The straight lines in the plot correspond almost exactly to curves at constant V_loop: moreover, it is evident that the larger the V_loop, the steeper the slope of n₀ versus n_G. This means that, with enough loop voltage, in the RFP it is possible to get n₀ = n_G with no disruption. This can be done at low current, since the ohmic input power is low despite the high voltage, as shown in figure 1(a) as red points: the Greenwald limit is exceeded by a factor ∼2. Moreover, a linear regression of V_loop as a function of n₀/n_G in the RFX database, represented by the grey points of figure 1(a), gives as a result V ≈ 12 + 60 n₀/n_G. The same linear law was seen to hold between the effective charge Z_eff and V_loop in the old RFX, with a bias loop voltage of V_dyn ∼ 17 V, which was interpreted as a minimum value to sustain the RFP dynamo [22]. This means that the limit has a radiative nature, with a difficulty in performing high current and high density discharges, that require an exceedingly large input power. Nevertheless, no clear disruption is seen, so that it is impossible in the RFP to define a 'critical' density, e.g. on the basis of a soft x-ray crash and associated loss of the plasma thermal content. This is a clear advantage of the RFP configuration over the tokamak, being summarized in the empirical, linear law

$$\begin{equation} V \approx 12 + 60 \, \frac{n_0}{n_{\rm G}} = V_{\rm dyn} + V_{\rm G} \frac{n_0}{n_{\rm G}} , \end{equation} \tag{ 1 }$$

where V_dyn = 12 V is the bias dynamo voltage, and V_G = 60 V is the voltage required to reach the Greenwald limit. Since the RFP heating is ohmic, multiplying equation (1) by n_G, we get (V − V_dyn) I_p = V_Gπa² n₀, which is a linear dependence of the central density on heating power (input power minus the dynamo term). This dependence has indeed been reported for RFX [23]: with larger input power, larger densities can be accessed. This is a result reminiscent of the Sudo scaling, and it has been shown in the past also in the TEXTOR tokamak [24]. In the same paper [23], it was reported that, with lithization, the power required to get to the same density is half as much as in standard RFX discharges. According to equation (1), this means that also V_G can be greatly reduced with wall conditioning. Another means of increasing central density without impacting too much on the loop voltage is pellet injection, as demonstrated in the Madison Symmetric Torus (MST) RFP where n_G can be easily exceeded with this technique [25]. In MST too, no disruption is seen crossing the Greenwald limit.

Zoom In Zoom Out Reset image size

A clear threshold with a critical density can be defined instead on a MHD basis, as done by Granetz on Alcator [3]. In fact, by increasing density, in the RFP we destabilize the m/n = 0/1 mode, which is responsible for the development of a MARFE, as it will be shown in section 3. In this sense, the RFP density limit is more a non-disruptive, MHD limit: in fact, it has already been published that the m = 0, n < 7 modes increase sharply as a function of n₀/n_G (at constant current) [21], and that this is a pre-requisite for the development of the MARFE [26]. In figure 1(b) we report results from recent experiments (April 2014), where we tried to stabilize the 0/1 mode with the feedback controlled, 192 active coil system of RFX [27]: indeed, the normalized, perturbed radial field at r = a, $B^r_{0,1}/B_{\theta,a}$ remains low, until a threshold n₀ ∼ 0.35 n_G is reached, then it shows an explosive behaviour. The threshold follows a Greenwald-like scaling, which recently has been put into relationship with the Prandtl number [28] which in visco-resistive MHD simulations governs the stability of the m = 0 modes [29]. This result is confirmed in a range of plasma currents I_p = 0.6 ÷ 1 MA with shallow reversal, q_a = −0.015 (corresponding to a reversal parameter F = B_{φ, a}/〈B_φ〉 = −0.08).

Let us consider now the FTU tokamak, which is a circular tokamak running in L-mode [30]. The edge density (r/a = 0.8) at the disruption is shown in figure 2(b): it follows a reduced, Greenwald scaling of the form

$$\begin{equation} n_{\rm edge} \sim 0.35 \, n_{\rm G} . \end{equation} \tag{ 2 }$$

It is striking that the threshold of equation (2) coincides with the RFP MHD limit, taking into account that density profiles in the RFP are usually flat at intermediate–low density, with n₀ ∼ n_edge: this is shown in figure 2(b) as light blue stars. Contrary to the RFP, in FTU the maximum achievable central density essentially depends on the toroidal magnetic field only [31], with a Granetz-like scaling n₀ = 0.19 B^1.5 [units 10²⁰ m⁻³ and T, see figure 2(a)]. This behaviour can be explained in terms of density peaking, since the edge density (open symbols in figure 2(a)) does not depend on B, as shown in equation (2). In fact, in FTU with a stronger B field (and the same current level) disruption happens later in the discharge, and larger core densities can be accessed, i.e. with larger peaking. This is shown in figure 6 in [32]. The behaviour of core and edge density at the disruption can be summarized as a scaling n₀/n_edge ∝ q_a, shown in figure 8 of [32]. A similar scaling |∇n|/n∝∇q/q has been put into relationship with the theory of the 'curvature pinch' [33] in a smaller database of auxiliary heated discharges [34] (not considered in this paper), but no clear theoretical answer to the peaking phenomenon has yet been given in FTU. Combining the empirical law n₀/n_edge ∝ q_a with the Greenwald edge scaling, it is obtained n₀ ∼ 0.7 B/μ₀R. Curiously, this is the Murakami parameter [35] which was used in the early research on the density limit. In other words, in experiment, the edge density is set by the plasma current: n_edge ∼ 0.35 n_G = 0.35 I_p/π a². For a given n_edge, the maximum n₀ prior to disruption is determined by the profile peaking, namely, by the value of q_a. For a given plasma current, this means that the maximum central density n₀ is determined by B. Conversely, the smaller B, the smaller n₀ and the flatter the density profile: in these cases, n₀ follows also a Greenwald scaling. This is the case of the RFX device (green points in figure 2), which can be operated also as a low-B tokamak with a flexible feedback control system [36]. In the RFX tokamak, n₀ follows the scaling with the magnetic field, but also a Greenwald scaling, since density profiles are essentially flat: as a consequence, full and open symbols in figure 2(a) collapse one over the other.

Zoom In Zoom Out Reset image size

Historically, the dependence on $|\vec{B}|$ was lost in the original 1988 paper by Greenwald [1]: in the Hugill plane [2], which displays the inverse safety factor versus the Murakami parameter [35], a straight line corresponds to a critical density n_c such that

$$\begin{equation} \frac{1}{q_a} \propto n_{\rm c} \frac{\mu_0 R}{2 B} \rightarrow n_{\rm c} = \frac{I_p}{\pi a^2} , \end{equation} \tag{ 3 }$$

which is the definition of n_G. In this way the dependence on $|\vec{B}|$ is completely lost. On the contrary, from the point of view of basic plasma physics, recently B was found to be the single parameter limiting the maximum achievable core density in laboratory and astrophysical plasmas [37]. The results from FTU [31, 32] confirm and strongly re-propose the B-limit for the central density, plus a reduced-Greenwald limit (equation (2)) for the edge density. As a matter of fact, the only machines operating with a central density n₀ ≈ n_G are those possessing the strongest fields, namely, FTU and ALCATOR. This could be a concern for DEMO and/or ITER, since with a lower range of B a lower density peaking can be obtained, although transient density peaking can be induced via pellet injection and/or NBI [25, 38, 39]. Therefore, density in ITER and DEMO could be essentially limited by the edge value n₀ ∼ 0.35 n_G, which we will prove to have a more fundamental MHD origin (see sections 3, 4).

An important point to raise is the role of the thermal instabilities in setting the environment for the development of the density limit. In both RFX and FTU, the density limit is associated with the appearance of the multifaceted asymmetric radiation from the edge (MARFE) [40, 41]. In FTU, the MARFE is an annulus of radiative-unstable plasma, which appears at n₀ > 0.4 n_G [31] as a toroidal ring, poloidally localized, as shown in figure 3(a). The threshold 0.4 n_G is consistent with the Lispschultz scaling of the MARFE critical density with I_p/π a² [40]. This toroidal MARFE is caused by a reduction of parallel electron conductivity, which prevents the energy from redistributing along the flux surface, and lets the formation of local cold regions with strong line emission, in the high field side, where the heat flux from the core is lower [42]. In TEXTOR, Tokar' provided a nice, simple model of this phenomenon, interpreted as a local imbalance between parallel heat diffusion and ionization/radial particle transport [43]. In that model, the critical density scales as

$$\begin{equation} n_0^{(\rm cr)} = \frac{1}{q_a} \, \sqrt{\frac{\kappa_{\parallel} a}{60 \pi R^2 D_{\perp} \, \sigma_{\rm ion}}} , \end{equation} \tag{ 4 }$$

κ_∥ being the parallel heat diffusivity, D_⊥ the cross-field particle diffusivity, and σ_ion the cross-section for the ionization process. This scaling is consistent with the scaling found in FTU.

Zoom In Zoom Out Reset image size

Contemporary with the MARFE, in FTU a 2/1 mode starts in the high density regime during the density ramp-up and grows up to high amplitude in the final phase of discharge preceding the density limit disruption [44]. The mode onset is in agreement with linear stability calculations, confirming that its appearance is correlated with a peaking of the current profile associated to the cooling of the edge plasma, as we will discuss in detail in section 4.

In the RFP, the MARFE appears at n₀/n_G > 0.5 as a poloidal ring of high radiation, toroidally localized, as shown in figure 3(b). Radiation comes from He-like impurity lines (mainly O and C) excited by low temperatures in the annulus: there is no evidence of strong recombination, as deduced from the behaviour of the ratio of the H_γ/H_α lines along the toroidal angle [20].

The unusual symmetry of the RFP MARFE is easily explained, taking into account that the equilibrium field is toroidal in the FTU edge, poloidal in the RFP edge. Rather than being a simple problem of energy redistribution (as it is the case with the Tokar' model), the MARFE in RFX is instead linked with a well-defined $\vec{E} \times \vec{B}$ flow pattern. The trigger of the MARFE mechanism in the RFP is the 0/1 island, which resonates at q = 0 in the RFP edge [45] and it is destabilized at n₀/n_G ∼ 0.35, as already mentioned in section 2.1. This island grows with n/n_G, as shown in figure 4: frame (a) reproposes the explosive growth of the 0/1 mode at I_p = 800 kA, already discussed in figure 1(b); frames (b)–(d) show a Poincaré plot of the island, in three discharges with increasing n/n_G, corresponding to the coloured circles in figure 4(a). Islands are plotted using as input the eigenfunctions calculated by the code NCT in toroidal geometry [46]: the radial width of the island is largely determined by the value of $B^r_{0,1}$ at r = a [47]. The island grows to a critical size at n₀/n_G ≈ 0.5, when it touches the wall (figure 4(d)): this is a necessary condition for driving electron transport along the magnetic field in regions of very short connection length L_c, as already discussed in the past on RFX [26, 48] and TEXTOR [49–51] edge stochastic layers. We can estimate the critical amplitude of the island as $B^r_{0,1}/B_{\theta,a} \sim 0.1\%$, which corresponds to $B^r_{0,1} \approx 3$ gauss at I_p = 800 kA and q_a = −0.015. The critical size of the island depends on q_a: with lower (more negative) q_a = −0.03 (reversal parameter F = −0.2), the critical amplitude for wall wetting almost doubles, with $B^r_{0,1} \approx 6\,{\rm G}$ (normalized value 0.3%), as reported in a previous work [21], since the resonance is moved further inside the plasma. Despite this, the threshold in n/n_G remains unchanged, because m = 0 modes are destabilized at lower q_a [29], so that the initial size of the 0/1 island is also larger. In this sense, the threshold in n/n_G is the result of a competing process between island width and wall proximity, as pointed out in [48].

Zoom In Zoom Out Reset image size

In summary, this paper strongly supports the idea that edge islands are a fundamental aspect in the Greenwald phenomenon: it is curious to notice that a similar behaviour for the 2/1 island in the tokamak, follows from a 3D generalization of the Rutherford equation with an explosive radiative term [52]. In the RFP, experiments show that the threshold for 0/1 destabilization is slightly smaller with strong Neon puffing (n₀/n_G ∼ 0.3), but impurities seem not to be a fundamental ingredient in the island growth process.

The growth of the 0/1 island is associated with a pathological behaviour of the toroidal component of the plasma flow v_φ, which in the RFP corresponds to the $\vec{E} \times \vec{B}$ flow (v_φ ≈ E_r/B_θ). The flow v_φ is generally negative along the toroidal angle [21]. When the island grows to the critical size shown in figure 4(d), v_φ reverses direction, with the formation of two null points, source and stagnation, as shown in figure 5(b). A similar reversal of the perpendicular flow is seen also in the C-mod tokamak when n₀/n_G ≳ 0.35 [9, 10], and, more recently, also on ASDEX at a comparable normalized density [12]. A crucial point of the flow analysis in RFX was to map the toroidal angle into the helical angle of the perturbation [53]: in the 0/1 case it is trivially u_{m, n} = mθ − nφ + φ = −φ + φ(φ phase of the mode). In this way, one can clearly recognize the presence of a source and a stagnation point, with the latter corresponding to the toroidally localized MARFE, as shown in figure 5(b). The source corresponds instead to the maximum H_α signal, which stands in between the O-point (OP) and X-point(XP) of the main island, at u ≈ π (figure 5(a)). In fact, it is a well-known result in the RFP that the phase-locking of the m = 1 TMs at u = π (i.e. shifted 90° toroidally from the OP of the 0/1 island) is a preferential source of particles and localized PWI [54].

Zoom In Zoom Out Reset image size

On the contrary, the stagnation point corresponds to the XP of the 0/1 island, at u ∼ 3/2π: here the density radial profile becomes markedly hollow, with an edge-localized peak as large as ∼1.5 n_G, as shown in figure 5(c). The cold and overdense plasma in the region of the XP is ultimately responsible for the radiative behaviour and the MARFE, as shown in figure 3(b). The connection between MARFE and m = 0 islands is confirmed by recent RFX experiments, where by removing the q = 0 resonance, edge density peaking is avoided [55]. It is also consistent with experiments done in TEXTOR, where the critical density for the MARFE onset was seen to decrease with the amplitude of static, edge islands produced via resonant magnetic perturbations (RMPs): the larger the island, the sooner the MARFE would appear [56].

Our interpretation of the MARFE as due to an $\vec{E} \times \vec{B}$ convective cell is radically different from the traditional model of equation (4): in fact, while in the traditional approach the critical density is set by the ratio $\sqrt{\kappa_{\parallel}/D_{\perp}}$ (parallel heat transport to perpendicular diffusivity), in the RFP we have shown that a new, critical element is the ratio between the two perpendicular transport terms, radial diffusion and toroidal convection, D_r/v_φL_n (with L_n the density gradient characteristic length). The convective term v_φL_n in the RFP can be two orders of magnitude larger than the diffusive term, as it has been shown elsewhere [21]. In RFX low current discharges, it is possible to measure also the radial component of the flow, v_r, with insertable probes (even if in a different symmetry, the helical 1/7), thus confirming the presence of a convective cell, dragging plasma along the helical angle [57, 58]. $\vec{E} \times \vec{B}$ convective cells have been measured around edge-resonant islands in tokamaks and stellarators [59–62] and provide an interesting link between the density limit and experiments with RMPs. In this sense, the density limit can be seen as a particular instance of the more general problem of the plasma response to RMPs.

In the RFP, the convective cell has been modelled with the guiding-centre code ORBIT [63], showing that the 0/1 flow pattern corresponds to an ambipolar potential, balancing electron radial diffusion between the OP and XP of the island, which is embedded in the RFP stochastic edge [64]. In fact, the connection length L_c is larger, by more than 3 orders of magnitude, at the XP, with respect to the OP: in classical mechanics, orbits (=magnetic field lines) take infinite 'time' (i.e. infinite parallel length) to perform a complete excursion around the XP along the separatrix. These long excursions take the name of 'homoclinic tangles', and are commonly found e.g. in diverted tokamaks [65]. Electrons, having a smaller Larmor radius, are more sensitive to the resonant tangle present around the XP [21]: ions instead average out the toroidal asymmetry of L_c. This causes a small charge imbalance [66] (order 10¹³ ÷ 10¹⁴ m⁻³) along u, which is immediately counteracted by an electrostatic potential with the same symmetry as the 'parent' island. This is shown in figure 6(a), where the 0/1 island is over-plotted to the contour map of the ambipolar potential: by definition, the flow is such that $\vec{v} \cdot \nabla \Phi = 0$, so the saddle point in the contour corresponds to the measured, convective cell. ORBIT predicts a correct amplitude, phase and geometry of the potential, with the potential well (positive bulge of v_φ) staying in proximity of the XP of the island, as it is evident by comparing figure 6(a) with figure 5(b). Moreover, if one searches for an algebraic solution to the ambipolar constraint Γ_e = Γ_i as a function of the potential phase $\tilde{\phi}$, one can find a second root, in addition to the known solution at the XP. This second root has the potential well at the OP (see figure 6(b)). The presence of two roots, instead of a single one as previously found [64], is a recent result which has been obtained after modifying ORBIT guiding centre equations [63] to correctly express electron drifts. Interestingly enough, recent simulations of the 4/1 islands generated with a RMP in TEXTOR, show a similar behaviour of the connection length L_c [51, 64], and also the presence of two roots for the model of ambipolar potential, one with the potential well at the OP, the other at the XP [67]. The solution matching TEXTOR experimental data is stable in the thermodynamic sense usually followed for stellarators [68], that opens interesting possibilities also for the RFX case: it shows that in principle it is possible, by acting on the T_e/T_i ratio, to make the system flip from one root to the other, e.g. with additional heating of the electron/ion channel (ECRH or ICRH). This is an interesting extension of the theory of ambipolar roots in the stellarator, applied to edge stochastic layers. Indeed, compelling experimental evidence in FTU with ECRH targeted on the 2/1 island will be shown in section 4.

Zoom In Zoom Out Reset image size

As a final remark, it is worth recalling that the presence of an XP, namely, of a resonance generated e.g. via magnetic field reconnection, seems to be a necessary condition for the reversal of the perpendicular component of the flow: a simple edge ripple is not sufficient to reverse the flow, as observed in the case of a 1/7 mode in RFX-mod, which resonates well inside the plasma [48, 53, 58, 66]. In fact, a ripple of the toroidal flux ψ along φ can be canonically transformed in a constant helical flux χ along the helical angle u [47]: electrons, as a first approximation, follow the helical flux χ, and thus only a simple modulation of density and temperature along φ is seen.

The role of edge islands in the stellarator (a.k.a. 'low-order' rationals) has been studied in relationship with turbulence: it has been reported that configurations with ι = 1/q close to rational values experience more easily confinement transitions to the H-mode, accompanied by a strong reduction of broadband turbulence [69]. The connection between these islands and radiation-driven collapses has not been investigated in detail, although edge islands could play a role by increasing edge temperature [70].