RF current condensation in magnetic islands and associated hysteresis phenomena

TL;DR

This paper investigates the nonlinear RF current condensation effect in magnetic islands, revealing enhanced stabilization, hysteresis phenomena, and potential for broader deposition profiles in plasma control.

Contribution

It introduces a simplified nonlinear model demonstrating how RF current condensation can improve magnetic island stabilization and describes associated hysteresis effects.

Findings

Nonlinear effects cause larger temperature variations.

Deposition widths become narrower with nonlinear effects.

Enhanced robustness in island stabilization due to RF condensation.

Abstract

The nonlinear RF current condensation effect suggests that magnetic islands might be well controlled with broader deposition profiles than previously thought possible. To assess this possibility, a simplified energy deposition model in a symmetrised 1D slab geometry is constructed. By limiting the RF wave power that can be absorbed through damping, this model describes also the predicted hysteresis phenomena. Compared to the linear model, the nonlinear effects lead to larger temperature variations, narrower deposition widths, and more robust island stabilisation. Although, in certain regimes, the island centre can be disadvantageously shaded because of the nonlinear effects, in general, the RF condensation effect can take place, with current preferentially generated, advantageously, close to the island centre.

Figures

Tables

Table 1. Table 1: Summary of the various scales relevant to obtain the form of Eqs. ( 8a ) and ( 8b ) in typical tokamak parameters of T = 10 𝑇 10 T=10~{} keV, a = 1 𝑎 1 a=1~{} m, R = 5 𝑅 5 R=5~{} m, W i = 0.05 subscript 𝑊 𝑖 0.05 W_{i}=0.05~{} m, n = 10 14 𝑛 superscript 10 14 n=10^{14}~{} cm -3 , χ e ⟂ = 0.5 subscript 𝜒 perpendicular-to 𝑒 absent 0.5 \chi_{e\perp}=0.5~{} m 2 s -1 , Z = 1 𝑍 1 Z=1 , q = 2 𝑞 2 q=2 .

Time scale	s
${\tau }_{\mathrm{se}}$	$1\times {10}^{-3}$
${\tau }_{⟂}$	$1\times {10}^{-3}$
${\tau }_D$	$3\times {10}^{-3}$
${\tau }_E$	$8\times {10}^{-3}$
${\tau }_{\mathrm{eq}}$	$4\times {10}^{-1}$
${\tau }_i,{\tau }_{\eta }$	$2\times {10}^2$

Table 2. Table 2: Initial conditions that specify the problem that is to be solved.

Var	$\tilde{x}$ limit	Var limit
$u$	$-{\alpha }_0$	0
$u^{\prime }$	$0$	$0$
$V$	$-{\alpha }_0$	$V_X$
$V^{\prime }$	$-{\alpha }_0$	$-V_X$

Table 3. Table 3: Reference description of variables in Sections III and IV

$\tilde{x}$	Dimensionless position
$u$	Dimensionless temperature variations
$T_0$	Background temperature
$u\left(0\right)$	$u$ at island centre
$V$	Dimensionless wave energy density
$V_X$	$V$ at left X-point
$V_0$	Scaled $V_X$: $V_0={\alpha }_0{V}_X$
$V_f$	$V$ leaving island
${\alpha }_0$	Half island to deposition width
	(also island edges)
${\alpha }_0/W_i$	Deposition strength
$W_i$	Island width
$x_{\mathrm{cent}}$	Deposition mid-point location
$w$	Ratio phase to thermal speeds
$\omega$	RF frequency

Equations

$$\frac{3}{2}n{k}_B{\partial }_{t}T-\nabla \cdot (\kappa \nabla T)=P$$

$$\frac{3}{2}n{k}_B{\partial }_{t}T-{\partial }_x({\kappa }_{\perp }{\partial }_{x}T)=P$$

$$\frac{3}{2}n{k}_B{\partial }_t\tilde{T}-{\kappa }_{\perp }{\partial }_x^2\tilde{T}=P$$

$$-{\kappa }_{\perp }{\partial }_x^2\tilde{T}=P$$

$$d_t\bar{V}=[-(\nabla \cdot v_g)+\frac{{\omega }_t}{\omega }+2\gamma ]\bar{V}\approx 2\gamma \bar{V},$$

$$d_x\bar{V}=2\frac{\gamma }{v_g}\bar{V},$$

$$\gamma \propto e^{-(v_{\parallel }/v_{T{e}_0}{)}^2}\exp \left(\frac{v_{\parallel }^2\tilde{T}}{v_{Te,0}^2{T}_0}\right)\to \frac{2\gamma }{v_g}=-\alpha e^u.$$

$$d_x\bar{V}=-\alpha (x)e^u\bar{V}$$

$$-{\kappa }_{\perp }{\partial }_x^2\tilde{T}(x)=-v_g\frac{{\bar{V}}^{\prime }(x)+{\bar{V}}^{\prime }(-x)}{2}$$

$${\bar{V}}^{\prime }(x)=-\alpha (x)e^u\bar{V}(x)$$

$$V^{\prime }(\tilde{x})$$

$$u^{\prime \prime }$$

$$u(\tilde{x})=2\log \gamma -\log \left[\sqrt{(\lambda +1{)}^2-{\gamma }^2}\cosh \gamma \tilde{x}+(\lambda +1)\right]$$

$${\left[{\gamma }^2-(1+\lambda )\right]}^2={\cosh }^2({\alpha }_0\gamma )\left[(\lambda +1{)}^2-{\gamma }^2\right]$$

