1. Introduction
In tokamak experiments, electron-cyclotron (EC) waves are launched into the plasma in the form of spatially narrow beams, and exchange energy with electrons passing the beam region when the specific resonance condition is satisfied (Brambilla Reference Brambilla1998),
\begin{equation} \omega-\frac{l_h\omega_c}{\gamma}-k_\parallel v_\parallel\approx 0, \quad l_h=\pm 1,\pm 2,\ldots, \end{equation}
with 
$\omega$, 
$\omega _c$ the wave and cyclotron frequencies, respectively, 
$\gamma$ the Lorentz factor and 
$k_\parallel$, 
$v_\parallel$ the wavevector and particle velocity components parallel to the magnetic field. When 
$k_\parallel \approx 0$, the result of the wave–particle interaction is plasma heating, whereas in cases where 
$k_{\parallel }\neq 0$ non-inductive current is generated (Erckmann & Gasparino Reference Erckmann and Gasparino1994). Apart from the wave frequency and geometry, the power absorption and current drive efficiency depend on the magnetic equilibrium and the spatial distribution of the plasma density and temperature. Since the tokamak magnetic field is non-uniform, under specific choices for the beam launching conditions, the EC resonance may be realized in a very narrow spatial region. In this respect, using EC waves gives the advantage of good localization in the power deposition, which is exploited in applications relevant to auxiliary plasma heating (Prater Reference Prater2004), plasma ramp-up and breakdown assistance (Granucci et al. Reference Granucci, Ricci, Farina, Figini, Iraji, Tudisco, Ramponi and Bin2012), sawtooth and neoclassical tearing mode (NTM) instability control (LaHaye Reference LaHaye2006), as well as plasma diagnostics based on electromagnetic (EM) waves like, for example, reflectometry and interferometry (Hutchinson Reference Hutchinson2005).
From the EC-relevant applications, NTM stabilization is of central interest in this work. Neoclassical tearing modes are resistive magnetohydrodynamic (MHD) instabilities that limit the performance of tokamaks, because the magnetic islands they generate reduce the plasma energy and angular momentum, leading to significant loss of confinement (LaHaye Reference LaHaye2006). An NTM appears in a combination of the finite plasma resistivity and magnetic perturbations, which leads to the appearance of helical currents and tearing of magnetic flux surfaces of rational safety factor, i.e. 
$q=m/n$ where 
$m$, 
$n$ are the poloidal and toroidal mode numbers (Wilson Reference Wilson2012). As a result of energy conservation, what follows is magnetic reconnection and the formation of island structures, which in turn provoke the loss of axial symmetry and plasma pressure flattening. The NTM process starts from an initially small island (also known as seed island) and can grow to a large amplitude due to the lack of pressure-gradient-driven bootstrap current, leading to a reduction of the plasma current and, possibly, device disruption. Supported also by theory, various techniques are under development for the control of these unstable modes, among which the application of localized EC current drive (ECCD) has been demonstrated as very effective in several experiments (see e.g. Maraschek (Reference Maraschek2012) for work done in ASDEX Upgrade).
In principle, theory models for EC propagation should be based on coupling the wave equation with the medium response function. The corresponding full-wave methods are the most complete, but also most resource-demanding option. For high-frequency waves in fusion devices, in most cases of interest the wavelength is much shorter than the plasma inhomogeneity scale, and the complexity is reduced by applying frequency-domain asymptotic methods. These techniques are based on geometric optics (GO), where the solution is obtained via canonical differential equations, using as Hamiltonian the dispersion function (Bernstein & Friedland Reference Bernstein and Friedland1980). The most popular of these tools is ray tracing, which provides a physical picture of the propagation in terms of wave rays; however, its validity becomes questionable when diffraction effects are present (Tsironis Reference Tsironis2013). More efficient asymptotic techniques, which retain the GO description along propagation but also solve for the transverse beam profile including diffraction, are quasi-optical (Nowak & Orefice Reference Nowak and Orefice1993) and paraxial beam tracing (Pereverzev Reference Pereverzev1998). Whenever the short-wavelength limit breaks down, the tools coming from asymptotic methods should be replaced by full-wave solvers. Such tools span from simpler techniques such as the transfer matrix (Born & Wolf Reference Born and Wolf1999), which is based on the electric field continuity conditions across physical boundaries, to more comprehensive methods such as time-domain finite differences (Taflove & Hagness Reference Taflove and Hagness2005), which discretize differential equations by approximating the derivatives as finite differences on a properly defined grid.
For the closure of the wave–plasma description, the plasma response is required in terms of the wave-induced electric current. Most models assume cold plasma dispersion for the propagation plus cyclotron damping only within the resonance area, motivated by the fact that, in modern fusion devices, the wave intensity is small and the wave–plasma coupling far from EC resonance is very weak (see e.g. Tsironis (Reference Tsironis2013) for details). The GO-based solvers implement this style by solving the real cold-plasma dispersion relation for the real part of the wavenumber and calculate its damping-related imaginary part from the complex dispersion relation only in the resonance area (Bernstein & Friedland Reference Bernstein and Friedland1980), whereas full-wave solvers employ a cold-plasma fluid equation or an effective hot-plasma dielectric tensor based on the linearized Vlasov equation (Tsironis, Samaras & Vlahos Reference Tsironis, Samaras and Vlahos2008). In all cases, turbulent structures may be incorporated as a perturbation function to the background plasma density. In case the incorporation of such structures without involving the associated steep density gradients is desirable, one may resort to EM homogenization methods (Mackay & Lakhtakia Reference Mackay and Lakhtakia2015) for deriving a plasma dielectric tensor that includes the density perturbations as a spatially averaged effect within the area of interest.
With reference to ECCD, theoretical modelling of the NTM stabilization process focuses on the effect of EC power deposition and current drive on the nonlinear mode dynamics. The magnetic island growth and decay may be computed using simple tools, based on the magnetic field diffusion equation (Wilmot-Smith, Priest & Hornig Reference Wilmot-Smith, Priest and Hornig2005), or more sophisticated models which solve the full MHD problem in the presence of kinetic effects. In the former category, the most established model is the modified Rutherford equation (MRE), which is a generalization of the Rutherford equation for classical tearing modes (Rutherford Reference Rutherford1973) by including the neoclassical physics term of the bootstrap current and the external electric current effect (the reader is pointed to Urso (Reference Urso2009) for details). Other effects which may play a crucial role in the EC deposition have been included in the context of the MRE, such as the significant changes introduced by the island topology in the magnetic field equilibrium and the plasma density and temperature profiles as compared with axisymmetric configurations, for example due to island asymmetry (Lazzaro & Nowak Reference Lazzaro and Nowak2009), particular flux surface nesting (Isliker et al. Reference Isliker, Chatziantonaki, Tsironis and Vlahos2012) or plasma rotation (Ayten & Westerhof Reference Ayten and Westerhof2012). More advanced models treat the wave–island interaction self-consistently by employing a closed set of EM wave and MHD equations (an example is the numerical code NIMROD (Jenkins & Kruger Reference Jenkins and Kruger2012)).
Dedicated experiments and theoretical studies with the tools mentioned above have established that, in order to achieve effective mode stabilization, the driven current should be highly localized around the island's 
$O$-point (in terms of both the peak value position and the profile width of the current density) and also its direction should be aligned with the equilibrium bootstrap current (Urso Reference Urso2009). However, the EC beam before reaching the targeted region crosses the plasma edge, where density fluctuations connected to turbulence may appear and induce additional wave distortion. This situation has a different consequence compared with the normal refraction in inhomogeneous plasma, which always generates an effect in the same direction to the propagation angle; turbulent structures scatter the wave spectrum components randomly, resulting more to diffractive-like broadening of the beam than a shift in its position. For larger tokamaks, like for example ITER and DEMO, where the distance from the launch points to the flux surfaces of rational 
$q$ may be several metres, even a small modification of the beam properties (also known as broadening and/or axial shift) may lead to a considerable deviation in the intended EC deposition location and extent (Tsironis et al. Reference Tsironis, Peeters, Isliker, Strintzi, Chatziantonaki and Vlahos2009).
Numerous articles have been devoted to the effect of electron density perturbations on the propagation of tokamak-relevant EC beams, and, for the most part, verified the aspects described above. We refer to a cross-section of these studies (see references therein for the complete amount of work), which involve the majority of asymptotic and full-wave models available in the literature: ray tracing (Hizanidis et al. Reference Hizanidis, Ram, Kominis and Tsironis2010; Peysson et al. Reference Peysson, Decker, Morini and Coda2011), quasi-optical beam tracing (Balakin, Bertelli & Westerhof Reference Balakin, Bertelli and Westerhof2011; Sysoeva et al. Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2015), Helmholtz equation solver (Ram & Hizanidis Reference Ram and Hizanidis2016), finite differences in time domain (Williams et al. Reference Williams, Koehn, O'Brien and Vann2014; Köhn et al. Reference Köhn, Holzhauer, Leddy, Thomas and Vann2016) and frequency domain (Papadopoulos et al. Reference Papadopoulos, Glytsis, Ram, Valvis, Papagiannis, Hizanidis and Zisis2019), finite elements (Chellai et al. Reference Chellai, Alberti, Baquero-Ruiz, Furno, Goodman, Labit, Maj, Ricci, Riva and Guidi2018) and wave kinetic solver (Weber, Maj & Poli Reference Weber, Maj and Poli2015; Snicker et al. Reference Snicker, Poli, Maj, Guidi, Koehn, Weber, Conway, Henderson and Saibene2017). A series of wave effects that may affect the intended beam quality have been identified, like for example partial reflection of the beam, radiation emission by the density perturbations and conversion of the wave polarization. Moreover, dedicated experiments for measuring the edge turbulent structures in the presence of EC wave beams have been conducted in TCV (Chellai et al. Reference Chellai, Alberti, Furno, Goodman, Koehn, Figini, Ricci, Hizanidis, Papagiannis and Tsironis2017), and the experimental results agree adequately with the theoretical predictions. Overall, assuming characteristic values within 
$[0.05, 0.5]$ of the background density for the fluctuations amplitudes and 
$[0.01, 0.2]$ minor radii for their spatial scales, the typical effects to the beam span within [
$10^{-3}, 10^{-1}$] degrees for the angular shift and 
$[0.05, 0.25]$ propagation lengths for the spot size broadening.
Something that has not been addressed extensively in the literature is the quantification of the effect of the beam alteration on the instability control effort (see e.g. Poli et al. (Reference Poli, Angioni, Casson, Farina, Figini, Goodman, Maj, Sauter, Weber and Zohm2015)). In this work, our main objective is to investigate the undesired effects of edge-resident density fluctuations to the EC wave propagation aimed at plasma heating and current drive, in connection with the efficiency of NTM stabilization. We use a modelling loop that incorporates all the system physics: the propagation of the wave beam through the perturbed layer and up to the magnetic island, the power absorption and current drive on the island's flux surfaces and, finally, the ECCD-driven mode dynamic evolution. Our tools include a full-wave solver for the propagation part prior to and inside the region of density perturbations, a ray tracing code (equipped with ECCD computation) for the propagation past the turbulent region, and an MRE solver for the computation of the island growth and decay. By assessing the input and output at each step of the modelling chain, an estimation of the consequences to the NTM control is made.
The structure of the paper is as follows. In § 2 we refer to the models used in our study (full-wave propagator in § 2.1, GO-based solution in § 2.2 and island dynamics in § 2.3). In § 3 we describe the set-up of the numerical simulations, giving special attention to the input–output connection of the three models and explaining the selection of the physics parameters, and in § 4 we present and discuss all the results, focusing on the implications of the density fluctuation intensity to the NTM decay time. Finally, in § 5 we draw the conclusions of this work and comment on the limitations of our modelling.
2. Description of the physics models
In this section we present the physics modelling chain for the beam propagation in the inhomogeneous plasma, the absorption and current drive on the magnetic island, and the mode response to the external current. The EC propagation through the tokamak edge is modelled using the transfer matrix technique for computing the wave electric field, in conjunction with a homogenized description of the plasma medium within the turbulent layer. The propagation after the edge until the island region, together with the current deposition on the flux surfaces of interest, is computed with a ray tracing code including the EC resonance and island geometry features. In the last part, the unstable mode dynamics is modelled in terms of the MRE in the presence of the external current.
2.1. Transfer matrix technique and EM homogenization method
In a plasma that contains density fluctuations, the use of asymptotic methods is not recommended because the short-wavelength limit breaks down when the size of these perturbations is comparable to the wavelength. There, the toolbox for modelling propagation requires a full-wave method for computing the beam field as well as a description of the plasma medium that includes the density perturbations. Our choice has been to combine an implementation of the transfer matrix technique for the wave propagation with an EM homogenization method for the perturbed plasma region.
We assume a Gaussian EC beam in a Cartesian coordinate system (
$x,y,z$) aligned with the tokamak geometry (
$\,\hat {\boldsymbol {y}}$ is the toroidal direction and 
$\hat {\boldsymbol {z}}$ the vertical direction). The beam propagates along the direction 
$\hat {\boldsymbol {k}}$ indicated by the wavevector, whereas its electric field amplitude spans across 
$\hat {\boldsymbol {k}}$ over a half-width 
$w$. For our calculation, the beam is analysed into 
$l$ plane-wave modes, where each mode propagates in the direction 
$\hat {\boldsymbol {k}}_{\boldsymbol {l}}$ implied by its specific wavevector, and their wavevectors are distributed symmetrically around the one of a central mode 
$l=0$. The lengths in the directions of beam propagation and amplitude profile are mapped to the Cartesian coordinate system via the projections of 
$\hat {\boldsymbol {k}}_{\boldsymbol {l}}$ on the position radius vector 
$\boldsymbol {r}$, which are correspondingly 
$\hat {\boldsymbol {k}}_{\boldsymbol {l}}\boldsymbol {\cdot }\boldsymbol {r}$ and 
$r[1-(\hat {\boldsymbol {k}}_{\boldsymbol {l}}\boldsymbol {\cdot }\hat {\boldsymbol {r}})^{2}]^{1/2}$.
Furthermore, each mode's propagation angle and refraction index with respect to the direction parallel to the magnetic field is defined accordingly as 
$\phi _l=\cos ^{-1}(\hat {\boldsymbol {k}}_{\boldsymbol {l}}\boldsymbol {\cdot }\boldsymbol {B}/|B|)$ and 
$k_{\parallel l}= \hat {\boldsymbol {k}}_{\boldsymbol {l}}\boldsymbol {\cdot }\boldsymbol {B}/|B|$. In this context, the electric field 
$E_l$ per mode is given by
\begin{equation} \frac{E_l\left (r\right )}{E_l^{\mathrm{max}}}=\exp\left\{-\displaystyle{\left [1-\left (\hat{\boldsymbol{k}}_{\boldsymbol{l}} \boldsymbol{\cdot}\boldsymbol{\hat{r}}\right )^{2}\right ]\frac{r^{2}}{w^{2}}}\right\}\exp\left (\mathrm{i}\boldsymbol{k}_{\boldsymbol{l}}\boldsymbol{\cdot}\boldsymbol{r}\right ), \end{equation}and the total beam electric field by summing all modes over their propagation angles and parallel refraction indexes (Valvis Reference Valvis2019),
\begin{align} \displaystyle{\frac{E\left (x,y,z\right )}{E^{\mathrm{max}}}}&= \sum_{\phi_l}\sum_{k_{\parallel l}}\exp\left \{\mathrm{i}\left [x k_{x}\left (k_{\parallel l},\phi_l\right )+yk_{y}\left (k_{\parallel l},\phi_l\right )+zk_{z}\left (k_{\parallel l},\phi_l\right )\right ]\right\} \nonumber\\ & \quad\times \exp\left\{-\left ( \displaystyle{\frac{w}{2}}\right )^{2}\left [k_x^{2}\left (k_{\parallel l},\phi_l\right )+k_y^{2}\left (k_{\parallel l},\phi_l\right )+k_z^{2}\left (k_{\parallel l},\phi_l\right )\right ]\right \}. \end{align} The transfer matrix technique is based on a phasor representation of the Faraday and Ampere laws (also known as Maxwell's curl equations). Denoting the electric and magnetic field normalized complex vector phasors as 
$\hat {\boldsymbol {e}}$ and 
$\hat {\boldsymbol {h}}$, these equations read
\begin{equation} \hat{\boldsymbol{h}}=\boldsymbol{n}\times\hat{\boldsymbol{e}}, \quad \boldsymbol{\mathsf{K}}\boldsymbol{\cdot}\hat{\boldsymbol{e}}=\boldsymbol{n}\times\hat{\boldsymbol{h}}, \end{equation}
where 
$\boldsymbol {N}=c\boldsymbol {k}/\omega$ is the refraction index vector and 
$\boldsymbol{\mathsf{K}}$ is the plasma dielectric permittivity tensor. If one formulates (2.3) in Cartesian coordinates (
$x,y,z$) and then eliminates the components 
$\hat {e}_x$, 
$\hat {h}_x$ using two from these equations, thus considering solving only for the field components along the toroidal and vertical directions 
$y$ and 
$z$, the four equations of interest may be cast in the form (Born & Wolf Reference Born and Wolf1999)
\begin{equation} N_xK_{xx}\boldsymbol{V}=\boldsymbol{\mathsf{M}}(\boldsymbol{N},\boldsymbol{\mathsf{K}})\boldsymbol{\cdot}\boldsymbol{V}, \end{equation}
where 
$\boldsymbol {V}=[\hat {e}_y\ \hat {h}_z\ \hat {e}_z\ \hat {h}_y]^{\mathrm {T}}$ is the solution vector and the 
$4\times 4$ matrix 
$\boldsymbol{\mathsf{M}}$ is given by
\begin{align} \boldsymbol{\mathsf{M}}=\begin{bmatrix} -N_yK_{xy} & K_{xx} - N_y^{2} & -N_yK_{xz} & N_yN_z \\ K_{xx}(K_{yy} - N_z^{2}) - K_{yx}K_{xy} & -N_yK_{yx} & K_{xx}(N_yN_z + K_{yz}) - K_{yx}K_{xz} & N_zK_{yx}\\ -N_zK_{xy} & -N_zN_y & -N_zK_{xz} & N_z^{2} - K_{xx} \\ K_{zx}K_{xy} - K_{xx}(N_yN_z + K_{zy}) & N_yK_{zx} & K_{xx}(N_y^{2} - K_{zz}) + K_{zx}K_{xz} & -N_z K_{zx} \end{bmatrix}. \end{align}
In general, the solution consists of two forward (
$+$) and two backward (
$-$) propagating waves, which may be determined by the real part of the 
$x$-component of the associated Poynting vector phasor 
$\hat {s}_x=\frac {1}{2}(\hat {e}_y\hat {h}_z^{*}-\hat {e}_z \hat {h}_y^{*})$: the forward waves correspond to 
$\textrm {Re}(\hat {s}_x)>0$, while negative values indicate backward waves and zero values indicate evanescence.
The problem may be solved in terms of the associated eigenvalue/eigenvector problem, where the respective eigenvalues specify the refraction indexes in the direction tangential to the problem boundaries. In physical space, the solution leads to (Valvis Reference Valvis2019)
\begin{equation} \boldsymbol{V}\left (x,y,z\right )=\boldsymbol{{\varPi}}\boldsymbol{\cdot} \boldsymbol{{\varPhi}}\left (x\right )\boldsymbol{\cdot} \boldsymbol{\mathcal{A}} \exp\left [\mathrm{i}\left (yk_y+zk_z\right )\right ], \end{equation}
where 
$\boldsymbol {\mathcal {A}}=[\mathcal {A}_1^{+}\ \mathcal {A}_2^{+}\ \mathcal {A}_1^{-} \ \mathcal {A}_2^{-}]^{\mathrm {T}}$ is the amplitude vector, 
$\boldsymbol {{\varPi }}=[\boldsymbol {\rm \pi }_{\boldsymbol {1}}^{\boldsymbol {+}}\ \boldsymbol {\rm \pi }_{\boldsymbol {2}}^{\boldsymbol {+}}\ \boldsymbol {\rm \pi }_{\boldsymbol {1}}^{\boldsymbol {-}}\ \boldsymbol {\rm \pi }_{\boldsymbol {2}}^{\boldsymbol {-}}]$ is the polarization tensor with the normalized vectors 
$\boldsymbol {\rm \pi }_{\boldsymbol {1},\boldsymbol {2}}^{\boldsymbol {\pm }}=[\rm \pi _{\hat {e}_y,1,2}^{\pm } \ {\rm \pi} _{\hat {h}_z,1,2}^{\pm }\ {\rm \pi} _{\hat {e}_z,1,2}^{\pm }\ {\rm \pi} _{\hat {h}_y,1,2}^{\pm }]^{\mathrm {T}}$ and 
$\boldsymbol {{\varPhi }}(x)= \textrm {diag} [\exp (\mathrm {i}k_{x,1}^{+}x)\ \exp (\mathrm {i}k_{x,2}^{+}x)\ \exp (\mathrm {i}k_{x,1}^{-}x)\ \exp (\mathrm {i}k_{x,2}^{-}x)]$ is the phase tensor. In accordance with the above mentioned for 
$\hat {s}_x$, the choice between 
$+$, 
$-$ and 1, 2 is determined by the quantity 
${\rm \pi} _{\hat {e}_y,1,2}^{\pm } ({\rm \pi} _{\hat {h}_z,1,2}^{\pm } )^{*}-{\rm \pi} _{\hat {e}_z,1,2}^{\pm } ({\rm \pi} _{\hat {h}_y,1,2}^{\pm } )^{*}$. The complex amplitude coefficients 
$\mathcal {A}_{1,2}^{\pm }$ are evaluated from the boundary conditions, which, for each transition through a boundary (denoted as 
$j\rightarrow j+1$) have the form 
$\boldsymbol {V}^{j} (x_j,y,z )=\boldsymbol {V}^{j+1} (x_{j+1},y,z)$.
The calculation of the wave's EM field still requires the determination of the plasma dielectric response, in order to provide the elements of tensor 
$\boldsymbol{\mathsf{K}}$ appearing in (2.4) and (2.5). For our studies, we consider a cold plasma (i.e. far from EC resonance) which is divided into three adjacent regions separated by two interfaces: the first and the third region are of the background plasma, while there is a second region in-between which contains turbulent blob-like structures. In the plasma, the background electron density varies according to a given radial profile 
$n_e(r)$, whereas the blob density is specified by the fluctuation amplitude 
$\delta n_e$. For the background plasma regions, the elements of 
$\boldsymbol{\mathsf{K}}$ are derived on the basis of the standard cold-plasma permittivity tensor, after rotational transformations over the propagation angles 
$\phi$ and 
$\theta$ (Brambilla Reference Brambilla1998):
\begin{gather} K_{xx}=\mathcal{S}-(\mathcal{S}-\mathcal{P})\cos^{2}\theta\sin^{2}\phi, \end{gather}
\begin{gather}K_{xy}=K_{yx}=\mp \mathrm{i}\mathcal{D}\cos\phi+(\mathcal{S}-\mathcal{P})\sin\theta\cos\theta\cos^{2}\phi, \end{gather}
\begin{gather}K_{xz}=K_{zx}=\pm \mathrm{i}\mathcal{D}\sin\theta\sin\phi+(\mathcal{S}-\mathcal{P})\cos\theta\sin\phi\cos\phi, \end{gather}
\begin{gather}K_{yy}=\mathcal{S}-(\mathcal{S}-\mathcal{P})\sin^{2}\theta\sin^{2}\phi, \end{gather}
\begin{gather}K_{yz}=K_{zy}=\pm \mathrm{i}\mathcal{D}\cos\theta\sin\phi-(\mathcal{S}-\mathcal{P})\sin\theta\sin\phi\cos\phi, \end{gather}
\begin{gather}K_{zz}=\mathcal{S}^{2}\sin^{2}\phi+\mathcal{P}^{2}\cos^{2}\phi. \end{gather}
In the above, 
$\mathcal {S}$, 
$\mathcal {D}$ and 
$\mathcal {P}$ are the Stix elements for cold magnetized plasma waves, which are known functions of the cyclotron, plasma and wave frequencies.
The intermediate region consists of a mixture of field-aligned density blobs, which are described by an arbitrary permittivity tensor 
$\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{B}}}$. The blobs are located in a cold magnetized plasma with dielectric tensor 
$\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{P}}}$ of the form (2.7), and are elliptical on the direction transverse to 
$\boldsymbol {B}$. The apparent method to tackle the problem is to define 
$\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{B}}}$ based on existing know-how (see § 1). However, the result would provide the dielectric response for only one realization of the blobs in the plasma, whereas many such realizations are required in order to have a sufficient statistical base towards the calculation of the effects on the propagation. To avoid this drawback, the problem is formulated in terms of EM homogenization methods, according to which the blob region is replaced by an equivalent homogenized dielectric region and the plasma is described by an effective permittivity tensor 
$\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{eff}}}$ similar to 
$\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{P}}}$. There are a variety of homogenization formalisms, most of which employ depolarization dyadics arising from associated Green's functions (the reader is pointed to Mackay & Lakhtakia (Reference Mackay and Lakhtakia2015) for a survey).
Concerning our effort, the popular methods have the shortcoming that the implemented depolarization dyadic does not take into account the blob size, assuming that it is much smaller than the wavelength. This issue may be overcome by adapting accordingly the core of the previous methods in order to include the wavelengths and blob sizes relevant to our problem. In this direction, starting from a technique proposed by Sihvola (Reference Sihvola1996), a generalized EM homogenization (gEMH) method has been developed (Bairaktaris Reference Bairaktaris2019). For the derivation, it has been assumed that 
$\boldsymbol {G}(\boldsymbol {r}-\boldsymbol {r'})$ is the dyadic Green's function of the wave equation's differential operator and 
$\boldsymbol {\mathcal {G}}(\boldsymbol {k}- \boldsymbol {k'})$ is its spatial Fourier transform. The integral equation that the EM field must satisfy is
\begin{equation} \boldsymbol{E}\left (\boldsymbol{r}\right )=\boldsymbol{E}^{\boldsymbol{P}}\left (\boldsymbol{r}\right )+\mathrm{i}\int_{V_B} \boldsymbol{G}(\boldsymbol{r}-\boldsymbol{r'})\boldsymbol{\cdot}\left (\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{P}}}-\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{B}}}\right )\boldsymbol{\cdot}\boldsymbol{E}\left (\boldsymbol{r'}\right )\mathrm{d}^{3}r', \end{equation}
with 
$\boldsymbol {E}^{\boldsymbol {P}}$ the solution in the absence of blobs in the plasma and 
$V_B$ the volume occupied by an ellipsoidal blob inclusion. By a Fourier transformation of the above equation, we obtain the useful relation 
$\boldsymbol {\mathcal {E}}= \boldsymbol {\mathcal {E}^{P}}+\textrm {i}\mathcal {U}_{B}\boldsymbol {\mathcal {G}}\boldsymbol {\cdot }(\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{P}}}-\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{B}}}) \boldsymbol {\cdot } \boldsymbol {\mathcal {E}}$, where 
$\boldsymbol {\mathcal {E}}$, 
$\boldsymbol {\mathcal {E}}^{\boldsymbol {P}}$ and 
$\mathcal {U}_{B}$ are the Fourier transforms of 
$\boldsymbol {E}$, 
$\boldsymbol {E}^{\boldsymbol {P}}$ and of the Heaviside step function inside the blob.
Using the Fourier transform 
$\boldsymbol {\mathcal {G}}$, as found above, (2.8) can be handled within the frame of Liouville–Neumann infinite series (Bairaktaris Reference Bairaktaris2019). Ultimately, one is led to a differential equation for the unknown tensor 
$\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{eff}}}$ as follows:
\begin{equation} {\frac{\mathrm{d}\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{eff}}}\left (\nu\right )}{ \mathrm{d}\nu}}\simeq {\frac{\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{B}}}- \boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{eff}}}\left (\nu\right )}{1-\nu}} \left\{1+ \epsilon\boldsymbol{\mathsf{Q}}_{\boldsymbol{\mathsf{1}}}\left [\boldsymbol{N},\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{eff}}} \left(\nu\right )\right ]+ \epsilon^{2}\boldsymbol{\mathsf{Q}}_{\boldsymbol{\mathsf{2}}} \left[\boldsymbol{N},\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{eff}}}\left (\nu\right )\right ]\right\}, \end{equation}
where 
$\nu$ is the fraction of the total blob volume versus the volume of the turbulent region, 
$\boldsymbol{\mathsf{Q}}_{\boldsymbol{\mathsf{1}},\boldsymbol{\mathsf{2}}}$ are functions arising from the Fourier integration in wavevector space, and the scaling parameter 
$\epsilon \approx (L_b/\lambda _\omega ) (|\delta n_e|/n_e) \leqslant 1$ is proportional to the ratio of the blob size over the wavelength multiplied by the relative density contrast of the blobs versus background plasma. Equation (2.9) is further simplified by decomposing the effective tensor into orders of 
$\epsilon$ as 
$\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{eff}}}=\boldsymbol{\mathsf{K}}_{\boldsymbol{\mathsf{0}}}^{ \boldsymbol{\mathsf{eff}}}+\epsilon \boldsymbol{\mathsf{K}}_{\boldsymbol{\mathsf{1}}}^{ \boldsymbol{\mathsf{eff}}}+ \epsilon ^{2}\boldsymbol{\mathsf{K}}_{\boldsymbol{\mathsf{2}}}^{\boldsymbol{\mathsf{eff}}}+\cdots$. In this fashion, separate equations emerge per order of 
$\epsilon$, which can be solved independently. In zero order, an analytical solution is feasible, which yields 
$\boldsymbol{\mathsf{K}}_{\boldsymbol{\mathsf{0}}}^{\boldsymbol{\mathsf{eff}}}= \nu \boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{B}}}+ (1-\nu )\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{P}}}$. The complexity of the higher-order differential equations rises proportionally to their order: the first-order one involves integrals of six variables, the second-order one integrals of nine variables, etc. Solving up to second order, and also using the zero-order solution for the dielectric tensor inside the integrand, is a simplification sufficient for approximating the behaviour of the effective permittivity.
2.2. Ray tracing algorithm
The approach of frequency-domain asymptotic techniques amounts to modelling the wave propagation in terms of rays or beams which are continuously refracted by a slowly space-varying medium, in the same way as the trajectory of a particle is deflected by a scalar potential (Bernstein & Friedland Reference Bernstein and Friedland1980). In this context, our problem is formulated as the propagation of a constant-frequency, linear-amplitude EM wave in a stationary, anisotropic and weakly inhomogeneous magnetized plasma, and the wave physics is described by the vector Helmholtz equation (Brambilla Reference Brambilla1998),
\begin{equation} \boldsymbol{\nabla}\times\left [\boldsymbol{\nabla}\times\boldsymbol{E}(\boldsymbol{r})\right ]-\frac{\omega^{2}}{{c}^{2}}\boldsymbol{\mathsf{K}} (\boldsymbol{r},\omega)\boldsymbol{\cdot}\boldsymbol{E}(\boldsymbol{r})=0. \end{equation}
The conditions implied by the slow spatial variation of the medium are summarized in the relation 
$\kappa =\omega L_p/{c} \gg 1$, where 
$\kappa$ is known as the short-wavelength-limit parameter and 
$L_p$ is the inhomogeneity scale length of the plasma, defined as the maximum value of the relative gradients 
$\boldsymbol {\nabla }\mathcal {X}/\mathcal {X}$ among the parameters 
$\mathcal {X}$ affecting propagation (i.e. entering the dispersion relation). Assuming the validity of this condition, the vector Helmholtz equation allows solutions of a generalized plane-wave ansatz 
$\boldsymbol {E}(\boldsymbol {r})=\boldsymbol {A}(\boldsymbol {r})\exp [\textrm {i}\kappa \sigma (\boldsymbol {r}) ]$, where the phase 
$\boldsymbol {k}\boldsymbol {\cdot }\boldsymbol {r}$ has been replaced by an eikonal function defined as 
$\boldsymbol {k}=\kappa \boldsymbol {\nabla }\sigma$.
For each ray, one can determine the ‘backbone’ of the wave field by means of a set of ordinary differential equations that give the variation of the wave phase and the electric field amplitude along the ray. This is achieved with the expansion of the amplitude in an asymptotic series over the large parameter 
$\kappa$ (
$\boldsymbol {A}= \boldsymbol {A}_{\boldsymbol {0}}+\kappa ^{-1} \boldsymbol {A}_{\boldsymbol {1}}+\kappa ^{-2}\boldsymbol {A}_{\boldsymbol {2}}+\cdots$), followed by the insertion of the amplitude series in (2.10), the separation of terms of different order in different equations, and the solution of these equations. In this process, an additional assumption is that the effect of the resonant absorption on the wave dispersion is not significant; this is quantified by setting the antihermitian part of the dielectric tensor to the next asymptotic order with respect to its Hermitian part, thus implying the relation 
$\boldsymbol{\mathsf{K}}=\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{H}}}+ \textrm {i}\kappa ^{-1}\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{A}}}$. The sequence reveals that in zero-order one follows the ray position and wavenumber, and in first order the amplitude transport, whereas the inclusion of higher-order equations partly describes wave interference and diffraction (for details the reader is pointed to Tsironis (Reference Tsironis2013)).
To zero order, the asymptotic expansion terms build up the following relation:
\begin{equation} \left [\frac{{c}^{2}}{\omega^{2}}\left (-k^{2}\boldsymbol{I}+\boldsymbol{k}\boldsymbol{k}\right )+\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{H}}}\right ]\boldsymbol{\cdot} \boldsymbol{A}_{\boldsymbol{0}}=\boldsymbol{\varLambda}\boldsymbol{\cdot} \boldsymbol{A}_{\boldsymbol{0}}=0, \end{equation}
where 
$\boldsymbol {I}$ is the unitary tensor and 
$\boldsymbol {k}\boldsymbol {k}$ is a dual product. The solvability condition of (2.11) reads 
$\det [\boldsymbol {\varLambda }(\boldsymbol {r},\omega )]=0$, and apart from being the EC wave dispersion relation it is also a Hamilton–Jacobi equation for the eikonal function. The corresponding solution yields Hamiltonian equations for 
$\boldsymbol {r}$ and 
$\boldsymbol {k}$ of each ray,
\begin{equation} \frac{\mathrm{d}\boldsymbol{r}}{\mathrm{d}\tau}=\frac{\partial H}{\partial\boldsymbol{k}}, \quad \frac{\mathrm{d}\boldsymbol{k}}{\mathrm{d}\tau}=-\frac{\partial H}{\partial\boldsymbol{r}}, \end{equation}
with 
$H=\det (\boldsymbol {\varLambda })$ playing the role of the Hamiltonian function, and 
$\tau$ being a normalized coordinate measuring the length along propagation. Going now to the first-order terms, these provide the following equation for the wave amplitude damping:
\begin{equation} \frac{\mathrm{d}|\boldsymbol{A}_{\boldsymbol{0}}|^{2}}{\mathrm{d}\tau}+\left (\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{v}_{\boldsymbol{g}}+ 2\hat{\boldsymbol{e}^{*}}\boldsymbol{\cdot} \boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{A}}} \boldsymbol{\cdot}\hat{\boldsymbol{e}}\right)|\boldsymbol{A}_{\boldsymbol{0}}|^{2}=0. \end{equation}
where 
$\boldsymbol {v}_{\boldsymbol {g}}=\partial \omega /\partial {\boldsymbol {k}}$ is the group velocity and the other addend stands for the damping coefficient. Equation (2.13) implies that the wave energy propagates in the direction of the group velocity and that absorption is proportional to the projection of the dielectric tensor's antihermitian part onto the polarization vector. It may be transformed to an equation for the absorption of wave power 
$P$, reading 
$\mathrm {d}P/\mathrm {d}\tau =-2\,\textrm {Im} (\boldsymbol {k})\boldsymbol {\cdot } \boldsymbol {v}_{\boldsymbol {g}}P$. The driven current is then given by 
$I_{\textrm {CD}}=\eta _{\textrm {CD}}P_{\textrm {abs}}/2{\rm \pi} R_{\textrm {maj}}$, with 
$P_{\textrm {abs}}$ the absorbed power (computed from the previous equation), 
$R_{\textrm {maj}}$ the tokamak major radius, and the current drive efficiency 
$\eta _{\textrm {CD}}$ is calculated by solving the associated drift kinetic problem, including the effects of trapped particles and ion–electron collisions (see Prater (Reference Prater2004) for details).
As seen from (2.11) and (2.13), ray tracing requires knowledge of the Hermitian part of the dielectric tensor, the group velocity and the imaginary part of the wavenumber. Adopting cold plasma propagation, and given radial profiles for the magnetic field and the plasma density and temperature, for 
$\boldsymbol{\mathsf{K}}^{\boldsymbol{\mathsf{H}}}$ we make use of the permittivity tensor of (2.7). Then, 
$\boldsymbol {v_g}$ and 
$\textrm {Im} (\boldsymbol {k})$ may be found from the real and imaginary part of the dispersion relation, respectively. In this context, the modifications to the system geometry due to the magnetic islands are added from the viewpoint of the associated effects to wave propagation. The first modification amounts to a poloidal perturbation in the magnetic field strength, which may be added to the equilibrium field within its Clebsch representation in the toroidal and poloidal angle coordinates 
$\varphi$ and 
$\vartheta$ (Morrison Reference Morrison2006),
\begin{equation} \boldsymbol{B}=\boldsymbol{\nabla}\varPsi_\varphi\boldsymbol{\cdot}\boldsymbol{\nabla}\vartheta- \boldsymbol{\nabla}\varPsi_\vartheta\boldsymbol{\cdot}\boldsymbol{\nabla}\varphi. \end{equation}
The toroidal flux is given by 
$\partial \varPsi _\varphi /\partial r=rB_\varphi$ and the poloidal flux is perturbed on the rational surfaces 
$q(r_s^{m,n})=m/n$, with the zero-order term given by 
$\partial \varPsi _{\vartheta 0}/\partial r=RB_\vartheta$ and the perturbation term having the form (Isliker et al. Reference Isliker, Chatziantonaki, Tsironis and Vlahos2012)
\begin{equation} \varPsi_{\vartheta 1}\left (r,\vartheta,\varphi\right )=-\sum_{m,n}\frac{1}{m}\delta B_\vartheta^{m,n}r(r-r_s^{m,n}) \cos (m\vartheta-n\varphi). \end{equation}
The second modification is the fixation of the plasma pressure inside the island, which, according to the equation of state, corresponds to a flattening of the electron density and temperature profiles within the island region. Denoting 
$n_e$, 
$T_e$ as 
$\mathcal {X}$, their profiles outside the island as 
$\mathcal {X}_E(r)$ and the island half-width as 
$W$, the radial profiles take the form
\begin{equation} \mathcal{X}(r)=\mathcal{X}_E(r)\left [{U}\left (r_s-W-r\right )+{U}\left (r-r_s-W\right )\right ]+ \mathcal{X}_E(r_s){U}\left (W-|r-r_s|\right ). \end{equation}2.3. MRE
The evolution of the magnetic islands resulting from NTM instability growth may be calculated with the MRE model, which is based on a nonlinear equation for the island's width. This equation stems by applying a modified Ohm's law on the plasma circuit around the mode's flux surface, and has the form (Hegna & Callen Reference Hegna and Callen1997)
\begin{equation} \frac{T_r}{r_s}\frac{\mathrm{d}W}{\mathrm{d}t}=r_s\left (\varDelta'_{N}+\varDelta'_{\textrm{BS}}+\varDelta'_{\textrm{CD}}\right ). \end{equation}
In the above, 
$T_r=0.82{\mu _0}r_s^{2}/\varrho _p$ is the resistive time scale for the mode growth, which is proportional to the square of the radial coordinate 
$r_s$ of the island's 
$O$-point and inversely proportional to the plasma resistivity 
$\varrho _p$ (
$\mu _0$ is the vacuum magnetic permeability), and the 
$\varDelta '$ indexes determine altogether the NTM stability: if their sum is positive then the mode is unstable, whereas if their sum is negative then the mode is stable. There are a variety of physics effects impinging on NTM dynamics for which a 
$\varDelta '$-term has been matched (the reader is pointed to Urso (Reference Urso2009) for a survey). In our description, we include the three most important indexes regarding ECCD-driven modes: neoclassical stability (
$\varDelta '_{N}$), bootstrap current (
$\varDelta '_{\textrm {BC}}$) and external current drive (
$\varDelta '_{\textrm {CD}}$).
The neoclassical stability index incorporates the (neoclassical) nonlinear saturation effect to the (classical) linear tearing-mode stability
\begin{equation} \varDelta'_{N}(W)=-\frac{m}{r_s}-\frac{W}{2.44r_{\min}^{2}}, \end{equation}and the bootstrap current index includes the effect of the associated helical current
\begin{equation} \varDelta'_{\textrm{BS}}(W)=\frac{\sqrt{\varepsilon_A}\beta_P}{W}\frac{L_q}{L_p}\left (\frac{W^{2}}{W^{2}+W_{\textrm{cs}}^{2}}+\frac{W^{2}}{W^{2}+28W_b^{2}}- \frac{W_{\textrm{pol}}^{2}}{W^{2}}\right ), \end{equation}
where 
$\varepsilon _A=r_{\min }/R_{\textrm {maj}}$ is the tokamak aspect ratio (
$r_{\min }$ the minor radius), 
$\beta _P$ is the plasma beta, 
$L_q$ and 
$L_p$ the shear lengths of the safety factor and the plasma pressure and 
$W_{\textrm {cs}}$, 
$W_b$, 
$W_{\textrm {pol}}$ are threshold values for the island width: 
$W_d$ is the critical width for classical destabilization, 
$W_b$ the width below which banana orbits weighs in to the bootstrap current and 
$W_{\textrm {pol}}$ the width below which the current response to diamagnetic island rotation is significant; a detailed physics analysis may be found in LaHaye (Reference LaHaye2006) and references therein. Finally, the current drive index describes the effect of ECCD,
\begin{equation} \varDelta'_{\textrm{CD}}(W)=-{U}\left (t-t_{\textrm{CD}}^{\textrm{on}}\right )\varpi_{\textrm{CD}} \frac{16\mu_0}{\sqrt{\rm \pi}}\frac{r_sL_q}{B_\vartheta} \frac{d_{\textrm{CD}}j_{\textrm{CD}0}}{W^{2}} \exp\left (-\frac{\delta r_{\textrm{CD}}^{2}}{d_{\textrm{CD}}^{2}}\right ), \end{equation}
with the EC turned on at 
$t=t^{\textrm {on}}_{\textrm {CD}}$, 
$j_{\textrm {CD}0}$ the current density per unit area, 
$d_{\textrm {CD}}$ the deposition width, 
$\delta r_{\textrm {CD}}=r_{\textrm {CD}}-r_s$ the radial misalignment of the deposition with respect to the 
$O$-point, and the optimization factor 
$\varpi _{\textrm {CD}}$ (there we denote 
$u_{\textrm {CD}}=W/d_{\textrm {CD}}$) is
\begin{align} \varpi_{\textrm{CD}}=0.07u_{\textrm{CD}}^{2}+\left (0.34- 0.07u_{\textrm{CD}}^{2}\right )\left [0.3u_{\textrm{CD}}{U} \left (2 - u_{\textrm{CD}}\right ) +\exp\left (-u_{\textrm{CD}}\right ){U} \left (u_{\textrm{CD}} - 2\right ) \right ]. \end{align} The MRE reproduces well the main aspects of the NTM dynamics. First, by looking at the inequality 
$\mathrm {d}W/\mathrm {d}t \geqslant 0$, one sees that a ‘seed’ island (matched to a minimum initial width) is required for the mode onset. The bootstrap current term is destabilizing for normal tokamak operation, since to have 
$\varDelta '_{\textrm {BS}}<0$ requires reversed shear. On the contrary, ECCD appears a net stabilizing effect provided that 
$\delta r_{\textrm {CD}}\ll W\ll 3d_{\textrm {CD}}$. Finally, an unstable mode that goes beyond the limit 
$W_d=r_{\min }-r_s$ leads to disruption, since the island size exceeds the tokamak boundaries and, as a consequence, particles trapped in the island may follow trajectories that intersect the tokamak wall.
3. Set-up of the numerical simulations
The scenario simulated here is the stabilization of low-order NTMs in ITER-like plasma and in the presence of edge density fluctuations, using an EC wave beam for generating current drive near the 
$O$-point of the associated magnetic island. The beam is launched into the plasma at an angle to the toroidal magnetic field with a Gaussian electric field amplitude profile and, before reaching the island, encounters the electron density perturbations in the form of curvilinear-shaped blobs. The modified beam is then absorbed at the EC resonance layer, which is aligned to the island's flux surfaces, adding up its stabilizing effect to the mode evolution. In figure 1, an indicative case for the geometry of the processes composing the aforementioned scenario is visualized. Below we define the parameters for the EC wave, the tokamak device (magnetic field and electron plasma), the blob structures and the NTM instability, and also describe the linking of the theoretical tools presented in § 2 while building the overall modelling chain.
Figure 1. Example geometry for the problem of NTM control with ECCD in tokamak plasma: (a) three-dimensional view of the magnetic field configuration, (b) poloidal projection of the beam propagation and the plasma equilibrium in the presence of magnetic islands and density blobs.
The physics parameters of the background plasma and the EC wave are relevant to the ITER design (Holtkamp Reference Holtkamp2009). The magnetic equilibrium has a circular cross-section and its magnetic field is represented as 
$\boldsymbol {B}=B_\varphi \hat {\boldsymbol {\varphi }}+B_\vartheta \hat {\boldsymbol {\vartheta }}$, where the toroidal component is axisymmetric and the poloidal component is defined in terms of the safety factor
\begin{equation} B_\varphi(r,\vartheta)=\frac{R_{\textrm{maj}}}{R_{\textrm{maj}}+r\cos\vartheta}B_0, \quad B_\vartheta(r,\vartheta)=\frac{r}{R_{\textrm{maj}}}\frac{B_\varphi(r,\vartheta)}{q(r)}. \end{equation}The radial profiles of the electron density and temperature, as well as the safety factor profile, are parabolic functions of the form
\begin{equation} \frac{\mathcal{X}_E(r)-\mathcal{X}_E(0)}{\mathcal{X}_E(r_{\min})-\mathcal{X}_E(0)}=\frac{r^{2}}{r_{\min}^{2}}. \end{equation}
The parameter values which are required in order to configure (3.1) and (3.2) are the following: 
$R_{\textrm {maj}}=6.2$ m, 
$r_{\min }=1.9\ \textrm {m}$, 
$B_0=5.1\ \textrm {T}$, 
$q(0)=1$ and 
$q(r_{\min })=4$, 
$n_e(0)=10^{20}\ \textrm {m}^{-3}$ and 
$n_e(r_{\min })=10^{19}\ \textrm {m}^{-3}$, 
$T_e(0)=10\ \textrm {KeV}$ and 
$T_e(r_{\min })=1\ \textrm {KeV}$. As far as the wave launching conditions are concerned, the frequency is 
$\omega /2{\rm \pi} =170\ \textrm {GHz}$, which corresponds to the first EC harmonic for the specific value of 
$B_0$, and the mode of electric field polarization is ordinary (O). The Gaussian beam is assumed to be injected from the point 
$(r_0,\vartheta _0,\varphi _0)=(r_{\min },0,0)$ at angles 
$\theta _0=0^{\circ }$ in the poloidal direction and 
$\phi _0=-5^{\circ }$ toroidally, with power 
$P_0=1\ \textrm {MW}$ and half-width 
$w_0=3\ \textrm {cm}$. The beam is launched at its waist position, orthogonally and without initial focusing, which practically corresponds to a wavefront curvature radius much larger than 
$R_{\textrm {maj}}$.
Regarding the MHD instability properties, here we consider NTMs of order equal to 2/1 and 3/2. These modes, according to the 
$q$-profile defined above, appear and develop on the flux surfaces 
$r_s^{2,1}=0.58r_{\min }$ and 
$r_s^{3,2}=0.4r_{\min }$, respectively. The plasma resistivity is set equal to 
$\varrho _p=1.1\times 10^{-9}\ \Omega \ \textrm {m}$, which, taking into account the specific value of 
$r_s$ for each mode, yields 
$T_r^{2,1}\approx 1251\ \textrm {s}$ and 
$T_r^{3,2}\approx 595\ \textrm {s}$ for the growth time scale in each case. The plasma beta is set to 
$\beta _p=0.6$, and the ratio of the pressure and safety factor shear lengths is assumed to be 
$L_p/L_q=1$, where the value of 
$L_q$ per mode is computed from its definition 
$L_q=q(r_s)/[\mathrm {d}q(r_s)/ \mathrm {d}r]$ on the basis of the 
$q$-profile. The last parameters to be specified are the characteristic threshold values defined along (2.19), which are 
$W_{\textrm {cs}}=0.01r_s$, 
$W_b=0.02r_s$ and 
$W_{\textrm {pol}}=0.015r_s$. As initial conditions, we choose that the island is detected at 
$t=0$ with a size equal to 6 cm, and that the EC wave source is turned on immediately after detection (
$t^{\textrm {on}}_{\textrm {CD}}=0$). For the synchronization of the current drive with the island rotation, in our model we have assumed that the EC power source is modulated very close to the island rotation frequency and the power-on phase is centred around the 
$O$-point passage through the beam.
Finally, we specify the parameters related to the electron density fluctuations. These appear close to the tokamak edge in the form of ellipsoidal blobs, and occupy a thin plasma layer containing a few flux surfaces (see figure 1b). The blob layer is centred around normalized radius 
$\rho _b=r_b/r_{\min }=0.8$, and its full-width 
$2\delta \rho _b$ is equal to 0.04. We note here that the restriction 
$\epsilon \leqslant 1$ imposes limitations to the selection of the values for 
$L_b$ and 
$\delta n_e$; e.g. since for frequency 170 GHz the wavelength is around 2 mm, the blob size cannot exceed 4 mm for density contrast 50 % or 2 mm for density contrast 100 %. In this frame, we assume for each blob a characteristic length 
$L_b=0.2\ \textrm {cm}$ and ellipticity factors 
$\alpha _{b1}$, 
$\alpha _{b2}$ that distribute uniformly within 0.8 and 1.2; in terms of these, the single blob volume may be calculated from the formula 
$V_B=4{\rm \pi} \alpha _{b1}\alpha _{b2}L_b^{3}/3$. The inhomogeneous plasma content of the blob region (as described up to this point) is replaced with an equivalent homogeneous plasma by applying the gEMH method shown in § 2.1; for the values of the parameters involved, the density perturbation 
$\delta n_e$ varies from 50 % to 100 % of the local background density, and that the blob filling ratio 
$\nu$ is fixed to 0.3.
The simulation procedure constitutes of three parts in sequential order, summarized in figure 2. In the first part (propagation through density fluctuations), we provide at input the launched beam profile (in the form (2.2)) and the plasma configuration, and then a MATLAB algorithm decomposes the beam into plane wave modes, computes the effective permittivity tensor of the turbulent region (solving (2.9)) and applies the transfer matrix for advancing each mode in space. The initial refraction index for each mode is found via the dispersion relation (solvability condition of (2.11)), and then decomposed per propagation direction according to the angles 
$\theta _0$ and 
$\phi _0$. Then, one can solve the first eigenvalue problem in (2.6) and determine the polarization in the region of incidence, as well as the four eigenvalues for the 
$z$-component. For the incident wave with 
$O$-mode polarization, one eigenvalue coincides with the one of the four solutions above. Then, the same values for the perpendicular refraction indexes, which are preserved within the boundary conditions, apply to the next problem that determines the polarizations and the associated eigenvalues for the refraction index along 
$x$. The algorithm continues the propagation up to the leftmost boundary of the turbulent layer, located at 
$\rho _b-\delta \rho _b$; in the output, the Cartesian components of the beam field are provided at each point.
Figure 2. Block diagram showing the various stages of the simulation approach.
The second part involves the beam trajectory from the stopping point of the full-wave calculation to the magnetic island. The computation is performed with a FORTRAN 95 code which solves the Hamiltonian equations (2.12) and the power damping law that stems from (2.13). As input, apart from the magnetic equilibrium and the plasma profiles (including NTM islands), the code requires initial conditions per ray for the position, the wavevector and the power fraction. Within our approach, the beam is described as the synthesis of three different rays: one ‘central’ ray, which represents the beam axis and along which the power absorption is computed, and two ‘peripheral’ rays for the description of the beam width in each direction transverse to the propagation. By processing the output from the preceding computation, all the required initial values may be expressed in terms of the local Poynting vector and the moments of the electric field distribution. After the initial condition assignment, the code propagates the rays and computes the transmitted power and ECCD along their trajectories. The computation terminates on the high-field side of the island's outer flux surface (located at 
$r=r_s-W$), and yields the radial profiles of the absorbed power and driven current.
The sequence for defining the initial conditions is the following: the initial position of the central ray 
$\boldsymbol {r}_{\boldsymbol {1}}^{\boldsymbol {c}}$ coincides with the maximum value of 
$E$ on the surface 
$r=r_b-\delta r_b$; then, the wavenumber 
$k_1^{c}$ is found from the dispersion relation at this position, and the direction of propagation 
$\hat {\boldsymbol {k}}_{\boldsymbol {1}}^{\boldsymbol {c}}$ emanates from Poynting's vector and yields the toroidal and poloidal propagation angles 
$\phi _1^{c}=\sin ^{-1}(\hat {\boldsymbol {k}}_{\boldsymbol {1}}^{\boldsymbol {c}}\boldsymbol {\cdot }\hat {\boldsymbol {y}})$ and 
$\theta _1^{c}=\cos ^{-1} (\hat {\boldsymbol {k}}_{\boldsymbol {1}}^{\boldsymbol {c}} \boldsymbol {\cdot }\hat {\boldsymbol {x}}/\cos \phi _1^{c})$. Next, the moments of 
$E^{2}$ are defined in a Cartesian coordinate system (
$\chi ,\psi ,\zeta$) attached to the beam, where, at each propagation point, the wavevector lies on the 
$\chi$-axis and the field profile spans over the 
$\psi$–
$\zeta$ plane (Tsironis & Papadopoulos Reference Tsironis and Papadopoulos2014),
\begin{equation} \langle\psi^{j_\psi}\zeta^{j_\zeta}\rangle(\chi)=\frac{\displaystyle{\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \psi^{j_\psi}\zeta^{j_\zeta}E^{2}(\chi,\psi,\zeta)\,\mathrm{d}\psi\,\mathrm{d}\zeta}}{ \displaystyle{\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} E^{2}(\chi,\psi,\zeta)\,\mathrm{d}\psi\,\mathrm{d}\zeta}}. \end{equation}
At 
$\boldsymbol {r}=\boldsymbol {r}_{\boldsymbol {1}}^{\boldsymbol {c}}$, the first-order moments give the transverse coordinates of the beam centre as 
$\psi _1^{c}=\langle \psi ^{1}\zeta ^{0}\rangle (\chi _1^{c})$ and 
$\zeta _1^{c}=\langle \psi ^{0}\zeta ^{1}\rangle (\chi _1^{c})$, while the second-order moments provide the associated beam widths as 
$w_{\psi ,1}^{2}=\langle \psi ^{2}\zeta ^{0}\rangle -(\psi _1^{c})^{2}$ and 
$w_{\zeta ,1}^{2}= \langle \psi ^{0}\zeta ^{2}\rangle -(\zeta _1^{c})^{2}$. With this information available, one may define the initial positions for the peripheral rays as 
$\boldsymbol {r}_{\boldsymbol {1}}^{\boldsymbol {fp}}=\boldsymbol {r}_{\boldsymbol {1}}^{\boldsymbol {c}}+w_{\psi ,1}\hat {\boldsymbol {\psi }}$ and 
$\boldsymbol {r}_{\boldsymbol {1}}^{\boldsymbol {sp}}= \boldsymbol {r}_{\boldsymbol {1}}^{\boldsymbol {c}}+w_{\zeta ,1}\hat {\boldsymbol {\zeta }}$. The axes of the local coordinate system are mapped to the ones of the global Cartesian system using the transformation relation
\begin{equation} \begin{bmatrix} \hat{\boldsymbol{\chi}} \\ \hat{\boldsymbol{\psi}} \\ \hat{\boldsymbol{\zeta}} \end{bmatrix}= \begin{bmatrix} \cos\phi_1^{c} & 0 & -\sin\phi_1^{c} \\ 0 & 1 & 0 \\ \sin\phi_1^{c} & 0 & \cos\phi_1^{c} \end{bmatrix} \cdot \begin{bmatrix} \cos\theta_1^{c} & -\sin\theta_1^{c} & 0 \\ \sin\theta_1^{c} & \cos\theta_1^{c} & 0 \\ 0 & 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} \hat{\boldsymbol{x}} \\ \hat{\boldsymbol{y}} \\ \hat{\boldsymbol{z}} \end{bmatrix}, \end{equation}
which emerges from the fact that these two coordinate systems may be aligned after two successive rotations (one per propagation angle). The following step is to define the initial wavevectors of the peripheral rays; along propagation, these match up to the wavevector on the beam axis, but, in the transverse directions, they also have components with their signs reflecting the beam focusing status (
$-$ for focused, 
$+$ for defocused). In this sense, the relations providing the associated refraction indexes are 
$\boldsymbol {N}_{\boldsymbol {1}}^{\boldsymbol {fp}}=\boldsymbol {N}_{\boldsymbol {1}}^{\boldsymbol {c}}+ w_{\psi ,1}\mathcal {R}_{\psi ,1}^{-1}\hat {\boldsymbol {\psi }}$ and 
$\boldsymbol {N}_{\boldsymbol {1}}^{\boldsymbol {sp}}=\boldsymbol {N}_{\boldsymbol {1}}^{\boldsymbol {c}}+w_{\zeta ,1}\mathcal {R}_{\zeta ,1}^{-1}\hat {\boldsymbol {\zeta }}$, where the curvature radii are given by 
$w_{\psi ,1}^{2}\mathcal {R}_\psi ^{-1}=\langle \psi ^{1}\zeta ^{0}\rangle \digamma (\langle \psi ^{1}\zeta ^{0}\rangle )$ and 
$w_{\zeta ,1}^{2}\mathcal {R}_\zeta ^{-1}=\langle \psi ^{0}\zeta ^{1}\rangle \digamma (\langle \psi ^{0} \zeta ^{1}\rangle )$; here, 
$\digamma$ denotes that the moments are calculated in Fourier space for the transformed electric field. Finally, the beam power is found from Poynting's vector (referring to the beam spot area size) and assigned solely to the central ray; the physics formula reads 
$P_1={\rm \pi} w_{\psi ,1}w_{\zeta ,1}S$.
In the third part, we compute the evolution of the magnetic island growth as influenced by the ECCD. We use a MATLAB script for solving (2.17), which gets as input the equilibrium, profiles and all the other parameters appearing in (2.18), (2.19) and (2.20) (see the third paragraph in this section), and gives as output the time series of the island width. Since each equilibrium geometry corresponds to a certain size of the island, the ray computation must be performed for every value of 
$W$ yielded by the MRE. This feedback-loop situation can be simplified by employing the linearity of current drive: we run the ray tracing code for different islands and apply regression analysis to derive a relation connecting the ECCD parameters with 
$W$. The data related to current drive, as delivered from the ray tracing code, require further processing in order to make up consistent input for the MRE solver. First, for obtaining the proper numbers for 
$r_{\textrm {CD}}$, 
$j_{\textrm {CD}0}$ and 
$d_{\textrm {CD}}$ from the ECCD profile, the driven current should be translated to area current density. Instead of dividing the values of 
$I_{\textrm {CD}}$ with their equivalent values of flux surface area, we perform this more easily with the relation 
$j_{\textrm {CD}}=\eta _{\textrm {CD}}\,\mathrm {d}P_{\textrm {abs}}/\mathrm {d}V_s$. The flux surface volumes are calculated from (Isliker et al. Reference Isliker, Chatziantonaki, Tsironis and Vlahos2012)
\begin{equation} V_s=-\frac{1}{n}\int_0^{2{\rm \pi}n}\int_{\xi_1}^{\xi_2}\int_{r_1}^{r_2}\left (R_{\textrm{maj}}+r\cos\vartheta\right )r\,\mathrm{d}r\,\mathrm{d}\xi \,\mathrm{d}\vartheta, \end{equation}
where 
$\xi =m\vartheta -n\varphi$ is the helical angle, and the integration limits for 
$\xi$ and 
$r$ are defined on the basis of the following normalized parameter:
\begin{equation} \varOmega (r,\xi)=-\frac{n}{2}\frac{B_0}{\delta B_\theta^{m,n}}\frac{1}{R_{\textrm{maj}}}\left (\frac{1}{r_s}-\frac{1}{L_q}\right )\left (r-r_s\right )^{2}+\frac{r}{r_s}\left (1+\frac{r-r_s}{r_{\min}}\right )\cos\xi, \end{equation}
that labels the flux surfaces as follows: 
$|\varOmega |>1$ yields surfaces outside the island, whereas 
$|\varOmega | \leqslant 1$ refers to the surface nest of the island (
$\varOmega =1$ corresponds to the separatrix).
4. Results and discussion
This section is devoted to the numerical results of each physics model in our sequence (transfer matrix in gEMH-configured plasma, ray tracing in non-axisymmetric geometry, MRE solution with ECCD term), focusing on the assessment of the effect of the beam distortion by the density blobs to the NTM stabilization quality. To this end, the EC beam spectrum that propagates through the density blobs, as well as the absorbed power and driven current, are compared with the ones of a beam computed in unperturbed plasma with ray tracing. The important quantities to be tracked are the propagation direction, the beam width, the EC deposition radius (especially with respect to the NTM surfaces), the peak value and radial extent of the ECCD profile as well as the magnetic island width.
We start with results coming from the full-wave propagation data and its post-processing. In figure 3 we visualize the beam path in its course of crossing the region populated by density fluctuations (contained within the grey planes), in terms of plotting the spatial distribution of the real part of Poynting vector's magnitude. In the four-dimensional plot, the Cartesian axes refer to the local beam coordinate system (
$\chi ,\psi ,\zeta$) introduced in § 3, but with an additional shift of the origin up to the geometric centre of the blob layer, whereas 
$\textrm {Re}(S)$ is represented by a colour code and normalized to its maximum value on the incident wave. The numerical computations were performed with the dedicated model described previously, and for the parameter values of the EC wave, background plasma and density perturbations mentioned there. The following cases are considered for 
$\delta n_e$ as percentage of the local background density: (a) 50 % and (b) 100 %.
Figure 3. Visualization of the EC beam propagation through the homogenized plasma region containing electron density fluctuations, in terms of the isosurface of the real part of the normalized Poynting vector. System parameters correspond to the case 
$\phi _0=-5^{\circ }$, 
$\rho _b=0.8$, 
$\delta \rho _b= 0.02$, 
$\nu = 0.3$ and 
$\delta n_e$ equal to (a) 50 % and (b) 100 % of the local background density.
By following the shape and trajectory of the beam prior to and past the perturbed region, it is observed that the wave characteristics are altered by the density fluctuations. Specifically, the beam distortion mainly amounts to a shift of the geometrical axis centre, an effect owed to refraction, and to the broadening of the transverse electric field profile, which is owed to diffraction and being quantified in terms of the beam width. Such behaviour has been recently computed in similar scenarios by many authors (see e.g. Ram & Hizanidis (Reference Ram and Hizanidis2016), Köhn et al. (Reference Köhn, Holzhauer, Leddy, Thomas and Vann2016) and Snicker et al. (Reference Snicker, Poli, Maj, Guidi, Koehn, Weber, Conway, Henderson and Saibene2017)). The axial shift and the broadening of the beam are computed using the values of the beam axis coordinates and the width at the entry and exit points on the perturbed layer, which may be found by applying the relevant equations given in § 3 at these specific points. Regarding the beam broadening definition, we should mention here that we have defined the broadening as the difference between the beam width when propagating in the presence of perturbations from the width of unperturbed plasma propagation. In case (a), the axial shift is found equal to 0.14 cm and the beam broadening is 0.09 cm, whereas in case (b) these numbers are 0.33 cm and 0.21 cm, respectively. In both cases, the computed reflection of the wave power (visualized with the dark curve in figure 3) remains below 2 %, and thus it does not constitute an important effect.
We move now to the ray tracing analysis and calculations that follow towards the quantification of NTM dynamics. As described in § 3, the EC beam is modelled via the synthesis of three rays, which are input into the code in order to compute each ray's position, the beam width per transverse direction and the current drive profile. For including the island geometry to the ECCD parameters that model the stabilizing effect in the MRE, it is necessary to derive relations connecting those parameters with the island size. This has been performed as described in Chatziantonaki et al. (Reference Chatziantonaki, Tsironis, Isliker and Vlahos2013); in each case for the NTM order and the blob amplitude, we have run ray tracing simulations for different island sizes and, assuming the Gaussian form 
$j_{\textrm {CD}}(r)=j_{\textrm {CD}0}\exp {[-(r-r_{\textrm {CD}})^{2}]/d_{\textrm {CD}}^{2}]}$ for the ECCD density radial profile, we have applied regression analysis on the output data in order to calculate 
$j_{\textrm {CD}0}$ and 
$d_{\textrm {CD}}$. The results are given in figure 4: in panels (a) and (b) we visualize the poloidal projection of the beam and the ECCD profile indicatively for the case of a 3/2 NTM with 
$\delta n_e/n_e(r_b)=1$, whereas in panels (c) and (d) we present the calculated functions 
$j_{\textrm {CD}0}(W)$ and 
$d_{\textrm {CD}}(W)$ for all cases of 
$m/n$ and 
$\delta n_e$.
Figure 4. Ray tracing analysis of the EC propagation past the density fluctuations and up to the magnetic islands: (a) poloidal projection of the beam path and (b) radial ECCD density profile, in example for the case of a 3/2 NTM, and the calculated regression of (c) 
$j_{\textrm {CD}0}$ and (d) 
$d_{\textrm {CD}}$ versus the island width in all cases for the NTM order.
In figure 4(b), the radial profile of the ECCD is visualized with and without density fluctuations present, so that comparisons are made regarding the effect of the blob structure on the wave deposition. The latter result is based on the additional ray computation mentioned above; this was performed by launching a beam similar to the one considered in the full-wave computations but modelling it using ray tracing, since, in the absence of density fluctuations, the asymptotic approximation is valid. By examining the profiles, the deposition radius for each case is found: for 
$m/n=3/2$, if 
$\delta n_e/n_e(r_b)=0.5$ it is 
$r_{\textrm {CD}}=0.3985r_{\min }$ and if 
$\delta n_e/n_e(r_b)=1$ it is 
$r_{\textrm {CD}}=0.3972r_{\min }$, and for 
$m/n=2/1$, if 
$\delta n_e/n_e(r_b)=0.5$ it is 
$r_{\textrm {CD}}=0.5029r_{\min }$ and if 
$\delta n_e/n_e(r_b)=1$ it is 
$r_{\textrm {CD}}=0.5037r_{\min }$. Under the set-up that, in the absence of blobs, the deposition occurs exactly on the flux surface of interest, these values correspond, respectively, to misalignments of 0.29 cm and 0.53 cm for the cases of 
$\delta n_e/n_e(r_b)$ equal to 0.5 and 1 when 
$m/n=3/2$, and 0.55 cm and 1.08 cm for the same cases when 
$m/n=2/1$.
From previous works, it is well known that the dependence of the EC current density amplitude and radial profile width on the NTM-driven magnetic island size stems from the fact that, in the presence of the island, the beam power is deposited into volumes increasingly smaller than those in its absence, which in turn leads to increasingly larger values of the absorbed power density and the driven current (the reader is referred to Isliker et al. (Reference Isliker, Chatziantonaki, Tsironis and Vlahos2012) and Chatziantonaki et al. (Reference Chatziantonaki, Tsironis, Isliker and Vlahos2013) for details on these issues). The dependence of 
$d_{\textrm {CD}}$ on 
$W$ is much weaker than the one of 
$j_{\textrm {CD}0}$, however, it has been included in the modelling in order to ensure consistency. For practical reasons, we have introduced appropriate functional forms 
$j_{\textrm {CD}0}=f_j(W)$ and 
$d_{\textrm {CD}}=f_d(W)$ for the self-consistent determination of the ECCD parameters, which are quantified by fitting over the exact values computed with ray tracing. Since the ECCD computation is done in terms of the linear adjoint method, it is expected that the scaling of 
$j_{\textrm {CD}0}$ and 
$d_{\textrm {CD}}$ with 
$W$ will be, to a good approximation, linear. In figures 4(c) and 4(d), we plot the functions 
$j_{\textrm {CD}0}(W)$ and 
$d_{\textrm {CD}}(W)$, respectively, for all the four cases that emerge from the combinations of the selected values of 
$m/n$ and 
$\delta n_e/n_e(r_b)$, including also (for comparisons) the computations where the blobs are absent.
As a first comment, regarding the functional dependencies, we can clearly see that the result of linear fitting is indeed very successful in all cases, with the regression error for 
$f_j(W)$ being less than 
$1\,\%$ and for 
$f_d(W)$ less than 
$3\,\%$. By comparing the curves for the different cases, for both 
$j_{\textrm {CD}0}$ and 
$d_{\textrm {CD}}$ it is observed that the slopes of the linear functions do not depend explicitly on the blob amplitude, but there is a weak dependence on the NTM order. To be exact, the slope corresponding to 
$j_{\textrm {CD}0}$ is slightly larger for 
$m/n=2/1$, whereas the slope of 
$d_{\textrm {CD}}$ it is slightly larger for 
$m/n=3/2$; these effects are related to the differences in the island geometry between the two cases, as explained in Isliker et al. (Reference Isliker, Chatziantonaki, Tsironis and Vlahos2012). Nevertheless, the ranges of values for 
$j_{\textrm {CD}0}$ and 
$d_{\textrm {CD}}$ depend both on 
$m/n$ and 
$\delta n_e/n_e(r_b)$. The dependence on 
$m/n$ is deduced through the comparison of curves with the same value of blob amplitude, and has been analysed in Tsironis (Reference Tsironis2013) for the mode 
$2/1$ and in Isliker et al. (Reference Isliker, Chatziantonaki, Tsironis and Vlahos2012) for the mode 
$3/2$. Regarding the dependence on the density fluctuation strength, by inspecting the curves referring to the same mode order it becomes apparent that, as 
$\delta n_e/n_e(r_b)$ increases, the values of 
$j_{\textrm {CD}0}$ reduce and the ones of 
$d_{\textrm {CD}}$ increase accordingly. This was expected due to the effects of beam axis shifting and electric field profile broadening by the blob structures on the wave propagation, which are deleterious with reference to the control requirement of accurate ECCD deposition on the magnetic island's flux surfaces.
In the last part of this section, we study the properties of the MRE in the presence of the effects and restrictions mentioned above. The central point in our analysis is to solve the MRE self-consistently (i.e. including the modification of the island size at each time step) by using the data output for the current drive from the previous modelling step, and to compare the results between the scenarios where the density fluctuations are present and absent for both cases of the NTM instability order. Within our study, the parameters which are bound to affect the form of the self-consistent solution are the initial width 
$W_0$, because for larger islands the undesired effects on the ECCD stabilizing term are more important, and the density fluctuation amplitude because, as seen previously, it affects the shape of the ECCD profile via the parameters 
$j_{\textrm {CD}0}$ and 
$d_{\textrm {CD}}$. Since the effect of the value of 
$W_0$ has been studied extensively in Chatziantonaki et al. (Reference Chatziantonaki, Tsironis, Isliker and Vlahos2013), here we only consider all the aforementioned cases for 
$m/n$ and 
$\delta n_e/n_e(r_b)$ at an ITER-specific initial condition of the island width (6 cm). In figure 5 we visualize the calculated solution of the MRE, which is expressed in terms of the time series of the island size 
$W(t)$, and in our model the applied ECCD is computed in the exact (non-axisymmetric) geometry of the magnetic islands.
Figure 5. Magnetic island size, as computed from the self-consistent solution of the MRE with the ECCD stabilizing term computed in non-axisymmetric geometry, for all the aforementioned values of the initial island size and density fluctuation strength in the two different cases for the mode order: (a) 
$m/n=2/1$, (b) 
$m/n=3/2$.
From the figures given above, it is found that the full stabilization of the mode is delayed in the presence of density fluctuations. As expected, the specific time lag of delay does not depend strongly on the mode order, but it has a significant dependence on the fluctuation amplitude. To be specific, and with respect to the scenario where density fluctuations are not present in the plasma, for the mode 2/1 the stabilization occurs with a delay of approximately 0.8 s when 
$\delta n_e/n_e(r_b)=0.5$ and of nearly 3.5 s when 
$\delta n_e/n_e(r_b)=1$, whereas, in the same manner, for the mode 3/2 the stabilization occurs roughly 0.6 s later when 
$\delta n_e/n_e(r_b)=0.1$ and almost 2.9 s later when 
$\delta n_e/n_e(r_b)=1$. To summarize, this means that the time needed for the full mode stabilization in the presence of blobs is an increasing function of both the mode order as well as the density perturbation strength. The behaviour described here is strongly connected to the displacement of the ECCD deposition with respect to the magnetic island's 
$O$-point and to the broadening of the current drive profile, which are induced by the beam axis shifting and width broadening brought up by the density fluctuations.
5. Conclusion
In this paper, we have presented a first-principles approach to the development of an integrated model for NTM dynamics in the presence of stabilizing ECCD and plasma density fluctuations. The toolbox consists of three models, associated with the different physics processes, which are applied in sequential order: an advanced model for the beam propagation through the density blobs, which is based on a transfer matrix applied to an electromagnetically homogenized medium; a simpler model for the propagation outside the blob, which is based on asymptotic ray tracing; and a solver for the nonlinear evolution of the magnetic island growth. The simulation results indicate that the density fluctuations affect the mode stabilization through the alteration of the intended current drive, namely the shifting of the deposition radius with respect to the island's 
$O$-point and the broadening of the deposition profile. These effects impose an increase in the time needed for mode stabilization, which scales proportionally with the perturbation strength.
The aforementioned effects may be of importance in some cases, and therefore should be taken into account by the control system commissioned with the NTM stabilization (involving steerable wave launchers). For efficient mode control in ITER, it is expected that the beam power should be deposited very close to the island centre with an allowed misalignment of 1–2 cm in time scales smaller than 1 s (Henderson et al. Reference Henderson, Heidinger, Strauss, Bertizzolo, Bruschi, Chavan, Ciattaglia, Cirant, Collazos and Danilov2008). Both the displacement and the broadening of the ECCD deposition are equivalently undesirable effects, with the former one being slightly more necessary to control. In this context, and within the parameter range simulated in this work, it becomes apparent that for large density fluctuation amplitudes the increase in the needed time for stabilization is critical and may not be manageable by the control system, thus annihilating the stabilizing effect of the ECCD. Our current study is in a limited range of values for the density fluctuation parameters that are effective, namely the perturbation amplitude, the blob size and the turbulent layer width, which we will make an effort to extend in future work.
The computational resources (CPU time and memory) needed by our model are not extreme; a full-wave simulation, in the range of parameters explored here, takes less than 10 minutes and has normal memory demands, whereas the ray tracing code and MRE solver run in less than a minute. The total CPU time required also reduces because we use ray tracing for the propagation outside the blob region, instead of a full-wave simulation from beam launch up to the magnetic island. The latter would increase significantly the overall computational burden not only due to the need for full-wave computations on a longer propagation path, but also due to the required evaluation of the hot-plasma dielectric response near the EC resonance, where asymptotic codes possess an advantage in treating more efficiently the wave absorption and current drive (see e.g. Tsironis & Papadopoulos (Reference Tsironis and Papadopoulos2014) for details). The difficulty in applying our model sequence lies in the different steps required in order to supply each model part with the proper input: as seen in § 3, the assignment of proper beam input to the ray tracing code requires the computation of the moments of the full-wave EM field output, whereas the coupling of the ECCD to the MRE solver requires a regression analysis over the output data of several ray tracing computations. These actions require an important amount of post-processing effort which, for the moment, cannot be performed automatically.
A discussion on the limitations of our models is necessary. First, the full-wave solver and the plasma homogenization technique may be improved by including more realistic magnetic equilibria and plasma profiles (coming for example from experimental data), a beam description based on non-plane waves, as well as the reflection effect at the vacuum–plasma crossing (which will allow a more realistic simulation of the beam launching conditions). Second, the ray tracing code, apart from accommodating realistic magnetic field and plasma profiles, requires also an increase of the consistency in the computation of the ECCD profile (e.g. by using a larger number of rays), as well as the development of a more realistic description for the wave beam (such as, for example, in paraxial beam tracing). Finally, improvements to the MRE solver are identified in the inclusion of the description of rotating magnetic islands, gyrotron power modulation and EC plasma heating, via the proper formulation of additional 
$\varDelta ^{\prime }$ terms. All these extensions are considered in terms of current work on a more thorough analysis of the effects presented in this paper.
Acknowledgements
This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
Editor Paolo Ricci thanks the referees for their advice in evaluating this article.
Declaration of interests
The authors report no conflict of interest.
Open access
