1 Introduction

The beam tracing code GRAY performs the computation of the quasi-optical propagation of a Gaussian beam of electron cyclotron waves in a general tokamak equilibrium, and of the power absorption and driven current [1]. The propagation of a general astigmatic Gaussian beam is described within the framework of the complex eikonal approach in terms of a set of “extended” rays that allow for diffraction effects. The absorbed power and the driven current density are computed along each ray solving the fully relativistic dispersion relation for electron cyclotron wave and by means of the neoclassical response function for the current.

The aim of the present note is to document the current version of the GRAY code, by describing the main code features and listing the inputs and outputs. To this goal in sec. 2 the main equations and models used in the GRAY code to compute Gaussian beam propagation, absorption and current drive are summarized shortly. In sec. 3, the details on the code inputs and outputs are presented.

2 Physics

2.1 Coordinate Reference systems

A few sets of coordinate systems are used in the code. The reference system is the right handed Cartesian orthogonal system $(x, y, z)$ with $z$ axis being the tokamak symmetry axis. For the purpose of the physics analysis this coordinate system may be rotated around the $z-$axis so that the $x z$ plane contains the launching point, i.e., $z$ vertical, $x$ radially outward through the port center, and $y$ pointing in the counter clockwise direction when viewed from above.

In addition to the right handed Cartesian orthogonal system specified above, we introduce also a right-handed cylindrical system $(R,φ,Z)$ with transformation from the cylindrical to the Cartesian system given by $x= R\cos φ$, $y=R\sin φ$, $z=Z$.

2.2 Quasi-optical approximation

In the complex eikonal framework, a solution of the wave equation for the electric field is looked for in the form

such that it allows for Gaussian beam descriptions. In eq. 1, $ω$ is the real frequency, $k_0 = ω/c$ the wavevector amplitude in vacuum, ${\mathbf e}({\mathbf x})$ the normalised polarisation (Jones) vector and $E_0({\mathbf x})$ the slowly varying wave amplitude.

The function $S({\mathbf x})$ is the complex eikonal, $S = S_R({\mathbf x}) + i S_I ({\mathbf x})$, in which the real part $S_R({\mathbf x})$ is related to the beam propagation as in the geometric optics (GO), and the imaginary part $S_I({\mathbf x}) (<0)$ to the beam intensity profile shape, as it is apparent writing eq. 1 as

2.3 Gaussian beams

We introduce a reference system $({\bar x},{\bar y},{\bar z})$, in which the $\bar z$ axis is directed along the direction of propagation of the beam and the $\bar x$ axis lies in the horizontal plane (i.e., $z=\text{const}$), and two additional coordinate systems, $(\xi_w,\eta_w)$ and $(\xi_R,\eta_R)$ in the $(\bar x, \bar y)$ plane, rotated by the angles $φ_w$ and $φ_R$, respectively,

In the $(\xi_w,\eta_w)$ and $(\xi_R,\eta_R)$ systems, the axes are aligned with the major and minor axes of the intensity and phase ellipses respectively, and the general astigmatic Gaussian beam in vacuum takes the simple form [2]

Note that a general astigmatic Gaussian beam is described in terms of six parameters: the beam widths $w_{\xi,\eta}$, the phase front curvature radii $R_{c\xi,\eta}$ and the intensity and phase ellipses rotation angles $φ_{w,R}$.

Simple astigmatic beams can be described in terms of 5 parameters only, because the phase and intensity ellipses are aligned, i.e., $φ_w=φ_R\equiv φ$: $w_{\xi,\eta}$, $R_{c\xi,\eta}$, $φ$ or alternatively by the beam waists $w_{0\xi,\eta}$, the waists $\bar z$ coordinates $d_{0\xi,\eta}$, and $φ$, where $R_{c\xi,\eta}$, $w_{\xi,\eta}$ are related to $d_{0\xi,\eta}$, $w_{0\xi,\eta}$ by the following equations:

and $z_{Rj}= k_0 w_{0j}^2/2$ is the Raylegh length. According to eq. 5, a convergent beam (${\bar z} < d_{0j}$) has $R_{cj}<0$, while a divergent beam has $R_{cj}>0$.

2.4 QO beam tracing equations

The beam tracing equations and the algorithm for their solution are described in detail in [1].