$${\gamma }^2=(2\lambda +1)+(V_X-\lambda {)}^2$$

$${\alpha }_0\sim 0.1-3V_0\sim 10^1$$

$$V^{\prime }(\tilde{x})$$

$$u^{\prime \prime }$$

$$S=\frac{V(\tilde{x})+V(-\tilde{x})}{2}$$

$$A=\frac{V(\tilde{x})-V(-\tilde{x})}{2}$$

$$u^{\prime \prime }={\left[\frac{V(\tilde{x})-V(-\tilde{x})}{2}\right]}^{\prime }=A^{\prime }$$

$$(A+S{)}^{\prime }=-e^u(A+S)$$

$$S^{\prime }=-e^{u}A$$

$$A^{\prime }=-e^{u}S$$

$$u^{\prime }=A$$

$$u^{\prime \prime }=e^u[e^u-(1+\lambda )]$$

$$u^{\prime 2}=e^u\left[e^u-2(\lambda +1)\right]+\mathbf{C}$$

$$z^{\prime 2}=1-2(\lambda +1)z+\mathbf{C}{z}^2$$

$$\mathbf{C}={\gamma }^2,B=\frac{\lambda +1}{{\gamma }^2},A^2=\frac{(\lambda +1{)}^2-{\gamma }^2}{4{\gamma }^4}$$

$${\left[{\gamma }^2-(1+\lambda )\right]}^2={\cosh }^2{\alpha }_0\gamma \left[(\lambda +1{)}^2-{\gamma }^2\right]$$

$$u(\tilde{x})=2\log \gamma -\log \left[\sqrt{(\lambda +1{)}^2-{\gamma }^2}\cosh \gamma \tilde{x}+(\lambda +1)\right]$$

$${\gamma }^2=(2\lambda +1)+(V_X-\lambda {)}^2$$

$$u^{\prime \prime }=\frac{V^{\prime }}{2}\to u^{\prime }=\frac{V}{2}-\frac{V_X}{2}$$

Full text

††preprint: AIP/123-QED

I Introduction:

Magnetic confinement approaches to fusion rely on the ordered topology of nested magnetic surfaces to prevent the plasma from escaping. Devices such as tokamaks and stellerators are designed towards this; however, in reality, magnetic fields are not perfect and are subject to error fields. These change the magnetic topology and result in magnetic islands appearing at rational surfaces.

Magnetic islands are characterised by flat density and temperature profiles due to enhanced transport through them. This reduction of pressure gradients suppresses bootstrap current within the island, generally making the island grow until saturationChang et al. (1995); LaHaye et al. (1997); Zohm et al. (1997); Gates et al. (1997). As a result, the confinement ability of the system decreases, paired with the occurrence of so called neoclassical tearing modes (NTMs)Buttery et al. (2000).

NTMs were recognised as a source of major disruptions in experiments such as JETde Vries et al. (2011, 2014), and thus their stabilisation is central. Amongst proposed stabilisation approaches, driving currentReiman (1983); Yoshioka, Kinoshita, and Kobayashi (1984) at the islands with RF waves has stimulated a long list of added efforts,Sauter (2004); Kamendje et al. (2005); La Haye (2006); La Haye et al. (2008); Henderson et al. (2008); De Lazzari and Westerhof (2009); Volpe et al. (2009); Sauter et al. (2010); Bertelli, De Lazzari, and Westerhof (2011); Hennen et al. (2012); Smolyakov et al. (2013); Ayten et al. (2014); Borgogno et al. (2014); Volpe et al. (2015a); Fevrier et al. (2016); Wang and Ma (2015); Li et al. (2017); Grasso et al. (2016, 2018); Poli et al. (2015) including many experimental demonstrationsBernabei et al. (1998); Warrick et al. (2000); Gantenbein et al. (2000); Zohm et al. (2001); Isayama et al. (2000); La Haye et al. (2002); Petty et al. (2004). By driving current at the centre of the island using ECCDFisch and Boozer (1980); Karney and Fisch (1981) or LHCDFisch (1978); Karney and Fisch (1979) one may balance the lack of bootstrap current, and prevent the island growth.De Lazzari and Westerhof (2009) This technique is, however, limited to its application to islands of smaller size due to available power constraints. This makes driving current precisely at their centre difficult, as the deposition width is comparable to the island size.

It has been recently suggestedReiman and Fisch (2018) that some of these stringent requirements may be relaxed due to the so called ”RF current condensation effect”. This effect takes into consideration the non-linear feedback of temperature variationsWesterhof et al. (2007) resultant from RF wave heating onto the deposition itself. It was found that condensation could improve mitigation as well as reduce radial sensitivity.

Formally, in Ref. [Reiman and Fisch, 2018] this non-linear feedback was modeled using a simplified diffusion energy balance equation that included resonant power deposition for a prescribed profile. Yet, the lack of a dissipation mechanism (eg. radiation) or the unlimited absorbable power from RF waves, gave, as they observed, an nonphysical temperature blow-up beyond a bifurcation point. Here that model is extended to include the damping of the RF wave.

In what follows, this extended model is first introduced in detail. Next, the equation is analytically and numerically solved and a hysteresis effect related to the island heating is described. Having introduced this phenomena, the effects of the non-linearities on different RF deposition schemes are explored. In order to evaluate these effects fairly, comparisons are made to the analytic linear solution, which is taken as representative of current approaches.

