gray/doc/2.physics.md

# Physics {#sec:physics}

## 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$.


## Quasi-optical approximation

In the complex eikonal framework, the a solution of the wave equation for the
electric field is looked for in the form
$$
  {\bf E}({\bf x},t) =
    {\bf e}({\bf x}) E_0({\bf x})
    e^{-i k_0 S({\bf x}) + iωt}
$$ {#eq:eikonal-ansatz}
such that it allows for Gaussian beam descriptions.
In [@eq:eikonal-ansatz], $ω$ is the real frequency, $k_0 = ω/c$ the
wavevector amplitude in vacuum, ${\bf e}({\bf x})$ the polarisation versor and
$E_0({\bf x})$ the slowly varying wave amplitude.

The function $S({\bf x})$ is the complex eikonal, $S = S_R({\bf x}) + i S_I
({\bf x})$, in which the real part $S_R({\bf x})$ is related to the beam
propagation as in the geometric optics (GO), and the imaginary part $S_I({\bf
x}) (<0)$ to the beam intensity profile shape, as it is apparent writing
[@eq:eikonal-ansatz] as
$$
  {\bf E}({\bf x},t) =
    {\bf e}({\bf x}) E_0({\bf x})
    e^{k_0 S_I({\bf x})} e^{-i k_0 S_R({\bf x})+i ωt}
$$ {#eq:efri}


## 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,
$$
\begin{aligned}
  \bar x &= \xi_w \cos φ_w - \eta_w \sin φ_w
          = \xi_R \cos φ_R - \eta_R \sin φ_R \\
  \bar y &= \xi_w \sin φ_w + \eta_w \cos φ_w
          = \xi_R \sin φ_R + \eta_R \cos φ_R
\end{aligned}
$$ {#eq:phiwr}
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
[@gaussian-beam]
$$
  E ({\bf x}) \propto
    \exp{\left[- \left(\frac{{\xi}_w^2}{w_\xi^2}
          +\frac{{\eta}_w^2}{w_\eta^2}\right)
         \right]}
    \exp{\left [-i k_0 \left({\bar z} +\frac{\xi_R^2}{2 R_{c\xi}}
               +\frac{\eta_R^2}{2 R_{c\eta}}\right )
         \right]}.
$$ {#eq:gb}

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:
$$
  \begin{aligned}
    R_{cj} &= [({\bar z}- d_{0j})^2+z_{Rj}^2]/({\bar z}- d_{0j}) \\
    w_j    &= w_{0j} \sqrt{1+({\bar z}- d_{0j})^2/z_{Rj}^2},
  \end{aligned}
$$ {#eq:rciw}

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


## QO beam tracing equations

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

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 {\bf x}}{dσ} &=
    +{∂ Λ \over ∂ {\bf N}} \biggr |_{Λ=0} \\
  \frac{d {\bf N}}{dσ} &=
    -{∂ Λ \over ∂ {\bf x}} \biggr |_{Λ=0} \\
  \frac{∂ Λ}{∂ {\bf N}} &\cdot ∇ S_I = 0
\end{aligned}
$$
where the function $Λ ({\bf x},{\bf k},ω)$ is the QO dispersion
relation, which reads
$$
  Λ  = N² -  N_c²({\bf x}, N_\parallel, ω)
          - |∇ S_I|² + \frac{1}{2}(\mathbf{b} ⋅ ∇ S_I )^2
                       \frac{∂² N_s²}{∂{N_\parallel}²} = 0
$$ {#eq:eqlam}
being $\mathbf{b}=\mathbf{B}/B$, $N_\parallel = {\mathbf N} \cdot \mathbf{b}$,
and $N_c({\bf 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:
$$
\begin{aligned}
  \frac{d{\bf x}}{ds} &=
    +\frac{∂ Λ /∂ {\bf N}}
    {|∂ Λ /∂ {\bf N}|} \biggr |_{Λ=0}  \\
  \frac{d{\bf N}}{ds} &=
    -\frac{∂ Λ /∂ {\bf x}}
    {|∂ Λ /∂ {\bf N}|} \biggr |_{Λ=0} \\
  \frac{∂ Λ}{∂ {\bf N}} &\cdot ∇ S_I = 0
\end{aligned}
$$ {#eq:qort}

## Ray initial conditions

The QO ray equations [@eq:qort] 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 \label{eq:rhomax} $\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 [@gray].

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:qort]. 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 [@eq:eqlam;@eq:qort]. 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.


## 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 $\bf 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}
$$


## EC Launching angles ($α,β$)

The poloidal and toroidal angles $α, β$ are defined in terms of the
cylindrical components of the wavevector
$$
  \begin{aligned}
    N_{R0}  &= -\cosβ \cosα, \\
    N_{φ0}  &= +\sinβ,       \\
    N_{Z0}  &= -\cosβ \sinα
  \end{aligned}
$$ {#eq:ncyl}
with $-180° ≤ α ≤ 180°$, and $-90° ≤ β ≤ 90°$, so that
$$
  \begin{aligned}
    \tanα &= N_{Z0}/N_{R0}, \\
    \sinβ &= N_{φ 0}
  \end{aligned}
$$ {#eq:albt}

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
[@angles].


## ECRH and absorption models

The EC power $P$ is assumed to evolve along the ray trajectory obeying to the
following equation
$$
  \frac{dP}{ds} = -α P,
$$ {#eq:pincta}
where here $α$ is the absorption coefficient
$$
  α = 2 \frac{ω}{c} \frac {{\text{Im}}(Λ_w)}
        {|∂ Λ /∂ {\bf{N}}|} \biggr|_{Λ=0}
    ≈ 4 \frac{ω}{c} {{\text {Im}}(N_{\perp w})}
        \frac {N_{\perp}} {|{∂ Λ}/{∂ {\bf N}|}}  \biggr|_{Λ=0}
    = 2{{\text{Im}}(k_{\perp w})} \frac{v_{g\perp}} v_{g}.
$$ {#eq:alpha}
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 [@dispersion].

Integration of [@eq:pincta] 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
$$
  p = P₀ α(s) e^{-τ(s)} \frac{δs}{δV}
$$ {#eq:pav}
$δs$ being the ray length between two adjacent magnetic surfaces, and $δV$ the
associated volume.


## EC Current Drive

Within the framework of the linear adjoint formulation, the flux surface
averaged EC driven current density is given by
$$
  \langle J_{\parallel}\rangle = {\mathcal R}^* \, p
$$ {#eq:jav}
where
${\mathcal R}^*$ is a current drive efficiency, which can be expressed as a ratio
between two integrals in momentum space
$$
  {\mathcal R}^*= \frac{e}{m c \nu_c} \frac{\langle B \rangle}{B_m}
  \frac{\int{d{\bf u} {\mathcal P}({\bf u}) \,
  \eta_{\bf u}({\bf u})}}{\int{d{\bf u} {\mathcal P}({\bf u}) }}
$$ {#eq:effr}
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}({\bf u})$ and $\eta_{\bf u}({\bf u})$ are the
normalized absorbed power density and current drive efficiency per unit
momentum ${\bf u}={\bf p}/mc$ [@gray].
Note that the warm wave polarisation is used to compute ${\mathcal P}({\bf u})$.
In the adjoint formulation adopted here, the function $\eta_{\bf u}({\bf 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:jav] can be written as
[@gray]
$$
  \langle J_{\parallel}\rangle =
    P_0 α(s) e^{-τ(s)} {\mathcal R}^*(s) \frac{δs}{δV}
$$ {#eq:jrtav}
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}
    = \langle \frac{1}{R^2} \rangle \frac{f(\psi)}{ \langle B\rangle}
    = \frac{ \langle {B_φ}/{R} \rangle}{ \langle B\rangle}
$$
being $f(\psi) =B_φ R$ the poloidal flux function.


## ECCD Models

Two models for $\eta_{\bf u}({\bf u})$ efficiency in [@eq:effr] 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
[@gray]. The Marushchenko module [@marushchenko] 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.


### 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{{\bf J}_{cd} \cdot {\bf B}}{B} \right \rangle
    = \frac{\langle {{\bf J}_{cd} \cdot {\bf B}}\rangle}
           {{\langle B^2 \rangle/}{\langle B \rangle}}.
$$
a *toroidal* driven current density $J_φ$ defined as
\begin{equation}
  J_φ =\frac{δ I_{cd}} {δ A}
  \label{eq:jphia}
\end{equation}
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
$$
  J_{cd} = \frac{\langle {\bf J} \cdot {\bf B} \rangle} {B_{ref}}
$$ {#eq:jcd}
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.,
$$
  \begin{aligned}
    J_φ &= \frac{f(\psi)}{\langle B \rangle}
           \frac{\langle {1/R^2} \rangle}{\langle{1/R} \rangle}
               {\langle J_\parallel \rangle } \\
    J_{cd} &= \frac{\langle B^2 \rangle }{\langle B\rangle B_{ref}}
              \langle J_\parallel \rangle \\
    J_φ &= \frac{B_{ref} f(\psi)}{\langle B^2 \rangle}
           \frac{\langle {1/R^2} \rangle}{\langle{1/R} \rangle} J_{cd}
  \end{aligned}
$$ {#eq:ratj}


## 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:ratj] 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)$
$$
  \langle ρ \rangle_p = \frac{\int dV ρ p(ρ)}{\int dV  p(ρ)} , \qquad
  \langle ρ \rangle_j = \frac{\int dA ρ | J_{φ}(ρ)|}  {\int dA |J_{φ}(ρ)|}
$$ {#eq:rav}
and average profile width ${δρ}_a$ defined in terms of the variance as
$$
  δ ρ_a = 2 \sqrt{2} \langle δ ρ \rangle_a
  \qquad  \mathrm {with } \qquad
  \langle δ ρ \rangle_a^2 = \langle ρ^2 \rangle_a-(\langle ρ \rangle_a)^2
$$ {#eq:drav}
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 [@eq:rav;@eq:drav] and same total absorbed power $P_{abs}$ and
driven current $I_{cd}$
$$
  p_0 = \frac{2}{\sqrt{\pi}} \frac{P_{abs}}{{ δ ρ}_p
        \left ({dV}/{d ρ}\right)_{\langle ρ \rangle_p}},
  \qquad
  J_{φ0} = \frac{2}{\sqrt{\pi}} \frac{I_{cd}}{{ δ ρ}_j
           \left ({dA}/{d ρ}\right)_{\langle ρ \rangle_j}}.
$$ {#eq:pjgauss}


## Reflection at inner wall and polarisation

A model for wave reflection on a smooth surface is included in GRAY. This is
used to describe 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
${\bf N}_{\rm{refl}}$ and the unit electric field $\hat {\bf e}_{\rm{refl}}$ of
the reflected wave are
\begin{equation}
  {\bf N}_{\rm{refl}} =
    {\bf N}_{\rm{in}} - 2 ({\bf N}_{\rm{in}}
    \cdot \hat {\bf n}) \hat {\bf n}, \qquad
  \hat {\bf e}_{\rm{refl}} =
    -\hat {\bf e}_{\rm{in}}
    + 2 (\hat {\bf e}_{\rm{in}} \cdot \hat {\bf n}) \hat {\bf n},
\end{equation}
being ${\bf N}_{\rm{in}}$ and $\hat {\bf e}_{\rm{in}}$ the vector refractive
index and the unit electric field of the incoming wave, and $\hat {\bf n}$ the
normal unit vector to the wall at the beam incidence point.

The Stokes parameter for the unit electric vector $\hat {\bf e}$ in vacuum are
defined in the beam reference system $({\bar x},{\bar y},{\bar z})$ as
$$
  \begin{aligned}
    I &= \vert \hat e_{\bar x} \vert^2 + \vert \hat e_{\bar y} \vert^2 = 1 \\
    Q &= \vert \hat e_{\bar x} \vert^2 - \vert \hat e_{\bar y} \vert^2 \\
    U &= 2 \cdot {\rm Re} (\hat e_{\bar x} \hat e_{\bar y}^*) \\
    V &= 2 \cdot {\rm Im} (\hat e_{\bar x} \hat e_{\bar y}^*).
  \end{aligned}
$$ {#eq:stokes}
Alternatively, the two angles $\psi_p$ and $\chi_p$ can be used:
$$
  \begin{aligned}
    Q &= \cos {2 \psi_p} \cos {2 \chi_p} \\
    U &= \sin {2 \psi_p} \cos {2 \chi_p} \\
    V &= \sin {2 \chi_p}
  \end{aligned}
$$
which define respectively the major axis orientation and the ellipticity of the
polarisation ellipse. The polarisation parameters of the reflected wave are
used to compute the coupling with the Ordinary (OM) and Extraordinary (XM)
modes at the vacuum-plasma interface before the calculation of the second pass
in the plasma. At the second pass both modes are traced, taking into account
that the power fraction coupled to each mode is
$$
  P_{\rm O,X} =
    \frac{P_{\rm in}}{2}
    (1 + Q_{\rm in} Q_{\rm O,X}
       + U_{\rm in} U_{\rm O,X}
       + V_{\rm in} V_{\rm O,X}).
$$

Note that the polarisation vectors of OM and XM form an orthogonal base:
$\psi_{p{\rm O}}=\psi_{p{\rm X}}+\pi/2$, $\chi_{p{\rm O}}=-\chi_{p{\rm X}}$ and
as a consequence $Q_{\rm O}=-Q_{\rm X}$,  $U_{\rm O}=-U_{\rm X}$,  and $V_{\rm
O}=-V_{\rm X}$, so that $P_{\rm O} + P_{\rm X} = P_{\rm in}$, i.e. all the
incoming power is coupled to the plasma.