The “extended” rays obey to the following quasi-optical ray-tracing equations that are coupled together through an additional constraint in the form of a partial differential equation: $$\begin{aligned} \frac{d {\mathbf x}}{dσ} &= +\frac{∂Λ}{∂\mathbf N} \biggr |_{Λ=0} \\ \frac{d {\mathbf N}}{dσ} &= -\frac{∂Λ}{∂\mathbf x} \biggr |_{Λ=0} \\ \frac{∂ Λ}{∂ {\mathbf N}} &\cdot ∇ S_I = 0 \end{aligned}$$ where the function $Λ ({\mathbf x},{\mathbf k},ω)$ is the QO dispersion relation, which reads

being $\mathbf{b}=\mathbf{B}/B$, $N_\parallel = {\mathbf N} \cdot \mathbf{b}$, and $N_c({\mathbf x}, N_\parallel, ω)$ the solution of the cold dispersion relation for the considered mode.

In GRAY three choices for the integration variable $σ$ are available, i.e.: 1. the arclength along the trajectory $s$, 2. the time $τ=ct$, and 3. the real part of the eikonal function $S_R$.

The default option is the variable $s$ and the QO ray equations become:

2.5 Ray initial conditions

The QO ray equations eq. 7 are solved for $N_T= N_r \times N_\vartheta +1$ rays distributed in order to simulate the Gaussian pattern of an actual antenna, with initial position on a suitable surface at the antenna centered on the beam axis.

The $N_r$ rays are distributed radially up to a “cut-off” radius $\tilde ρ_{max}$ defined as $\tilde ρ_{max}^2=-k_0 S_{I,max}$ such that the beam carries a fraction of the input power equal to $[1-\exp({-2 \tilde ρ_{max}^2)]}$. The $N_\vartheta$ angular rays are distributed at constant electric field amplitude (i.e. at $S_I = \text{const}$). Details are given in [1].

Care must be taken in the proper choice of the integration step to avoid the occurrence of numerical instabilities due to the last equation in the set eq. 7. The value must be tuned with respect to the number of rays (i.e., to the distance between rays).

The code can be run also as a “standard” ray-tracing code, simply imposing $S_I = 0$ in eqns. 6, 7. In this case the initial conditions are given to asymptotically match, for $\vert {\bar z}- d_{0j} \vert \gg z_{Rj}$, the ray distribution used for the QO ray-tracing.

2.6 Launching coordinates and wave vector

The launching coordinates of the central ray of the EC beam will be denoted either as $(x_0, y_0, z_0)$, or $(R_0, φ_0, Z_0)$, depending on the coordinate system used (cartesian or cylindrical) $$\begin{aligned} x_0 &= R_0\cos φ_0 \\ y_0 &= R_0\sin φ_0 \\ z_0 &= Z_0. \end{aligned}$$ and the launched wavevector $\mathbf N$ will have components $(N_{x0}, N_{y0}, N_{z0})$, and $(N_{R0}, N_{φ 0}, N_{Z0})$, related by $$\begin{aligned} N_{x0} &= N_{R0} \cos φ_0 - N_{φ 0} \sin φ_0, \\ N_{y0} &= N_{R0} \sin φ_0 + N_{φ 0} \cos φ_0, \\ N_{z0} &= N_{Z0} \end{aligned}$$

2.7 EC Launching angles ($α,β$)

The poloidal and toroidal angles $α, β$ are defined in terms of the cylindrical components of the wavevector

with $-180° ≤ α ≤ 180°$, and $-90° ≤ β ≤ 90°$, so that

A 1-D scan of launch angle with constant toroidal component at launch ($N_{φ 0}$) is achieved by varying only $α$, keeping $β$ fixed. Injection at $β=0, α=0$ results in a ray launched horizontally and in a poloidal plane towards the machine centre. The above choice is quite convenient to perform physics simulations, since EC results are invariant under toroidal rotation, due to axisymmetry. This convention is the same used for the EC injection angles in ITER [3].

2.8 ECRH and absorption models

The EC power $P$ is assumed to evolve along the ray trajectory obeying to the following equation

where here $α$ is the absorption coefficient

being $N_{\perp w}$ (and $k_{\perp w}$) the perpendicular refractive index (and wave vector) solution of the relativistic dispersion relation for EC waves $$Λ_w = N^2-N_{\parallel}^2-N_{\perp w}^2=0$$

The warm dispersion relation $Λ_w$ is solved up to the desired Larmor radius order either in the weakly or the fully relativistic approximation as described in [4].