II Fundamental theoretical model

The model describes both the temperature variations of a magnetic island and the RF wave deposition. The latter may be described as a wave that is being damped along a ray trajectory, so as to provide the plasma with energy. The temperature of the plasma, driven by the RF waves, is described as part of an energy balance model with thermal diffusion. To construct such a model, various approximations are introduced.

Take, as a starting point, the transport equationBraginskii (1965), representing the second moment of the Boltzmann equation, to describe temperature, $T$, dynamics. Under the assumption of no significant flows, and considering fast equilibration between electrons and ions, a single equation may be written combining the two-fluid Braginskii equations

[TABLE]

where $n$ is the plasma density, $\kappa$ is the generalised heat conductivity tensor, and $P$ a volumetric power deposition which will be later related to the RF power input.

It has been stated as an assumption for Eq. (1) that electrons and ions are effectively equilibrated. Thus, the equation should only apply to those time scales larger than the typical equilibration time $\tau_{\mathrm{eq}}$; this constitutes the first temporal constraint: $t\gg\tau_{\mathrm{eq}}$, where $t$ denotes the time scales that the model is suited to describe.

A second point is related to the use of temperature as a measure of plasma energy. The concept of temperature customarily applies only to thermalised systems, in which the populations in $v$-space are Maxwellian distributed. However, the continuous injection of RF waves distorts the distribution function locally so as to deposit energy resonantly onto a small fraction of faster electrons (eg. for LH waves $v\approx 4.5v_{Te}$ and for EC $v\approx 3v_{Te}$)Fisch (1987); Karney and Fisch (1981). Hence, the plasma is made up of a Maxwellian bulk with a well defined $T$, but also a resonant minority population. If the bulk $T$ is to represent the total internal energy of the plasma, then the energy drawn locally in $v$-space needs to be redistributed quickly. Here a second time ordering is introduced: collisional thermalisation ($\tau_{\mathrm{se}}$) and isotropisation ($\tau_{\perp}$) processes must be faster than the time scales of interest ($t\gg\tau_{\mathrm{se}},\tau_{\perp}$). In that case non-Maxwellian features that affect a minute fraction of the total population may be ignored to leading order. Kinetic details will still prove important for $P$.

As it stands, Eq. (1) remains a three dimensional problem, but it may be cast into an approximate reduced 1D problem assuming the following. Consider as it happens (see Table 1)Spitzer and Härm (1953); Westerhof et al. (2007), transport over a given magnetic flux surface to be much faster than perpendicular to it; that is, $f_{\kappa}\equiv\kappa_{\parallel}/\kappa_{\perp}\gg 1$. In that case, magnetic surfaces will be approximately isothermal, simplifying the derivatives in the diffusive term (i.e. the second term in Eq. (1)). For analytic simplicity, the geometrical factor associated with the particular shape of island flux surfaces will be ignored, using instead a single slab coordinate $x$ (one may think of making a cut to an elongated, narrow island). This reduced form is,

[TABLE]

As a result of this slab adaptation, areal weighting is made equal for all points, though in reality this should be larger for the edges. A more complete treatment considering the flux coordinate is left for future work, though it was shown itn [Reiman and Fisch, 2018] that the slab geometry shared the qualitative physics with the more realistic geometry.

Nonetheless, the model retains an important feature of the magnetic island geometry: the closed nature of magnetic surfaces about the island centre. Consequently, in the reduced 1D model, temperature solutions are required to be even about the centre $x=0$.

Proceed now to linearise Eq. (2). Assume that the temperature changes ($\widetilde{T}$) in the island due to directing RF heating ($P$) to it are small; ie. $\epsilon=\widetilde{T}/T_{0}\ll 1$, where $T_{0}$ is the equilibrium temperature. Seeking precision in this definition, $T_{0}$ is defined to be the island temperature when the RF power is not aimed directly at the island, but is rather part of the total power budget that heats the centre of the tokamak (see Fig. 1). Because $T_{0}$ is constant over the island (as is the density), any island inhomogeneity that develops will be at least $O(\epsilon)$; and with this ordering we drop derivatives with respect to $\kappa_{\perp}$ in Eq. (2),

[TABLE]

The consideration of stiff temperature profiles that could modify $\kappa_{\perp}$ non-smoothly is left for future work.

To drop the time dependence of Eq. (3), the time scales of concern should exceed those of energy diffusion ($\kappa_{\perp}\partial^{2}_{x}\widetilde{T}$ term) and the driving times ($P$ term). Consider the former; clearly for the dominance of the diffusion term $t\gg\tau_{D}=W_{i}^{2}/\chi_{\perp}$. Here $W_{i}$ represents the island width and $\chi_{\perp}\approx\kappa_{\perp}/nk_{B}$ is the heat diffusion coefficient.

For the latter, $nk_{B}\partial_{t}T\sim P$ suggests that the time derivative may be dropped provided the power density of the wave has had enough time to deposit all the needed thermal energy; i.e., $t\gg\tau_{E}=nk_{B}\Delta T/P$, where $\Delta T\sim\epsilon T_{0}$ is the characteristic variation of the island temperature. Then, with $t$ in this regime, and being consistent with all of previous requirements,

[TABLE]

