1. General considerations on wave propagation for plasma wave heating and current drive
It is useful to start by citing some reviews and books that provide the background on the different kinds of waves used to heat fusion plasmas. General information about waves in plasmas can be found in the books by Swanson (Reference Swanson2003) and Stix (Reference Stix1992); and the required physics background on waves for heating plasmas and the early history of wave heating are available in Cairns (Reference Cairns1991). Reviews of electron cyclotron waves can be found in Bornatici et al. (Reference Bornatici, Cano, de Barbeiri and Englemann1983), of lower-hybrid wave in Bonoli (Reference Bonoli1985), and of ion cyclotron wave in Jaeger, Batchelor & Stallings (Reference Jaeger, Batchelor and Stallings1993), Fuchs et al. (Reference Fuchs, Ram, Schultz, Bers and Lashmore-Davies1995) and Wilson & Bonoli (Reference Wilson and Bonoli2015). Bearing in mind the main elements of wave propagation and damping mechanisms, as well as the requirements for ITER (International Tokamak Experimental Reactor) heating (Poli et al. Reference Poli, Kessel, Bonoli, Batchelor, Harvey and Snyder2014), the method for modelling heating systems and the latest physical phenomena to take into account following new experimental evidence will be presented.
In order to model or simulate plasma heating in fusion devices, the first question to be addressed is this: What are the minimal physical processes that should be taken into account? In fact, this question is not at all simple, since, as the plasma is heated, its parameters change. Therefore, a self-consistent description becomes complex, owing to the fact that transport is, at least, a function of the local temperature gradients, which are modified by the wave energy deposition, which is itself a function of the local values of the plasma parameters. In previous work, we have only considered smooth plasma parameters, which used to be the standard approach. Recently, the role of density fluctuations started to be seen as an important parameter, especially at the plasma edge, in a passive way, as well as for the electron cyclotron heating (Tsironis et al. Reference Tsironis, Peteers, Isliker, Strinzi and Chatziantonaki2009; Decker, Peysson & Coda Reference Decker, Peysson and Coda2012) or the lower-hybrid current drive (Peysson & Decker Reference Peysson and Decker2014). It also affects diagnoses that use waves, since it changes the properties of the probing beam at a higher order. At the next order, wave trapping effects, multi-scattering and nonlinear reflection are introduced (Heuraux et al. Reference Heuraux, Gusakov, Popov, da Silva and Irzak2010; da Silva et al. Reference da Silva, Heuraux, Gusakov and Popov2010). Possibilities to generate parametric instabilities during heating (Gusakov, Popov & Saveliev Reference Gusakov, Popov and Saveliev2014) or mode conversion to remove accessibility conditions associated with cut-offs (Irzak & Popov Reference Irzak and Popov2008) are reconsidered. In these last cases, the code should describe multi-mode propagation, and the plasma domain can be reduced to a poloidal cross-section, owing to the very small parallel wavenumber values for the turbulence. Considering these remarks is enough (i) to deduce that the minimal set of equations driving the wave propagation should be written at least in two dimensions, (ii) to be able to integrate fluctuating plasma parameters, and (iii) to describe several modes simultaneously. That is to say, a relevant code should compute all electric field components in two dimensions for a given parallel wavenumber. The radio-frequency (RF) sheath effects have to be integrated in the simulation of ion cyclotron resonant heating (ICRH) simulations through the particular boundary conditions to describe the spurious effects seen in experiments (Wukitch et al. Reference Wukitch, LaBombard, Lin, Lipschultz, Marmar and Reinke2009; Jacquet et al. Reference Jacquet, Colas, Mayoral and Arnoux2011), which requires the computation of slow and fast wave fields at the same time. In general, the available codes offer a single-mode or multi-mode description with some sophistication, such as a non-Maxwellian distribution function (Brambilla & Bilato Reference Brambilla and Bilato2009; Dumont & Zarzoso Reference Dumont and Zarzoso2013; Peysson & Decker Reference Peysson and Decker2014) and ad hoc tools to cross a wave resonance. Concerning this last point, new mathematical developments allow for solutions without the introduction of artificial dissipation (Després, Imbert-Gérard & Weder Reference Després, Imbert-Gérard and Weder2014).
2. Wave heating simulations for fusion plasmas
2.1. Introduction
An initial idea on the approximations made to describe wave propagation using Maxwell’s equations can be the understanding of the idealized model analysis written in the book by Jackson (Reference Jackson1999). With that, it is possible to create a clear picture of the different approximation levels. At the same time, a second question arises: How does one describe the behaviour of the medium versus the studied wave? A third question is: How does one do this without missing a relevant mechanism? To address them, more and more physics needs to be taken into account. This increase of the physics content induces an increase of the computation time, which can exceed the possibilities of the actual computer capabilities. So a compromise has to be found to work at same computation time between the physical inputs and numerical accuracy, that is to say, the number of grid points or finite elements. That has an impact on the exactness of the simulation carried out. This point raises questions: How does one achieve the best numerical description of a given heating experiment whilst meeting the numerical constraints (computation time, stability, accuracy)? How does one know if there are missing constraints (for example, voltage limitation to avoid arcing) or not? How does one introduce the technical constraints when one wants to optimize a heating system or attain a given objective such as the induced plasma current in a given ITER scenario?
An example to illustrate this is to find the optimal efficiency as a function of the incident angle of the mirror, knowing the energy deposition into the plasma and the depolarization rate, which depends on the aging of the material of the mirror. The latter remains unknown even if, to prevent this problem, a cleaning mirror system is designed with the fusion plasma exposure in mind. In fact, we are far from the possibility of optimizing a heating system due to the lack of information. For example, the ITER plasma density profile in front of the ICRH antenna is not available, although this input is essential to determine the properties of the plasma–antenna coupling. Another distant possibility is to find the set of parameters in which a given heating system is able to work, following the specifications of an experiment. It is still not possible to determine the margins in which this experiment can be performed, knowing that heating induces changes to the plasma parameters, including the turbulent transport and the magnetohydrodynamic (MHD) activities.
As a partial conclusion, the description of the wave propagation, including absorption (damping), requires the use of Maxwell’s equations and adapted partial differential equations (PDEs) that describe the response of the plasma. These equations could be linear or nonlinear, depending on the wanted requirements and the working conditions. In any case, the geometry of the system should be preserved, in order to keep the exact boundary conditions regarding the possibility of having oblique RF sheaths. These RF sheaths strongly change the behaviour of the slow wave in front of an ion cyclotron antenna. The wave polarization’s role on the damping mechanisms is essential. Therefore, the next generation of codes describing a wave heating system needs to include a multi-mode description able to follow the wave polarization changes in order to optimize energy deposition into the plasma.
As a preliminary to that future development, the role of turbulence in the linear mode conversion efficiency in the O-X-B heating scheme was recently published in Popov (Reference Popov2015) and further developments of the SWITCH code for ICRH heating in Jacquot et al. (Reference Jacquet, Colas, Mayoral, Arnoux, Bobkov, Brix, Coad, Czarnecka, Dodt, Durodie, Ekedahl, Frigione, Fursdon, Gauthier, Goniche, Graham, Joffrin, Korotkov, Lerche, Mailloux, Monakhov, Noble, Ongena, Petrzilka, Portafaix, Rimini, Sirinelli, Riccardo, Vizvary, Widdowson and Zastrow2014). Other points also have to be considered in the time-dependent cases, where Maxwell’s equation solvers are coupled with J-solvers, numerical stability and energy conservation. Both points become problems at a large number of iterations, although they are still below what is required for ITER (da Silva et al. Reference da Silva, Campos Pinto, Després and Heuraux2015). This kind of numerical problem disappears when the time-dependent part is removed under the following assumptions: the launch wave is considered as monochromatic (continuous emission at the same frequency) and the plasma parameters present time-scale evolution much larger than the time of flight of the heating wave. Extra requirements should be considered for the solvers, which should be able to describe wave propagation in inhomogeneities present in experimental configurations as well as a high level of plasma fluctuations and to be able to cross wave resonances. The previous assumptions are generally used in the cases of wave heating simulations, which require also the description of source terms and boundary conditions in a relevant way. To finish providing relevant information about the behaviour of the heating system and its impact on the plasma parameters, these solvers have to be included in an integrated tokamak modelling framework, but they have to be optimized before their integration into such platform.
2.2. Basic equations
In the most general situation, Maxwell’s equations are used and can be written as follows:
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{E}=\frac{{\it\rho}(\boldsymbol{E},\boldsymbol{B})}{{\it\epsilon}_{0}}\quad \text{(Poisson's law),} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0\quad \text{(magnetic flux conservation law),} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\times \boldsymbol{E}=-\frac{\partial \boldsymbol{B}}{\partial t}\quad \text{(Faraday's law),} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\rm\nabla}}\times \boldsymbol{B}={\it\mu}_{0}\boldsymbol{J}(\boldsymbol{E},\boldsymbol{B})+\frac{1}{c^{2}}\frac{\partial \boldsymbol{E}}{\partial t}\quad \text{(Ampère's law),} & \displaystyle\end{eqnarray}$$
$\boldsymbol{J}$
and charge density
${\it\rho}$
evolutions as functions of (
$\boldsymbol{E},\boldsymbol{B}$
), where
${\it\epsilon}_{0}$
is the dielectric constant of vacuum,
${\it\mu}_{0}$
is the permeability of vacuum and
$c$
is the speed of light in vacuum. Here, we consider only the plasma parameter set of magnetic fusion, meaning the plasma quantum effects are neglected and first-order relativistic corrections are taken into account. In this framework, the most general description of plasma responses to wave excitation is driven by Vlasov’s equation. This includes kinetic effects able to describe precisely all the wave–particle interactions such as Landau damping or cyclotron resonances. Since the plasma is composed of different particle species, for each kind of particle, one Vlasov’s equation should be written and coupled through Maxwell’s equations or the introduction of a Fokker–Planck equation as described in van Eester (Reference van Eester2012). Therefore, the set of PDEs necessary to define (
$\boldsymbol{J},{\it\rho}$
) can be written as follows: 
$$\begin{eqnarray}\frac{\partial f_{s}}{\partial t}+\boldsymbol{v}_{\boldsymbol{s}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{\boldsymbol{r}}f_{s}+\frac{q_{s}}{m_{s}}[\boldsymbol{E}+\boldsymbol{v}_{\boldsymbol{s}}\times (\boldsymbol{B})]\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{\boldsymbol{v}}\,f_{s}=0,\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle {\it\rho}(\boldsymbol{r},t)=\sum q_{s}\int \text{d}\boldsymbol{v}\,f_{s}(\boldsymbol{r},\boldsymbol{v},t),\quad \boldsymbol{J}(\boldsymbol{r},t)=\sum q_{s}\int \text{d}\boldsymbol{v}\,\boldsymbol{v}f_{s}(\boldsymbol{r},\boldsymbol{v},t), & & \displaystyle\end{eqnarray}$$
where
$q_{s}$
is the charge,
$m_{s}$
is the mass,
$\boldsymbol{v}_{\boldsymbol{s}}$
is the particle velocity and
$f_{s}$
is the distribution function in phase space of the species
$s$
.
In the three-dimensional (3D) case, this PDE set (2.5) is too demanding in terms of computational resources. Thus, one needs to remove part of the physics by adding relevant assumptions, whilst preserving the main effects of wave injection for heating or current drive, even if you can use near- to far-field transformation (Taflove & Hagness Reference Taflove and Hagness2000). That means that one has to focus on specific studies and forget, momentarily, having a universal tool able to describe the wave heating system in its entirety. To take into account all the possible synergy available or the spurious effects generated by the wave itself, such as plasma density fluctuations, one needs a suite of integrated codes assumed to converge, which is not at all guaranteed. To be practical, we now focus on possible assumptions to reduce the complexity of the model, first in general, and then for each wave heating system.
The most common assumption within a good framework for wave propagation for a polarization of the launched wave is the cold plasma approximation. Generally, it is mostly because the hot plasma corrections are small. The relativistic corrections have to be taken into account earlier than the hot plasma corrections when the electron temperature becomes greater than a few kiloelectronvolts (Bindslev Reference Bindslev1991), and these relativistic effects can be introduced through an electron mass correction. The relativistic effects induce an increase of the effective mass of the electron, which drives a frequency down-shift of the electron cyclotron frequency. The damping processes in tokamak plasma are essentially a result of a resonant interaction between waves and particles, which requires knowledge of the velocity distribution function in an accurate manner. That means the cold plasma approximation is inappropriate to describe the wave absorption in the tokamak plasma core. To compute the damping rate accurately, good knowledge is required of the velocity distribution function, which is generally provided by a Fokker–Planck equation, which permits one to introduce the trapped particle contributions. A detailed discussion on this point can be found in Gnesin (Reference Gnesin2011) for electron cyclotron heating, in Peysson & Decker (Reference Peysson and Decker2014) for the lower-hybrid heating and current drive, and in Brambilla & Bilato (Reference Brambilla and Bilato2009) for the ICRH.
Usually in wave heating simulation the monochromatic wave assumption is used. In this case, the most used solver to describe wave propagation allows interference, scattering and multi-reflection effects to be taken into account. This solver can describe a single wave equation for a given polarization or multi-modes depending on the computed components of the electric field. This solver type is called full wave, and corresponds to a Helmholtz-like equation, which comes from Maxwell’s equations or their equivalent in terms of
$\boldsymbol{A}$
and
${\it\phi}$
and reduces to
$$\begin{eqnarray}\left\{\begin{array}{@{}l@{}}\displaystyle -\boldsymbol{{\rm\nabla}}\times \boldsymbol{{\rm\nabla}}\times \boldsymbol{E}(\boldsymbol{r},t)+\frac{{\it\omega}^{2}}{c^{2}}(\boldsymbol{E}(\boldsymbol{r},t))+\frac{\text{i}}{{\it\omega}{\it\epsilon}_{0}}\boldsymbol{J}_{\boldsymbol{p}}(\boldsymbol{E})=-\text{i}{\it\omega}{\it\mu}_{0}\boldsymbol{J}_{ant}(\boldsymbol{r}_{ant}),\quad \\ \displaystyle \boldsymbol{J}_{\boldsymbol{p}}(\boldsymbol{E})=\overline{\overline{{\bf\sigma}}}\boldsymbol{\cdot }\boldsymbol{E}\quad (\text{linear case})\quad \text{or}\quad \boldsymbol{J}_{\boldsymbol{p}}(\boldsymbol{r},t)=\sum q_{s}\int \text{d}\boldsymbol{v}\,\boldsymbol{v}f_{s}(\boldsymbol{r},\boldsymbol{v},t),\quad \end{array}\right.\end{eqnarray}$$
where
$\overline{\overline{{\bf\sigma}}}$
is the cold or hot plasma conductivity tensor corresponding to the heating frequency
${\it\omega}$
; and
$\boldsymbol{J}_{ant}$
is the current density source term corresponding to the antenna emission, for which the current density expression in space should represent the effective wavenumber spectrum of the antenna, taking into account the plasma properties in front of it. It is clear that, for a hot plasma description, the evolution of the plasma conductivity is taken into account and it is generally done using the Vlasov equation or the Fokker–Planck equation as mentioned before, which can also include other heating contributions through adapted quasi-linear coefficients. Another possibility also explored is using a particle-in-cell (PIC) code coupled to Maxwell’s equation, and use of an implicit scheme overcomes the stability problem (Smithe Reference Smithe2007). However, the requirements in terms of computational resources can become intractable in ITER cases.
2.3. Numerical tools, source terms and boundary conditions
A recent review describes briefly the numerical techniques and schemes developed to compute the properties of the different heating schemes (van Eester Reference van Eester2012). Brief descriptions of the different kinds of wave solvers can be found in Heuraux & da Silva (Reference Heuraux and da Silva2012) and Heuraux et al. (Reference Heuraux, Faudot, da Silva, Jacquot, Colas, Hacquin, Teplova, Sysoeva and Gusakov2014), in order to clarify what are their relevant domains in wave heating. However, several important points still have to be discussed in order to complete the landscape of wave heating simulations in magnetized plasmas: the source terms and the boundary conditions. Usually, the boundary conditions are specific to a kind of heating, which can be very different from the standard electromagnetic boundary conditions, such as perfect electric conductor (PEC) or others described elsewhere (Taflove & Hagness Reference Taflove and Hagness2000). A complex boundary condition exists in the case of ICRH due to the difference between electron and ion mobilities. Indeed, a charge separation can be induced and modifies the parallel electric field component (parallel to the DC magnetic field). That changes the behaviour of the wave against the metallic wall, and at the same time generates a rectified DC potential, owing to a lack of electrons induced by the RF electric field. Therefore, it is necessary to introduce sheath boundary conditions for the ion cyclotron heating to take into account the fact that the parallel component of the incident wave can be more or less screened at the wall. This depends on the incident angle of the magnetic field line on the wall and on the local values of the plasma parameters (Kohno, Myra & D’ipolitto Reference Kohno, Myra and D’ipolitto2012; Jacquot et al. Reference Jacquet, Colas, Mayoral, Arnoux, Bobkov, Brix, Coad, Czarnecka, Dodt, Durodie, Ekedahl, Frigione, Fursdon, Gauthier, Goniche, Graham, Joffrin, Korotkov, Lerche, Mailloux, Monakhov, Noble, Ongena, Petrzilka, Portafaix, Rimini, Sirinelli, Riccardo, Vizvary, Widdowson and Zastrow2014). We will further develop this concept of boundary conditions later on in § 3.
Another boundary condition is often used to reduce the mesh size, which corresponds to free wave propagation through the mesh edge without reflection, which is the so-called perfect matching layer (PML) boundary condition (Bérenger Reference Bérenger1994). Although very useful, the PML boundary condition does not work in anisotropic media, with propagating modes having opposite phase velocities or negative value of the components of the dielectric tensor (Jacquot et al. Reference Jacquot, Colas, Clairet, Goniche, Heuraux, Hillairet, Lombard and Milanesio2013), as is the case for a meta-material. A refined computation of the near fields in the case of wave heating can be possible using PML layers, assuming that the single-pass absorption is valid. The spatial domain is then reduced to the region where the near fields are of interest (Jacquot et al. Reference Jacquet, Colas, Mayoral, Arnoux, Bobkov, Brix, Coad, Czarnecka, Dodt, Durodie, Ekedahl, Frigione, Fursdon, Gauthier, Goniche, Graham, Joffrin, Korotkov, Lerche, Mailloux, Monakhov, Noble, Ongena, Petrzilka, Portafaix, Rimini, Sirinelli, Riccardo, Vizvary, Widdowson and Zastrow2014). The implementation of the PML layers depends on the working domain, frequency or time domain (Taflove & Hagness Reference Taflove and Hagness2000).
In order to define the boundary conditions appropriately (PEC,
$\text{PML},\dots$
), an exact knowledge of the geometry of the launching structure is required. In the case of the sheath boundary condition (SBC), knowledge of the angle between the magnetic field line and the wall is essential, which requires a good description of the antenna design and the magnetic configuration. This is particularly important at grazing incidence, when the Debye sheath disappears. As a consequence, there is no charge separation on average (Stangeby Reference Stangeby2012) and rectification of the RF potential is reduced or disappears as shown using two-dimensional (2D) PIC simulations depending of the RF potential amplitude (works under way in the framework of a EUROfusion Enabling Research contract). The complete knowledge of the antenna structure and surroundings is also needed in order to compute the near fields in front of it. This computation is essentially devoted to the so-called plasma–antenna coupling codes.
Then, its map can be used as an input in the wave equation solver as a source term. Generally speaking, this kind of source term is provided by antenna–plasma coupling codes to stay as close as possible to the experiments. Depending on the studies, synthetic source terms can be implemented as a Gaussian beam or anything else one wants. However, some caution has to be taken to avoid bidirectional emission (
$-\boldsymbol{k}$
and
$+\boldsymbol{k}$
at the same time). Having unidirectional emission of the launching wave is less simple to implement or requires sophisticated tools if one wants to introduce a transparent source (Taflove & Hagness Reference Taflove and Hagness2000; da Silva et al.
Reference da Silva, Heuraux, Hacquin and Manso2005). To simulate an antenna in interaction with its surroundings, the source term should be as close as possible to the hardware components to simulate accurately the radiation wave pattern of the wave launcher. This would allow one to inject power into the antenna, for example, to study the resonances between a horn antenna put in a cavity between two ITER blanket modules (da Silva, Heuraux & Manso Reference da Silva, Heuraux and Manso2006a
). At the moment, the wave launcher will be equipped with diagnostics to measure the changes in the density profile in front of the antenna to study the interplay between the edge plasma parameters and the injected power in the cases of ICRH and lower-hybrid current drive (LHCD), which is not actually well included in the present wave heating simulations, and changes are required to understand the plasma–antenna coupling.
2.4. Heating simulation
In general, as the launched frequency is fixed for a given plasma shot, the wave description can move to the frequency domain, which is generally more comfortable than the time domain, which normally has to fulfil a CFL condition (Courant, Friedrichs & Lewy Reference Courant, Friedrichs and Lewy1928) to be numerically stable. A first problem on the exactness can appear after this choice, especially when two or more modes exist in the computation domain at the same time for this chosen frequency, as is the case for the ICRH simulations, where slow and fast waves exist at the plasma edge (Colas et al. Reference Colas, Jacquot, Heuraux, Faudot, Crombé, Kyrytsya, Hillairet and Goniche2012; Jacquot et al. Reference Jacquet, Colas, Mayoral, Arnoux, Bobkov, Brix, Coad, Czarnecka, Dodt, Durodie, Ekedahl, Frigione, Fursdon, Gauthier, Goniche, Graham, Joffrin, Korotkov, Lerche, Mailloux, Monakhov, Noble, Ongena, Petrzilka, Portafaix, Rimini, Sirinelli, Riccardo, Vizvary, Widdowson and Zastrow2014). A question then arises: Are these modes decoupled or not, and, if not, how are they translated in terms of mode coupling or linear mode conversion? For example, the ordinary mode (O-mode) and extraordinary mode (X-mode) can be coupled linearly through density parallel to the magnetic field or its fluctuations (Colas et al. Reference Colas, Zou, Paume, Chareau, Guiziou, Hoang, Michelot and Gresillon1998), and a wanted (or not) depolarization of the launching wave can occur. This phenomenon can possibly take place during electron cyclotron heating, and it can be used after a good choice of the oblique incidence of the launching wave on the O-mode, as done in the O-X-B double-mode conversion heating scheme (Laqua et al. Reference Laqua, Erckmann, Hartfuss and Laqua1997). In such a heating scheme, two questions are still open: the role of the density fluctuations in the conversion regions (Popov Reference Popov2015); and the effect of the beam widening induced by the fluctuations at the edge, as shown for ITER electron cyclotron heating cases (Peysson & Decker Reference Peysson and Decker2014; Sysoeva et al. Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2014). Following this example and trying to answer these questions, we need to look at the different wave heating systems independently, since each frequency range is indeed separated from the others.
Now, we focus on heating in tokamaks. Three main types of solvers can be identified, ray tracing, Helmholtz’s equation and Maxwell’s equations, corresponding to different orders of approximations.
The simulations using the highest level of approximation associated with a ray tracing code correspond to a single-mode description in the Wentzel, Kramers, and Brillouin (WKB) approximation which includes no cut-off (local index equal to zero), no resonance (local index value going to infinity) in the computation, and where wave propagation is described without interference, scattering or mode conversion. An additional equation can be added to provide the spatial evolution of the amplitude, as for the ray tracing code in toroidal geometry with adapted coordinate system, as done in Peysson, Decker & Morini (Reference Peysson, Decker and Morini2012) and Marushchenko, Turkin & Maassberg (Reference Marushchenko, Turkin and Maassberg2014) or for Gaussian beam (Bertelli et al. Reference Bertelli, Maj, Poli, Harvey, Wright, Bonoli, Phillips, Smirnov, Valeo and Wilson2012). While the dispersion relation is known, it is possible to apply this tool if the initial conditions are also known. The vacuum–plasma border crossing should be treated very carefully. A stopping condition based on wave absorption can be added, which permits, at the same time, the determination of where the energy is deposited. This does not take into account scattering, cavity effects, mode conversion, or the interference pattern. Surprisingly, the results obtained using the ray tracing code provide an overall good behaviour of the wave.
2.5. Ion cyclotron, lower-hybrid and electron cyclotron heating codes
ICRH codes. A brief description of the full-wave codes used for the ion cyclotron frequency heating simulation is given in Budny et al. (Reference Budny, Berry, Bilato, Bonoli, Brambilla, Dumont, Fukuyama, Harvey, Jaeger, Indireshkumar, Lerche, McCune, Phillips, Vdovin and Wright2012). In a variational formalism used in the EVE code (Dumont & Zarzoso Reference Dumont and Zarzoso2013) there is no possibility to model the antenna–plasma coupling self-consistently, because the current density on the antenna is fixed and cannot be modified by the current induced on the antenna structure, which is limited to the straps description. To obtain relevant near-field computations, other ways are used to build antenna–plasma coupling codes, for example, boundary element method in the ICANT code, and integral method in TOPICA (see Pécoul et al. Reference Pécoul, Heuraux, Koch and Leclert2002; Lancellotti et al. Reference Lancellotti, Milanesio, Maggiora, Vecchi and Kyrytsya2006). Another one, the TORIC code, using spectral decomposition, was also developed to study the ICRH scenarios (Brambilla Reference Brambilla1999). The SCENIC code is described in Jucker et al. (Reference Jucker, Graves, Cooper, Mellet, Johnson and Brunner2011), where the full-wave solver part should be improved to become relevant to ITER and stellarator applications. These codes follow a suite of home-made codes, previously mentioned for ion cyclotron, plus the TOMCAT code and its upgrades, in which a coupling code is included (van Eester & Koch Reference van Eester and Koch2001). Recently, commercial software packages are being used to solve the Helmholtz equation applying the finite element method (FEM). For instance, COMSOL is used in many cases: the ion cyclotron in the SSWICH code (Jacquot et al. Reference Jacquot, Colas, Clairet, Goniche, Heuraux, Hillairet, Lombard and Milanesio2013), and the lower-hybrid (LH) cases (a) and (b).
LH codes. Two codes based on COMSOL software deal with the lower-hybrid heating, (a) the most recent the LHEAF code (Shiraiwa et al. Reference Shiraiwa, Meneghini, Parker and Bonoli2010) and (b) the coupling code PICCOLO2D including nonlinear effects (Preynas et al. 2013), which is based on the same principle as the one-dimensional (1D) code ALOHA (Hillairet et al. Reference Hillairet, Voyer, Ekedahl, Goniche, Kazda, Meneghini, Milanesio and Preynas2010). The lower-hybrid wave propagation code TORLH (Wright et al. Reference Wright, Valeo, Phillips, Bonoli and Brambilla2008) follows the same numerical method as the TORIC code. The coupling code TOPLHA (Milanesio et al. Reference Milanesio, Meneghini, Maggiora, Guadamuz, Hillairet, Lancellotti and Vecchi2012) uses the same tools as developed in the TOPICA code. To study the LHCD, a ray-tracing code coupled to a Fokker–Planck solver has been used, including the effect of density fluctuations, which can also be applied for studying electron cyclotron current drive (Peysson & Decker Reference Peysson and Decker2014). An up-to-date review of the lower-hybrid heating and LHCD can be found in Bonoli (Reference Bonoli2014).
EC codes. In the frequency range of electron cyclotron heating, the Helmholtz solvers were used for diagnostics (Fanack et al. Reference Fanack, Boucher, Heuraux, Leclert, Clairet and Zou1996; Mazzucato Reference Mazzucato1998) a long time ago. Until recently, the ray tracing method was estimated to be enough to study beam propagation in ITER plasmas (Prater et al. Reference Prater, Farina, Gribov, Harvey, Ram, Lin-Liu, Poli, Smirnov, Volpe, Westerhof and Zvonkov2008), but it is clear now that the density fluctuations have to be included (Decker et al. Reference Decker, Peysson and Coda2012), and, for these studies, full-wave codes are required for the frequency range of electron cyclotron heating or current drive (Sysoeva et al. Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2014, Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2015).
Before going further, a remark has to be made on the commercial software available, since the tendency to use multi-physics software is becoming stronger and stronger. This can induce problems if specific difficulties associated with new types of numerical schemes are needed that are not implemented in the commercial software, such as the asymptotic preserving scheme for very high anisotropy (Crouseilles, Lemou & Méhats Reference Crouseilles, Lemou and Méhats2013), just like the dielectric tensor of a magnetized plasma for ICRH. Furthermore, new developments of no interest for industry can lead to a lack of numerical tools in commercial software. For example, PML in magnetized plasma or unsolved numerical instabilities found in it (Heuraux et al. Reference Heuraux, Faudot, da Silva, Jacquot, Colas, Hacquin, Teplova, Sysoeva and Gusakov2014). Their use on a high-performance computer (HPC) can be another source of difficulties due to the commercial policy of ‘one licence per node’.
The next degree of sophistication is to obtain the direct response of current density,
$\boldsymbol{J}$
, and the charge separation,
${\it\rho}$
, using a set of differential equations driving their time evolution during the heating, to insert as the source terms in Maxwell’s equations. Depending on the degree of approximation, one is able to describe cold plasma linear response in the case of the equations of motion of particles (Després et al.
Reference Després, Imbert-Gérard and Weder2014), or kinetic responses for the charge separation
${\it\rho}$
and the current density
$\boldsymbol{J}$
as defined (2.5) in the case of the Poisson–Vlasov equations, with one equation for each species (Kuley et al.
Reference Kuley, Wang, Lin and Wessel2013), or PIC code in Smithe (Reference Smithe2007) and Smithe, Myra & D’Ippolito (Reference Smithe, Myra and D’Ippolito2014). These last descriptions based on the mean-field theory allow one to describe, in principle, nonlinear and boundary condition effects, but they require HPC machines. There are numerical stability problems due to the anisotropy introduced by the external DC magnetic field (Heuraux et al.
Reference Heuraux, Faudot, da Silva, Jacquot, Colas, Hacquin, Teplova, Sysoeva and Gusakov2014). One is connected to the very large difference between responses parallel to and perpendicular to the magnetic field, and another relates to the accumulation of round-off error in the long computation time (more than
$10^{6}$
time steps) that drives numerical instabilities (NIs) as shown in da Silva et al. (Reference da Silva, Campos Pinto, Després and Heuraux2015). On such long runs, energy conservation should be fulfilled to preserve the statistical properties of the physical system. This is especially true if you want to study, after averaging, the macroscopic effects of turbulence on the wave propagation. This last part is an on-going project showing that new methods improve the wave propagation descriptions (Tierens & de Zutter Reference Tierens and de Zutter2012) and the PIC code responses (Campos Pinto et al.
Reference Campos Pinto, Sonnendrücker, Friedman, Grote and Lund2014) or gyrokinetic responses using an asymptotic preserving scheme (Crouseilles et al.
Reference Crouseilles, Lemou and Méhats2013).
Up to now, we have only introduced the elements needed to describe the wave propagation in a given plasma, considering, to a greater or lesser extent, the local modifications of the plasma properties. To see the full impact of the heating on the target plasma, the equilibrium and transport equations have to be solved in parallel. This is necessary in order to provide new values of the plasma parameters for the wave heating codes, including an equation of current density transport (in the case of current drive (Peysson & Decker Reference Peysson and Decker2014)). Owing to the role of density fluctuations at the edge on the wave heating performance, a turbulence code also has to be used. This can be done iteratively using an integrated tokamak modelling, as in CRONOS (Artaud et al. Reference Artaud, Basiuk, Imbeaux, Schneider, Garcia, Giruzzi, Huynh, Aniel, Albajar, Ané, Bécoulet, Bourdelle, Casati, Colas, Decker, Dumont, Eriksson, Garbet, Guirlet, Hertout, Hoang, Houlberg, Huysmans, Joffrin, Kim, Kochl, Lister, Litaudon, Maget, Masset, Pégourié, Peysson, Thomas, Tsitrone and Turco2010) or in the European ITM platform (Falchetto et al. Reference Falchetto, Coster, Coelho, Scott, Figini, Kalupin, Nardon, Nowak, Alves, Artaud, Basiuk, Bizarro, Boulbe, Dinklage, Farina, Faugeras, Ferreira, Figueiredo, Huynh, Imbeaux, Ivanova-Stanik, Jonsson, Klingshirn, Konz, Kus, Marushchenko, Pereverzev, Owsiak, Poli, Peysson, Reimer, Signoret, Sauter, Stankiewicz, Strand, Voitsekhovitch, Westerhof, Zok and Zwingmann2014). At present, we are rather far from this situation. The focus of the effort should be on improving the depolarization mechanisms and the performance of each module (particularly on wave solvers), in which specific boundary conditions have to be implemented. Now, we present two points that are under development: the SBCs in magnetized plasma, and how the role of the density fluctuations can be simulated.
3. Sheath boundary conditions (SBC)
3.1. Context
The wave propagation of RF waves, radiated by complex antennae, and their energy deposition in the core of magnetized plasmas, have been described, for a long time, using sophisticated first-principles models in realistic geometry, as mentioned in the previous section. Nevertheless, these simulations do not explain the experimental facts seen at the plasma edge during ion cyclotron heating (Jacquet et al.
Reference Jacquet, Colas, Mayoral and Arnoux2011; Jacquot et al.
Reference Jacquet, Colas, Mayoral, Arnoux, Bobkov, Brix, Coad, Czarnecka, Dodt, Durodie, Ekedahl, Frigione, Fursdon, Gauthier, Goniche, Graham, Joffrin, Korotkov, Lerche, Mailloux, Monakhov, Noble, Ongena, Petrzilka, Portafaix, Rimini, Sirinelli, Riccardo, Vizvary, Widdowson and Zastrow2014). Comparatively, the simulation of anomalous RF power losses in the plasma edge is still less advanced, although the nonlinear wave–plasma interactions in the plasma edge often set the operational limits of RF heating systems. Peripheral wave damping in the ion cyclotron range of frequencies (30–100 MHz in present fusion devices) is attributed to a DC biasing of the edge plasma by RF sheath rectification induced by the different mobilities of ions and electrons (Noterdaeme & van Oost Reference Noterdaeme and van Oost1993). The first modelling of this nonlinear process was provided by Perkins (Reference Perkins1989), in an analogy of a double Langmuir probe driven by an oscillating RF voltage, estimated as the field-line-integrated RF field
$E_{\Vert }$
parallel to the confinement magnetic field
$B_{0}$
. However, in this kind of model, each open flux tube in the scrape-off layer (SOL) is considered independent of its neighbours. The nonlinear coupling between flux tubes can be found in Faudot, Heuraux & Colas (Reference Faudot, Heuraux and Colas2006). Its effects for an ICRH ITER antenna are available in Faudot et al. (Reference Faudot, Heuraux, Colas and Gunn2010). The main ingredient to determine the RF potential along a magnetic field line
$E_{\Vert }$
is generally computed from conventional antenna codes in the absence of sheaths, i.e. PECs are assumed to be in direct contact with the plasma. This procedure, although clearly not self-consistent, was widely implemented as the only tool able to model realistic wave launching structures (Pécoul et al.
Reference Pécoul, Heuraux, Koch and Leclert2002; Lancellotti et al.
Reference Lancellotti, Milanesio, Maggiora, Vecchi and Kyrytsya2006). The theoretical background of the SSWICH code that considers SBC is given in Colas et al. (Reference Colas, Jacquot, Heuraux, Faudot, Crombé, Kyrytsya, Hillairet and Goniche2012). However, only magnetic field lines perpendicular to the wall are considered in that paper. A crude model of oblique SBCs is described in Kohno et al. (Reference Kohno, Myra and D’ipolitto2012). A recent work on sheaths in magnetized plasmas mentions that the Debye sheath disappears at grazing angle (Stangeby Reference Stangeby2012). Consequently, the entire DC potential drop is associated with the magnetic pre-sheath. To evaluate the impact of the oblique SBCs, PIC simulation using VORPAL (commercial software) has been reviewed in the paper of Smithe et al. (Reference Smithe, Myra and D’Ippolito2014). This simulation, when done on a tiny domain, gives the feel of this phenomenon and provides accurate physics, but without the possibility of being relevant for ITER cases in the near future due to the computation requirements, in the knowledge that a
$10~\text{cm}\times 10~\text{cm}$
area needs
$10^{7}$
particles per species to obtain a significant result above the noise level. Recent experiments introduce new facts in contradiction with the existing models, e.g. the scaling of rectified DC plasma potential with RF power (Wukitch et al.
Reference Wukitch, LaBombard, Lin, Lipschultz, Marmar and Reinke2009), the radial penetration of the plasma bias (Faudot et al.
Reference Faudot, Heuraux, Colas and Gunn2010) or the nonlinear generation of edge DC currents by RF waves (Bobkov et al.
Reference Bobkov, Braun, Dux, Herrmann, Giannone, Kallenbach, Rohde, Schweinzer, Sips and Zammuto2010), for which a fluid model has been built (Faudot et al.
Reference Faudot, Heuraux, Kubic, Gunn and Colas2013). Although no consensus presently exists over an alternative approach, RF sheath physics needs improvements using first principles.
Figure 1. RF model to describe the emitted wave taken into RF SBCs.
Figure 2. DC model part associated to RF SBCs.
3.2. Modelling and results
The proposed SSWICH model can be thought of as a minimal self-consistent
$\text{RF}+\text{DC}$
approach, able to capture the experimental phenomenology. It was motivated by a closer proximity to first principles and the allowance for DC current circulation in the SOL. The model involves DC plasma quantities as well as harmonic RF waves oscillating at the ion cyclotron wave frequency
${\it\omega}_{0}$
. RF and DC quantities are coupled nonlinearly via two processes representing the sheath physics and implemented as boundary conditions. The plasma medium is still highly anisotropic, with high conductivity in the parallel direction. However, the neighbouring field lines are coupled via RF and DC current exchanges and, if they are sufficiently long, their two extremities can behave independently. Being minimal, the model is not fully complete in its present simplified formulation, with a single-mode description associated with the slow wave that is assumed to be the source of the RF potential
$V_{RF}$
driving the rectification process, following the expression
$$\begin{eqnarray}V_{RF}=\int E_{\Vert }\,\text{d}s,\end{eqnarray}$$
where
$s$
is the curvilinear coordinate associated with a magnetic field line and
$E_{\Vert }$
is the parallel electric field of the slow wave. Owing to the difference between the ion and electron mobilities, a DC potential is generated, on average, over a period of the wave heating frequency. This DC potential modifies the sheath width, which changes the capacitance and the resistance of the sheath, e.g. the electric properties of the sheath against the RF electric field of the slow wave. With this scenario, it is clear that a self-consistent treatment is required to describe the slow wave interaction with the wall, assuming that the magnetic field lines are perpendicular to the wall. This last assumption is not always valid, and the future description with oblique incidence to the wall should be done in more complete versions than those given in Kohno et al. (Reference Kohno, Myra and D’ipolitto2012), where the sheath is associated with a dielectric constant corresponding to vacuum. Both figures 1 and 2 describe a crude model of an antenna, for which we want to introduce sheath boundary conditions. The model equation driving the slow wave is shown in figure 1, and the charge exchange in figure 2. The RF part of the model is shown in figure 1, where
$x,y,z$
denote, respectively, the local radial, poloidal and parallel directions of a flattened tokamak. The simulation is a 3D collection of straight open field lines aligned along
$z$
in the SOL. For the sake of simplicity, Cartesian geometry will be considered here, with parallel (respectively, radial) extension
$L_{\Vert }~(L_{\bot })$
. The equations for the three physical quantities are coupled together by SBC applied on all lateral boundaries, as indicated on figure 1, and can be written as follows:
$$\begin{eqnarray}\displaystyle & {\it\varepsilon}_{\Vert }{\it\Delta}_{\Vert }E_{\Vert }+{\it\varepsilon}_{\bot }{\it\Delta}_{\bot }E_{\Vert }+{\it\varepsilon}_{\Vert }{\it\varepsilon}_{\bot }({\it\omega}_{0}/c)^{2}E_{\Vert }=0\quad \text{(slow wave propagation)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\varepsilon}_{\bot }{\it\Delta}_{\bot }V_{RF}=\pm {\it\varepsilon}_{\Vert }\partial _{\Vert }E_{\Vert }\quad \text{(RF voltage excitation from RF fields)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\sigma}_{\Vert DC}{\it\Delta}_{\Vert }V_{DC}+{\it\sigma}_{\bot DC}{\it\Delta}_{\bot }V_{DC}=0\quad \text{(DC charge conservation w/o sources)}, & \displaystyle\end{eqnarray}$$
${\it\varepsilon}_{\Vert }$
and
${\it\varepsilon}_{\bot }$
are the diagonal elements of the plasma dielectric tensor at the RF pulsation
${\it\omega}_{0}$
,
${\it\sigma}_{\Vert DC}$
is the DC Spitzer parallel conductivity,
${\it\sigma}_{\bot DC}$
is a small effective DC perpendicular conductivity, and
${\it\Delta}_{\Vert }$
is the parallel part of the Laplacian (parallel means aligned to the magnetic field). The SBC expression comes from the fact that, at the RF frequency, the RF current through the sheath is supposed to be mainly a displacement current. In this framework, the sheath is assimilated to a parallel-plate capacitor of width
${\it\delta}$
, filled with a dielectric material of dielectric constant
${\it\varepsilon}_{sh}$
. This description is motivated by the fact that the sheath is a region depleted of electrons (similar to vacuum), as was done in D’Ippolito & Myra (Reference D’Ippolito and Myra2006). The dielectric constant is of the order of
${\it\varepsilon}_{sh}\sim 1$
, while
${\it\varepsilon}_{\Vert }$
in the SOL is highly negative. However, in the presence of high-power waves, the real sheath width is subject to large-amplitude RF oscillations around its mean value. This may affect the effective sheath capacitance. Within these simple assumptions, the RF SBC is linear and reads as (D’Ippolito & Myra Reference D’Ippolito and Myra2006) 
$$\begin{eqnarray}E_{t}^{pl}=\boldsymbol{{\rm\nabla}}_{t}V_{RF}=\boldsymbol{{\rm\nabla}}_{t}({\it\delta}D_{n}^{pl}/{\it\varepsilon}_{sh}).\end{eqnarray}$$
Here,
$E$
is the RF electric field and
$D=\bar{{\bf\varepsilon}}E$
is the RF electric displacement. RF quantities are evaluated at the plasma side of the sheath–plasma interface (superscript
$pl$
), where
$D$
is continuous and
$E$
exhibits a jump, because the dielectric tensor
$\bar{{\bf\varepsilon}}$
changes abruptly. The subscripts
$n$
and
$t$
refer, respectively, to the direction that is locally normal (towards plasma) and transverse to the wall. The sheath width
${\it\delta}$
is allowed to vary spatially, and its radial/poloidal distribution needs to be subsequently determined self-consistently from the DC sheath potential. Assuming that the magnetic field line follows the normal incidence, the SBC expression can be simplified, and, for the sheath width provided by the Child-Langmuir law, assuming a plane geometry and high RF field amplitude, can be written as
$$\begin{eqnarray}\displaystyle & E_{\Vert }={\it\varepsilon}_{sh}V_{RF}/{\it\varepsilon}_{\Vert }{\it\delta}\quad \text{with }{\it\delta}={\it\lambda}_{e}(V_{DC}/T_{e})^{3/4}\quad \text{(sheath capacitance)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle i^{+}\left[1-\exp \left(\frac{V_{b}-V_{DC}}{T_{e}}\right)\right]=\bar{\bar{{\bf\sigma}}}_{DC}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}_{n}V_{DC},\quad & \displaystyle \nonumber\\ \displaystyle & V_{b}=V_{f}+\ln [I_{0}(|V_{RF}/T_{e}|)]\quad \text{(rectification)}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & V_{RF}=0\quad \text{(at both radial ends of lateral boundaries)}. & \displaystyle\end{eqnarray}$$
The main results are that the electric field of the slow wave can be drastically changed, and the SBC behaviour can go from metallic conditions to transparent boundary conditions, depending on the plasma parameter and on the amplitude of the RF potential. This means that the perpendicular DC conductivity is a crucial parameter to determine the radial extent of the DC potential, thus (impacting/explaining/leading to) the power losses and spurious effects induced by the energetic ion flux on the wall (impurity injection, hot spots). The experimental results obtained with the new Faraday screen of Tore Supra ICRH antenna (Jacquot et al. Reference Jacquet, Colas, Mayoral, Arnoux, Bobkov, Brix, Coad, Czarnecka, Dodt, Durodie, Ekedahl, Frigione, Fursdon, Gauthier, Goniche, Graham, Joffrin, Korotkov, Lerche, Mailloux, Monakhov, Noble, Ongena, Petrzilka, Portafaix, Rimini, Sirinelli, Riccardo, Vizvary, Widdowson and Zastrow2014) show that further work has to be done to obtain a relevant modelling for any ITER scenarios, for example, the improvement of the RF sheath at oblique incidence developed in Kohno et al. (Reference Kohno, Myra and D’ipolitto2012).
To reach the goal with a predictive tool that is relevant for ITER and a fusion reactor, some physical issues have to be studied: (i) the role of radial particle turbulent transport on the RF sheath properties needs to be determined (this can be important in the case of an oblique RF sheath near grazing angle to determine the DC potential drop as a function of the RF potential amplitude and particle turbulence flux level); (ii) the secondary emission is also an important parameter to determine the RF sheath properties; and (iii) since the RF sheath takes place in the neighbourhood of the wall or antenna structure, the role of recycling has to be studied, which also allows one to refine the knowledge on the contributions of a ion species mixture or dust to the RF sheath properties. These tasks will have an answer in the framework of the EUROfusion Enabling Research project devoted to the simulation of the ICRH heating including the RF sheath physics.
The knowledge of the total flux, coming from the integration of all the fluxes reaching the wall at any tilting angle of incidence, associated with the rectified potential, is a key point for the safety of the ICRH antenna design. Simulations can provide crucial information on this aspect, if we are able to generate particle turbulent flux with the relevant properties, which are included in the model. There are a few possibilities: (i) to implement an RF antenna in a turbulent code (seems too complex); (ii) to have a kinetic simulation (including all the system within the condition) to simulate the turbulence; or (iii) to introduce a term in the equation of motion able to produce a random walk with an average velocity not equal to zero, which can be associated with turbulent transport. Including dissipation, the natural model corresponds to the Langevin’s equations, instead of the motion equations in PIC codes. Such an operator induces a global particle flux, keeping the energy of the system constant. However, the number of particles required is so high that accurate parallel computing is needed especially for ITER or further fusion reactors, that is, if we want to have a predictive tool, since we have to combine different phenomena (including macroscopic changes of the plasma target) during heating, until we reach an asymptotic state. For such studies, the computation should definitely cover the overall plasma (including the interaction of a turbulent plasma with the wall (Kuhn, Tskhakaya & Tskhakaya Reference Kuhn, Tskhakaya and Tskhakaya2007)), because the RF sheath properties depend on the species contained and, eventually, on accelerated dust in the RF sheath (Ticos, Stoica & Delzanno Reference Ticos, Stoica and Delzanno2012) if we want to estimate the lifetime of an antenna structure in a fusion reactor.
4. Role of density fluctuations on wave heating
4.1. Context
Nowadays, the role of density fluctuations is reconsidered, because the numerical tools, such as full-wave codes, are now able to investigate their effects as a result of improvement in the stability of numerical schemes and their performance using HPCs (da Silva et al. Reference da Silva, Campos Pinto, Després and Heuraux2015). We can also add the displacement of the cut-off layer on the ICRH antenna–plasma coupling (Clairet et al. Reference Clairet, Colas, Heuraux and Lombard2004; Jacquot et al. Reference Jacquot, Colas, Clairet, Goniche, Heuraux, Hillairet, Lombard and Milanesio2013) (which has not been taken into account until now), as well as the depolarization induced by the density fluctuations in fusion devices. In general, the depolarization mechanism should be studied in detail, by itself, applied to fusion plasmas, in order to evaluate the impact on the wave heating or on current drive efficiency. In addition, as shown in Peysson & Decker (Reference Peysson and Decker2014), energy deposition is affected by the density fluctuations at the plasma edge changing the radiation pattern of the antenna. It can also explain the spectral gap if the wavenumber spectrum contains the specific range of parallel wavenumber values. This type of study can have a big impact especially on the neoclassical tearing mode control by electron cyclotron current drive where the beam property requirements are strong (Comisso & Lazzaro Reference Comisso and Lazzaro2010) and particularly when a narrow beam is needed (knowing that density fluctuations at the plasma edge induce a broadening of the launched Gaussian beam as in ITER cases (Peysson & Decker Reference Peysson and Decker2014; Sysoeva et al. Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2014, Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2015)). In these works, the depolarization effect is not taken into account and has not been studied in detail, even though it has been shown that the magnetic shear in a tokamak plays no role (Mazzucato Reference Mazzucato1998). This phenomenon is known in other domains of telecommunication (Tenouxa & Lostanlen Reference Tenouxa and Lostanlen2012), space plasma physics (Jandieri, Ishimaru & Jandieri Reference Jandieri, Ishimaru, Jandieri and Zhukova2011) or laser plasma physics (Javan & Adli Reference Javan and Adli2013), and only recently in fusion in Mirnov et al. (Reference Mirnov, Brower, Hartog and Ding2014). The polarization changes have been studied in astrophysics (Ledenev et al. Reference Ledenev, Tirsky, Tomozov and Zlobec2002). The tools developed in these domains can be useful to study the full impact of turbulence on the wave launched to heat the fusion plasmas. The simulation of the impact of the turbulence on the wave propagation requires full-wave codes, in order to take into account multi-scattering inducing beam broadening, multi-reflection able to modify the polarization, and wave trapping able transiently to enhance the heating power due to local enhancement of the wave amplitude (Heuraux et al. Reference Heuraux, Gusakov, Popov, da Silva and Irzak2010). However, to integrate all the phenomena, a full-wave code has to describe all the components of the wave electric fields at least in two directions (radial, poloidal or toroidal). A 3D simulation (no mode decomposition) over the realistic space domain of ITER or DEMO is too demanding in terms of computation resources at present. We now detail how to proceed in order to build a full-wave simulation to compare theory and numerical results for beam broadening induced by density fluctuations at the plasma edge (in the case of electron cyclotron resonant heating or current drive using O-mode or X-mode), and look at the main results for ITER cases.
4.2. Beam broadening simulation, from numerical scheme to results in the electron cyclotron frequency range
Before simulating such Gaussian beam broadening induced by density fluctuations, we have to choose between different options, bearing in mind the problem to solve the study of the beam shape evolution as a function of the radial coordinate for a given turbulence with known characteristics (wavenumber spectrum, density fluctuations profile) for a given set of plasma parameters (magnetic field, density profile, temperature profile). The choices are explained by answering questions concerning numerical schemes, parallel computing and stability, before presenting the results on beam widening. The inputs can come from measurements in experiments, provided by equilibrium and turbulence codes, or built from analytical expressions. The properties of the launched wave beam, including the polarization (O-mode or X-mode here), can be idealized through the expression of a Gaussian beam, or computed using a complete description of the launching system (waveguide, mirror, horn). In possession of all the required data, we have to define the computational domain and the properties of the border. That is to say, the boundary conditions (PEC, PML, partial absorbing boundary condition (ABC), unidirectional transparent source (UTS) or not for the wave launcher) as a function of the assumptions (single-pass absorption here) or imposed by the knowledge of the plasma properties (with or without mode conversion). At this stage, one needs to choose the numerical scheme, which depends on different options and physical phenomena that must be considered. Here, the options are as follows: (i) ergodic treatment associated with an average of the wave intensity over all samples of turbulence matrices (Helmholtz, frequency domain (Heuraux & da Silva Reference Heuraux and da Silva2012)); or (ii) time averaging of the wave intensity (FDTD Taflove & Hagness Reference Taflove and Hagness2000) for a time-dependent turbulence matrix. It is clear that numerical diagnosis has to be defined at the same time as a function of the studied parameters (wave intensity and poloidal correlation of the wave intensity) and as the control parameters (wave energy conversion, stopping criterion, etc.). For this, a minimal knowledge of the computation parameters is necessary, as are the required number of time steps, computation time per time step, memory requirements, among others. A parallel computing strategy has to be defined for the implementation (determining the algorithm), assuming this numerical scheme is stable and meets the imposed criteria on accuracy, energy conservation, and other constraints. At the end, if it is necessary in the case of a long run, a possible strategy is to dump the data according to the storage possibilities.
Now, we go into the detail to answer the following question: Why do we have to choose a time-dependent algorithm for this study? The choice is driven first by the physical constraints. We want to describe the O- and X-modes in the same plasma conditions. In the frequency domain, the Helmholtz solver for the X-mode is more problematic than the O-mode, because a set of coupled PDEs containing cross-derivatives has to be solved, which introduces numerical difficulties. Therefore, in the frequency domain, we have two codes to implement. One of the codes (for the X-mode) contains numerical difficulties, with no possibility to introduce mode conversion without writing another, even more complex, code. A Maxwell’s equation finite-difference time-domain (FDTD) solver allows both polarizations to be treated with the same algorithm, and only the source term has to be defined, by changing the excited wave fields. The second point concerns the choice of the time averaging. Knowing that the turbulence has its own correlation time
${\it\tau}_{c}$
, the macroscopic behaviour of the wave is exhibited when the time averaging is greater than
${\it\tau}_{c}$
, in practice, at least
$4{\it\tau}_{c}$
. Maxwell’s equations solver also has the advantage of extending to other applications without changing its kernel, even if the computation time is longer than for a dedicated code. After choosing Maxwell’s equation solver, a second question quickly arises: How does one generate the turbulent matrix? The standard method corresponds to a sum of trigonometric functions fulfilling the prescribed correlation time and wavenumber spectrum of the turbulence. The direct summation of modes is irrelevant, owing to the time required to compute a turbulence matrix. The use of a fast Fourier transform is the fastest way to address this, associated with the random phase between each mode (da Silva et al.
Reference da Silva, Heuraux, Gusakov and Popov2010). With such a method, it is possible to introduce a phase relationship in time that can mimic a dispersion relation or a specific turbulence behaviour. Here the goal is to obtain an average value of the wave intensity at the asymptotic state. A way to speed up the computation is to generate a bigger turbulence matrix extended in the poloidal direction (keeping the imposed wavenumber spectrum), and to move it at a given constant speed (less than 10 % of the speed of light) in the poloidal direction (da Silva, Heuraux & Manso Reference da Silva, Heuraux and Manso2006b
). The price to pay is the storage of a huge turbulence matrix, which introduces limitations, depending on the computation resources. Another constraint can appear for the X-mode: the number of time steps needed before the birth of an NI associated with a high level of turbulence. This NI is described in da Silva et al. (Reference da Silva, Campos Pinto, Després and Heuraux2015), as are the possibilities to remove it. An important point to mention is that the numerical scheme used should preserve the energy of the system, in order to conserve the statistical properties we want to study. This kind of algorithm is presented in da Silva et al. (Reference da Silva, Campos Pinto, Després and Heuraux2015). To avoid huge storage, the diagnostic is defined on lines at constant radial values (assuming that the beam spatial evolution is small between each line). This can be optimized to speed up the computation and reduce data storage. As a Gaussian beam is used, it allows evaluation of the radial distance between each of the diagnosis lines. An unexpected problem arises if the number of time steps necessary to obtain an asymptotic state for the wave intensity (more than
$1.5\times 10^{7}$
time steps for an O-mode case and twice as many for X-mode) is too great to compute it in a single run. Thus, the choice was to average over time more than one million time steps, on 20 samples of the turbulence matrix. The total computation time is approximately one week on a single processor. It is clear in such cases that parallel computing is required.
Another question is this: What is the best choice for the parallel computing, MPI and/or OpenMP communication tools? In our case, a parallel version of the full-wave FDTD code REFMULX was developed within the framework of a high-level support team (HLST) project. One of the dimensions of the 2D domain was decomposed, using the MPI standard in such a manner as to maximize the contiguousness of the code’s data representation in memory. The solver routines were modified accordingly, to allow each MPI task to work on a different subdomain. The coupling between the neighbouring subdomains (fixed by the scheme’s numerical stencil) was ensured by the use of ghost cells, which store the data that need to be exchanged via standard point-to-point MPI calls. Following the trend of modern hierarchical architectures on which multiple cores share memory, a threaded-based parallelism was also implemented in the same direction of REFMULX’s domain. Such a hybrid implementation allows several threads to exist within each subdomain (i.e. MPI task). Since the threads share the same memory space, there is no need for inter-node communication nor guard cells, typically leading to better memory usage at no performance cost. This expectation was confirmed with the parallel REFMULX code, whose scaling measurements on the HELIOS-CSC Supercomputer in Japan showed a speedup of two orders of magnitude on a couple of hundred cores for a typical
$1500\times 1000$
grid-count case, compared to the original serial version. Beyond this point, there is not enough work per core to justify the usage of further resources. However, bigger problems (e.g. an ITER case) provide that extra computational work. In other words, new physics studies are potentially enabled by the parallel version of REFMULX. Additionally, in light of the results obtained, extending the modelling to three spatial dimensions now becomes feasible, since all the parallelization techniques developed here will be directly deployed in that case. Since the grid count increases by two or even three orders of magnitude, further extending the domain decomposition might be required. However, it is already clear that some effort needs to be put into implementing parallel input/output (I/O) in the code, which was already measured to be a bottleneck in the biggest 2D test cases ran so far.
Figure 3. Inputs for an ITER case: (a) density and density fluctuation profiles, and (b) wavenumber spectrum of density fluctuations identical in
$x$
and
$y$
directions.
Figure 4. (a) Contour plot superimposed on a snapshot of the positive part of the wave electric field over the averaged intensity in an ITER case and the density profiles in the turbulent zone. (b) Spatial evolution in the plasma zone of the quasi-Gaussian beam width in vacuum (black), in plasma without turbulence (blue), in turbulence
${\it\delta}n=3\,\%$
(green dashed) and in turbulence
${\it\delta}n=5\,\%$
(red).
To improve the stability of our code, a smoothing of the turbulence matrix is carried out just at the plasma edge to avoid stiff index gradients. It is also amazing to see that the optimization done during the parallelization of the code has also improved the code stability. The code remains stable above the turbulence level expected in the ITER case. The code here has to run in unknown conditions, as for DEMO or for data coming from a given turbulence code. Cases in which the turbulence level is higher than the threshold of the NI should trigger a warning. In addition, the turbulence is rescaled below the threshold and the wanted results are extrapolated knowing the nonlinear dependences. The NI threshold value is high enough to evaluate all the nonlinear phenomena found previously, e.g. nonlinear Bragg backscattering (Gusakov, Heuraux & Popov Reference Gusakov, Heuraux and Popov2009) or wave trapping (Heuraux et al. Reference Heuraux, Gusakov, Popov, da Silva and Irzak2010; Gusakov et al. Reference Gusakov, Heuraux, Irzak and Popov2011). The possibility to have nonlinear effects is experimental evidence in the stellarators through measurements of Bragg backscattering (Batanov et al. Reference Batanov, Borzosekov, Kovrizhnykh, Kolik, Konchekov, Malakhov, Petrov, Sarksyan, Skvortsova, Stepakhin and Kharchev2013).
To illustrate beam broadening, ITER-like plasma parameters are used as inputs, including all the possible inhomogeneities described by the theoretical models presented in Sysoeva et al. (Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2015), and those associated with electron cyclotron heating or current drive. The frequency is taken equal to 170 GHz. The density fluctuations (the root mean square value), the density profiles and the turbulence perpendicular wavenumber spectrum used in the analysis are presented in figure 3. The spectrum density, which corresponds to a rescaled wavenumber spectrum coming from experimental data (Gerbaud et al.
Reference Gerbaud, Clairet, Sirinelli, Sabot, Heuraux, Vermare and Leclert2006), is in agreement with GYRO simulations (Casati et al.
Reference Casati, Grangirard, Bourdelle, Hennequin, Gerbaud, Heuraux, Vermare, Aniel, Clairet, Garbet, Imbeaux, Candy and Waltz2009). The radial evolution of the beam width shows the role of turbulence on the beam divergence changes (see figure 4
b), and the enhanced divergence after the crossing of the turbulence zone. The slight inflection of the beam width for the highest level of turbulence corresponds to the fact that part of the beam enters into the PML, showing that we have to be careful when numerical tools are applied automatically. The snapshot of the wave electric field shows that turbulence completely breaks the wavefront, gives a filamentary structure to the wave beam and spreads the wave intensity (see figure 4). After averaging over a very long time (
${\sim}3~{\rm\mu}\text{s}$
) compared to the time of flight (few nanoseconds), it is amazing to see that the wave intensity recovers a Gaussian beam shape, but with a wider angular divergence. This lack of directivity can be reduced if the magnetic field is increased for the same density profile or the density at the centre is reduced. In this case, for a fusion reactor, it will be better to work at higher magnetic field intensities. Looking at the computation time (rescaled into real time (few microseconds)) that is necessary to obtain a Gaussian beam (on average), it is much less than the time needed to measure the wave intensity experimentally (some milliseconds). This fact can explain why the modelling seems to work, but not with the adequate or expected parameters. This is simply due to the omission of the effect on wave propagation of plasma edge turbulence that drives the use of the Gaussian beam parameters evaluated in vacuum, instead of those that consider the density fluctuation effects. The role played by the wavenumber spectrum characteristics is important (see da Silva et al.
Reference da Silva, Heuraux, Gusakov and Popov2010). The dependence on poloidal wavenumber is detailed in Sysoeva et al. (Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2015). Thus, knowledge of the density fluctuation characteristics is required to evaluate these effects, which is not really the case for ITER and even less so for DEMO. Therefore, in future, we need to have a measurement of the turbulence characteristics in real time, to be able to work with the wave heating or current drive systems in an optimal way.
In conclusion, the control of the growing magnetic islands associated with neoclassical tearing modes can be less easy due to beam broadening (Comisso & Lazzaro Reference Comisso and Lazzaro2010). However, tools are now available to evaluate the corrections that need to be introduced in a ray tracing code after the crossing of a high level of turbulence at the plasma edge. The full-wave code was used in Sysoeva et al. (Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2015) to validate the theoretical models. One can note that the full-wave simulation is able to make precise the upper limit of the turbulence level for which the theoretical models are applicable. It is easy to see that, for such studies, the interplay between different codes is important. A turbulence code will be interesting for providing the turbulence characteristics, especially the poloidal wavenumber spectrum and the density fluctuation profile. However, the real-time use of turbulence code data in the full-wave code is risky. The risk stems from the fact that, usually, an interpolation is necessary or a smoothing should be applied. To have a universal tool to adapt the synthetic data coming from a turbulence code, without changing the properties of the turbulence, is not an easy task. This kind of method is, at present, turbulence-code-dependent, probably due to the different boundary conditions of these codes (flux or gradient). The use of a full-wave code can be limited to the plasma edge area, in order to determine the beam properties inside the plasma at a given position for which the density fluctuation level is low enough to keep the beam characteristics constant. Above this position, beam ray tracing can be used (Figini et al. Reference Figini, Decker, Farina, Marushchenko, Peysson, Poli and Westerhof2012) to determine the energy deposition volume.
The role of turbulence at the plasma edge has to be taken into account for any wave heating system, because of changes induced on: (i) the RF SBCs surrounding the antenna structure in the ion cyclotron frequency range due to the sensitivity of coupling as a function of the cut-off position (Clairet et al. Reference Clairet, Colas, Heuraux and Lombard2004); (ii) the effective wavenumber spectrum launched into the plasma by the lower-hybrid grid (Peysson & Decker Reference Peysson and Decker2014); and (iii) the beam properties in the electron cyclotron frequency range (Sysoeva et al. Reference Sysoeva, da Silva, Gusakov, Heuraux and Popov2015).
5. Conclusions and new developments
Over the past year, simulations on wave heating in fusion plasmas have been widely improved. They started being integrated in a platform (e.g. ITM platform; Falchetto et al. Reference Falchetto, Coster, Coelho, Scott, Figini, Kalupin, Nardon, Nowak, Alves, Artaud, Basiuk, Bizarro, Boulbe, Dinklage, Farina, Faugeras, Ferreira, Figueiredo, Huynh, Imbeaux, Ivanova-Stanik, Jonsson, Klingshirn, Konz, Kus, Marushchenko, Pereverzev, Owsiak, Poli, Peysson, Reimer, Signoret, Sauter, Stankiewicz, Strand, Voitsekhovitch, Westerhof, Zok and Zwingmann2014) able to provide the changes of the macroscopic plasma parameters induced by wave heating or current drive. Following an iterative process, we can generally have access to an asymptotic state that allows the behaviour of an experiment to be predicted. Even if the computational resources increase (allowing one to reach a 3D description), important parts of the physical mechanisms are missing, crudely described or treated independently. Knowing the start-up phase of the stellarator W7X, the description of the wave propagation in the stellarator configurations should be rapidly improved, which is not obvious without making crude assumptions (Jucker et al. Reference Jucker, Graves, Cooper, Mellet, Johnson and Brunner2011). Recently, the role of density fluctuations has been reconsidered (as mentioned in Peysson & Decker Reference Peysson and Decker2014), and the impact of turbulence is being actively studied. However, these works ignore the depolarization processes because they are negligible in quiet magnetized plasmas (even if there is magnetic shear (Mazzucato Reference Mazzucato1998)); or the wave reflection is always considered as perfect without changing the polarization, which can be far from real conditions, as shown in another domain Bai, Li & Liu (Reference Bai, Li, Liu, Xu, Shi and Xie2014). Nevertheless, in turbulent plasmas, the picture can be different from space plasmas (Jandieri et al. Reference Jandieri, Ishimaru, Jandieri and Zhukova2011) or laser plasmas (Javan & Adli Reference Javan and Adli2013). Therefore, many questions have to be addressed, such as these: What is the full role of turbulence on the coupling and on the depolarization, especially when the plasma flows intermittently? What is the role of single dense events? How does the time averaging introduced by the coupling system act on it? Is it relevant to consider only the averaged values of the plasma parameters integrated over the time integration of the coupling system? In other words, where is the effective position of the cut-off layer in a turbulent plasma (considering the second-order correction of the effective index associated with the density fluctuations (da Silva et al. Reference da Silva, Heuraux, Manso and Varela2003)) to consider for the evaluation of the coupling resistance? Is a model based on an effective index introducing a constant correction at the second order relevant for ICRH? A similar question on the accessibility condition and conversion layer has just now had a response on the impact of the turbulence on O-X mode conversion (Popov Reference Popov2015), which is part of the O-X-B heating scheme (Laqua et al. Reference Laqua, Erckmann, Hartfuss and Laqua1997). At present, for ICRH, the SBC condition description becomes more and more sophisticated. However, it remains crude, especially in the oblique case, where the magnetic field intercepts the wall far from the normal incidence. The 1D PIC simulations at grazing incidence angles show that the time scale to reach equilibrium is too long, and perpendicular transport as well as recycling and ionization have to be taken into account to determine the RF sheath properties that are useful for SBCs. The ponderomotive effects have been evaluated considering the power density injected into the plasma. However, the crude expressions of the ponderomotive force, the disregard of the role of wave polarizations in its expression, and the essential role of quasi-neutrality compensate its effects (Heuraux, Leclert & Hadjoudj Reference Heuraux, Leclert and Hadjoudj1994). Other nonlinear effects are present, e.g. the parametric instability shown in Gusakov et al. (Reference Gusakov, Popov and Saveliev2014), which need to be considered for a relevant description of a wave heating system.
Figure 5. Contour plots of the modulus of the electric field difference between perturbed and unperturbed density profiles averaged over 200 turbulence matrix samples (scattered wave in Born approximation): (a) in the case of a wavenumber spectrum with a wide range of
$k$
-values covering all the possible scattering processes; (b) same as (a) but for an intermediate range of
$k$
, for which the radiation pattern presents a hole; and (c) same as (a) but with a very narrow spectrum around zero, able to induce forward scattering only for
$k<k_{Airy}$
. In all cases, the density fluctuation amplitude is
${\it\delta}dn/n_{c}=0.1\,\%$
.
A numerical effort should be made to build a fast unconditional stable scheme that is able to preserve the energy of the system for all possible modes that one can have in a time-dependent magnetized turbulent plasma. That is to say, we must have an improved PML or ABCs for the relevant modes, which is not the case at present (Jacquot et al. Reference Jacquot, Colas, Clairet, Goniche, Heuraux, Hillairet, Lombard and Milanesio2013). The SBC description should be improved and valid, at any angle. A way to cross the wave resonance (the real part of the index going to infinity) should be found without artificial damping. Right now, we are far from this goal. In addition, this wave propagation code should integrate a coupling code, or have a full geometrical description of the antenna structure able to support the induced currents on the antenna structure, in order to describe self-consistently the wave emission and the propagation of the emitted wave. No development has been made so far on the self-consistent numerical description of a heating system (ICRH, or LHCD) going from the back side of the antenna to the plasma core, mainly due to too large computational resource needs. That is to say, no integration of a coupling antenna–plasma code into a wave propagation code describing the plasma heating (which should also integrate the macroscopic plasma parameter modifications) has taken place.
The use of the subdomain decomposition should be generalized, since the true need of a full-wave code takes place in a limited area. To finalize, this code should be part of an integrated modelling framework, as planned for the ERCC 3D code that is under development (Coelho et al. Reference Coelho, Äkäslompolo, Dinklage, Kus, Sundén, Blanco, Conway, Hacquin, Heuraux, Lechte, Silva and Sirinelli2013); the role of transport through the evolution of the plasma parameters should be considered; and access to turbulence properties in a self-consistent description must be available. These can be solved properly only if the problem of data exchange is resolved adequately; knowing that, for numerical constraints, the meshes in each module are different, and thus a coherent interpolation should be done. However, it is necessary to be aware that the spectra can be truncated or distorted in this interpolation, or intrinsic, due to the mesh choice. For instance, spatial steps that are too wide are used in turbulence codes to describe properly the small spatial scales involved in the scattering processes really present in an experiment, e.g. in backscattering (Batanov et al. Reference Batanov, Borzosekov, Kovrizhnykh, Kolik, Konchekov, Malakhov, Petrov, Sarksyan, Skvortsova, Stepakhin and Kharchev2013). That is to say, for the turbulence code to be relevant, it should include such possibilities, by choosing a spatial grid able to precisely describe these phenomena, which could be incompatible with the numerical constraints. Three different cases are shown in figure 5 just to illustrate this purpose of the dependences on the scattered field as a function of the turbulence spectrum. In the first case, all spatial scales are present, and the scattered one is more or less isotropic after averaging. Then, part of the spectrum is truncated and, as a consequence, an anisotropy appears in the scattered field pattern, according to the Bragg rule (see Fanack et al. Reference Fanack, Boucher, Heuraux, Leclert, Clairet and Zou1996). In the last one, only forward scattering is possible, showing a narrow lobe for the scattered wave emission. Therefore, this shows that the simulation, in order to be relevant, has to take care of the scales in time and space present in an experiment, which is not always possible, mainly due to the limitations of available numerical resources. As a conclusion, the simulations on wave heating help us to improve the heating systems when we add more physics to it, as this was shown, partially.
Acknowledgements
The authors acknowledge the support of ANR under contracts ANR-12-BS01-0006-01, ANR-12-BS09-0028-02 and ANR-12-JS09-0013-01. Moreover, this work was carried out within the framework of the European Fusion Development Agreement and the French Research Federation for Fusion Studies. Part of this work was carried out using the HELIOS Supercomputer system at the Computational Simulation Centre of the International Fusion Energy Research Centre (IFERC-CSC), Aomori, Japan, under the Broader Approach collaboration between Euratom and Japan, implemented by Fusion for Energy and JAEA. This project has received funding from the European Unions Horizon 2020 research and innovation programme, under grant agreement no. 633053. IST, Lisbon, activities also received financial support from Fundação para a Ciência e Tecnologia through project Pest-OE/SADG/LA0010/2013. It is supported by the European Communities under the Contract of Association between Euratom and CEA. The views and opinions expressed herein do not necessarily represent those of the European Commission.