Integration of eq. 10 yields the local transmitted and deposited power in terms of the optical depth $τ= \int_0^{s}{α(s') d s'}$ as $$P(s)=P_0 e^{-τ(s)}, \quad \mathrm{and} \quad P_{abs} (s)=P_0 [1-e^{-τ}] ,$$ respectively, being $P_0$ the injected power.

The flux surface averaged absorbed power density $p(ρ)=dP_{abs}/dV$ is computed as the the ratio between the power deposited within the volume $dV$ between two adjacent flux surfaces and the volume itself. At each position along the ray trajectory (parametrized by $s$), the absorbed power density can be written in terms of the absorption coefficient as

$δs$ being the ray length between two adjacent magnetic surfaces, and $δV$ the associated volume.

2.9 EC Current Drive

Within the framework of the linear adjoint formulation, the flux surface averaged EC driven current density is given by

where ${\mathcal R}^*$ is a current drive efficiency, which can be expressed as a ratio between two integrals in momentum space

where $\nu_c=4 \pi n e^4 Λ_c/(m^2 c^3)$ is the collision frequency, with $Λ_c$ the Coulomb logarithm, and $B_m$, $\langle B \rangle$ are the minimum value and the flux surface averaged value of the magnetic field on the given magnetic surface, respectively.

The functions ${\mathcal P}({\mathbf u})$ and $\eta_{\mathbf u}({\mathbf u})$ are the normalized absorbed power density and current drive efficiency per unit momentum ${\mathbf u}={\mathbf p}/mc$ [1]. Note that the warm wave polarisation is used to compute ${\mathcal P}({\mathbf u})$. In the adjoint formulation adopted here, the function $\eta_{\mathbf u}({\mathbf u})$ is written in terms of the response function for the current, and its explicit expression is related to the chosen ECCD model.

The flux surface average driven current density eq. 13 can be written as [1]

and the equation for the current evolution $I_{cd}$ along the ray trajectory as $$\frac{dI_{cd}}{ds} = -{\mathcal R}^*(s)\frac{1}{2 \pi R_J } \frac{dP}{ds},$$ where $R_J(\psi)$ is an effective radius for the computation of the driven current $$\frac{1}{R_J} = \left\langle\frac{1}{R^2}\right\rangle \frac{f(\psi)}{ \langle B\rangle} = \frac{ \langle {B_φ}/{R} \rangle}{ \langle B\rangle}$$ being $f(\psi) =B_φ R$ the poloidal flux function.

2.10 ECCD Models

Two models for $\eta_{\mathbf u}({\mathbf u})$ efficiency in eq. 14 are implemented for ECCD calculations, a Cohen-like module in the high-velocity limit and the momentum conserving model developed by Marushenko. The used Cohen-like module, developed explicitly for GRAY, is described in [1]. The Marushchenko module [5] has been incorporated into GRAY as far as the energy part is concerned, while the pitch-angle part on trapping is based on a local development.

2.10.1 Current density definitions

In GRAY, three outputs for the EC driven current density are given. The EC flux surface averaged driven parallel current density $\langle J_{\parallel}\rangle$, that is the output of the ECCD theory, defined as $$\langle J_{\parallel}\rangle = \left \langle\frac{{\mathbf J}_{cd} \cdot {\mathbf B}}{B} \right \rangle = \frac{\langle {{\mathbf J}_{cd} \cdot {\mathbf B}}\rangle} {{\langle B^2 \rangle/}{\langle B \rangle}}.$$ a toroidal driven current density $J_φ$ defined as

being $δ I_{cd}$ the current driven within the volume $δ V$ between two adjacent flux surfaces, and $δ A$ the poloidal area between the two adjacent flux surfaces, such that the total driven current is computed as $I_{cd}= \int J_φ dA$.

Finally, an EC flux surface averaged driven current density $J_{cd}$ to be compared with transport code outputs

with the $B_{ref}$ value dependent on the transport code, i.e, $B_{ref}=B_0$ for ASTRA and CRONOS, and $B_{ref}={\langle B \rangle}$ for JINTRAC.

The above definitions are related to each other in terms of flux surface averaged quantities, dependent on the equilibrium, i.e.,

2.11 ECRH & CD location and profile characterization

Driven current and absorbed power density profiles, $J_{cd}(ρ)$, $p(ρ)$, can be characterized in term of suitable quantities. In GRAY, two approaches are followed, both available at each computation, that yields the same results in case of almost Gaussian profiles. Here, the flux label $ρ$ denotes the normalized toroidal radius defined as the square root of the toroidal flux normalized to its edge value.

In the first case, the profiles are characterized in terms of three quantities: the peak value of the toroidal current density $J_φ$, the radius $ρ$ corresponding to the peak, and the full profile width at 1/e of the peak value. In addition, the ratio between $J_{cd}/ J_φ$ is computed at the peak radius via eq. 18 for the two $B_{ref}$ choices.

The second approach applies also to non monotonic profiles. Two average quantities are computed for both power and current density profiles, namely, the average radius $\langle ρ \rangle_a$ $(a=p,j)$

and average profile width ${δρ}_a$ defined in terms of the variance as

Factor $\sqrt{8}$ is introduced to match with the definition of the full profile width in case of Gaussian profiles. Consistently with the above average definitions, we introduce suitable peak values $p_{0}$ and $J_{φ 0}$, corresponding to those of a Gaussian profile characterized by eqns. 19, 20 and same total absorbed power $P_{abs}$ and driven current $I_{cd}$

2.12 Mode coupling and reflection at inner wall

The polarisation of the beam is used to compute the coupling to the Ordinary (O) and Extraordinary (X) plasma modes when the beam crosses the vacuum-plasma interface. The fraction of power converted into a mode is given by the coupling coefficient $$c_\text{mode} = (\hat{\mathbf e}_\text{mode}⋅\hat{\mathbf e})²,$$ where $\hat{\mathbf e}_\text{mode},\hat{\mathbf e}$ are the plasma mode and beam Jones vectors, respectively. The mode vector is defined as the eigenvector of the cold plasma dielectric tensor in the low density limit. The beam vector at launch is computed from the polarisation ellipse parameters using the formula:

The following convention is assumed (illustrated in fig. 1):

• $ψ$ is the angle between the $x$ axis and the major axis.

• $χ{=}\tan(b/a)$ where $a,b$ are ellipse major and minor semi-axes, respectively.

• A positive $ψ$ corresponds to an ellipse rotated counterclockwise in the $x∧y$ plane.

• A positive $χ$ corresponds to an ellipse traced clockwise in the $x∧y$ plane with the $z$ axis in the direction of the wave propagation. In other words: negative helicity (projection of spin angular momentum unto wavevector) or left handed wave (IEEE convention).

If the initial polarisation is not specified, 100% coupling to a given mode is assumed.

A model for wave reflection on a smooth surface is included in GRAY. This is used to describe the beam reflection on the inner wall of the tokamak in the cases where only partial absorption occurs at the first pass in the plasma. An ideal conductor is assumed for the reflecting surface, so that the full power of the incident beam is transferred to the reflected one. The vector refractive index ${\mathbf N}_\text{refl}$ and the Jones vector $\hat {\mathbf e}_\text{refl}$ of the reflected beam are $${\mathbf N}_\text{refl} = {\mathbf N}_\text{in} - 2 ({\mathbf N}_\text{in} \cdot \hat {\mathbf n}) \hat {\mathbf n}, \qquad \hat {\mathbf e}_\text{refl} = -\hat {\mathbf e}_\text{in} + 2 (\hat {\mathbf e}_\text{in} \cdot \hat {\mathbf n}) \hat {\mathbf n},$$ being ${\mathbf N}_\text{in}$ and $\hat {\mathbf e}_\text{in}$ the vector refractive index and the Jones vector of the incoming wave, and $\hat {\mathbf n}$ the normal unit vector to the wall at the beam incidence point.

The reflected beam Jones vector is again used to compute the coupling to the plasma modes at the second and successive pass, with potentially $2^{n-1}$ independent modes being traced after n reflections.

Note that the Jones vectors of the ordinary and extraordinary modes are orthogonal w.r.t. the standard Hermitian product: $\hat{\mathbf{e}}_\mathrm{O}⋅\hat{\mathbf{e}}_\mathrm{X}^* = 0$. From eq. 22 it then follows that these relations hold: $$\begin{aligned} ψ_\textrm{O} &= ψ_\textrm{X} + \frac{π}{2} \\ χ_\textrm{O} &= -χ_\textrm{X} \\ 1 &= c_\textrm{O} + c_\textrm{X} \end{aligned}$$ with the latter meaning all the incoming power is coupled to the plasma.

3 Input and output files

This section describes the input data required by GRAY and the output files that are saved and stored for the GRT-161 analysis. In tbl. 1 the assigned unit numbers used for various input and output files in GRAY are listed. Shell scripts allow to run the code for various input data that are varied by means of suitable loops.

3.1 Input Files

To run GRAY, the user must supply a gray_params.data file and two files, containing information about the MHD equilibrium and the kinetic profiles, respectively.

The variables and quantities in gray_params.data are listed in tbls. 6-8, suitably grouped together. The equilibrium information is provided via a G-EQDSK file (with extension .eqdsk), with conventions specified as in [6]. The kinetic profiles are provided in an ASCII file (with extension .prf). The format of the G-EQDSK file is described in sec. 3.2, while the format of the .prf file for the profiles is the following quantities defined in tbl. 9:

  read (98,*) npp
  do i=1,npp
    read(98,*) psin(i),Te(i),ne(i),Zeff(i)
  end do

The beam parameters are read either from the gray_params.data file, or from an ASCII file (with extension .txt), depending on the value of the ibeam parameter as specified in tbl. 6. If the ASCII file is used, the user can supply multiple launching conditions, and/or use a general astigmatic Gaussian beam. The format of the .txt file is the following, for ibeam=1 and ibeam=2 respectively (quantities defined in tbl. 10):

  ! +++ IBEAM=1 +++
  read(97,*) nsteer
  do i=1,nsteer
    read(97,*) gamma(i),alpha0(i),beta0(i),x0mm(i),y0mm(i),z0mm(i), &
                w0xi(i),d0eta(i),w0xi(i),d0eta(i),phiw(i)
  end do

  ! +++ IBEAM=2 +++
  read(97,*) nsteer
  do i=1,nisteer
    read(97,*) gamma(i),alpha0(i),beta0(i),x0mm(i),y0mm(i),z0mm(i), &
               wxi(i),weta(i),rcixi(i),rcieta(i),phiw(i),phir(i)
  end do

3.2 G-EQDSK format

The G EQDSK file provides information on:

1. pressure,
2. poloidal current function,
3. $q$ profile,
4. plasma boundary,
5. limiter contour

All quantities are defined on a uniform flux grid from the magnetic axis to the plasma boundary and the poloidal flux function on the rectangular computation grid. A right-handed cylindrical coordinate system $(R, φ, Ζ)$ is used.

3.2.1 Variables

In order of appearance in the file:

3.2.2 Toroidal Current Density

The toroidal current $J_T$ (A/m²) is related to $P'(ψ)$ and $FF'(ψ)$ through $$J_T = R P'(ψ) + FF'(ψ) / R$$

3.2.3 Example Fortran 77 code

The following snippet can be used to load a G-EQDSK file:

      character*10 case(6)
      dimension psirz(nw,nh),fpol(1),pres(1),ffprim(1),
     .                    pprime(1),qpsi(1),rbbbs(1),zbbbs(1),
     .                    rlim(1),zlim(1)
c
      read (neqdsk,2000) (case(i),i=1,6),idum,nw,nh
      read (neqdsk,2020) rdim,zdim,rcentr,rleft,zmid
      read (neqdsk,2020) rmaxis,zmaxis,simag,sibry,bcentr
      read (neqdsk,2020) current,simag,xdum,rmaxis,xdum
      read (neqdsk,2020) zmaxis,xdum,sibry,xdum,xdum
      read (neqdsk,2020) (fpol(i),i=1,nw)
      read (neqdsk,2020) (pres(i),i=1,nw)
      read (neqdsk,2020) (ffprim(i),i=1,nw)
      read (neqdsk,2020) (pprime(i),i=1,nw)
      read (neqdsk,2020) ((psirz(i,j),i=1,nw),j=1,nh)
      read (neqdsk,2020) (qpsi(i),i=1,nw)
      read (neqdsk,2022) nbbbs,limitr
      read (neqdsk,2020) (rbbbs(i),zbbbs(i),i=1,nbbbs)
      read (neqdsk,2020) (rlim(i),zlim(i),i=1,limitr)
c
 2000 format (6a8,3i4)
 2020 format (5e16.9)
 2022 format (2i5)

3.3 Output Files

GRAY outputs are given in ASCII files at each GRAY execution. tbls. 11, 12 describe the outputs reported in units 7 and 48, respectively.

4 Implementation

4.1 index_rt

The index_rt is a unique index assigned to each combination of beam propagation mode and number of passes into the plasma. Initially index_rt is 1 for the ordinary mode or 2 for the extraordinary mode. Due to the mode mixing, on subsequent passes each beam splits into two modes and the index_rt is updated as:

  index_rt = 2*index_rt + 1  ! for the O mode
  index_rt = 2*index_rt + 2  ! for the X mode

It follows that ordinary(extraordinary) modes respectively have odd(even) indices and the number of passes is given by $\lfloor \log₂(1 + \texttt{index\_rt}) \rfloor$. For example, an index_rt=19 denotes the following chain: $$\begin{aligned} \text{mode:} && O &→ X → O → O \\ \texttt{index\_rt:} && 1 &→ 4 → 9 → 19 \end{aligned}$$

5 Changelog

The code GRAY is maintained under Git revision control in a local repository at ISTP-CNR. The changes in input/output variables and parameters names are listed in tbl. 13, below.

Bibliography