Prior to detailing the form of $P$, the problem should be closed by both defining the spatial domain and setting appropriate boundary conditions. Naturally, one defines the last closed surface of the island including the X-points as boundaries of the domain of $x$, ie. $|x|\leq W_{i}/2$. The island width, $W_{i}$, will be kept constant. That means that $t$ must be shorter than the typical island growth $\tau_{i}=\left(\partial\ln W_{i}/\partial t\right)^{-1}\approx\tau_{\eta}$, where $\tau_{\eta}$ is the global resistive time scale.Rutherford (1973) Table 1 shows this last requirement is consistent with previous time orderings. Width changes may then be treated adiabatically; ie. the steady state Eq. (4) may be taken to be satisfied at all times as $W_{i}$ is changed artificially.

The boundary condition on temperature encodes the influence of the island on the remaining of the plasma (and vice versa). To specify it, a simplified treatment of the energy dynamics of the rest of the tokamak is done, using a steady state diffusion model like Eq. (1). It is convenient to apply Gauß’ theorem to magnetic flux surfaces, so that $\int_{\psi}PdV=-\partial_{\psi}T\int_{\partial\psi}\kappa_{t}\hat{n}\cdot\nabla\psi~{}dS$, where $\psi$ is the magnetic flux coordinate. This shows that the slope of the temperature profile at a particular flux surface is determined by the power deposited inside it. For simplicity, let us associate a spatial 1D coordinate $R$, in the absence of islands monotonic with $\psi$, and take the heat conductivity, $\kappa_{t}$, to be constant; then, the slope of temperature at some $R_{0}$ is determined by the power deposited at $R<R_{0}$. With this in mind, let there be some heating in the tokamak centre: $P_{c}$ (fusion power, Ohmic heating, etc.) and $P_{\mathrm{RF}}$ (RF heating). The temperature profile is then determined by the heat flux and the fixed plasma edge temperature (see Fig. 1a).

Now, let there be an island of size $W_{i}$ at a distance $R_{i}$ from the core over which the temperature profile is flat (see Fig. 1b). Because heat sources have not changed, the temperature slope remains unchanged elsewhere. Let then $P_{\mathrm{RF}}$ be redirected to the island (ie. the case of interest). For those magnetic flux surfaces at $R>R_{i}+W_{i}/2$, the enclosed total power does not change, and thus the slope of $T$ should neither (see Fig. 1c). Given that the tokamak plasma edge temperature is fixed, the temperature of at the edge of the island, $T_{0}$, remains unchanged. The boundary condition for our island temperature may then be taken to be $\widetilde{T}(x=\pm W_{i}/2)=0$.

It is now the turn of specifying $P$ in Eq. (4) to represent the energy deposition from RF waves. Adopting a geometrical optics (GO) description of the wave envelopeTracy et al. (2014), the evolution of the energy of the wave may be written as

[TABLE]

where $\bar{V}$ represents the wave energy density, $v_{g}$ is the group velocity of the wave, $\omega_{t}$ represents the time derivative of frequency due to a time dependent medium, $d_{t}$ represents the total time derivative following a wave along a ray and $\gamma$ represents collisionless damping rate. Assuming the medium to be stationary in the wave damping time scale ($\omega_{t}/\omega\ll\gamma$) and the spatial inhomogeneity to be much smaller than the variation resultant from the damping ($\nabla\cdot v_{g}\ll\gamma$), the last approximated equality follows. This condition is not difficult to satisfy considering only small variations are created within the island.

Expressing Eq. (5) in terms of $x$, the distance along the ray $d_{t}x=v_{g}$,

[TABLE]

which has the form of damped propagation. The factor $\gamma$ may be obtained under the assumption of a Maxwellian magnetised backgroundStix (1992); Karney and Fisch (1981), one may show that for EC and LH, $\gamma\propto\exp(-\chi^{2})$ where $\chi=(\omega-n\Omega_{e})/kv_{Te}=v_{\parallel}/v_{Te}$ where $v_{\parallel}$ is the phase velocity of the wave, $v_{Te}$ is the electron thermal speed, $\Omega_{e}$ is the electron cyclotron frequency, and $n=0$ corresponds to LH waves and $n=-1$ to EC waves. The power deposition both for electron cyclotron and lower hybrid waves occurs on the tail of the Maxwellian velocity distribution, with damping exponentially small in the lowest resonant velocity.

It is this exponential factor which makes deposition highly sensitive to variations in temperature. Indeed, considering the phase velocity of the wave to remain constant over the extent of the island,

[TABLE]

It is convenient here to define the dimensionless variable $u\equiv v_{\parallel}^{2}\widetilde{T}/v_{Te,0}^{2}T_{0}=w^{2}\widetilde{T}/T_{0}$, where $w=v_{\parallel}/v_{T_{e},0}$. The location of the wave damping within the whole plasma is generally dependent on $T$, $B$, $\omega$ and $k$. This location can be determined using ray tracingPrater et al. (2008); Bonoli and Englade (1986). Here the picture is simplified by artificially restricting the damping to a particular defined region within the island , while keeping $v_{\parallel}$ constant, and that way allowing for the $\widetilde{T}$ expansion. A fully self consistent, full GO analysis is left for future work.

Let us express Eq. (6) as

[TABLE]

