\documentclass[a4paper]{article}

\usepackage{amsfonts}
\usepackage{amsmath}
\usepackage[english]{babel}
\usepackage[squaren]{SIunits}
\usepackage[latin1]{inputenc}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{epsfig}
\usepackage{subfig}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Renew commands
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newcommand{\de}{\partial}
\newcommand{\dd}{\mathrm{d}}

\newcommand{\er}{\mathbf{e_r}}
\newcommand{\ez}{\mathbf{e_z}}
\newcommand{\et}{\mathbf{e_\theta}}
\newcommand{\ep}{\mathbf{e_\phi}}
\newcommand{\B}{\mathbf{B}}
\newcommand{\E}{\mathbf{E}}
\renewcommand{\j}{\mathbf{j}}
\newcommand{\psih}{\hat{\psi}}

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}

\newcommand{\ba}{\begin{align}}
\newcommand{\ea}{\end{align}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\didi}[2]{\ensuremath{\frac{\dd{#1}}{\dd{#2}}}}
\newcommand{\dede}[2]{\ensuremath{\frac{\partial{#1}}{\partial{#2}}}}
\newcommand{\ddedde}[2]{\ensuremath{\frac{\partial^2{#1}}{\partial{#2}^2}}}
\newcommand{\psihat}{\hat{\psi}}
\newcommand{\Tehat}{\hat{T_e}}
\newcommand{\rhohat}{\hat{\rho}}
\newcommand{\psihb}{\hat{\boldsymbol{\psi}}}

\newcommand{\cmd}[1]{\textsf{#1}}

\newcommand{\Simulink}{\textsc{Simulink}}
\newcommand{\Matlab}{\textsc{Matlab}}

% units
\newcommand{\eV}{\electronvolt}
\newcommand{\keV}{\kilo\electronvolt}

% fonts
\newcommand{\bs}[1]{\boldsymbol{#1}}


%%%%%%%%%%%%%%%%%%%%%%

\title{RAPTOR equation summary}
\author{F. Felici, f.felici@tue.nl, \\A. Teplukhina anna.teplukhina@epfl.ch\\J. Citrin J.Citrin@differ.nl}

\bibliographystyle{unsrt}

\begin{document}
\maketitle


\section{Definitions \& useful equations}
RAPTOR uses right-handed $(R,\phi,Z)$ coordinate system and COCOS=11 with $\sigma_{I_p}=\sigma_{B_0} = 1$, since $I_p$ and $B_0$ are always assumed positive \cite{Sauter2012}.


\noindent Furthermore:\\
$\rho = \sqrt{\Phi / \pi B_0}, \quad  \rho_b = \sqrt{\Phi_b / \pi B_0}, \quad \hat{\rho} = \sqrt{\Phi/\Phi_b}$\\
$V'_{\rho} = \dede{V}{\rho}$, $V'_{\hat \rho} = \dede{V}{\hat \rho}$ \\
$g_0 = \left\langle \nabla V \right\rangle$,
$g_1 = \left\langle (\nabla V)^2 \right\rangle$,
$g_2 = \left\langle \frac{(\nabla V)^2}{R^2}\right\rangle$,
$g_3 = \left\langle \frac{1}{R^2}\right\rangle$.\\
$F = RB_\phi$ \\
$j_{ni} = \frac{\langle \mathbf{j_{ni}} \cdot \mathbf{B}\rangle}{B_0}$

From the chain rule, we can write the useful relation, for a quantity $A=A(t,\psi(\Phi))$
\be
\left.\dede{}{t}\right|_\psi A(t,\psi) = \left.\dede{A(\Phi,t)}{t}\right|_{\Phi} + \left.\dede{A(\Phi,t)}{\Phi}\right|_{t}\left.\dede{\Phi}{t}\right|_{\psi} 
\ee

\section{Flux transport}
Flux surface averaged ohm's law in general form:

\begin{align}
 \sigma_\parallel \left. \dede{\Phi}{t} \right|_\psi
 &= \frac{F^2}{2\pi \mu_0} \dede{}{\psi}\left[ \frac{g_2 }{F} \dede{\psi}{V} \right]
 -  \dede{V}{\psi} \langle \j_{ni} \cdot \B \rangle
\end{align}

We can then re-write this with $\psi$ or $\Phi$ or $\rho$ as independent variable:

\subsection{$\Phi$}
$\hat\Phi = \Phi/\Phi_b$, always [0,1]
\begin{align}
\sigma_\parallel \left( \left. \dede{\psi}{t} \right|_{\hat{\Phi}} - \frac{\hat\Phi \dot{\Phi_b}}{\Phi_b}\dede{\psi}{\hat\Phi} \right)
&= \frac{F^2}{4\pi^2 \mu_0\Phi_b^2} \dede{}{\hat\Phi}\left[ g_2 g_3 \dede{\psi}{\hat\Phi} \right]
 -  \frac{1}{\Phi_b}\dede{V}{\hat\Phi} \langle \j_{ni} \cdot \B \rangle
\end{align}
Which is equal to E. Fable's notes $\cite{Fable2013}$.

\subsection{$\rho$}
$\rho = \sqrt{\Phi/\pi B_0}, d\Phi = 2\pi B_0 \rho d\rho$, and $\hat\rho = \sqrt{\hat\Phi} = \sqrt{\Phi/\Phi_b}, d\Phi = 2\Phi_b \hat\rho d\hat\rho$
\begin{align}
\sigma_\parallel \left( \left. \dede{\psi}{t} \right|_\rho - \frac{\rho \dot B_0}{2B_0}\dede{\psi}{\rho}\right)
&= \frac{F^2}{16\pi^4 \mu_0 B_0^2 \rho} \dede{}{\rho}\left[ \frac{g_2 g_3}{\rho} \dede{\psi}{\rho} \right]
 - \frac{1}{2\pi \rho} \dede{V}{\rho} \frac{\langle \j_{ni} \cdot \B \rangle}{B_0}
 \end{align}
 Which is equal to $\cite{Pereverzev2002}$ since $\frac{G_2}{J} = \frac{R_0}{16\pi^4} \frac{g_2g_3}{\rho}$

 \begin{align}
 \sigma_\parallel \left( \left. \dede{\psi}{t} \right|_{\hat\rho} - \frac{\hat \rho \dot \Phi_b}{2\Phi_b}\dede{\psi}{\hat \rho}\right)
&= \frac{F^2}{16\pi^2 \mu_0 \Phi_b^2 \hat\rho} \dede{}{\hat \rho}\left[ \frac{g_2 g_3}{\hat \rho} \dede{\psi}{\hat \rho} \right]
 - \frac{B_0}{2\Phi_b \hat \rho} \dede{V}{\hat \rho} \frac{\langle \j_{ni} \cdot \B \rangle}{B_0}
\end{align}

This can be rewritten to have all terms correspond to a flux-surface-averaged $\j\cdot\B$ term:
\begin{align}\underbrace{
 \sigma_\parallel \left( \frac{2\Phi_b\hat\rho}{V'_{\hat \rho}}
 \left. \dede{\psi}{t} \right|_{\hat\rho} - \frac{\hat \rho^2 \dot
   \Phi_b}{V'_{\hat \rho}}\dede{\psi}{\hat \rho}\right)}_{\langle \mathbf{j_{\Omega}} \cdot \mathbf{B} \rangle}
&= \underbrace{\frac{F^2}{8\pi^2 \mu_0 \Phi_b V'_{\hat \rho}} \dede{}{\hat \rho}\left[ \frac{g_2 g_3}{\hat \rho} \dede{\psi}{\hat \rho} \right]}_{\langle \mathbf{j}\cdot \mathbf{B} \rangle}
 - \langle \j_{ni} \cdot \B \rangle
\label{eqn:jdotB}
\end{align}

\subsection{Finite Element forms}
\be
 \sigma_\parallel \left( \left. \dede{\psi}{t} \right|_{\hat\rho} - \frac{\hat \rho \dot \Phi_b}{2\Phi_b}\dede{\psi}{\hat \rho}\right)
= \frac{F^2}{16\pi^2 \mu_0 \Phi_b^2 \hat\rho} \dede{}{\hat \rho}\left[ \frac{g_2 g_3}{\hat \rho} \dede{\psi}{\hat \rho} \right]
 - \frac{B_0}{2\Phi_b \hat \rho} V'_{\hat \rho} j_{ni}
\ee
is rewritten to eliminate the term in front of the second order derivative, yielding
\be
 \frac{16\pi^2 \mu_0 \Phi_b^2 \hat\rho \sigma_{\parallel}}{F^2} \left. \dede{\psi}{t} \right|_{\hat\rho}
=
  \dede{}{\hat \rho}\left[ \frac{g_2 g_3}{\hat \rho} \dede{\psi}{\hat \rho} \right]
  +  \dot \Phi_b\frac{8\pi^2 \mu_0 \Phi_b \hat\rho^2 \sigma_{\parallel}}{F^2}\dede{\psi}{\hat \rho}
 -   \frac{8\pi^2 \mu_0 \Phi_bV'_{\hat \rho}}{F^2} \langle \j_{ni} \cdot \B \rangle
\ee
  or
  \be
  m_\psi \dede{\psi}{t}=  a_\psi \dede{\psi}{\hat\rho} + \dede{}{\hat\rho} d_{\psi} \dede{\psi}{\hat\rho}  + s_\psi
  \ee
  with
\begin{align}
  m_\psi &=16\pi^2 \mu_0 \hat\rho  \frac{\Phi_b^2 \sigma_{\parallel}}{F^2} \\
  a_\psi &= 8\pi^2 \mu_0 \dot \Phi_b \Phi_b \frac{\sigma_{\parallel}\hat\rho^2}{F^2}\\
  d_\psi &=  \frac{g_2 g_3}{\hat \rho}  \\
  s_\psi &=  -8\pi^2 \mu_0 \Phi_b\frac{V'_{\hat \rho}}{F^2}  \langle\j_{ni} \cdot \B \rangle \\
\end{align}
Now write $\psi$ as a sum of spatial basis functions
\be
\psi(\rho,t) = \sum_{\alpha = 1}^{n_{sp}} \Lambda_\alpha(\hat\rho) \hat{y}_\alpha(t)
\ee
then the weak form, after projection on $\Lambda_b$ and integration by parts is
\begin{align}
\sum_{\alpha=1}^{n_{sp}} \didi{\hat y_\alpha (t)}{t} \int_0^1 m\Lambda_\beta \Lambda_\alpha \dd \hat \rho &=
\sum_{\alpha=1}^{n_{sp}} \hat y_\alpha \int_0^1 a_\psi \Lambda_\beta \dede{\Lambda_\alpha}{\hat\rho} \dd \hat \rho\\
- \sum_{\alpha=1}^{n_{sp}} \hat y_\alpha \int_0^1 d_\psi \dede{\Lambda_\beta}{\hat\rho} \dede{\Lambda_\alpha}{\hat\rho} \dd \hat \rho&
 +\left[ d_\psi \Lambda_\beta \dede{\psi}{\hat\rho} \right]_0^1
+\int_0^1 \Lambda_\beta s_\psi \dd \hat\rho
\end{align}
which gives the matrix form
\be
\mathbf{M_\psi}\didi{\hat\psi}{t} = (-\mathbf{D}_\psi+ \mathbf{A}_\psi) \hat{\psi} + \mathbf{l} + \mathbf{s}
\ee

The boundary term $\mathbf{l}$ contains only the last element
\be
\left. d_\psi  \Lambda_\beta \dede{\psi}{\hat\rho} \right]_{\rho=1} = \left.  \frac{g_2g_3}{\hat\rho}\dede{\psi}{\hat\rho} \right]_{\rho=1} = \left.\frac{16\pi^3\mu_0\Phi_b}{ F}\right]_{\rho=1} I_p \label{eqn:psi_BC_Ipl}
\ee

\section{Thermal transport}
Equation for thermal transport
\be
\frac{3}{2}(V'_{\rho})^{-5/3}(\left.\dede{}{t}\right|_{\rho}-\frac{\dot B_0}{2B_0}\dede{}{\rho}\rho)[(V'_{\rho})^{5/3}n_eT_e]+\frac{1}{V'_{\rho}}\dede{}{\rho}(-\frac{g_1}{V'_{\rho}}n_e\chi_e\dede{T_e}{\rho}+\frac{5}{2}T_e{\Gamma_e}g_0)=P_e
\ee
or
\be
\frac{3}{2}(V'_{\hat\rho})^{-5/3}(\left.\dede{}{t}\right|_{\hat\rho}-\frac{\dot\Phi_b}{2\Phi_b}\dede{}{\hat\rho}\hat\rho)[(V'_{\hat\rho})^{5/3}n_eT_e]+\frac{1}{V'_{\hat\rho}}\dede{}{\hat\rho}(-\frac{g_1}{V'_{\hat\rho}}n_e\chi_e\dede{T_e}{\hat\rho}+\frac{5}{2}T_e{\Gamma_e}g_0)=P_e
\ee
Where $\Gamma_e$ is the electron particle flux defined in \eqref{eqn:Gamma}, in units of $\meter^{-2}\second^{-1}$. The equation is rewritten as
\begin{align}
\begin{split}
\frac{3}{2}V'_{\hat \rho}  n_e\dede{T_e}{t} = \dede{}{\hat \rho} \left( \left[\frac{3}{2} \frac{\dot \Phi_b}{2 \Phi_b}  V'_{\hat \rho} \hat \rho n_e - \frac{5}{2} \Gamma_e g_0 \right] T_e \right)+\dede{}{\hat \rho} \left(\frac{g_1}{V'_{\hat \rho}} n_e\chi_e \dede{T_e}{\hat \rho}\right) \\ + \left( \frac{\dot \Phi_b}{2 \Phi_b} \hat \rho n_e \dede{V'_{\hat \rho}}{\hat \rho}- \frac{5}{2}n_e \dede{V'_{\hat \rho}}{t} - \frac{3}{2} V'_{\hat \rho} \dede{n_e}{t} \right) T_e+V'_{\hat \rho} P_e
\end{split}
\end{align}

\subsection{Finite Element forms}
Equation in matrix form
\be
 m_{T_e} \dede{T_e}{t}=  \dede{}{\hat\rho}(a_{T_e} T_e) + \dede{}{\hat\rho} d_{T_e}\dede{T_e}{\hat\rho}  + h_{T_e} T_e + s_{T_e}
\ee
with
\begin{align}
m_{T_e} &= \frac{3}{2}V'_{\hat \rho} n_e\\
a_{T_e} &= \frac{3}{2} \hat \rho n_e V'_{\hat \rho} \frac{\dot \Phi_b}{2 \Phi_b} - \frac{5}{2} \Gamma_e g_0\\
d_{T_e} &= \frac{g_1}{V'_{\hat \rho}}n_e \chi_e\\
h_{T_e} &= \frac{\dot \Phi_b}{2 \Phi_b} \hat \rho n_e \dede{V'_{\hat \rho}}{\hat \rho}-\frac{5}{2} n_e \dede{V'_{\hat \rho}}{t} - \frac{3}{2} V'_{\hat \rho} \dede{n_e}{t}\\
s_{T_e} &= V'_{\hat \rho} P_e
\end{align}

Now write $T_e$ as a sum of spatial basis functions
\be
T_e(\hat \rho,t) = \sum_{\alpha = 1}^{n_{sp}} \Lambda_\alpha(\hat\rho) \hat{z}_\alpha(t)
\ee
then the weak form, after projection on $\Lambda_b$ and integration by parts is
\begin{align}
\sum_{\alpha=1}^{n_{sp}} \didi{\hat z_\alpha (t)}{t} \int_0^1 m_{T_e}\Lambda_\beta \Lambda_\alpha \dd \hat \rho &=
-\sum_{\alpha=1}^{n_{sp}} \hat z_\alpha \int_0^1 a_{T_e} \dede{\Lambda_\beta}{\hat\rho} \Lambda_\alpha \dd \hat \rho +\left[ a_{T_e} \Lambda_\beta T_e \right]_0^1\\
- \sum_{\alpha=1}^{n_{sp}} \hat z_\alpha \int_0^1 d_{T_e} \dede{\Lambda_\beta}{\hat\rho} \dede{\Lambda_\alpha}{\hat\rho} \dd \hat \rho
 &+\left[ d_{T_e} \Lambda_\beta \dede{T_e}{\hat\rho} \right]_0^1 \\
&+ \sum_{\alpha=1}^{n_{sp}} \hat z_\alpha \int_0^1 h_{T_e} \Lambda_\beta \Lambda_\alpha \de \hat \rho + \int_0^1 \Lambda_\beta s_{T_e} \dd \hat\rho
\end{align}
which gives the matrix form
\be
\mathbf{M}_{T_e}\didi{\hat T_e}{t} = (-\mathbf{A}_{T_e} - \mathbf{D}_{T_e} + \mathbf{H}_{T_e})
\hat{T_e} + \mathbf{l} + \mathbf{s}
\ee
 with the boundary term
\be
l = \left. a_{T_e} \Lambda_\beta T_e \right]_{\hat \rho = 1} + \left. d_{T_e} \Lambda_\beta \dede{T_e}{\hat\rho} \right]_{\hat \rho = 1} = \\ \left. \left( \left(\frac{3}{2} \hat \rho n_e V'_{\hat \rho} \frac{\dot \Phi_b}{2 \Phi_b} - \frac{5}{2} \Gamma_e g_0 \right) T_e + \frac{g_1}{V'_{\hat \rho}}n_e \chi_e \dede{T_e}{\hat\rho} \right) \right]_{\hat \rho = 1}
\ee

\section{Particle transport}
Starting from Hinton\&Hazeltine \cite{Hinton1976} equation (7.30) we find the particle transport equation (for species $s$), where $\dede{}{t}$ is intended at fixed $\Phi$ surface: 
\be
\left.\dede{}{t}\right|_{\Phi}(V_\rho'n_s) + \dede{}{\rho}\left(V_\rho'\Gamma_s \right) = V'S_s
\ee
Using the definition of $\rho$ and allowing a time-varying $B_0$, and expressing $\Gamma_s$ (in $\meter^{-2}\second^{-1}$) as

\be
\Gamma_s = -\frac{g_1}{{V'_\rho}^2} D_s \dede{n_s}{\rho} + \frac{g_0}{V'_\rho} V_s n_s,
\label{eqn:Gamma}
\ee
we arrive at
\be
\frac{1}{V'_{\rho}}(\left.\dede{}{t}\right|_\rho-\frac{\dot B_0}{2B_0}\dede{}{\rho}\rho)[(V'_{\rho})n_s]+\frac{1}{V'_{\rho}}\dede{}{\rho}(-\frac{g_1}{V'_{\rho}}D_s\dede{n_s}{\rho}+g_0V_sn_s)=S_s
\ee
or
\be
\frac{1}{V'_{\hat\rho}}(\left.\dede{}{t}\right|_{\hat\rho}-\frac{\dot \Phi_b}{2\Phi_b}\dede{}{\hat\rho}\hat\rho)[(V'_{\hat\rho})n_s]+\frac{1}{V'_{\hat\rho}}\dede{}{\hat\rho}(-\frac{g_1}{V'_{\hat\rho}}D_s\dede{n_s}{\hat\rho}+g_0V_sn_s)=S_s
\ee

Where $D_s$ is the particle diffusion, and $V_s$ the particle pinch (positive $V_s$ is outward pinch). 
\begin{align}
\begin{split}
V'_{\hat \rho}\dede{n_s}{t} & = \dede{}{\hat \rho} \left( \left[\frac{\dot \Phi_b}{2 \Phi_b}  V'_{\hat \rho} \hat \rho-V_sg_0\right] n_s\right) +\dede{}{\hat \rho}\left(\frac{g_1}{V'_{\hat \rho}} D_s \dede{n_s}{\hat \rho}\right) - \dede{V'_{\hat \rho}}{t} n_s +V'_{\hat \rho} S_s
\end{split}
\end{align}

\subsection{Finite Element forms}
Equation in matrix form
\be
m_{n_s} \dede{n_s}{t}=  \dede{}{\hat\rho} (a_{n_s} n_s) + \dede{}{\hat\rho} d_{n_s}\dede{n_s}{\hat\rho}  + h_{n_s} n_s + s_{n_s}
\ee
with
\begin{align}
m_{n_s} &= V'_{\hat \rho} \\
a_{n_s} &= \frac{\dot \Phi_b}{2 \Phi_b}  V'_{\hat \rho} \hat \rho-V_sg_0\\
d_{n_s} &= \frac{g_1}{V'_{\hat \rho}} D_s\\
h_{n_s} &= - \dede{V'_{\hat \rho}}{t}\\
s_{n_s} &= V'_{\hat \rho} S_s
\end{align}

Now write $n_s$ as a sum of spatial basis functions
\be
n_s(\hat \rho,t) = \sum_{\alpha = 1}^{n_{sp}} \Lambda_\alpha(\hat\rho) \hat{y}_\alpha(t)
\ee
then the weak form, after projection on $\Lambda_b$ and integration by parts is
\begin{align}
\sum_{\alpha=1}^{n_{sp}} \didi{\hat y_\alpha (t)}{t} \int_0^1 m_{n_s}\Lambda_\beta \Lambda_\alpha \dd \hat \rho &=
-\sum_{\alpha=1}^{n_{sp}} \hat y_\alpha \int_0^1 a_{n_s} \dede{\Lambda_\beta}{\hat\rho} \Lambda_\alpha \dd \hat \rho +\left[ a_{n_s} \Lambda_\beta n_s \right]_0^1\\
- \sum_{\alpha=1}^{n_{sp}} \hat y_\alpha \int_0^1 d_{n_s} \dede{\Lambda_\beta}{\hat\rho} \dede{\Lambda_\alpha}{\hat\rho} \dd \hat \rho
&+\left[ d_{n_s} \Lambda_\beta \dede{n_s}{\hat\rho} \right]_0^1 \\
&+ \sum_{\alpha=1}^{n_{sp}} \hat y_\alpha \int_0^1 h_{n_s} \Lambda_\beta \Lambda_\alpha \de \hat \rho + \int_0^1 \Lambda_\beta s_{n_s} \dd \hat\rho
\end{align}
which gives the matrix form
\be
\mathbf{M}_{n_s}\didi{\hat n_s}{t} = (- \mathbf{A}_{n_s}-\mathbf{D}_{n_s} + \mathbf{H}_{n_s})
\hat{n_s} + \mathbf{l} + \mathbf{s}
\ee
with the boundary term
\be
l =  \left. a_{n_s} \Lambda_\beta n_s \right]_{\hat \rho = 1} + \left. d_{n_s} \Lambda_\beta \dede{n_s}{\hat\rho} \right]_{\hat \rho = 1} = \left. \left( \left(\frac{\dot \Phi_b}{2 \Phi_b}  V'_{\hat \rho} \hat \rho-V_sg_0 \right) n_s + \frac{g_1}{V'_{\hat \rho}} D_s \dede{n_s}{\hat\rho} \right) \right]_{\hat \rho = 1}
\ee

\subsection{Choice of boundary conditions and constraints}
\begin{enumerate}
\item Edge source from input.
\item Change edge source to match input line average density.
\item Internal boundary condition density from input.
\item Change internal boundary condition to match reference line average density.
\end{enumerate}
In all cases:  transport coefficients from model.

For full discharge simulations: allow to switch BC type as a function of time.



\subsection{Choice of particle species simulations}
There is freedom to choose between predicting or prescribing either evolution of either the electron density or ion density.

Particle densities and $Z_{eff}$ are linked by the relations:
\begin{equation}\label{eqn:quasin}
n_e=\sum_jn_jZ_j
\end{equation}
and
\begin{equation}
\label{eqn:zeff}
Z_{eff}=\sum_j\frac{n_jZ_j^2}{n_e}
\end{equation}

If both $n_e$ and $n_i$ are known (either having been prescribed or predicted by solving a transport equation), then they must either exactly satisfy quasineutrality (with $Z_{eff}$ computed trivially from \eqref{eqn:zeff}), or one must introduce an additional impurity species.

If one introduces one impurity species with density $n_1$ and known $Z_1$, then one can prescribe two out of these three: $n_e, n_i, n_1$ and the third is calculated by \eqref{eqn:quasin}. $Z_{eff}$ is trivially calculated from \eqref{eqn:zeff}.

If one also prescribes $Z_{eff}$, then one can prescribe two out of these four: $n_e, n_i, n_1, Z_{eff}$ and the other two are computed from \eqref{eqn:quasin}, \eqref{eqn:zeff}.
In practice, one would usually prescibe in any case $Z_{eff}$ or $n_i/n_1$ plus either $n_e$ or $n_i$. These main options are described in the next section.

\subsubsection{Typical choices}
We list here the typical use cases we want to implement.
In the list below we indicate the main ion species by $n_i$ and any further species as $n_1,n_2$. $Z_1$, $Z_2$ are charges of the non-main ion species, which may be profiles to allow heavy impurities with various charge states. When there are two main ion species (e.g. D-T mixtures), then the ion with the lowest atomic number is considered as $n_i$ and subsequent species are numbered $n_1,n_2$ etc.

\begin{enumerate}
\item $n_e$ predicted, $n_i$ prescribed (or) $n_i=k n_e$, prescribe $Z_{eff}$, hence cheating (inconsistent with \eqref{eqn:zeff}), and inconsistent with \eqref{eqn:quasin} unless $k=1$
\item $n_e$ predicted or prescribed, and $Z_{eff}$ prescribed, and $Z_i$, $Z_1$ known. $\rightarrow$
Calculate $n_i$ and $n_1$ from \eqref{eqn:quasin}, \eqref{eqn:zeff}.
\item $n_1/n_i$ prescribed, $n_e$ predicted or prescribed, and $Z_i$, $Z_1$ known. $\rightarrow$ Calculate $n_i$, $n_1$ from \eqref{eqn:quasin}, \eqref{eqn:zeff}.
 \item $n_1$ and $n_i$ predicted or prescribed, and $Z_i$, $Z_1$ known. $\rightarrow$ Calculate $n_e$, $Z_{eff}$ from \eqref{eqn:quasin}, \eqref{eqn:zeff}. (D-T with $A_i=2.5$)
 \item $n_1$ and $n_i$ predicted or prescribed, and $Z_i$, $Z_1$, $Z_2$ and $Z_{eff}$ known $\rightarrow$ calculate $n_e$, $n_2$ from \eqref{eqn:quasin}, \eqref{eqn:zeff}. (Possibility for D-T)
\item N-dimensional generalizations of the above.
\end{enumerate}


Further details and general possibilities are discussed in Appendix \ref{sec:Particle species simulation options}


%%%%%%%%%%%%%%
\section{Geometric coefficients and quantities related to magnetic flux}
\subsection {Geometric coefficients from equilibrium code quantities}
Equilibrium codes (e.g. LIUQE, CHEASE) return contour integrals over flux surface quantities.

Definition of the $C_i$ coefficients in SI units:
\be
\{C_0,C_1,C_2,C_3,C_4\}= \oint\{\frac{1}{R}, 1, \frac{1}{R^2}, B_p^2,
R^2B_p^2\} \frac{\dd l_p}{B_p}
\ee
Definitions for $\dede{V}{\psi}$ and $B_p$ in terms of COCOS $\cite{Sauter2013}$:
\be
B_p = \frac{1}{(2\pi)^{e_{Bp}}} \cdot \frac{|\nabla \psi|}{R} = \frac{\sigma_{Bp}
  \sigma_{Ip}}{(2 \pi)^{e_{Bp}}} \cdot \frac{|\nabla \psi|}{R}
\ee
\be
\dede{V}{\psi} = \oint \sigma_{Bp} \sigma_{Ip} \cdot \frac{2 \pi R}{|\nabla
  \psi|} \dd l_p  = (2\pi)^{1-e_{Bp}} \sigma_{Bp} \sigma_{Ip} \oint \frac{\dd l_p}{B_p}
\ee
$C_i$ coefficients in terms of COCOS:
\begin{align}
C_{0} &= \oint \frac{1}{R} \frac{\dd l_p}{B_p} = (2 \pi)^{e_{Bp}} \oint \frac{\dd l_p}{|\nabla \psi|} \\
C_{1} &= \oint \frac{\dd l_p}{B_p} = (2 \pi)^{e_{Bp}} \oint \frac{R\dd l_p}{|\nabla \psi|} \\
C_{2} &= \oint \frac{1}{R^2} \frac{\dd l_p}{B_p} = (2 \pi)^{e_{Bp}} \oint\frac{\dd l_p}{R|\nabla \psi|} \\
C_{3} &= \oint B_p^2 \frac{\dd l_p}{B_p} = \oint B_p \dd l_p = \frac{1}{(2 \pi)^{e_{Bp}}} \oint
\frac{|\nabla \psi|}{R} \; \dd l_p \\
C_{4} &= \oint R^2 B_p^2 \frac{\dd l_p}{B_p} = \oint R^2 B_p \dd l_p= \frac{1}{(2 \pi)^{e_{Bp}}}
\oint R |\nabla \psi|  \; \dd l_p
\end{align}
If one wants to normalize with the coefficients $l_d$ for distance and
$l_B$ for magnetic field, one would get
\begin{align}
C_{0}^{norm} &= l_B\oint \frac{1}{R} \frac{\dd l_p}{B_p}= l_B C_0 \label {eqn:C0norm} \\
C_{1}^{norm} &= \frac{l_B}{l_d} \oint \frac{\dd l_p}{B_p}= \frac{l_B}{l_d} C_1 \\
C_{2}^{norm} &=l_dl_B \oint \frac{1}{R^2} \frac{\dd l_p}{B_p}=l_dl_B C_2 \\
C_{3}^{norm} &=\frac{1}{l_dl_B} \oint B_p^2 \frac{\dd l_p}{B_p}= \frac{1}{l_dl_B(2 \pi)^{2e_{Bp}}} \oint
\frac{(\nabla \psi)^2}{R^2}\frac{\dd l_p}{B_p}= \frac{1}{l_dl_B} C_3 \\
C_{4}^{norm} &= \frac{1}{l_d^3 l_B}\oint R^2 B_p^2 \frac{\dd l_p}{B_p}= \frac{1}{l_d^3 l_B(2 \pi)^{2e_{Bp}}}
\oint (\nabla \psi)^2 \frac{\dd l_p}{B_p}= \frac{1}{l_d^3 l_B}
C_4 \label{eqn:C4norm}
\end{align}
The system \eqref{eqn:C0norm}-\eqref{eqn:C4norm} provides the general definition of $C_0^{norm}-C_4^{norm}$ for any
coordinate conventions and any normalization. The $C_i$ coefficients are
obtained with $l_d=l_B=1$.

\section{Flux surface averaged quantities from contour integrals}
\subsection{General definitions in SI and COCOS}
We define $\dede{V}{\psi}, g_1, g_2, g_3, \dede{V}{\hat \rho}, q$ in terms of $C_0-C_4$ defined in
(42)-(46) and in terms of $C_0^{norm}-C_4^{norm}$ defined in (47)-(51).
\be
\dede{V}{\psi}=(2\pi)^{1-e_{Bp}} \sigma_{Bp} \sigma_{Ip} C_1=(2\pi)^{1-e_{Bp}}
\sigma_{Bp} \sigma_{Ip} \frac{l_d}{l_B} C_1^{norm}
\ee
\begin{align}
g_1= \left\langle (\nabla V)^2 \right\rangle= \left(\dede {V}{\psi}\right)^2 \cdot
\left\langle (\nabla \psi)^2 \right\rangle=(2 \pi)^{2(1-e_{Bp})} \oint \frac{\dd l_p}{B_p} \cdot
\oint (\nabla \psi)^2 \frac{\dd l_p}{B_p} \\ = (2 \pi)^{2(1-e_{Bp})} C_1 \cdot (2 \pi)^{2e_{Bp}} C_4= 4 \pi^2 C_1 C_4= 4 \pi^2 l_d^4  \cdot C_1^{norm} C_4^{norm}
\end{align}
\begin{align}
g_2= \left\langle \frac{(\nabla V)^2}{R^2} \right\rangle= \left(\dede {V}{\psi}\right)^2 \cdot
\left\langle \frac{(\nabla \psi)^2}{R^2} \right\rangle= (2 \pi)^{2(1-e_{Bp})} \oint \frac{\dd
  l_p}{B_p} \cdot \oint \frac{(\nabla \psi)^2}{R^2} \frac{\dd l_p}{B_p} \\
  = (2 \pi)^{2(1-e_{Bp})} C_1 \cdot (2 \pi)^{2e_{Bp}} C_3= 4 \pi^2 C_1C_3= 4 \pi^2 l_d^2 \cdot C_1^{norm}C_3^{norm}
\end{align}
\begin{align}
g_3= \left\langle \frac{1}{R^2} \right\rangle={\oint \frac{1}{R^2} \frac{\dd l_p}{B_p}}/{\oint \frac{\dd l_p}{B_p}} = \frac{C_2}{C_1}
=\frac{1}{l_d^2} \frac{C_2^{norm}}{C_1^{norm}}
\end{align}
\begin{align}
q=\frac{\sigma_{Bp} \sigma_{\rho \theta \phi}}{(2 \pi)^{1-e_{Bp}}}
\dede{\Phi}{\psi}=\frac{\sigma_{Bp} \sigma_{\rho \theta \phi}}{(2 \pi)^{1-e_{Bp}}} \dede{\Phi}{V} \dede{V}{\psi}= \sigma_{Ip} \sigma_{\rho \theta \phi} \cdot C_1 \cdot \frac{1}{2 \pi} \dede{}{V} \int
\frac{F}{R^2} \dd V \\ = \sigma_{Ip} \sigma_{\rho \theta \phi} \cdot C_1 \cdot
\frac{F}{2 \pi} \left\langle \frac{1}{R^2} \right\rangle = \sigma_{Ip} \sigma_{\rho \theta \phi}
\frac{1}{2 \pi} F C_2 =  \sigma_{Ip} \sigma_{\rho \theta \phi} \sigma_{B_0} \frac{1}{2 \pi} F^{norm} C_2^{norm}
\end{align}
\begin{align}
\dede{V}{\hat \rho}= \dede{V}{\psi} \dede{\psi}{\Phi} \dede{\Phi}{\hat \rho}=
(2 \pi)^{1-e_{Bp}} \sigma_{Bp} \sigma_{Ip} C_1 \cdot \frac{(2
  \pi)^{e_{Bp}}}{\sigma_{Ip} \sigma_{Bp} F C_2} \cdot 2 \Phi_b \hat \rho= 4 \pi
\frac{C_1}{F C_2} \Phi_b \hat \rho \\= 4 \pi \frac{l_d}{l_B} C_1^{norm}
\frac{1}{F^{norm} C_2^{norm} \sigma_{B_0}} \cdot \Phi_b^{norm} l_d^2 l_B
\sigma_{B_0} \hat \rho= 4 \pi l_d^3 \frac{C_1^{norm}}{F^{norm} C_2^{norm}}
\Phi_b^{norm} \hat \rho
\end{align}
Let us rewrite the find general relations
\\*
\bgroup
\def\arraystretch{1.5}
\begin{tabular} {|c|| c | c|}
\hline
$\dede{V}{\psi}$ & $(2\pi)^{1-e_{Bp}} \sigma_{Bp} \sigma_{Ip} C_1$ & $(2\pi)^{1-e_{Bp}} \sigma_{Bp} \sigma_{Ip} \frac{l_d}{l_B} C_1^{norm}$\\ \hline
$g_1$ & $4 \pi^2 C_1 C_4$ & $4 \pi^2 l_d^4  C_1^{norm} C_4^{norm}$\\ \hline
$g_2$ & $4 \pi^2 C_1C_3$ & $4 \pi^2 l_d^2 C_1^{norm}C_3^{norm}$ \\ \hline
$g_3$ & $C_2/C_1$ & $\frac{1}{l_d^2} C_2^{norm}/C_1^{norm}$ \\  \hline
$q$ & $\sigma_{Ip} \sigma_{\rho \theta \phi} \frac{1}{2 \pi} F C_2$ &
$\sigma_{Ip} \sigma_{\rho \theta \phi} \sigma_{B_0} \frac{1}{2 \pi} F^{norm} C_2^{norm}$ \\  \hline
$\dede{V}{\hat \rho}$ & $4 \pi \frac{C_1}{F C_2} \Phi_b \hat \rho$ & $4 \pi l_d^3 \frac{C_1^{norm}}{F^{norm} C_2^{norm}}
\Phi_b^{norm} \hat \rho$ \\ \hline
\end{tabular}
\egroup

\subsection{Behavior at LCFS for diverted plasmas}
Special care needs to be taken when evaluating these quantities at the last closed flux surface for diverted plasmas. In these cases, $|\nabla \psi|$ goes to 0 along the contour (at the x-point), and the integrals may diverge. In practice, $C_0$,$C_1$,$C_2$ are infinite, while $C_3$ and $C_4$, which do not have $|\nabla \psi|$ in the denominator of the integrand, remain finite. For a general flux surface averaged quantity, and for a single null, it holds that
\be
<Q> = Q_{x-point} \label{eqn:fse_xpoint}
\ee
This can be understood by the fact that the $1/|\nabla\psi|$ term effectively weighs the value at the x-point infinitely more than at other points, so only that value matters for the average.
For the geometric quantities described above we have:
\begin{align}
g_1(\psi_b) &\rightarrow \infty \quad \hbox {since $C_1=\dede{V}{\psi} \rightarrow \infty$}  \\
g_2(\psi_b) &\rightarrow \infty \quad \hbox {since $C_1=\dede{V}{\psi} \rightarrow \infty$} \\
g_3(\psi_b) &= \frac{1}{R_x^2} \\
q(\psi_b) &\rightarrow \infty \quad \hbox {since $C_2 \rightarrow \infty$} \\
\dede{V}{\hat \rho} &= \frac{4 \pi}{F} R_{x}^2 \Phi_b \hat \rho \quad\hbox{since $C_1/C_2 =1/g_3 $}
\end{align}
[later we should probably compute proper limits to write how these go to infinity]


\subsection{Other quantities related to magnetic flux}
q/safety factor = Rotational transform / $2\pi$. Actual rotational trasform is $\iota = 2\pi/q$, strictly speaking $1/q$ is called $\bar\iota$.\\
\be
\iota = \frac{1}{q}  = \dede{\psi}{\Phi} = \dede{\rho}{\Phi}\dede{\hat\rho}{\rho}\dede{\psi}{\hat\rho} =
\frac{1}{2\Phi_b\hat\rho}\dede{\psi}{\hat\rho}
\ee
\begin{align}
I_{pl} = \int j_{\phi} \dd A_{\phi} = \int \oint j_{\phi} \sigma_{B_p} \sigma_{I_p} \frac{\dd \psi}{|\nabla \psi|} \dd l_p = \int \left\langle
\frac{j_{\phi}}{R} \right\rangle \frac{\sigma_{B_p} \sigma_{I_p}}{(2 \pi)^{e_{Bp}}} \oint \frac{\dd l_p}{B_p} \dd \psi = \\
\int \frac{\sigma_{B_p}}{(2 \pi)^{e_{Bp}} \mu_0} \frac{4 \pi^2}{V'_{\rho}} \dede{}{\rho} \left[G_2 \dede{\psi}{\rho} \right] \frac{\sigma_{B_p} \sigma_{I_p}}{(2
  \pi)^{e_{Bp}}} \oint \frac{\dd l_p}{B_p} \dd \psi = \frac{\sigma_{B_p}}{(2
  \pi)^{e_{Bp} -1} \mu_0} G_2 \dede{\psi}{\rho}
\end{align}
One could use eq. for $G_2$, $g_2$, $g_3$ and $\dede{V}{\hat \rho}$ to get
the following equation for $I_p$
\be
 I_{pl} = \frac{\sigma_{B_p}}{(2 \pi)^{e_{Bp}+2}}\frac{1}{2 \mu_0 \Phi_b} F \frac{g_2
   g_3}{\hat \rho} \dede{\psi}{\hat \rho} = \frac{\sigma_{B_p}}{(2 \pi)^{e_{Bp}+2} }\frac{1}{\mu_0} \frac{F}{q} g_2  g_3 \label{eqn:Ipl}
\ee
\be
j_{tor} =  2\pi R_0 \dede{I_{pl}}{V}=\frac{\sigma_{B_p}}{(2 \pi)^{e_{Bp}+2}}
\frac{2 \pi R_0 }{2\mu_0 \Phi_b V'_{\hat \rho}} \dede{}{\hat \rho} \left[F \frac{g_2
    g_3}{\hat \rho} \dede{\psi}{\hat \rho} \right]
\ee
$j_{\parallel}$ (starting from \cite{Pereverzev2002} eq. 34), also \eqref{eqn:jdotB}
\be
j_{\parallel} \equiv \frac{\langle \mathbf{j}\cdot \mathbf{B}\rangle}{B_0} = \frac{2\pi}{B_0} F^2 \dede{}{V} \left( \frac{I_{pl}}{F} \right) =
\frac{2\pi}{B_0} \frac{\sigma_{B_p}}{(2 \pi)^{e_{Bp}+2}} \frac{1}{2\mu_0
  \Phi_b} \frac{F^2}{V'_{\hat \rho}} \dede{}{\hat\rho} \left[ \frac{g_2g_3}{\hat\rho} \dede{\psi}{\hat\rho}\right]
\ee

\subsection{Boundary conditions at LCFS for diverted plasmas}
Also here care must be taken when evaluating products $g_2g_3\dede{\psi}{\rho} $ since $\dede{\psi}{\rho} = 0$ but $g_2g_3 \rightarrow \infty$.
Using \eqref{eqn:Ipl}:
\be
\frac{g_2g_3}{2\Phi_b\rho}\dede{\psi}{\rho} = \frac{g_2g_3}{q} =  \frac{(2 \pi)^{e_{Bp}+2}}{\sigma_{B_p}}  \frac{\mu_0I_{pl}}{F}
\ee
Which is well-defined and finite at the edge. This is indeed what is used as the boundary condition \eqref{eqn:psi_BC_Ipl}.

Also, following \eqref{eqn:fse_xpoint}:
\be
j_{\parallel}(\psi_b) = \frac{\langle \mathbf{j} \cdot \mathbf{B}\rangle_{\psi_b}}{B_0} = \frac{j_xB_x}{B_0}
\ee
which is generally nonzero, but can be zero in special cases zero if $\sigma_{\parallel} = 0$ (e.g. if zero edge temperature is imposed) and $j_{ni}$ current vanishes simultaneously.



\section{RAPTOR quantities in SI units}
Lets define $\dede{V}{\psi}, g_1, g_2, g_3, \dede{V}{\hat \rho}, q$ for RAPTOR
in terms of $C_0-C_4$, F, $\Phi_b$.\\*
RAPTOR: $l_d=1$, $l_B=1$, $e_{Bp}=1$, $\sigma_{B_p}=1$,
$\sigma_{I_p}=1$, $\sigma_{R \phi Z}=1$ $\sigma_{\rho \theta \phi}=1$, $\sigma_{B_0}=1$,
$COCOS=11$. \\*
\\*
\bgroup
\def\arraystretch{1.5}
\begin{tabular} {|c|| c |}
\hline
$\dede{V}{\psi}$ & $C_1$\\ \hline
$g_1$ & $4 \pi^2 C_1 C_4$\\ \hline
$g_2$ & $4 \pi^2  C_1C_3$ \\ \hline
$g_3$ & $C_2/C_1$\\  \hline
$q$ & $\frac{1}{2 \pi} F C_2$\\  \hline
$\dede{V}{\hat \rho}$ & $4 \pi \frac{C_1}{F C_2} \Phi_b \hat \rho$\\ \hline
\end{tabular}
\egroup
\\*
\\*
The same way one could define $I_{pl}$, $j_{tor}$, $j_{\parallel}$ for RAPTOR.\\*
\\*
\bgroup
\def\arraystretch{2}
\begin{tabular} {|c||c|}
\hline
$I_{pl}$ & $  \frac{1}{8 \pi^3\mu_0}\frac{1}{ 2\Phi_b} F \frac{g_2
   g_3}{\hat \rho} \dede{\psi}{\hat \rho}  = \frac{g_2 g_3}{8 \pi^3 \mu_0} \frac{F}{q}$\\ \hline
$j_{tor}$ & $ \frac{2 \pi R_0}{16 \pi^3 \mu_0 \Phi_b V'_{\hat \rho}} \dede{}{\hat \rho}\left[F \frac{g_2g_3}{\hat \rho} \dede{\psi}{\hat \rho} \right]$\\ \hline
$j_{\parallel}$ & $\frac{2\pi}{B_0} \frac{1}{16 \pi^3 \mu_0 \Phi_b}
\frac{F^2}{V'_{\hat \rho}} \dede{}{\hat\rho} \left[ \frac{g_2g_3}{\hat\rho}  \dede{\psi}{\hat\rho}\right]$\\ \hline
\end{tabular}
\egroup

One could calculate the $C_0-C_4$ coefficients, F, $\Phi_b$ in terms of CHEASE outputs for example.\\*
CHEASE: $l_d=R_0$, $l_B=B_0$, $e_{Bp}=0$, $\sigma_{B_p}=1$,
$\sigma_{I_p}=1$, $\sigma_{\rho \theta \phi}=1$, $\sigma_{B_0}=1$, $COCOS=2$.
\begin{align}
C_0= \frac{1}{B_0} C_0^{CH} \\
C_1= \frac{R_0}{B_0} C_1^{CH} \\
C_2=\frac{1}{R_0B_0} C_2^{CH} \\
C_3= R_0B_0 C_3^{CH} \\
C_4= R_0^3 B_0 C_4^{CH}\\
F= R_0B_0 F^{CH}\\
\Phi_b= R_0^2B_0 \Phi_b^{CH}
\end{align}

\section{Plasma pressures}
Many definitions exist for the plasma pressures that make use of the volume integrated thermal pressure defined by $\langle p\rangle_V=\frac{\int p dV}{\int dV}=\frac{\int p dV}{V}$. RAPTOR includes in the output function three formulations:
\begin{itemize}
\item Experimental: $\beta_{exp}=\frac{\langle p\rangle_V}{B_0^2/(2\mu_0)}$, (multiplied by 100 in RAPTOR to get a value close to 1)
\item Normalized: $\beta_N=100\frac{\beta_{exp}}{I_p [\rm{MA}] B_0[\rm{T}] a[\rm{m}]}$
\item Poloidal: $\beta_p=\frac{8}{3}\frac{W_{\rm{th}}}{\mu_0 R_0 I_p^2}$. \\A (brief) derivation:\begin{itemize}
\item Start: $\beta_p=\frac{\langle p\rangle_V}{\bar{B}_p^2/(2\mu_0)}$, with $\bar{B}_p$ an averaged poloidal magnetic field to be defined.
\item We use the definition: $\bar{B}_p=\frac{\int B_p dl_p}{\int dl_p}=\frac{\mu_0 I_p}{\int dl_p}$, where $\int dl_p$ can be either calculated by the equilibrium code or approximated.
\item Gives: $\beta_p=\frac{ \langle p\rangle_V}{\frac{\mu_0 I_p^2}{2(\int dl_p)^2}}$
\item We use large aspect ratio approximation and elongated cylindrical cross section for the volume $V=2\pi^2 R_0 a^2 \kappa$ and for the perimeter we use $(\int dl_p)=\sqrt{\frac{2V}{R_0}}=\int dl_p=2\pi a \sqrt{\kappa}$
\item Hence we can approximate the averaged poloidal magnetic field by: $\bar{B}_p^2=\frac{\mu_0^2 I_p^2}{\frac{2V}{R_0}}$
\item Thus inserting $\bar{B}_p^2$ in $\beta_p$: $\beta_p=\frac{\langle p\rangle_V}{\frac{\mu_0^2 I_p^2}{\frac{2V}{R_0}}/(2\mu_0)}$
    \item Rewriting gives: $\beta_p=\frac{4 \langle p\rangle_V V}{\mu_0 R_0 I_p^2}$
    \item Using the relation: $\langle p\rangle_V =\frac{2}{3}\frac{W_{th}}{V}$
    \item Finally we get: $\beta_p=\frac{8}{3}\frac{W_{\rm{th}}}{\mu_0 R_0 I_p^2}$
\end{itemize}
\end{itemize}

\section{Neoclassical Tearing Modes}
RAPTOR contains an NTM \emph{prescribed} module for which entries in $v_k$ hold the island width $w$ for four modes (default: $q=\left\{\frac{2}{1},\frac{3}{2},\frac{4}{3},\frac{3}{1}\right\}$). Additionally there is the NTM \emph{mre} module which
include states for the NTM width $w$ (default: $q=\left\{\frac{2}{1},\frac{3}{2}\right\}$). In both cases, the width $w$ affects the thermal
conductivity $\chi_e$ in the same way, following:
\be
\chi_{e}(\rho)=\chi_{e,0}(\rho)\prod_{i=1}^{N_{NTMs}}\left(1+A_i\exp\left[\frac{-4\left(\rho-\rho_i\right)^2}{\left(Cw_i*w_i/a\right)^2}\right]\right)
\ee
in which $i$ indicates the mode, $\rho_i$ is the position of the mode, $A_i$ is the amplitude of the perturbation (\emph{chiefact\_A} in \emph{ntm\_mre.m}), $Cw_i$
influences the width of the perturbation relative to the island width (\emph{chiefact\_w} in \emph{ntm\_mre.m}) , $w_i$ is the island width of the $i$-th mode, and $a$ is the minor radius of the plasma.

\subsection{Modified Rutherford Equation}
The modified rutherford equation (mre) is one of the options for including NTMs in the simulation.
Details about the mre terms can be found in \cite{Sauter2015}.
If this is used, a seed island width is read from $v_k$. This value is set in the state and is evaluated following the equation
\begin{align}
 \frac{dw}{dt} &= \frac{\rho_s}{\tau_R} \biggl[ \sum{\rho_s \,\Delta'_{ii} } \biggr],
\end{align}
where $\tau_R$ is the resistive time, $\Delta'_{ii}$ correspond to the various terms, and
\begin{align}
\rho_s=\bar{\rho} &= \hat{\rho} \, a \;[\rm m].
\end{align}

The resistive time is taken as
\begin{align}
  \label{tauR}
 \tau_R &= \frac{\mu_0\,\rho_s^2}{1.22\,\eta_{neo}},
\end{align}
where $\eta_{neo}$ is the neoclassical resistivity at the rational surface and a factor (\emph{resis\_fact}) 
is introduced to reduce or increase local resistivity with respect to the resistivity calculated by the \emph{sigma\_neo} module in RAPTOR.

The classical growth is implemented as equation~(27) in \cite{Sauter2015}, with $\alpha$ and
$\bar{\rho}_{mn}\Delta'_0$ as parameters. The default value for $\bar{\rho}_{mn}\Delta'_0$ is $-m$.
In the square root a small number is added such that in case $\bar{\rho}_{mn}\Delta'_0=-m$ and 
$\alpha=0$ a division by zero does not occur.

The bootstrap current term is implemented as equation~(20) in \cite{Sauter2015}. The value for
$f_{bse}$ is obtained by calling \emph{jBSdotB}, and the values $w_{d,e}$, $w_{d,i}$, and $a_2$ are parameters.
For $j_{bs,mn}$ an average value is used such that an unperturbed bootstrap current is taken.

The Glasser-Greene-Johnson term is implemented as equation~(35), which is multiplied by a
correction factor \emph{a\_ggj}. The factor \emph{a\_ggj} should normally be one, but setting
it to zero will switch off the GGJ term. The value of $L_{bs}$ is obtained from the
value of $L_{31}$ that is calculated in \emph{jBSdotB}.

The current drive term is included following
equation~(39) in \cite{Sauter2015}, using equations~(41) and (42) for the effect of finite island width compared
to current drive width and the effect of misalignment respectively. The current drive density $j_{cd0}$ is calculated
using the \emph{echcd\_gaussian} module of RAPTOR. NTM suppression requires a narrower current drive than what a typical
RAPTOR simulation uses. Therefore, the peak current drive is recalculated for a width $w_{cd,NTM}$ (input parameter)
under the assumption that the total driven current with the larger RAPTOR current drive width is the same.


\appendix
\section{Further notes}

\subsection{Notes on choices and settings for particle transport}
\label{sec:Particle species simulation options}
[JC notes, to be refined]

For a simulation with a single ion species, $Z_{eff}$ is not prescribed, but trivially calculated. Both $Z_{eff}$ and $Z_i$ are profiles. Two $Z_{eff}$ arrays are provided: one for prescribed $Z_{eff}$ (e.g. from experimental data), and the other for calculated $Z_{eff}$ and is written over during the simulation, necessary for when ion species are predicted or separately prescribed (see below). $Z_i$ is a profile to allow for heavy impurities, under the bundling assumption (at each radial location, average over all charge states).
We allow up to 10 ion species to be present (Note (JC): this is to allow, for example, D+T+He3+He4+FastD+FastHe4+FastHe3+W+Kr+Be, which is the most imaginable to ever include). 

For non-predicted (interpretative) ion species, we must constrain the densities by setting their ratios either with respect to $n_e$, or one of the other ion species. We restrict the possible ratios to one of the following:
\begin{equation}
C_i=n_i/n_e
\end{equation}
\begin{equation}
\tilde{C}_i=n_i/n_1
\end{equation}
Where $n_1$ is the first ion species. For example, $C_i$ is useful to constrain minority species used for ICRH heating, which is typically prescribed with respect to $n_e$. On the other hand, $\tilde{C}_i$ may be useful for prescribing $n_D=n_T$, for example. $C_i$ and $\tilde{C}_i$ cannot be set together, and the simulation will not start if this inconsistency is present. If neither $C_i$ or $\tilde{C}_i$ are set, then the $Z_{eff}$ constraint is used for that species. This is covered more in detail below.

\subsubsection{Predicted electron density}
The transport model provides the transport coefficients for electron particle transport. $n_e$ is thus known, and all ion species are interpretative. For simplicity, $n_1$ is always considered unknown (i.e. one cannot set $C_1$). There is then freedom to set one other ion species unknown and constrain it with $Z_{eff}$, or otherwise prescribe the ratios of all other ion species and then calculate $Z_{eff}$ based on that. 

If all ion ratios are prescribed, then the $n_1$ ion density is simply from quasineutrality:
\begin{equation}
n_1=n_e\frac{\left(1-\sum_iZ_iC_i\right)}{Z_1+\sum_j\tilde{C}_jZ_j}
\end{equation}
Where $C_i$ and $\tilde{C_j}$ are the two types of ion density ratio. Of course, one can restrict themselves to purely using either $C_i$ or $\tilde{C_j}$. If $C_i$ is too large, the result may be unphysical (negative $n_1$). A check is set up, and simulation aborts with the appropriate warning if this occurs. $Z_{eff}$ is calculated according to Eq.~\ref{eqn:zeff}.  

When an additional ion species apart from $n_1$ is not prescribed, i.e. no value for $C$ or $\tilde{C}$ is set, then we solve its density according to both quasineutrality and the prescribed $Z_eff$. We name this additional unknown ion species $n_u$, and obtain the following matrix equation.
\begin{equation}
\begin{bmatrix}
	a_{1}       & b_{1} \\
	a_{2}       & b_{2} \\
\end{bmatrix}
\begin{bmatrix}
	n_1 \\
	n_u 
\end{bmatrix}
=
\begin{bmatrix}
c_1 \\
c_2 
\end{bmatrix}
\end{equation}
With:
\begin{equation}
a_1=Z_1+\sum_jZ_j\tilde{C_j}
\end{equation}
\begin{equation}
b_1=Z_u
\end{equation}
\begin{equation}
a_2=Z_1^2+\sum_jZ_j^2\tilde{C_j}
\end{equation}
\begin{equation}
b_2=Z_u^2
\end{equation}
\begin{equation}
c_1=1-\sum_iZ_i\tilde{C_i}
\end{equation}
\begin{equation}
c_2=Z_{eff}-\sum_iZ_i^2\tilde{C_i}
\end{equation}
The standard solution is:
\begin{equation}
n_1=D_x/D
\end{equation}
\begin{equation}
n_u=D_y/D
\end{equation}
With:
\begin{equation}
D=a_1b_2-b_1a_2
\end{equation}
\begin{equation}
D_x=c_1b_2-b_1c_2
\end{equation}
\begin{equation}
D_y=a_1c_2-c_1a_2
\end{equation}
We then set $n_i=C_in_e$ and $n_j=\tilde{C_j}n_1$, for all the other species.

\subsubsection{Predicted ion density}
The transport model provides the transport coefficients for ion particle transport. Some (or all) of the ions are predicted, and the rest are interpretative. For simplicity, we assume that $n_1$ is always predicted in this mode. 

If all ion species are either predicted or prescribed with $C_i$ and $C_j$ ratios, then $n_e$ is the only unknown:
\begin{equation}
n_e=\frac{\sum_kZ_kn_k+\sum_jZ_jn_1\tilde{C_j}}{1-\sum_iZ_iC_i}
\end{equation}
Where $n_k$ are all the known (predicted) ion species, including $n_1$. After calculation of $n_e$, we then set $n_i=C_in_e$ and $n_j=\tilde{C_j}n_1$, for all the other ion species, as well as $Z_{eff}$.

When we include an unknown ion species, constrained by the prescribed $Z_{eff}$, we then obtain a similar matrix as before:
\begin{equation}
\begin{bmatrix}
a_{1}       & b_{1} \\
a_{2}       & b_{2} \\
\end{bmatrix}
\begin{bmatrix}
n_e \\
n_u 
\end{bmatrix}
=
\begin{bmatrix}
c_1 \\
c_2 
\end{bmatrix}
\end{equation}
Now with:
\begin{equation}
a_1=-\left(1-\sum_iZ_i\tilde{C_i}\right)
\end{equation}
\begin{equation}
a_2=-\left(Z_{eff}-\sum_iZ_i^2\tilde{C_i}\right)
\end{equation}
\begin{equation}
b_1=Z_u
\end{equation}
\begin{equation}
b_2=Z_u^2
\end{equation}
\begin{equation}
c_1=-\left(Z_1+\sum_jZ_j\tilde{C_j}\right)
\end{equation}
\begin{equation}
c_2=-\left(Z_1^2+\sum_jZ_j^2\tilde{C_j}\right)
\end{equation}
And again the standard solution is:
\begin{equation}
n_1=D_x/D
\end{equation}
\begin{equation}
n_u=D_y/D
\end{equation}
With:
\begin{equation}
D=a_1b_2-b_1a_2
\end{equation}
\begin{equation}
D_x=c_1b_2-b_1c_2
\end{equation}
\begin{equation}
D_y=a_1c_2-c_1a_2
\end{equation}
We then set $n_i=C_in_e$ and $n_j=\tilde{C_j}n_1$, for all the interpretative ion species.



\begin{thebibliography}{3}
\bibitem{Hinton1976}
  F. L. Hinton and R. D. Hazeltine, ÒTheory of plasma transport in toroidal confinement systems,Ó Rev. Mod. Phys., vol. 48, no. 2, pp. 239--308, Apr. 1976.
\bibitem{Fable2013}
  E. Fable, et al, Nucl. Fusion $\textbf{53}$, 033002 (2013).
\bibitem{Pereverzev2002}
  G.V. Pereverzev, P.N. Yusmanov,
  ASTRA Automated System for TRansport Analysis in a Tokamak. Technical Report
  5/98, IPP Report, February 2002.
\bibitem{Sauter2015}
  O. Sauter, et al., Details on the expressions used in the various
terms of the modified Rutherford equation (MRE), Internal Report July 2015
\bibitem{Sauter2012}
  O. Sauter and S. Medvedev, Computer Physics Communications 184 (2013) 293--302
\end{thebibliography}

  \end{document}