where the damping strength $\alpha(x)\equiv 2\alpha_{0}f(x)/W_{i}$. The factor $\alpha_{0}$ represents the strength of the damping, but it is also defined in a dimensionless way to include the island width $W_{i}$. For example, for EC wavesKarney and Fisch (1981) $\alpha_{0}\approx W_{i}\sqrt{\pi}\omega_{pe}^{2}\exp(-w^{2})/2ckv_{Te}$. It is helpful to introduce a more physically motivated interpretation of $\alpha_{0}$. If a linear limit is taken of Eq. (7), the power deposition profile takes the form $|V^{\prime}|\propto\exp\left(-2\alpha_{0}x/W_{i}\right)$; i.e., $\alpha_{0}$ is the ratio of the island half-width to the characteristic deposition width. Note that the deposition has an exponential shape, and not its usual Gaussian form generally considered for electron-cyclotron wavesPrater et al. (2008); however, both schemes are peaked and of finite width, and ultimately quite similar.

Now, going back to the original question: how is $P$ related to this wave energy $\bar{V}$? From the damping of the wave along a ray, it is easily seen that the volumetric power deposition at a given point is given by $P=-d_{t}\bar{V}\approx-v_{g}d_{x}\bar{V}\equiv-v_{g}\bar{V}^{\prime}(x)$. But, because in this particular geometry the points $\pm x$ are linked together (recall this is true due to them belonging to the same flux surface) these points share the total deposition at $x$ and $-x$. All things considered,

[TABLE]

These equations may be non-dimensionalised, reducing them to

[TABLE]

where the new $V(\widetilde{x})=\bar{V}(x)W_{i}v_{g}v_{p}^{2}/2\alpha_{0}\kappa_{\perp}T_{0}v_{T}^{2}=\bar{V}(x)W_{i}\Upsilon^{2}/\alpha_{0}$ and $\widetilde{x}=2\alpha_{0}x/W_{i}$. Note that the edges are now at $x=\pm\alpha_{0}$.

In order to complete the setting of the problem, an initial value must be taken for Eq. (9a). Let $V_{X}\equiv V(-\alpha_{0})\equiv V_{0}/\alpha_{0}$, where $V_{0}$ is a constant representing some wave energy density input. For interpreting solutions, it is important to bear in mind that $V_{X}$ is independent of island width, but will however scale as $1/\alpha_{0}$ with the deposition strength.

For clarity in the following Sections III and IV, one may refer to Appendix B as a quick reference for the variables employed.

III Hysteresis phenomena

To investigate the effect of the non-linear wave deposition, consider the tractable basic problem of wave damping occurring everywhere within the island. This case, represented by $f(\widetilde{x})=1$, allows for an analytic solution of Eqs. (9a) and (9b) (see Appendix A for a detailed derivation). Implementing the appropriate boundary and initial conditions,

[TABLE]

where the paramters $\lambda$ and $\gamma$ are determined by,

[TABLE]

The integration constant is $\lambda=(V_{X}+V_{f})/2$, where $V_{f}$ is the energy density when exiting the island. Eqs. (11a-11b) solve $\lambda$ implicitly, which ultimately determines the temperature of the island as a function of $\alpha_{0}$ and $V_{X}$ from Eq. (10).

The dependence of perturbed central island temperature on these two parameters will be represented as contour curves (see for example Fig. 2). Two main representations are of particular physical interest. First, contours of constant $\alpha_{0}$ (ie. fixed deposition strength and island width) in the $u(0)-V_{0}$ plane. These contours show the effects of the wave power on temperature (see Fig. 2a). The second interesting picture is related to how the heating of the island evolves as its width or the wave profile width changes. This is captured by curves of constant $V_{X}$ or $V_{0}$ at fixed deposition strength in the $u(0)-\alpha_{0}$ plane (see Fig. 2b for an example).

There are a number of general features in the solutions to Eqs. (10)-(11b) worth highlighting. The first of those is the existence of bifurcation points. As previously observedReiman and Fisch (2018), for sufficiently broad depositions (small $\alpha_{0}$), saddle-node bifurcations appear, at which the two lower temperature solutions disappear. Such points are marked in Fig. 2a. The appearance of these points may be linked to the action of a self-focusing mechanism affecting RF waves. Schematically, below the bifurcation, significant energy leakage takes place. As the bifurcation is approached, island temperature perturbations become larger, while the power deposited by the wave increases accordingly. This positive feedback eventually extracts all RF energy effectively, reaching a higher temperature steady state and thus jumping into an upper branch of the solution.

This hot stable solution may be seen in Fig. 2, along with the asymptotic form of the solution as $V_{f}\rightarrow 0$ (dotted line). The proximity of the two solutions demonstrates that the upper branch indeed corresponds to nearly complete deposition of the wave energy in the island. Such a solution branch is also, immediately after the bifurcation point, significantly larger than the linear prediction (see broken lines).

A consequence of the solution structure obtained is the hysteresis behaviour of island heating. To illustrate such a process, take as a starting point the system to be in equilibrium at the lower temperature branch in Fig. 2a, and increase the absorbed energy of the incoming wave ($V_{0}$) gradually. As a result, the temperature of the island will grow until the bifurcation point is reached. Once at this point, and driven by the self-focusing feedback, the temperature of the island will rapidly increase towards the upper branch, which is the only stable solution at high absorbed power.

The hot island exhibits, at this point, a large temperature difference between the centre and the separatrix. A priori, this would help to absorb RF power closer to the O-point, and thus also drive thte central currentFisch (1987). Here we ignore current associated with the DC electric field, which could be due to Spitzer conductivity or due to the hot electron conductivityFisch (1985); Karney, Jobes, and Fisch (1985). These currents are less important than the directly driven RF current.Reiman and Fisch (2018) With such presumed centred current drive, the island would tend to shrink and stabilise, as governed by the Rutherford equationRutherford (1973). This size reduction corresponds to a leftwards displacement towards smaller values of $\alpha_{0}$ in Fig. 2b. In such a case, and if the energy available to the island is maintained, the reduction does not imply a return back to the original low temperature, but instead remains in the more effective current driving upper branch for some time. Similarly, once in the upper branch, driving power requirements are relaxed, and lower $V_{0}$ would still keep the plasma hot. This constitutes the hysteresis effect. A more careful discussion on the usefulness, accessibility and consistency of this sketched simplified picture for particular deposition schemes is the concern of following sections.

IV Effects of non-linear feedback

IV.1 Typical parameters

Before proceeding further, a brief estimate and collection of typical values for both $V_{0}$ and $\alpha_{0}$ is presented. We emphasise that our slab model provides a physical qualitatively correct picture of the problem, but only a rough guide into the quantitative behaviour of the more realistic geometry, as previous calculations suggestReiman and Fisch (2018). Other simplified features, such as the exponential form of the linear deposition profile, are also different when compared to actual experimentsPrater et al. (2008), but do however share the fundamental characteristics. Thus, the linear case will be taken as reference in guiding conclusions in the following sections, as well as orientative comparison standard to existing experimental parameters.

Focus first on the values for the parameter $\alpha_{0}$. Two different routes are taken at this point. One possible method uses, given the definition of $\alpha_{0}$ as the size of the linear power deposition width, typical deposition widths in tokamak experiments could be used to obtain $\alpha_{0}\sim 0.5-3$.Volpe et al. (2015b) It has been recently reported that current drive profiles are in experiment subject to broadeningBrookman, M. W. et al. (2017) by factors of 2–3 due to effects unaccounted for in ray tracing routines, such as edge density fluctuations. This effective broadening could make typical $\alpha_{0}$ values even lower, down to $\sim 0.2$.

Alternatively, one could use the form for $\alpha_{0}$ given before and obtained in the context of GO. In the case of ECCD, for instance, using typical approximated hydrogen tokamak values (see caption of Table 1), with $w^{2}\sim 10$,Karney and Fisch (1981) wavenumberPrater (2004) $k\sim 2\pi/(5~{}\mathrm{mm})$, density $n\sim 10^{20}$ m*-3* and $W_{i}\sim 10~{}$cm; $\alpha_{0}\approx W_{i}\sqrt{\pi}\omega_{pe}^{2}\exp(-w^{2})/2ckv_{Te}\sim 10^{-1}$.

Now consider the wave power density $V_{0}$. First, we may compare existing literatureWesterhof et al. (2007) where island temperature variations are computed to our linear toy model. Comparing values of $u(0)$,orientative typical power parameters on the order of $V_{0}\sim 5$ for 20 MW RF power are found. In a more first principle approach, we might use the rescaled definition of $V_{0}$ introduced before. Using an RF power on the order of $P\sim 10$ MW, with a beam of cross section $A\sim 1~{}$m2, with $\chi_{\perp}\sim 1$m2/s, temperature $T_{0}\sim 10$ keV, density $n\sim 10^{20}$ m*-3* and $w^{2}\sim 10$, then $V_{0}\approx PW_{i}w^{2}/2A\chi_{\perp}nk_{B}T_{0}\sim 10^{1}$.

Summarising:

[TABLE]

IV.2 Central deposition

The spatial distribution of the RF deposition strongly affects the final temperature of the island, as well as the island mitigation efficiency. In this section, the best case scenario is first analysed; i.e., deposition starting from the island centre. To formally emulate this ideal case, $f(\widetilde{x})=H(\widetilde{x})$, where $H$ is the Heaviside step function. Given this newly introduced asymmetry, well defined parity is lost from the equations and the solution to the equation is only found numerically.

Consider first the occurrence of bifurcation points, and in particular, how they depend on $\alpha_{0}$. To illustrate these points, Figures 3 and 4 are presented.

In the broad linear deposition limit, with $3\alpha_{0}\ll 1$, there always exists a bifurcation point (see Fig. 3). The turning point, however, occurs at increasingly larger wave energy densities $V_{0}$. One may understand this result by referring to the analytic asymptotic form of the solution at low $V_{0}$. In that limit, the system takes the form of the linear problem, for which $u(0)\sim V_{0}\left(\alpha_{0}-1+e^{-\alpha_{0}}\right)/2\alpha_{0}$. This shows that the temperature of the island becomes decreasingly responsive as $\alpha_{0}\rightarrow 0$ (see the decreasing initial slopes of curves in Fig. 3 inset), which is ultimately related to there being a significant wave energy leakage ($V_{f}=V_{X}e^{-\alpha_{0}}$).

Only for those cases for which the initial energy leakage $V_{f}$ is significant will a bifurcation occur. The bifurcation is a result of the system being able to access all that previously lost energy when the damping $e^{u}$ factor becomes significant. This jump will be associated with a narrowing of the deposition and a current that is more efficiently utilised in stabilisation.

As the initial deposition is reduced by increasing $\alpha_{0}$, $V_{f}$ in the lower branch decreases, and it eventually becomes too small to sustain a bifurcation. Fig. 3) shows the boundary value $\alpha_{0}=0.32$ beyond which no bifurcation occurs.

Where no bifurcation occurs, the temperature of the island only undergoes a smooth transition in temperature between the linear solution and the high temperature asymptote (see Fig. 3). That limiting form of the nonlinear solution as $V_{f}\rightarrow 0$ is $u(0)\sim V_{0}/2$ (see dashed line in Fig. 3 inset and Appendix C), which is also the limit as $\alpha_{0}\rightarrow\infty$ of the linear deposition. That is, the non-linear response serves as a short cut via self-focusing to the linear ideal infinitely narrow deposition. Thus, for the case of centralized deposition, the nonlinear mechanism always leads to an enhanced temperature increase.

However, the wave power required to obtain a substantial improvement exhibits strong dependence on $\alpha_{0}$ as shown in Fig. 4a. The large energy leakage and the small difference between linear and non-linear solutions at low and high $\alpha_{0}$ respectively leave a most easily accessible (lower $V_{0}$) central region at values $\alpha_{0}\sim 1-1.5$. To emphasise the second of these limitations Fig. 4b shows the ratio of the analytic asymptotic forms of the non-linear and linear solutions. Evidently, the differences become marginal (ie. the ratio tends to one) for stronger depositions, which explains why Fig. 4a diverges at $\alpha_{0}\sim 2.8$.

Finally, we examine the extent to which self focusing narrows the RF deposition profile. From Fig. 4c the self-focusing effect is apparent, and will undoubtedly improve stabilisation by bringing current drive closer to the O point. This reduction by a factor of $\sim 2-3$ in width opens the possibility of previously disregarded regimes of island stabilisation. Note also that for $\alpha_{0}\approx 2.8$, the non-linear deposition profile becomes similar to that of a delta function (within discretisation).

This analysis suggests that the region of interest and current experimental relevance may in some subset of cases (for the broadest depositions) show some hysteresis behaviour, but most will just show significant temperature variation. In addition, as a result of deposition narrowing due to the non-linearity, islands could be stabilised when traditionally predicted not to. This opens the door to experimental verification of the non-linear effect, as well as extension of mitigation schemes.

IV.3 Edge deposition

The scenario adopted for analysis before was that of central deposition. This is the ideal case, and so it presumes that one is experimentally capable of aiming perfectly at the centre of the island without depositing any energy before that. But, what would happen if the deposition departs from this idealised case? The worst case scenario is now presented. To that end, we recover the analytic solution from Eqs. (10)-(11b) which represents wave deposition from the very edge of the magnetic island.

With this analytic result at hand, let us explore first the limit of complete wave deposition: $u(0)\sim\log\left(V_{0}/2\alpha_{0}+1\right)$. The linear solution gives $u(0)=V_{0}e^{-\alpha_{0}}(\cosh\alpha_{0}-1)/\alpha_{0}$. It is remarkable that the non-linear model gives a logarithmic growth of the island temperature as the power input is increased, while the linear case grows linearly. It necessarily follows that some non-linear inhibition mechanism must take part. Indeed, one may relate this to the deposition profile becoming localised ever closer to the edge of the island.

It is the same self-focusing that narrowed the deposition closer to the centre when central deposition was considered, which displaces deposition towards the island edge (see Fig. 6). Physically, the RF wave becomes so strongly damped that it runs out of energy very close to the edge. There, the temperature slope is large, and thus as heat gets to the centre of the island, it is quickly lost across the edges, $u(0)$ becoming limited. For typical values $\alpha_{0}\sim 1$ and $V_{0}\sim 5-10$, the peak of deposition $x_{\mathrm{peak}}\sim-0.7$; ie. power is deposited somewhere between the X- and O-points.

As a result of this detrimental displacement, only over a limited regime will the non-linear solution be hotter than the linear one. As shown in Fig. 5, this interval is larger for the broader depositions, but tends to disappear as $\alpha_{0}\rightarrow\infty$. When the deposition is broader, the non-linear focusing increases the amount of deposited power significantly, a benefitial addition that outweights the inhibiting displacement for a more extended range of powers.

In addition to temperature, this deposition shift will also bring the driven current closer to the X-point. This displacement can be catastrophic when trying to mitigate the growth of magnetic islands. The proximity of the deposition to the edge may be seen in Fig. 6. The plot shows that indeed, for large $V_{0}$ values, the non-linear self-focusing brings the deposition ever closer to the edge. This is the result of the large damping $e^{u}$ which drains the incoming wave faster than the linearly increasing $V_{0}$. Nevertheless, there is an initial region in which the focusing is starting to affect the system, and clearly becomes beneficial in terms of drawing power towards the centre.

To further explore the sensitivity of stability on $x_{\mathrm{peak}}$, the procedure in [De Lazzari and Westerhof, 2009] is followed. The relative Fourier weighting to $\Delta^{\prime}$ in the Rutherford equation due to driving current at a particular flux coordinate, $\psi$ ($\Omega$ in the reference), may be estimated taking the 1D spatial variable $x$ in our model to match the spatial $x$ in [De Lazzari and Westerhof, 2009], looking at $\xi=0$ (see Fig. 7). The calculation shows that within approximately 90% of the island the drive is stabilising. Therefore, looking back at Fig. 6. The current drive will still be central enough to be stabilising for a significant fraction of the cases, even when the non-linear solution is colder than the linear one. No definitive conclusion may however be drawn on the precise fraction of the island that is truly stabilising, as the treatment of the island geometry in [De Lazzari and Westerhof, 2009] is different from that of our model. A fully consistent treatment is left for future work.

We have shown that an initially broad RF profile may then be used to stabilise islands, so long as the input power remains below some upper bound. This result is promising and an idea to further explore.

More generally, if deposition was to start midway between the X- and O-points, then deposition would be driven and narrowed towards that location instead. The trends will then fit between the extreme centre and edge cases shown explicitly, with less constraining requirements as the centre is reached. In this context, RF overshooting scenarios will never suffer from the inhibition that takes place when undershooting with increasing $V_{0}$.

Yet another possibility to circumvent inhibition might be to look for means to amplify the wave power within the island, such as through an $\alpha$-channeling effectFisch and Rax (1992). This could give a power not peaked at the island periphery. Unfortunately, this volumetric amplification does not occur for electron cyclotron waves, but might be exploited when using lower hybrid wavesOchs, Bertelli, and Fisch (2015). Despite lacking this possible enhancement, EC waves do however have the benefit of having a $B$-dependent resonance that allows for deposition starting at a particular point in space.

V Conclusion

The possibility of hysteresis involving the heating of magnetic islands with RF waves is shown for a symmetrised, 1D slab model. The wave power deposition self focuses mediated by island temperature, leading to higher temperatures. Past a bifurcation point, the island may remain in this high temperature solution even as it shrinks or power is reduced. Exploiting the hysteresis effect could thus provide an easier and improved way to eliminate magnetic islands.

In typical parameter regimes of current experiments, though, bifurcation is likely to occur only for the broadest profiles. Yet, centred deposition scenarios show that important temperature increments (on the order of $\sim 20\%$) occur. These differences are most significant in broader deposition schemes, where the self-focusing mechanism of the non-linear model is most different from the linear model. Alongside these thermal variations, the non-linear narrowing of the RF deposition profile by a factor of up to $\sim 2-3$ for typical values will improve the utility of current drive for purposes of stabilisation.

Broad schemes that deviate from the centre still lead to not dissimilar self-narrowing and stabilisation under certain circumstances. These circumstances involve the form of EC power density profiles, and must be considered in designing deposition scenarios. In particular, for deposition profiles that peak before the O-point, there is a threshold power density above which a self-inhibition mechanism is encountered; beyond the O-point, this threshold does not exist. This opens the door to exploring previously disregarded broad RF deposition stabilisation schemes.

Acknowledgements.

Thanks to Suying Jin for fruitful discussions and help. This work was supported by US DOE DE-AC02-09CH11466 and DE-SC0016072.

Appendix A Analytic solution to constant $\alpha$

Consider the coupled set of equations,

[TABLE]

Define the following symmetric and antisymmetric parts of the wave energy density,

[TABLE]

Given these, Eq. (13b) may be cast in the form,

[TABLE]

Substituting definitions (14) and (15) into Eq. (13a),

[TABLE]

and realising that the spatial derivative $\frac{d}{dx}$ is an odd operator while $u$ is an even function, the equation may be separated into its symmetric and asymmetric parts,

[TABLE]

which with Eq. (16) form a set of three coupled non-linear differential equations. The set avoids the explicit appearance of $-x$ on the expense of an additional equation.

Because of even parity $u^{\prime}(\widetilde{x}=0)=0$, so Eq. (16) may be integrated

[TABLE]

Substituting (19) and (16) into equations (17) and (18), putting together and integrating once

[TABLE]

where $\lambda$ is an integration constant that must satisfy the initial conditions specified. From Eq. (16) it follows that $\lambda=(V_{X}+V_{f})/2$, which limits $V_{X}/2<\lambda<V_{X}$.

To solve Eq. (20), multiply both sides of the equation by $u^{\prime}$ and integrate with respect to $\widetilde{x}$. Thus,

[TABLE]

where $\mathbf{C}$ is another integration constant to be determined later on.

Let the substitution $u\equiv-\log z$ be implemented in Equation (21),

[TABLE]

At his point, one may try a symmetric solution of the form $z=A\left(e^{\gamma\widetilde{x}}+e^{-\gamma\widetilde{x}}\right)+B$, where $A$, $B$ and $\gamma$ are to be reduced to a single integration constant. If this is possible, uniqueness guarantees this to be the general solution. After some manipulation,

[TABLE]

which indeed leaves a single degree of freedom, $\gamma$, as expected for a first order ODE. It is then time to implement boundary conditions to determine $\lambda$ and $\gamma$. From $u=0$ at $\widetilde{x}=-\alpha_{0}$, $z(\widetilde{x}=-\alpha_{0})=1$ is evaluated, and eliminating $A$,

[TABLE]

Having eliminated $A$, the solution for $u$,

[TABLE]

The additional boundary or initial condition may be imposed requiring $u^{\prime}(-\alpha_{0})=A(-\alpha_{0})=(V_{X}-V_{f})/2=V_{X}-\lambda$ from Eq. (20), which will introduce explicitly the physically relevant parameter $V_{X}$. Then,

[TABLE]

where $\gamma^{2}>1+2\lambda$.

Appendix B Summary of variables and abbreviations

A collection of the variables used in Sections III and IV is presented as a reference in Table 3. A brief description is also provided where relevant.

Appendix C Asymptotic limit non-linear central deposition

Let us consider the limiting case for complete RF power deposition in island stabilisation starting deposition from the centre of the magnetic island.

Begin with,

[TABLE]

Now take the wave energy to be damped quickly, so that $V\approx 0$ for $x>0$. In that case, and as the edge is located at $x=\alpha_{0}$:

[TABLE]