#  Spectral broadening from turbulence in multiscale lower hybrid current
drive simulations

Bodhi Biswas Plasma Science and Fusion Center, Massachusetts Institute of
Technology  Present affiliation: University of York, UK. E-mail:
<EMAIL_ADDRESS>Plasma Science and Fusion Center, Massachusetts
Institute of Technology  Paul Bonoli Plasma Science and Fusion Center,
Massachusetts Institute of Technology  Abhay Ram Plasma Science and Fusion
Center, Massachusetts Institute of Technology  Anne White Plasma Science and
Fusion Center, Massachusetts Institute of Technology

###### Abstract

_The scattering of lower hybrid (LH) waves due to scrape-off layer (SOL)
filaments is investigated. It is revealed that scattering can account for the
LH spectral gap without any ad hoc modification to the wave-spectrum. This is
shown using a multiscale simulation approach which allows, for the first time,
the inclusion of full-wave scattering physics in ray-tracing/Fokker-Planck
calculations. In this approach, full-wave scattering probabilities are
calculated for a wave interacting with a statistical ensemble of filaments.
These probabilities are coupled to ray-tracing equations using radiative
transfer (RT) theory. This allows the modeling of scattering along the entire
ray-trajectory, which can be important in the multi-pass regime. Simulations
are conducted for lower hybrid current drive (LHCD) in Alcator C-Mod,
resulting in excellent agreement with experimental current and hard X-ray
(HXR) profiles. A region in filament parameter space is identified in which
the impact of scattering on LHCD is saturated. Such a state coincides with
experimental LHCD measurements, suggesting saturation indeed occurs in C-Mod,
and therefore the exact statistical properties of the filaments are not
important._

##  1 Introduction

Lower hybrid waves are an efficient means to non-inductively drive current in
a tokamak via electron Landau damping (ELD) [1]. It is an attractive actuator
for current profile shaping, which has been successfully demonstrated in
several tokamaks[2, 3, 4, 5, 6]. There is also interest in targeted LHCD for
neoclassical tearing mode suppression [7, 8].

The condition for strong linear ELD is [9]

$\absolutevalue{N_{||}}\geq N_{||,\text{ELD}}\equiv\frac{c}{3v_{te}}$ (1)

where $N_{||}\equiv\mathbb{N}\cdot\mathbb{B}/B$ is the parallel refractive
index (with respect to the background magnetic field $\mathbb{B}$), $c$ is the
speed of light, and $v_{te}=\sqrt{2T_{e}/m_{e}}$ is electron thermal velocity
(with $T_{e}$ and $m_{e}$ being the electron temperature and mass,
respectively). If the initial parallel refractive index
$\absolutevalue{N_{||0}}<N_{\text{ELD}}$, the wave is expected to first
undergo $|N_{||}|$-upshift due to toroidicity until the condition in Eq. (1)
is met [10]. For a sufficiently large _spectral-gap_ between
$\absolutevalue{N_{||0}}$ and $N_{\text{ELD}}$, the wave may require multiple
passes through the plasma before this gap is bridged. This scenario is typical
of present-day tokamaks due to low core $T_{e}$.

A large spectral-gap poses difficulties for theory and modeling of LHCD. Ray-
tracing/Fokker-Planck modeling of LHCD reveal large discrepancies between
simulations and experiments in several tokamaks [11, 12, 13, 14]. In high
aspect ratio devices ($\mathcal{A}\equiv R_{0}/a\gtrsim 5$, where $R_{0}$ is
the major radius, and $a$ is the minor radius), like WEST and TRIAM-1,
toroidicity is insufficient to bridge the spectral-gap via
$\absolutevalue{N_{||}}$-upshift. In moderate aspect ratio devices
($\mathcal{A}\approx 3$), like Alcator C-Mod and EAST, toroidicity can
accurately model the _total_ current drive efficiency, but not the radial
current or HXR profiles [11, 13]. For example, in C-Mod, the simulated LH
current profiles are not smooth and are peaked far off-axis, while
experimental motional Stark effect (MSE) and HXR measurements suggest smooth
profiles robustly peaked on-axis [11]. Consequently, $N_{||}$-upshift due to
toroidicity is not adequate for understanding experimental observations. There
exist other contributing factors that must be taken into account in
theoretical models.

Several mechanisms have been proposed as possible explanations for the
spectral gap. Ad hoc modifications to the launched LH wave-spectrum reveal
that either $N_{||}$ broadening [15] or angle-broadening of the perpendicular
refractive index ($\mathbb{N}_{\perp}$) [16] can explain experimental current
drive results in C-Mod.

The two theorized causes of $N_{||}$ broadening are parametric decay
instabilities (PDI) [17, 18] and scattering from parallel density fluctuations
in front of the LH antenna [19]. In FTU, a combination of modeling and
experiment shows PDI is likely responsible for low LHCD efficiency at high
densities [20]. However, there is little theoretical support for - or
experimental evidence of - strong PDI in low and moderate density C-Mod
discharges ($\bar{n}_{e}<10^{20}\,\text{m}^{-3}$) where the spectral gap
persists [21]. Likewise, there is little evidence of density fluctuations with
the large parallel gradients required to induce a significant
$N_{||}$-broadening.

The most likely mechanism for $\mathbb{N}_{\perp}$ angle-broadening is
scattering from turbulent scrape-off layer (SOL) fluctuations. Prior models
have employed either ray-tracing [22, 23] or wave-kinetic treatments [10, 24,
25] for scattering in drift-wave turbulence. These models have demonstrated
modified current profiles, but cannot match experimental observations. Gas-
puff imaging (GPI) [26] and statistical analysis of Langmuir probe
measurements [27] in the SOL have motivated the modeling of LH scattering from
intermittent, field-aligned filaments. The extent of $\mathbb{N}_{\perp}$
angle-broadening from filaments is greater than in “equivalent” drift-wave
turbulence [28]. A recent hybrid full-wave/statistical model for wave-filament
interactions was developed to model the modification to the LH wave-spectrum
in front of the antenna [29]. Multiple wave-filament interactions are
accounted for using the radiative transfer (RT) equation. This treatment
allows the modeling of realistic turbulence parameters without being
restricted to the validity constraints of the ray-tracing or random phase
approximation. The study finds a large angle-broadening of the incident wave-
spectrum, enough to robustly direct a fraction of LH power to damp on-axis on
first pass through the plasma. In turn, the LH current profile is monotonic
and peaked on-axis, in much closer agreement with experiment. In addition,
asymmetric scattering of the LH wave is observed in angle-space. This is a
full-wave effect only possible with intermittent density fluctuations, and
therefore is not accounted for in prior ray-tracing or wave-kinetic
treatments.

The hybrid full-wave/statistical scattering model discussed above is limited
to a slab geometry with homogeneous background and turbulence parameters. It
also only treats slow wave to slow wave ($S\rightarrow S$) scattering and
ignores the fast ($F$) wave. Therefore this model can only approximately treat
scattering directly in front of the LH antenna. In this paper, the more
general RT equation is solved using a multiscale full-wave/ray-tracing solver.
It allows the modeling of arbitrary geometry and both like-mode ($S\rightarrow
S$, $F\rightarrow F$) and unlike-mode ($S\rightarrow F$, $F\rightarrow S$)
scatter. An arbitrary geometry allows accounting for realistically tapered SOL
turbulence profiles in a tokamak, and models scattering along the entire ray-
trajectory (important in the multi-pass regime). The inclusion of all
scattering modes is especially important near the mode-conversion density
where like- and unlike-mode scattering probabilities are comparable.

This multiscale model is applied to Alcator C-Mod, allowing, for the first
time, the inclusion of self-consistent full-wave scattering physics in LHCD
simulations. The resulting current and HXR profiles provide excellent matches
to experimental measurements. Thus, the spectral-gap in C-Mod is resolved via
turbulent scattering, and without any ad hoc modification to the wave-
spectrum.

This paper is organized as follows. Section 2 will describe the multiscale
scattering model. Section 3 discusses its application to LHCD in Alcator C-Mod
discharges. Section 4 provides a discussion and summary of the findings.

##  2 Multiscale scattering model

From a modeling perspective, the LH wave is in an interesting and challenging
range of wavelengths. The wavelength is small enough to employ ray-tracing in
the quiescent core [30], but large enough that full-wave modeling is required
for common SOL turbulence parameters. On the other hand, the LH wavelength is
small enough such that whole-device full-wave modeling in a turbulent plasma
is prohibitively expensive. Therefore, a multiscale model that employs ray-
tracing in the core and a full-wave solver in the SOL is a promising concept.
The implementation is as follows. A Mie-scattering approach efficiently
calculates the scattering probabilities for an incident plane wave interacting
with a filament. The probabilities are averaged over a statistical ensemble of
filament parameters. Next, the scattering of RF power in phase-space is
modeled using the RT approximation. The RT equation is solved in a ray-tracing
solver, where the scattering terms act as stochastic kicks to the ray-
trajectory. Finally, the rays are used to calculate current drive using the
quasi-linear Fokker-Planck equation.

###  2.1 Review of single wave-filament interaction

The Mie-scattering model for a single wave-filament interaction is briefly
reviewed. An incident LH plane wave, either the slow ($S$) or fast ($F$) mode,
travels through a cold magnetized plasma with a homogeneous background density
$n_{0}$ and magnetic field $\mathbb{B}$. The magnetic field is aligned along
$\mathbb{\hat{e}}_{z}$, such that the parallel wave-number $k_{||}=k_{z}$. The
wave trajectory is aligned such that
$\mathbb{v}_{gr\perp}=v_{gr\perp}\mathbb{\hat{e}}_{x}$, where
$\mathbb{v}_{gr\perp}$ is the perpendicular group velocity and
$v_{gr\perp}>0$. An infinitely long, cylindrical, field-aligned filament
passes through the origin with density $n_{b}$. Given this is a poloidally
symmetric system (where $\theta=\text{tan}^{-1}(y/x)$), the Jacobi-Anger
expansion is used to write the electric field in cylindrical coordinates
($\rho,\theta,z$). Consequently, the field everywhere is a series solution in
poloidal mode-numbers ($m\equiv k_{\theta}\rho$),

$\displaystyle E_{j\beta}$ $\displaystyle=e^{i(k_{||}z-\omega
t)}\sum_{m=-\infty}^{+\infty}E_{jm}W_{j\beta
m}e^{im\theta};\quad\quad\beta=\rho,\theta,z$ (2a) $\displaystyle W_{j\rho m}$
$\displaystyle=\xi_{jx}J_{m}^{{}^{\prime}}(k_{j\perp}\rho)-i\xi_{jy}\frac{m}{k_{j\perp}\rho}J_{m}(k_{j\perp}\rho)$
(2b) $\displaystyle W_{j\theta m}$
$\displaystyle=i\xi_{jx}\frac{m}{k_{j\perp}\rho}J_{m}(k_{j\perp}\rho)+\xi_{jy}J_{m}^{{}^{\prime}}(k_{j\perp}\rho)$
(2c) $\displaystyle W_{jzm}$ $\displaystyle=i\xi_{jz}J_{m}(k_{j\perp}\rho)$
(2d)

where $j=0,...,4$ is the wave index, indicating (0) the incident wave, (1,2)
the slow/fast mode inside the filament, and (3,4) the slow/fast mode outside
the filament. $\bar{\xi}_{j}=\\{\xi_{jx},\xi_{jy},\xi_{jz}\\}$ is the plane-
wave polarization of wave $j$. $J_{m}$ is the Bessel function of the first
kind and order $m$. $J_{m}^{{}^{\prime}}$ is the first derivative of $J_{m}$
with respect to its argument. The coefficients $E_{jm}$ for the incident plane
wave ($j=0$) are known through the Jacobi-Anger expansion. The coefficients
for the remaining waves are found by imposing Maxwell’s boundary conditions at
the filament perimeter. This is detailed in Ram _et al._ , (2016) [31].

This technique can be extended to a radially varying filament density profile,
as detailed in Biswas _et al._ , (2021) [29]. A similar method was used to
model ion-cyclotron wave scattering [32]. For the purposes of this work, it is
assumed the filament has a Gaussian density profile such that

$n(\rho)-n_{0}=n_{0}\left(\frac{n_{b}}{n_{0}}-1\right)e^{-\left(\frac{2\sqrt{\text{ln}(2)}\rho}{a_{b}}\right)^{2}}$
(3)

where $n$ is the density, $n_{b}/n_{0}$ is the relative density of the
filament, and $a_{b}$ is the full-width half max of the filament radial
profile.

Following [33, 29], the differential scattering-width is calculated by
evaluating the radially scattered Poynting flux in the far-field limit

$\sigma_{j\rightarrow j^{\prime}}(\theta)=\frac{\mp
2}{\pi}\frac{|\xi_{j^{\prime}y}|^{2}+|\xi_{j^{\prime}z}|^{2}+\frac{k_{z}}{k_{j^{\prime}\perp}}\Re{\xi_{j^{\prime}x}\xi_{j^{\prime}z}^{*}}}{k_{jx}(|\xi_{jy}|^{2}+|\xi_{jz}|^{2})-k_{z}\Re{\xi_{jx}\xi_{jz}^{*}}}\absolutevalue{\sum_{m=-\infty}^{+\infty}i^{\pm
m}E_{j^{\prime}m}e^{im\theta}}^{2}$ (4)

where now $j=S,F$ is the incident mode, $j^{\prime}$ is the scattered mode,
$\xi_{j\ell}$ is the $\ell$-component of the electric field polarization of
mode $j$, and $k_{j\ell}$ is the $\ell$-component of the wave-number of mode
$j$. The sign $\mp$ depends on the scattered mode. The differential
scattering-width is analogous to a differential scattering cross-section.
Likewise, $\sigma=\int\sigma(\theta)d\theta$ is the scattering-width, which is
analogous to a scattering cross-section. Lastly,
$\hat{\sigma}(\theta)\equiv\sigma(\theta)/\sigma$ is the normalized
differential scattering-width and will be useful in writing the RT equation.

###  2.2 Scattering through statistical ensemble of filaments

To account for the statistical variation in filament relative density
($n_{b}/n_{0}$) and radial width ($a_{b}$), a joint probability distribution
function $p(n_{b}/n_{0},a_{b})$ is defined for an ensemble of filaments. The
quantities $\langle n_{b}/n_{0}\rangle$ and $\langle a_{b}\rangle$ denote the
average relative density and radial width of filaments. The average, or
“effective”, differential scattering-width is

$\sigma_{\text{eff},j\rightarrow
j^{\prime}}(\theta)=\langle\sigma_{j\rightarrow
j^{\prime}}(\theta)\rangle=\int_{0}^{\infty}\text{d}a_{b}\int_{0}^{\infty}\text{d}\left(n_{b}/n_{0}\right)\,\sigma_{j\rightarrow
j^{\prime}}(\theta;n_{b}/n_{0},a_{b})p(n_{b}/n_{0},a_{b})$ (5)

A wave packet traveling in a straight line will encounter, on average,
$\frac{f_{p}}{\pi\langle a_{b}\rangle^{2}}$ filaments per unit length in the
perpendicular plane. Thus, the inverse mean-free-path to scatter is

$\Sigma_{\text{eff},j\rightarrow j^{\prime}}=\frac{f_{p}}{\pi\langle
a_{b}\rangle^{2}}\sigma_{\text{eff},j\rightarrow j^{\prime}}$ (6)

where $f_{p}$ is the packing-fraction, which is defined as the fractional area
(in the perpendicular plane) inhabited by filaments. Note that
$\Sigma_{\text{eff},j\rightarrow j^{\prime}}$ and
$\sigma_{\text{eff},j\rightarrow j^{\prime}}(\theta)$ are dependent on the
local plasma parameters, the incident wave’s frequency and wavenumber, and the
filament PDF. These functional dependencies are only stated explicitly when
needed.

###  2.3 Radiative transfer approximation

Radiative transfer (RT) theory models the RF wave-spectrum as a wave-packet
distribution function in phase-space, akin to the Fokker-Planck equation for a
particle distribution function [24]. We define $P_{j}(\mathbb{r},\mathbb{k})$
as the distribution corresponding to the LH mode $j$. The RT equation
governing $P_{j}$ is

$\left(\frac{dP_{j}}{dt}\right)_{\text{r}}+2\gamma(\mathbb{k}_{j},\mathbb{r})P_{j}=\left(\frac{dP_{j}}{dt}\right)_{\text{sct}}$
(7)

The $(...)_{\text{r}}$ term is the convective term accounting for the
trajectory for a wave-packet/ray. The second term accounts for wave damping,
where $\gamma$ is the damping rate. The left hand side (LHS) is routinely
solved using ray-tracing/Fokker-Planck codes. The right hand term accounts for
the added effect of scattering,

$\begin{split}\left(\frac{dP_{j}}{dt}\right)_{\text{sct}}&={\sum}_{j^{\prime}=S,F}-\Sigma_{\text{eff},j\rightarrow
j^{\prime}}(k_{||},\mathbb{r})|\mathbb{v}_{gr\perp}(\mathbb{k}_{j},\mathbb{r})|P_{j}\\\
&+{\sum}_{j^{\prime}=S,F}\Sigma_{\text{eff},j^{\prime}\rightarrow
j}(k^{\prime}_{||},\mathbb{r})|\mathbb{v}_{gr\perp}(\mathbb{k}_{j^{\prime}},\mathbb{r})|\int_{-\pi}^{\pi}\hat{\sigma}_{\text{eff},j^{\prime}\rightarrow
j}(\chi-\chi^{\prime};k_{||},\mathbb{r})P_{j^{\prime}}(\chi^{\prime},k_{||},\mathbb{r})d\chi^{\prime}\end{split}$
(8)

The transformation of $\theta\rightarrow\chi$ must first be explained. The
variable $\theta$ describes the angle of scatter in the frame of a single
filament-wave interaction, as described in Section 2.1. In contrast, $\chi$ is
defined through
$\mathbb{b}\cdot(\hat{\mathbb{e}}_{\psi}\times\mathbb{k}_{\perp})=k_{\perp}\sin{\chi}$,
where $\hat{\mathbb{e}}_{\psi}$ is the unit vector normal to the flux surface
and directed outwards, and $\mathbb{b}\equiv\mathbb{B}/B$ is the unit vector
along the magnetic field. Thus, by using $\chi$, the angular orientation of
$\mathbb{k}_{\perp}$ is unambiguously defined in a tokamak geometry. The
functional dependence of $P_{j}(\mathbb{r},\mathbb{k})$ can be mapped to
$P_{j}(\chi,k_{||},\mathbb{r})$ and vice-versa using the dispersion relation
for the appropriate mode, and the local magnetic geometry.

The two RHS terms in Eq. (8) account for angular ($\chi$) rotation of the
perpendicular wave-vector $\mathbb{k}_{\perp}$ due to a wave-filament
interaction. The first term accounts for _out-scatter_ of $P_{j}$ from a
phase-space volume element centered at $\chi$. The second term accounts for
_in-scatter_ into this volume element from all other angles $\chi^{\prime}$.
Note that the scattering probabilities and background plasma parameters are
allowed to vary in real-space. In addition, the summation over the slow and
fast modes now ensure that both like-mode and unlike-mode scatter
($S\rightarrow F,F\rightarrow S$) are accounted for. In contrast, the system
studied in [29] was restricted to homogeneous scattering and background
parameters, and neglected mode-conversion.

It should be clarified that, in generating scattering probabilities using a
full-wave formalism, interference effects during a single wave-filament
interaction are accounted for. However, in using the RT approximation, the
interference effects of simultaneous _multi_ -filament scattering are ignored.
This is a reasonable approximation as long as $k_{\perp}d\gg 1$, where $d$ is
the average distance between filaments [34]. Code comparison with a fully
numeric full-wave solver finds that this multiscale method generally
overestimates the effects of scattering for $f_{p}\gtrsim 0.15$ [29].
Nevertheless, the multiscale model retains many important full-wave effects
and is therefore a significant improvement over prior reduced RF-turbulence
scattering models in the SOL.

###  2.4 Coupling to ray-tracing code GENRAY

The initial wave-spectrum, launched from the LH antenna, is discretized into
rays. Each ray-trajectory is evolved using the ray-tracing equations in GENRAY
[35]. The quasi-linear calculation of $\gamma$, accounting for Landau damping
and collisions, is accomplished with the Fokker-Planck solver CQL3D [36]. In
this way, the usual ray-tracing/Fokker-Planck calculation solves for the LHS
of Eq. (7).

The right hand side (RHS) of Eq. (7) introduces stochastic kicks to the ray-
trajectories via a Monte Carlo process. Previous ray-tracing Monte Carlo
scattering models exist for drift-wave turbulence [10, 25]. In GENRAY, after
one time-step of $\Delta t$, the ray-trajectory is evolved from $\mathbb{r}$
to $\mathbb{r}+\Delta\mathbb{r}$, where $\mathbb{r}$ is a spatial coordinate
and $\Delta\mathbb{r}=\mathbb{v}_{gr}\Delta t$. Given
$\Sigma_{\text{eff},j\rightarrow j}$ is the inverse mean free path for like-
mode scatter, the probability of this scattering event during this time-step
is $p_{\text{sct}}=1-e^{-\Sigma_{\text{eff},j\rightarrow
j}|\Delta\mathbb{r}_{\perp}|}\approx\Sigma_{\text{eff},j\rightarrow
j}|\Delta\mathbb{r}_{\perp}|$, where $|\Delta\mathbb{r}_{\perp}|$ is the
distance travelled perpendicular to $\mathbb{b}$. The approximation holds for
$\Sigma_{\text{eff},j\rightarrow j}|\Delta\mathbb{r}_{\perp}|\ll 1$. The
probability of unlike-mode scatter, and no scatter, are calculated in a
similar manner. After each time-step, a uniformly distributed random variable
$X_{1}\in[0,1]$ is sampled. Assuming the incident ray is a slow wave, there
are three possibilities,

  1. 1.

$X_{1}<\Sigma_{\text{eff},S\rightarrow S}|\Delta\mathbb{r}_{\perp}|$. Then a
$S\rightarrow S$ scatter occurs.

  2. 2.

$\Sigma_{\text{eff},S\rightarrow S}|\Delta\mathbb{r}_{\perp}|\leq
X_{1}<(\Sigma_{\text{eff},S\rightarrow S}+\Sigma_{\text{eff},S\rightarrow
F})|\Delta\mathbb{r}_{\perp}|$. Then a $S\rightarrow F$ scatter occurs.

  3. 3.

$X_{1}\geq(\Sigma_{\text{eff},S\rightarrow S}+\Sigma_{\text{eff},S\rightarrow
F})|\Delta\mathbb{r}_{\perp}|$. Then no scatter occurs.

The above stochastic scattering procedure is only valid when the scattering
probability between time-steps is small ($(\Sigma_{\text{eff},S\rightarrow
S}+\Sigma_{\text{eff},S\rightarrow F})|\Delta\mathbb{r}_{\perp}|\ll 1$). Once
a scattering event is initiated, another uniformly distributed random variable
$X_{2}\in[0,1]$ is sampled. A new variable is introduced:

$\tilde{\sigma}_{\text{eff},j\rightarrow
j^{\prime}}(\chi)=\frac{1}{\sigma_{\text{eff},j\rightarrow
j^{\prime}}}\int_{-\pi}^{\chi}\sigma_{\text{eff},j\rightarrow
j^{\prime}}(\chi^{\prime})d\chi^{\prime}$ (9)

where $\tilde{\sigma}_{\text{eff},j\rightarrow j^{\prime}}(\Delta\chi)$ is the
cumulative distribution function (CDF) of rotating the incident ray
$\mathbb{k}_{\perp}$ by angle $\Delta\chi$. Inversely,

$\Delta\chi=\tilde{\sigma}_{\text{eff},j\rightarrow j^{\prime}}^{-1}(X_{2})$
(10)

Therefore, the incident ray is converted to a scattered ray, by first
undergoing mode-conversion in the case of unlike-mode scatter, and then
rotating $\mathbb{v}_{gr\perp}$ by $\Delta\chi$ with respect to $\mathbb{b}$.
Note that, in the case of like-mode scatter, this is equivalent to rotating
$\mathbb{k}_{\perp}$ by $\Delta\chi$. In the case of unlike-mode scatter, the
direction of $\mathbb{k}_{\perp}$ is rotated by $\Delta\chi+\pi$ because the
fast wave is forward-propagating
($\mathbb{k}_{\perp}\cdot\mathbb{v}_{gr\perp}=|k_{\perp}v_{gr\perp}|$) while
the slow wave is backward-propagating
($\mathbb{k}_{\perp}\cdot\mathbb{v}_{gr\perp}=-|k_{\perp}v_{gr\perp}|$).

The values of $\Sigma_{\text{eff},j\rightarrow j^{\prime}}$ and
$\tilde{\sigma}_{\text{eff},j\rightarrow j^{\prime}}(\chi)$ depend on $k_{||}$
and $\mathbb{r}$. They are implemented in GENRAY as lookup tables
parameterized by $k_{||}$, $\psi$, and $\theta_{p}$ (where $\psi$ is the
normalized radial coordinate and $\theta_{p}$ is the poloidal coordinate
typically used in tokamaks). The filament statistics (and therefore the
scattering probabilities) are made to vary such that relative density
fluctuations increase with outward distance from the separatrix. Due to the
ballooning nature of SOL turbulence, fluctuations are smaller on the high
field side. This is accounted for via the $\theta_{p}$ dependence of the
scattering probabilities. For simplicity, $f_{p}$ is assumed constant (though
it can also be made to vary spatially). Details on the radial and poloidal
tapering of scattering probabilities can be found in Appendix A.

##  3 Application to Alactor C-Mod

A well-studied [11, 28, 16, 29] low-$\bar{n}_{e}$, fully non-inductive
discharge (#1101104011) is modeled in GENRAY/CQL3D. Both HXR and MSE
measurements are available, allowing for direct comparison to the simulated
current density HXR count profiles. In modelling a quasi-steady-state time-
slice at which $V_{\text{loop}}\approx 0$, it can be safely assumed that Ohmic
contribution to the current profile is negligible.

Smooth electron density and temperature profiles are fit to Thompson
scattering measurements. A kinetic EFIT profile is used (constrained by MSE
measurements and produced in [11]). To model the SOL geometry, the two-point
model is used to generate $n,T$ profiles outside the separatrix [37]. It is
assumed that the SOL e-folding width $\lambda_{\text{SOL}}=7\,$mm, and that
$Z_{\text{eff}}=1.5$ everywhere.

###  3.1 Scattering probabilities

The filament joint-PDF $p(n_{b}/n_{0},a_{b})$ is assumed to be a bivariate
normal distribution. Relevant parameters include the average
($\langle...\rangle$), standard deviation ($s$), and skewness ($\Gamma$) for
$n_{b}/n_{0}$ and $a_{b}$. (See Appendix A for a complete definition of these
parameters, and how they vary radially and poloidally. Figures 12 and 13
summarize the spatial variations of the filament PDF).

For the purposes of illustrating how scattering probabilities vary in the SOL,
it is assumed that $\langle n_{0}/n_{b}\rangle_{\text{grill}}=3$, $\langle
a_{b}\rangle=0.5\,$cm, $s_{a_{b}}=0.2\,$cm, and $\Gamma_{a_{b}}=7$. (The
parameter $\langle n_{0}/n_{b}\rangle_{\text{grill}}$ sets $\langle
n_{0}/n_{b}\rangle$ at the LH antenna.) Figures 1 plot
$\Sigma_{\text{eff},i\rightarrow j}(\psi,\theta_{p})$ for $N_{||}=1.6$ and
$f_{p}=0.2$ (though one can argue that so high a $f_{p}$ is not strictly valid
given the discussion in Section 2.3). As expected, scattering is strongest in
the far-SOL at the outer mid-plane, where $\Sigma_{\text{eff},S\rightarrow
F}\approx 40\,\text{m}^{-1}$. Conversely, scattering is weakest further into
the core and closer to the inner mid-plane. Striations and local maxima in
these plots, most notably for $\Sigma_{\text{eff},S\rightarrow S}$ and
$\Sigma_{\text{eff},F\rightarrow F}$, are evidence of resonances akin to those
found in Mie-scattering.

Figure 1: $\Sigma_{\text{eff},i\rightarrow j}(\psi,\theta_{p},N_{||}=1.6)$
plotted for $\langle n_{0}/n_{b}\rangle_{\text{grill}}=3$, $\langle
a_{b}\rangle=0.5\,$cm, $s_{a_{b}}=0.2\,$cm, $\Gamma_{a_{b}}=7$, and
$f_{p}$=0.2. The white regions correspond to an evanescent FW. The red line
indicates background density.

It should be noted that, while $S\rightarrow F$ and $F\rightarrow S$ are
“opposite” scattering events, their probabilities are not equal. Within the
range of parameters tested, it is most common that
$\Sigma_{\text{eff},F\rightarrow S}>\Sigma_{\text{eff},S\rightarrow F}$. In
other words, there is preferential scattering from $S\rightarrow F$.

Asymmetric scatter (in angle space) is a full-wave effect, and is found to
make a significant modification to the LH wave-spectrum. The physical
explanation for asymmetric scatter is discussed in [29]. Figures 2 plots the
effective asymmetric scatter metric $\alpha_{\text{eff},i\rightarrow j}$,
defined as:

$\alpha_{\text{eff},j\rightarrow
j^{\prime}}\equiv\int_{0}^{\pi}\langle\,\hat{\sigma}_{j\rightarrow
j^{\prime}}(\chi)\rangle\text{d}\chi-\frac{1}{2}$ (11)

An $\alpha_{\text{eff},j\rightarrow j^{\prime}}>0$ ($<0$) corresponds to
scatter from mode $j$ to $j^{\prime}$ preferentially in the $+\chi$ ($-\chi$)
direction. Since LH waves are generally launched with $N_{||}<0$ in order to
drive co-current, it can be shown that preferential $+\chi$ rotation equates
to preferential scatter _away_ from the plasma core. Conversely, $-\chi$
rotation leads to scattering _towards_ the core. Figures 1 and 2 show the
strongest asymmetry is from $S\rightarrow S$ scatter (in the sense that
$\Sigma_{\text{eff},S\rightarrow S}$ is rather large and
$|\alpha_{\text{eff},S\rightarrow F}|$ is robustly positive). Thus, there is a
net scatter of power away from the core everywhere in the SOL. Although not
shown here, this trend is insensitive to $N_{||}$ in the ranges expected in
C-Mod.

Figure 2: Asymmetric scattering metric $\alpha_{\text{eff},i\rightarrow
j}(\psi,\theta_{p},N_{||}=1.6)$ plotted for same filament joint-PDF’s as in
Figure 1. White regions correspond to $\Sigma_{\text{eff},i\rightarrow
j}<1\,\text{m}^{-1}$, where scattering is negligible. The green line indicates
background density.

###  3.2 Parametric scan of turbulence parameters

As a metric for simulation accuracy compared to experiment, we define the
_normalized mean squared error_ for an arbitrary simulated measurement $X$ as

$\bar{S}^{2}=\frac{||X_{\text{sim}}-X_{\text{exp}}||}{||X_{\text{simref}}-X_{\text{exp}}||}$
(12)

where “sim” denotes the simulated values, “exp” the experimental values, and
“simref” the simulated values in the reference (no turbulence) case. Further,
we define $\bar{S}^{2}_{J_{\phi}}$ as replacing the array $X$ in Eq. (12) with
toroidal current density ($J_{\phi}$) evaluated at the radial grid points
modeled in CQL3D.

A parametric scan of parameters $\langle n_{b}/n_{0}\rangle_{\text{grill}}$,
$\langle a_{b}\rangle$, and $f_{p}$ is conducted. Figure 3 plots contours of
$\bar{S}^{2}_{J_{\phi}}$. It is evident that $\bar{S}^{2}_{J_{\phi}}$
decreases as turbulence intensity increases. Here, we define “turbulence
intensity” to mean the spatially averaged inverse scale length of density

$\langle\langle L^{-1}\rangle\rangle\equiv\langle\langle|\nabla
n|/n_{0}\rangle\rangle\approx f_{p}\langle n_{b}/n_{0}\rangle/\langle
a_{b}\rangle$ (13)

where $\langle\langle...\rangle\rangle$ denotes a spatial average. A minimum
$\bar{S}^{2}_{J_{\phi}}\approx 0.15$ is achieved for $\langle\langle
L^{-1}\rangle\rangle\approx 30\,\text{m}^{-1}$, signifying a near seven-fold
decrease in error. Interestingly, $\bar{S}^{2}_{J_{\phi}}$ remains near 0.15
for $\langle\langle L^{-1}\rangle\rangle\gtrsim 30\,\text{m}^{-1}$. This
suggests that the impact of scattering on LHCD saturates for sufficiently high
turbulence intensity, and that this condition is met in experiment.

Figure 3: Parametric scan of turbulence parameters in GENRAY/CQL3D for
discharge # 1101104011. Plotted are contours of $\bar{S}^{2}_{J_{\phi}}$
(defined in Eq. (12) with $X$ replaced by the toroidal current drive profile).
Green box: low turbulence case. Blue box: high turbulence case.

This is further illustrated by focusing on two points in the turbulence
parameter space. The point $\left[\langle
n_{b}/n_{0}\rangle_{\text{grill}},\langle
a_{b}\rangle,f_{p}\right]=[3,0.5\,\text{cm},0.2]$ is termed “high” turbulence,
as it results in saturated state (see blue box in Fig. 3). The point
$\left[\langle n_{b}/n_{0}\rangle_{\text{grill}},\langle
a_{b}\rangle,f_{p}\right]=[1.5,1\,\text{cm},0.2]$ is termed “low” turbulence,
as $\bar{S}^{2}_{J_{\phi}}\approx 0.6$ and is therefore not fully saturated
(see green box in Fig. 3). Figure 4 plots the corresponding current density
and HXR profiles. The reference case predicts a largely hollow current profile
with peaks at $\psi\approx 0.1$ and 0.8. Correspondingly, the HXR profile is
flat since this is a line-integrated diagnostic. These profiles are a poor
match to the experiment. In contrast, the high turbulence case fully mitigates
the off-axis peak in the current density profiles, and fills in the current
“valley” in the $\psi\approx 0.2-0.7$ region. (Note that small fluctuations in
the simulated “high turb.” current profile are a numerical artefact. They can
be minimized by further increasing the number of rays simulated. In reality,
fluctuations at this scale would be smoothed over by radial diffusion of fast
electrons.) Consequently, the HXR profile is peaked in the middle channels and
smoothly decays outwards. These profiles exhibit greatly improved agreement
with experiment, although there remains a $\sim 20\%$ underestimation of
current density on-axis. The low turbulence case is qualitatively a mix of the
saturated and reference case. Off-axis peaks are somewhat mitigated, and near-
axis current is also somewhat increased compared to the reference case.

Figure 4: Current density (left) and HXR (right) profiles calculated in
GENRAY/CQL3D for discharge # 1101104011. Low and high turbulence correspond to
points in Fig. 3. Note: CQL3D top-down symmetrizes the plasma profile along
the mid-plane, resulting in simulated HXR profiles roughly symmetric with
respect to the middle channels (16 and 17). For proper comparison, the
experimental HXR profile is also symmetrized.

###  3.3 Diffusion of ray-trajectories

The scattering-induced diffusion of rays can be observed in both real- and
phase-space. Figure 5 plots ray-trajectories in the poloidal plane during the
first pass. In the reference case, ray damping in the core is negligible. In
the high turbulence case, the initial ray bundle has greatly “fanned” due to
scattering. It should be emphasized that in Fig. 5 is only plotting ray-
trajectories during _first-pass_ , so what is shown is purely a scattering
effect, and not due to ray-stochasticity following multiple
reflections/passes. A fraction of rays scatter inwards, such that they
strongly Landau damp in the hot near-axis region. This seeds a super-thermal
electron tail near-axis, on which rays can further damp on subsequent passes.
This explains the significant increase in near-axis current drive when
scattering is included in the simulations.

Figure 5: Ray-trajectories during first pass through the core. Plotting rays
launched with $N_{||0}=1.6\pm 0.1$. Ray color denotes $\text{log}_{10}$ of
normalized ray power. Green patch denotes region with $\psi\leq 0.2$.

Figure 6 plots ray-trajectories in phase-space for the reference and high
turbulence case, showing only rays with $N_{||,0}=1.6\pm 0.1$. These rays are
confined between the strong ELD contour (Eq. (1)) in dashed-red, and the Stix-
Golant [38] accessibility contour in dashed-blue. After multiple passes, they
eventually damp near the ELD contour. Two observations are made about the
reference case. (1) Even though a wide sample of rays in $N_{||0}$ space is
plotted, the allowable phase-space is only partially and non-uniformly filled
in. (2) There is a distinct bundle of rays that strongly Landau damps at
$\psi\approx 0.8$. This causes the off-axis peak in the current density
profile in Figure 4. Compare this to the high turbulence case, where phase-
space is more uniformly filled in due to a strong diffusion of ray-
trajectories caused by scattering. One particular result of this is the
complete disappearance of the off-axis ray bundle.

Figure 6: Ray-trajectories in phase-space for the reference (left) and high
turbulence (right) cases. Plotting rays with $N_{||0}=1.6\pm 0.1$. The thick
black bar denotes start of rays, and ray color denotes imaginary component of
wave-vector due to Landau damping. Red dashed like: strong linear ELD
condition as defined in Eq. (1). Blue dashed like: the Stix-Golant
accessibility condition. For completeness, we also plot the generalized
accessibility condition [39, 40] in solid-black, though it is not a
particularly relevant limit in C-Mod.

###  3.4 High density partially inductive discharges

Two higher density discharges, with partially non-inductive current drive, are
modeled in GENRAY/CQL3D.

In the first discharge (#1140411016), 460kW of LH power launched at
$N_{||0}=1.9$ is coupled to an L-mode plasma with $\bar{n}_{e}=0.8\times
10^{20}\,\text{m}^{-3}$ and $T_{e0}=2.2\,$keV. During the LH phase,
$V_{\text{loop}}$ decreases from 1V to 0.4V. In the second discharge
(#1140411015), 450kW of LH power launched at $N_{||0}=1.9$ is coupled to
plasma with $\bar{n}_{e}=1.1\times 10^{20}\,\text{m}^{-3}$ and
$T_{e0}=1.8\,$keV. Correspondingly, loop voltage drops from 0.95V to 0.75V.
Note that this plasma is above the LHCD density limit [41], which explains the
relatively small change in loop voltage. The DC electric field is neglected
during CQL3D modeling of both discharges. There remain large uncertainties
about the DC field radial profile, which can significantly modify the current
profile. This can be accurately accounted for using time-dependent simulations
as done by Poli _et al._ , (2016) [42].

Figures 7 and 8 plot the simulated current density and HXR profiles for the
$\bar{n}_{e}=0.8$ and $1.1\times 10^{20}\,\text{m}^{-3}$ discharges,
respectively. The cases with turbulence use the same filament parameters as
the high turbulence case identified in Section 3.2. No attempt is made to
compare to experimental current densities. MSE measurements were not available
for these discharges to reconstruct the current profile. Even if it were
available, it is difficult to calculate the Ohmic and non-inductive
contributions separately because the DC electric field profile is not known.

Notably, for the $\bar{n}_{e}=0.8\times 10^{20}\,\text{m}^{-3}$ discharge, the
current valley at $\psi\approx 0.3$ disappears when scattering is modelled.
For both discharges, turbulence results in greater current drive near-axis
($\psi<0.4$), ultimately leading to improved agreement with HXR profile
measurements.

Figure 7: Current density (left) and HXR (right) profiles calculated in
GENRAY/CQL3D for discharge #1140411016 for which $\bar{n}_{e}=0.8\times
10^{20}\,\text{m}^{-3}$, $T_{e0}=2.2\,$keV, and $\left[\langle
n_{b}/n_{0}\rangle_{\text{grill}},\langle
a_{b}\rangle,f_{p}\right]=[3,0.5\,\text{cm},0.2]$. Figure 8: Current density
(left) and HXR (right) profiles calculated in GENRAY/CQL3D for discharge
#1140411015 for which $\bar{n}_{e}=1.1\times 10^{20}\,\text{m}^{-3}$,
$T_{e0}=1.8\,$keV, and $\left[\langle
n_{b}/n_{0}\rangle_{\text{grill}},\langle
a_{b}\rangle,f_{p}\right]=[3,0.5\,\text{cm},0.2]$.

Similar to the low-density fully non-inductive case, a scan of filaments
parameters is conducted for the $\bar{n}_{e}=0.8\times 10^{20}\,\text{m}^{-3}$
discharge in Fig. 9. Given that the experimental current drive profile cannot
be isolated, we instead plot $\bar{S}^{2}_{\text{HXR}}$, for which the array
$X$ in Eq. (12) is replaced with count rates from the 32-channel HXR
diagnostic. It is evident that the saturation effect of SOL turbulence on LHCD
is also present at high densities, in this case for $\langle\langle
L^{-1}\rangle\rangle\gtrsim 15\text{m}^{-1}$. Thus, it is confirmed that the
high turbulence parameters identified as being in the saturated regime for the
low density ($\bar{n}_{e}=0.52\times 10^{20}\,\text{m}^{-3}$) case, is also in
the saturated regime for higher density discharges.

Figure 9: Parametric scan of turbulence parameters in GENRAY/CQL3D for
discharge #1140411015. Plotted are contours of $\bar{S}^{2}_{\text{HXR}}$
(defined in Eq. (12) with $X$ replaced by the HXR profile).

###  3.5 Impact of scattering asymmetry

As mentioned in Sec. 3.1, the dominant asymmetric scattering effect is
$S\rightarrow S$ scatter. On average, this leads to rays rotating away from
the core. In turn, these rays propagate longer distances in the SOL, where
parasitic collisional damping is non-negligible. One might then expect that
less current is driven when asymmetric scatter is accounted for. A simple
gedanken experiment is conducted to test this hypothesis. Three cases are run:
(1) the reference case, (2) a case with the usual scattering model, and (3) a
case with the scattering model assuming reversed parity. Here, a reversal of
parity means setting $\sigma_{\text{eff},j\rightarrow
j^{\prime}}(\chi)\rightarrow\sigma_{\text{eff},j\rightarrow
j^{\prime}}(-\chi)$, such that, on average, rays scatter _towards_ the core.

First, this gedanken experiment is conducted for the low density fully non-
inductive C-Mod discharge. Figure 10 plots resulting current density profiles.
It is revealed that flipping the scattering parity does indeed affect the
profile. In the reverse parity case, a large current density peak forms at
$\psi\approx 0.2-04$. Relative to the normal scattering case, the reverse
parity case results in $\sim 30\%$ greater current drive. This trend is as
expected. This experiment is repeated with the $\bar{n}_{e}=1.1\times
10^{20}\,\text{m}^{-3}$ discharge, with Figure 11 plotting current density
profiles. A reversed scattering parity results in significantly greater near-
axis current deposition, and $\sim 25\%$ greater total current drive compared
to the normal scattering case.

Figure 10: Current density profiles calculated in GENRAY/CQL3D for discharge
#1101104011. Testing a reversed scattering parity. Horizontal dashed line
denotes experimental value. Figure 11: Current density profiles calculated in
GENRAY/CQL3D for discharge #1140411015. Testing a reversed scattering parity.
Simulated current profile not compared to experimental measurements because
the Ohmic contribution is neglected in the model.

Contrary to the hypothesis, the normal parity case results in greater current
than in the reference case. An explanation for this can be found in examining
Fig. 5. While there is a net rotation of rays away from the core, there is
still a large fraction of rays rotated towards the core such that they damp
on-axis during first pass. In contrast, virtually no power damps on-axis
during first pass in the reference case. Thus, this effect is enough to
counteract the net outward scatter.

##  4 Summary and discussion

A multiscale model has been developed for LH wave propagation in a magnetized
plasma with filamentary turbulence. Wave propagation is computed using ray-
tracing. Relatively expensive full-wave calculations are reserved only for
wave-filament interactions, which then act as stochastic kicks to the ray-
trajectories. The scattering probabilities are spatially dependent, which
allows the prescription of a radially and poloidally tapered SOL turbulence
profile in a tokamak geometry. In addition, both like- and unlike-mode
scattering is modelled.

This multiscale scattering model is applied to C-Mod discharges, where the
impact of scattering on LHCD is shown to saturate at sufficiently high levels
of turbulence intensity. In the saturated phase, ray-trajectories are
significantly diffused in phase-space. In turn, LHCD predictions are
insensitive to the assumed turbulence parameters. At low density
($\bar{n}_{e}\approx 0.52\times 10^{20}\,\text{m}^{-3}$), the saturated case
results in excellent agreement with experimental MSE and HXR measurements.
Here, a turbulence intensity of $\langle\langle L^{-1}\rangle\rangle\gtrsim
30\,\text{m}^{-1}$ resulted in saturation with $\bar{S}^{2}_{J_{\phi}}<0.15$,
signifying a near seven-fold decrease in error for the current profile
(compared to the simulation with no scattering). This saturation is also
present in the high density case above the LHCD density limit. Here,
saturation at $\langle\langle L^{-1}\rangle\rangle\gtrsim 15\,\text{m}^{-1}$
results in significantly improved match to experiment
($\bar{S}^{2}_{\text{HXR}}<0.1$). This suggests that, experimentally, the
impact of wave-filament scatter is saturated in C-Mod L-mode discharges for
all background densities. In future work, it will be worthwhile to study how
the threshold $\langle\langle L^{-1}\rangle\rangle$ depends on plasma and LH
launch parameters.

While this modeling suggests that C-Mod LH launch is in this saturation
regime, it is unclear whether this is corroborated by edge turbulence
measurements. Analysis of Langmuir probe measurements in the SOL of C-Mod
measure $\langle n_{b}/n_{0}\rangle\approx 2.1$ [43]. Similarly, GPI
measurements indicate $\langle a_{b}\rangle\approx 1\,$cm [44] and $f_{p}\sim
0.05$ [45]. However, these burst statistic analyses require setting an
arbitrary threshold above r.m.s. values for identifying a filament. This
threshold affects the computed $\langle n_{b}/n_{0}\rangle$, $\langle
a_{b}\rangle$, and greatly affects $f_{p}$. Given the large experimental
uncertainties in measuring these parameters, no attempt is made here in
prescribing the likeliest filament PDF or packing fraction. Rather, it is
argued that SOL measurements in C-Mod do not rule out the possibility of a
saturated regime for LH wave scattering.

The impact of asymmetric scattering was tested by reversing the scattering
parity. The reversed parity leads to increased scatter of LH power into the
core. As expected, this results in greater current drive and HXR count
compared to simulations with the correct scattering parity. It is important to
note, however, that scatter with correct parity, which results in net outward
scatter, still leads to greater current drive than the case with no
scattering. This is attributed to the non-zero fraction of LH power that
scatters inwards and Landau damps near-axis on first pass.

In regards to the three C-Mod discharges modelled in this work, it is
interesting to note that the simulated HXR count, for both the reference and
high turbulence cases, agree within an order of magnitude with the
experimental values. This is in stark contrast to previous analysis, in which
standard ray-tracing/Fokker-Planck simulations could not replicate the
anomalous drop in total current driven and HXR count for $\bar{n}_{e}\gtrsim
10^{20}\,\text{m}^{-3}$ [41, 11]. However, later simulations including a more
realistic collisional SOL could successfully do so [37]. This present work
uses this improved SOL model. Overall, this suggests that proper modeling of
SOL collisional losses is needed to accurately calculate _total_ current and
HXR count, while proper modeling of scattering is needed to accurately
calculate the _radial profile_ of current and HXR count.

Lastly, it should be emphasized that these multiscale simulations allow, for
the first time, the coupling of wave-filament scattering physics, beyond ray-
tracing, to the down-stream physics of core damping and current drive. This is
of practical importance for the following reason. It has been shown, through
the use of this multiscale model, that accurate simulation of LHCD requires
accounting for wave-filament scatter from sufficiently intense SOL turbulence.
This corresponds to a situation in which ray-tracing is invalid, and so
higher-order scattering effects must be taken into account. This is evidenced
by the relatively poor performance of the filament refraction model which
employed ray-tracing alone [28]. While global full-wave simulations are
computationally in-feasible, this multiscale approach applies full-wave
calculations only locally where it is needed (at the wave-filament
interaction). In this way, fast parametric scans and intershot analysis is
made possible on modest, university-scale computing clusters. For example,
consider Fig. 3. This parametric scan consisted of 196 simulations and
required $\sim 3500$ total CPU-hours on the Engaging cluster. (Equivalently,
$\sim 17.5$ CPU-hours were required per simulation.)

## Acknowledgements

This work was supported by US DoE under contract numbers: DE-SC0018090
supporting the RF-SciDAC 4 project, DE-SC0014264 supporting PSFC MFE projects,
and DE-FC02-04ER54698 supporting DIII-D projects. The GENRAY/CQL3D and PETRA-M
simulations presented in this paper were performed on the MIT-PSFC partition
of the Engaging cluster at the MGHPCC facility (www.mghpcc.org) which was
funded by DoE grant number DE-FG02-91-ER54109.

## Declaration of interest

The authors report no conflict of interest.

##  Appendix A Radially and poloidally tapered turbulence

GENRAY uses $(\psi,\theta_{p})$-coordinates in the poloidal plane, where
$\psi\equiv\sqrt{\Psi_{\phi}/\Psi_{\phi,a}}$ is the normalized radial
coordinate (and $\Psi_{\phi}$ is toroidal magnetic flux). The evaluation of
coordinate $\psi$ is approximately extended outside the separatrix. It is
approximated that $n(\mathbb{r})\approx n(\psi,\theta_{p})$. Given the
relatively slow spatial variation of $B$, it is acceptable for our purposes of
calculating scattering probabilities in the SOL to assume

$B\approx B_{\phi}\approx\frac{B_{0}}{1+\frac{a}{R_{0}}\cos(\theta_{p})}$
(A.1)

where $B_{\phi}$ is the toroidal magnetic field and $a,R_{0}$ are the minor,
major plasma radius.

In assuming that the turbulence parameters vary radially and poloidally (i.e.
$p(n_{b}/n_{0}(\mathbb{r}),a_{b}(\mathbb{r}))=p(n_{b}/n_{0},a_{b};\psi,\theta_{p})$),
one can parametrize the scattering probabilities as
$\hat{\sigma}_{\text{eff},j\rightarrow
j^{\prime}}(\chi;\psi,\theta_{p},k_{||})$ and $\Sigma_{\text{eff},j\rightarrow
j^{\prime}}(\chi;\psi,\theta_{p},k_{||})$. These form the lookup tables in
GENRAY.

As in [29], it is assumed that the filament joint-PDF is a bivariate skewed
Gaussion[46], now with radial and poloidal dependence. Its analytic form is
written as a function of the 2D column vector
$\zeta=\left[n_{b}/n_{0};a_{b}\right]$, while parametrization depending on
radial and poloidal coordinates is made explicit.

$\displaystyle p(\zeta;\psi,\theta_{p})$
$\displaystyle=2\phi_{2}(\zeta-Q(\psi,\theta_{p})\Omega(\psi,\theta_{p}))\Phi(\Gamma^{\text{T}}(\psi,\theta_{p})\left[\zeta-Q(\psi,\theta_{p})\right])$
(A.2a) $\displaystyle\Omega(\psi,\theta_{p})$
$\displaystyle=\begin{bmatrix}s_{n_{b}/n_{0}}^{2}(\psi,\theta_{p})&\eta
s_{n_{b}/n_{0}}(\psi,\theta_{p})s_{a_{b}}\\\ \eta
s_{n_{b}/n_{0}}(\psi,\theta_{p})s_{a_{b}}&s_{a_{b}}^{2}\end{bmatrix}$ (A.2b)
$\displaystyle Q(\psi,\theta_{p})$
$\displaystyle=\langle\zeta\rangle(\psi,\theta_{p})-\sqrt{\frac{2}{\pi}}\delta(\psi,\theta_{p})$
(A.2c) $\displaystyle\delta(\psi,\theta_{p})$
$\displaystyle=\frac{1}{\sqrt{1+\Gamma^{\text{T}}\Omega(\psi,\theta_{p})\Gamma}}\Omega(\psi,\theta_{p})\Gamma$
(A.2d) $\displaystyle\Gamma(\psi,\theta_{p})$
$\displaystyle=\begin{bmatrix}\Gamma_{n_{b}/n_{0}}(\psi,\theta_{p})\\\
\Gamma_{a_{b}}\end{bmatrix}$ (A.2e)

where $\phi_{2}(\zeta,\Omega)$ is a 2D Gaussian PDF with zero mean and
correlation matrix $\Omega$, $\Phi(x)$ is the 1D Gaussian CDF for scalar input
$x$, $\Gamma$ is now a 2D column vector of skewness factors, and
$\langle\zeta\rangle\equiv\left[\langle n_{b}/n_{0}\rangle;\langle
a_{b}\rangle\right]$ is the 2D column vector of mean $n_{b}/n_{0}$ and
$a_{b}$. Similarly, $s\equiv[s_{n_{b}/n_{0}};s_{a_{b}}]$ are standard
deviations, and $\eta$ is the scalar correlation coefficient. It is therefore
possible to parameterize $p(\zeta)$ as a function of
$\langle\zeta\rangle,\Gamma,s$, $\eta$, $\psi$, and $\theta_{p}$. In contrast
to [29], the filament relative density mean ($\langle n_{b}/n_{0}\rangle$),
standard deviation ($s_{n_{b}/n_{0}}$), and skewness ($\Gamma_{n_{b}/n_{0}}$)
are allowed to vary with $\psi$ and $\theta_{p}$. For simplicity, the filament
width statistics ($\langle a_{b}\rangle$, $s_{a_{b}}$, and $\Gamma_{a_{b}}$)
are kept homogeneous. Further, it is assumed that $s_{a_{b}}=0.2\,$cm and
$\Gamma_{a_{b}}=7$. The correlation coefficient $\eta=0.9$, signifying strong
positive correlation between $n_{b}/n_{0}$ and $a_{b}$ [47].

The radial profile for fluctuations should smoothly increase from zero inside
the seperatrix to the user-specified
$\langle\frac{n_{b}}{n_{0}}\rangle_{\text{grill}}$ at the LH grill. The
following analytic forms are assumed for $\langle n_{b}/n_{0}\rangle$,
$s_{n_{b}/n_{0}}$ and $\Gamma_{n_{b}/n_{0}}$:

$\displaystyle\langle n_{b}/n_{0}\rangle$ $\displaystyle=\left[\langle
n_{b}/n_{0}\rangle_{\text{min}}+(\langle
n_{b}/n_{0}\rangle_{\text{grill}}-\langle
n_{b}/n_{0}\rangle_{\text{min}})\bar{\psi}-1\right]^{g(\theta_{p})}+1$ (A.3a)
$\displaystyle s_{n_{b}/n_{0}}$
$\displaystyle=\left[s_{n_{b}/n_{0},\text{min}}+(s_{n_{b}/n_{0},\text{grill}}-s_{n_{b}/n_{0},\text{min}})\bar{\psi}\right]^{g(\theta_{p})}$
(A.3b) $\displaystyle\Gamma_{n_{b}/n_{0}}$
$\displaystyle=\left[\Gamma_{n_{b}/n_{0},\text{min}}+(\Gamma_{n_{b}/n_{0},\text{grill}}-\Gamma_{n_{b}/n_{0},\text{min}})\bar{\psi}\right]^{g(\theta_{p})}$
(A.3c)

where

$\bar{\psi}(\psi)=\begin{cases}0&\psi<\psi_{\text{min}}\\\
(\psi-\psi_{\text{min}})/(\psi_{\text{grill}}-\psi_{\text{min}})&\psi_{\text{min}}<\psi<\psi_{\text{grill}}\\\
1&\psi>\psi_{\text{grill}}\end{cases}$ (A.4)

and $\psi_{\text{grill}}$ is the radial location of the LH launcher (a.k.a.
grill). In this study, the minimum radial location of non-zero turbulence
($\psi_{\text{min}}$) is set to 0.85. It is assumed that $\langle
n_{b}/n_{0}\rangle_{\text{min}}=1$, $s_{n_{b}/n_{0},\text{min}}=0.1$,
$\Gamma_{n_{b}/n_{0},\text{min}}=0.01$, and
$\Gamma_{n_{b}/n_{0},\text{max}}=7$. The following linear relation is used to
prescribe $s_{n_{b}/n_{0},\text{grill}}$:

$s_{n_{b}/n_{0},\text{grill}}=c_{1}\left\langle\frac{n_{b}}{n_{0}}\right\rangle_{\text{grill}}+c_{2}$
(A.5)

where $c_{1}$ and $c_{2}$ are set such that $s_{n_{b}/n_{0},\text{grill}}=1.1$
for $\langle\frac{n_{b}}{n_{0}}\rangle_{\text{grill}}=4$ and
$s_{n_{b}/n_{0},\text{grill}}=0.22$ for
$\langle\frac{n_{b}}{n_{0}}\rangle_{\text{grill}}=1.1$.

Lastly, the function $g(\theta_{p})$ introduces a poloidal variation and is
prescribed to be

$g(\theta_{p})=0.1+0.9\text{cos}^{2}(\theta_{p}/2)$ (A.6)

This accounts for experimental measurements that density fluctuations on the
HFS are greatly suppressed [48, 49].

Figure 12 plots background density and example turbulence profiles at the
outer mid-plane for the non-inductive C-Mod discharge. In this case, the LH
grill position corresponds to $\psi_{\text{grill}}\approx 1.02$. The free-
parameter $\langle n_{0}/n_{b}\rangle_{\text{grill}}=3$ prescribes mean
relative density of filaments at the LH grill. The mean relative density of
filaments monotonically decrease to unity at $\psi=0.85$, meaning they are
indistinguishable from the background, thereby smoothly transitioning into a
quiescent core plasma.

Figure 12: Background density and modeled turbulence profile at the outer mid-
plane in a C-Mod discharge. The only user-specified parameter is $\langle
n_{0}/n_{b}\rangle_{\text{grill}}=3$. The $\langle n_{0}/n_{b}\rangle$ and
$s_{n_{b}/n_{0}}$ profiles are prescribed according to Eq. (A.3-A.5).

Figure 13 shows an example of how the joint-PDF at $\psi=1.02$ (in the SOL) is
tapered poloidally such that filaments are, on average, less dense on the
high-field side.

Figure 13: Joint-PDF of filament parameters for $\psi=1.02$. The only user-
specified parameters are $\langle n_{0}/n_{b}\rangle_{\text{grill}}=3$ and
$\langle a_{b}\rangle=1\,$cm. Red cross-hairs denote $\pm$ 1 standard
deviation in $n_{b}/n_{0}$ and $a_{b}$.

## References

  * [1] N. J. Fisch, “Theory of current drive in plasmas,” Reviews of Modern Physics, vol. 59, no. 1, pp. 175–234, 1987.
  * [2] S. Ide, O. Naito, T. Oikawa, T. Fujita, T. Kondoh, M. Seki, K. Ushigusa, and J.-. Team, “LHCD current profile control experiments towards steady state improved confinement on JT-60U,” Nuclear Fusion, vol. 40, pp. 445–451, Mar. 2000.
  * [3] Y. Peysson and T. S. Team, “High power lower hybrid current drive experiments in the Tore Supra tokamak,” Nuclear Fusion, vol. 41, pp. 1703–1713, Nov. 2001.
  * [4] X. Litaudon, F. Crisanti, B. Alper, J. F. Artaud, Y. F. Baranov, E. Barbato, V. Basiuk, A. Bécoulet, M. Bécoulet, C. Castaldo, C. D. Challis, G. D. Conway, R. Dux, L. G. Eriksson, B. Esposito, C. Fourment, D. Frigione, X. Garbet, C. Giroud, N. C. Hawkes, P. Hennequin, G. T. A. Huysmans, F. Imbeaux, E. Joffrin, P. J. Lomas, P. Lotte, P. Maget, M. Mantsinen, J. Mailloux, D. Mazon, F. Milani, D. Moreau, V. Parail, E. Pohn, F. G. Rimini, Y. Sarazin, G. Tresset, K. D. Zastrow, M. Zerbini, and c. t. t. E.-J. Workprogramme, “Towards fully non-inductive current drive operation in JET,” Plasma Physics and Controlled Fusion, vol. 44, pp. 1057–1086, July 2002.
  * [5] V. Pericoli-Ridolfini, E. Barbato, S. Cirant, H. Kroegler, L. Panaccione, S. Podda, F. Alladio, B. Angelini, M. L. Apicella, G. Apruzzese, L. Bertalot, A. Bertocchi, M. Borra, G. Bracco, A. Bruschi, G. Buceti, P. Buratti, A. Cardinali, C. Centioli, R. Cesario, S. Ciattaglia, V. Cocilovo, F. Crisanti, R. De Angelis, F. De Marco, B. Esposito, D. Frigione, L. Gabellieri, G. Gatti, E. Giovannozzi, C. Gourlan, G. Granucci, F. Gravanti, M. Grolli, A. Imparato, M. Leigheb, L. Lovisetto, G. Maffia, G. Maddaluno, M. Marinucci, G. Mazzitelli, P. Micozzi, F. Mirizzi, S. Nowak, F. P. Orsitto, D. Pacella, M. Panella, L. Pieroni, G. B. Righetti, F. Romanelli, F. Santini, M. Sassi, S. E. Segre, A. Simonetto, C. Sozzi, E. Sternini, O. Tudisco, V. Vitale, G. Vlad, N. Tartoni, B. Tilia, A. A. Tuccillo, V. Zanza, M. Zerbini, and F. Zonca, “High Plasma Density Lower-Hybrid Current Drive in the FTU Tokamak,” Physical Review Letters, vol. 82, pp. 93–96, Jan. 1999.
  * [6] S. Shiraiwa, O. Meneghini, R. Parker, G. Wallace, J. Wilson, I. Faust, C. Lau, R. Mumgaard, S. Scott, S. Wukitch, W. Beck, J. Doody, J. Irby, P. MacGibbon, D. Johnson, A. Kanojia, P. Koert, D. Terry, R. Vieira, and the Alcator C-Mod team, “Design, and initial experiment results of a novel LH launcher on Alcator C-Mod,” Nuclear Fusion, vol. 51, p. 103024, Oct. 2011.
  * [7] R. J. La Haye, “Neoclassical tearing modes and their control,” Physics of Plasmas, vol. 13, no. 5, p. 055501, 2006.
  * [8] A. H. Reiman and N. J. Fisch, “Suppression of Tearing Modes by Radio Frequency Current Condensation,” Physical Review Letters, vol. 121, p. 225001, Nov. 2018.
  * [9] P. Bonoli, “Linear theory of lower hybrid heating,” IEEE transactions on plasma science, vol. 12, no. 2, pp. 95–107, 1984.
  * [10] P. T. Bonoli and E. Ott, “Toroidal and scattering effects on lower-hybrid wave propagation,” Physics of Fluids, vol. 25, p. 359, Mar. 1982.
  * [11] R. T. Mumgaard, Lower Hybrid Current Drive on Alcator C-Mod: Measurements with an Upgraded MSE Diagnostic and Comparisons to Simulation.  PhD thesis, MIT, 2015.
  * [12] Y. Peysson and T. S. Team, “Progress towards high-power lower hybrid current drive in TORE SUPRA,” Plasma Physics and Controlled Fusion, vol. 42, no. 12B, p. B87, 2000.
  * [13] A. Garofalo, X. Gong, J. Qian, J. Chen, G. Li, K. Li, M. Li, X. Zhai, P. Bonoli, D. Brower, L. Cao, L. Cui, S. Ding, W. Ding, W. Guo, C. Holcomb, J. Huang, A. Hyatt, M. Lanctot, L. Lao, H. Liu, B. Lyu, J. McClenaghan, Y. Peysson, Q. Ren, S. Shiraiwa, W. Solomon, Q. Zang, and B. Wan, “Development of high poloidal beta, steady-state scenario with ITER-like tungsten divertor on EAST,” Nuclear Fusion, vol. 57, p. 076037, July 2017.
  * [14] Y. Peysson, D. Mazon, J.-F. Artaud, A. Ekedahl, L. Delpech, J. Hillairet, T. Hoang, X. L. Zou, WEST Team, X. Y. Bai, Y. P. Zhang, HL-2A Team, K. Krol, J. Bielecki, A. Jardin, M. Scholz, D. Dworak, and J. Decker, “Lower Hybrid Current Drive in High Aspect Ratio Tokamaks,” Journal of Fusion Energy, vol. 39, pp. 270–291, Oct. 2020.
  * [15] Y. Peysson, J. Decker, E. Nilsson, J.-F. Artaud, A. Ekedahl, M. Goniche, J. Hillairet, B. Ding, M. Li, P. T. Bonoli, S. Shiraiwa, and M. Madi, “Advances in modeling of lower hybrid current drive,” Plasma Physics and Controlled Fusion, vol. 58, p. 044008, Apr. 2016.
  * [16] S. Baek, B. Biswas, G. Wallace, P. Bonoli, B. Ding, M. Li, Y. Li, Y. Wang, M. Wang, C. Wu, G. Yan, J. Chen, X. Zhai, A. Garofalo, W. Choi, F. Poli, and S. Shiraiwa, “A model investigation of the impact of lower hybrid wave scattering angle on current drive profile in EAST and Alcator C-Mod,” Nuclear Fusion, vol. 61, p. 106034, Oct. 2021.
  * [17] M. Porkolab, “Theory of parametric instability near the lower-hybrid frequency,” Physics of Fluids, vol. 17, no. 7, p. 1432, 1974.
  * [18] R. Cesario, A. Cardinali, C. Castaldo, F. Paoletti, and D. Mazon, “Modeling of a Lower-Hybrid Current Drive by Including Spectral Broadening Induced by Parametric Instability in Tokamak Plasmas,” Physical Review Letters, vol. 92, p. 175002, Apr. 2004.
  * [19] M. Madi, Y. Peysson, J. Decker, and K. Y. Kabalan, “Propagation of the lower hybrid wave in a density fluctuating scrape-off layer (SOL),” Plasma Physics and Controlled Fusion, vol. 57, no. 12, p. 125001, 2015.
  * [20] R. Cesario, L. Amicucci, A. Cardinali, C. Castaldo, M. Marinucci, F. Napoli, F. Paoletti, D. De Arcangelis, M. Ferrari, A. Galli, G. Gallo, E. Pullara, G. Schettini, and A. Tuccillo, “Spectral broadening of parametric instability in lower hybrid current drive at a high density,” Nuclear Fusion, vol. 54, p. 043002, Apr. 2014.
  * [21] S. Baek, R. Parker, P. Bonoli, S. Shiraiwa, G. Wallace, B. LaBombard, I. Faust, M. Porkolab, and D. Whyte, “High density LHRF experiments in Alcator C-Mod and implications for reactor scale devices,” Nuclear Fusion, vol. 55, p. 043009, Apr. 2015.
  * [22] P. M. Bellan and K. L. Wong, “Effect of density fluctuations on lower hybrid resonance cone propagation,” Physics of Fluids, vol. 21, no. 4, pp. 592–599, 1978.
  * [23] Y. Peysson, J. Decker, L. Morini, and S. Coda, “RF current drive and plasma fluctuations,” Plasma Physics and Controlled Fusion, vol. 53, no. 12, p. 124028, 2011.
  * [24] P. L. Andrews and F. W. Perkins, “Scattering of lower-hybrid waves by drift-wave density fluctuations: Solutions of the radiative transfer equation,” Physics of Fluids, vol. 26, no. 9, pp. 2537–2545, 1983.
  * [25] N. Bertelli, G. Wallace, P. T. Bonoli, R. W. Harvey, A. P. Smirnov, S. G. Baek, R. R. Parker, C. K. Phillips, E. J. Valeo, J. R. Wilson, and J. C. Wright, “The effects of the scattering by edge plasma density fluctuations on lower hybrid wave propagation,” Plasma Physics and Controlled Fusion, vol. 55, p. 074003, July 2013.
  * [26] S. J. Zweben, D. P. Stotler, J. L. Terry, B. Labombard, M. Greenwald, M. Muterspaugh, C. S. Pitcher, K. Hallatschek, R. J. Maqueda, B. Rogers, J. L. Lowrance, V. J. Mastrocola, and G. F. Renda, “Edge turbulence imaging in the Alcator C-Mod tokamak,” Physics of Plasmas, vol. 9, no. 5, pp. 1981–1989, 2002.
  * [27] J. P. Graves, J. Horacek, R. A. Pitts, and K. I. Hopcraft, “Self-similar density turbulence in the TCV tokamak scrape-off layer,” Plasma Phys. Control. Fusion, vol. 47, pp. 1–9, 2005.
  * [28] B. Biswas, S. Gyou Baek, P. Bonoli, S. Shiraiwa, G. Wallace, and A. White, “Study of turbulence-induced refraction of lower hybrid waves using synthetic scrape-off layer filaments,” Plasma Physics and Controlled Fusion, vol. 62, p. 115006, Nov. 2020.
  * [29] B. Biswas, S. Shiraiwa, S.-G. Baek, P. Bonoli, A. Ram, and A. E. White, “A hybrid full-wave Markov chain approach to calculating radio-frequency wave scattering from scrape-off layer filaments,” Journal of Plasma Physics, vol. 87, no. 5, p. 35, 2021.
  * [30] S. J. Frank, J. P. Lee, J. C. Wright, I. H. Hutchinson, and P. T. Bonoli, “Verifying raytracing/Fokker-Planck lower-hybrid current drive predictions with self-consistent full-wave/Fokker-Planck simulations,” arXiv, July 2022.
  * [31] A. K. Ram and K. Hizanidis, “Scattering of radio frequency waves by cylindrical density filaments in tokamak plasmas,” Physics of Plasmas, vol. 23, no. 2, p. 022504, 2016.
  * [32] W. Zhang, W. Tierens, V. Bobkov, A. Cathey, I. Cziegler, M. Griener, M. Hoelzl, O. Kardaun, and V. Bobkov, “Interaction between filaments and ICRF in the plasma edge,” Nuclear Materials and Energy, vol. 26, p. 100941, Mar. 2021.
  * [33] J. R. Myra and D. A. D’Ippolito, “Scattering of radio frequency waves by blob-filaments,” Physics of Plasmas, vol. 17, no. 10, p. 102510, 2010.
  * [34] M. I. Mishchenko, Electromagnetic Scattering by Particles and Particle Groups: An Introduction.  Cambridge: Cambridge University Press, 2014.
  * [35] A. Smirnov and R. Harvey, “The GENRAY Ray Tracing Code,” tech. rep., CompX, 2001.
  * [36] R. W. Harvey and M. G. Mccoy, “The CQL3D Fokker-Planck Code (reconstituted in editable form from General Atomics Report GA-A20978, 1992),” tech. rep., 2005.
  * [37] S. Shiraiwa, S. G. Baek, I. Faust, G. Wallace, P. Bonoli, O. Meneghini, R. Mumgaard, R. Parker, S. Scott, R. W. Harvey, B. J. Ding, M. H. Li, S. Y. Lin, and C. Yang, “Impact of SOL plasma profiles on lower hybrid current drive: Experimental evidence, mitigation and modeling approaches,” in RADIO FREQUENCY POWER IN PLASMAS: Proceedings of the 21st Topical Conference, (California, USA), p. 030016, 2015.
  * [38] T. H. Stix, Waves in Plasmas.  Springer Science & Business Media, 1992.
  * [39] H. Takahashi, “The generalized accessibility and spectral gap of lower hybrid waves in tokamaks,” Physics of Plasmas, vol. 1, no. 7, pp. 2254–2276, 1994\.
  * [40] F. Paoletti, D. W. Ignat, J. Kesner, S. Bernabei, R. Kaita, B. Leblanc, F. M. Levinton, and S. C. Luckhardt, “Lower hybrid current drive accessibility study with reconstructed magnetic equilibria in PBX-M,” Nuclear Fusion, vol. 34, no. 6, pp. 771–776, 1994.
  * [41] G. M. Wallace, A. E. Hubbard, P. T. Bonoli, I. C. Faust, R. W. Harvey, J. W. Hughes, B. L. Labombard, O. Meneghini, R. R. Parker, A. E. Schmidt, S. Shiraiwa, A. P. Smirnov, D. G. Whyte, J. R. Wilson, J. C. Wright, and S. J. Wukitch, “Lower hybrid current drive at high density in Alcator C-Mod,” Nuclear Fusion, vol. 51, no. 8, p. 083032, 2011.
  * [42] F. M. Poli, P. T. Bonoli, M. Chilenski, R. Mumgaard, S. Shiraiwa, G. M. Wallace, R. Andre, L. Delgado-Aparicio, S. Scott, J. R. Wilson, R. W. Harvey, Y. V. Petrov, M. Reinke, I. Faust, R. Granetz, J. Hughes, and J. Rice, “Experimental and modeling uncertainties in the validation of lower hybrid current drive,” Plasma Physics and Controlled Fusion, vol. 58, no. 9, p. 095001, 2016.
  * [43] R. Kube, A. Theodorsen, O. E. Garcia, B. Labombard, and J. L. Terry, “Fluctuation statistics in the scrape-off layer of Alcator C-Mod,” Plasma Physics and Controlled Fusion, vol. 58, no. 5, p. 054001, 2016.
  * [44] R. Kube, O. Garcia, B. LaBombard, J. Terry, and S. Zweben, “Blob sizes and velocities in the Alcator C-Mod scrape-off layer,” Journal of Nuclear Materials, vol. 438, pp. S505–S508, July 2013.
  * [45] S. J. Zweben, J. L. Terry, B. Labombard, M. Agostini, M. Greenwald, O. Grulke, J. W. Hughes, D. A. D’Ippolito, S. I. Krasheninnikov, J. R. Myra, D. A. Russell, D. P. Stotler, and M. Umansky, “Estimate of convective radial transport due to SOL turbulence as measured by GPI in Alcator C-Mod,” Journal of Nuclear Materials, vol. 415, no. 1, pp. S463–S466, 2011.
  * [46] A. Azzalini and A. Capitanio, “Statistical applications of the multivariate skew normal distribution,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 61, pp. 579–602, Aug. 1999.
  * [47] G. Decristoforo, F. Militello, T. Nicholas, J. Omotani, C. Marsden, N. Walkden, and O. E. Garcia, “Blob interactions in 2D scrape-off layer simulations,” Physics of Plasmas, vol. 27, p. 122301, Dec. 2020.
  * [48] J. L. Terry, S. J. Zweben, K. Hallatschek, B. Labombard, R. J. Maqueda, B. Bai, C. J. Boswell, M. Greenwald, D. Kopon, W. M. Nevins, C. S. Pitcher, B. N. Rogers, D. P. Stotler, X. Q. Xu, R. Boivin, S. M. Kaye, and I. Introducton, “Observations of the turbulence in the scrape-off-layer of Alcator C-Mod and comparisons with simulation,” Physics of Plasmas, vol. 10, no. 5 II, pp. 1739–1747, 2003.
  * [49] N. Smick, B. Labombard, and C. S. Pitcher, “Plasma profiles and flows in the high-field side scrape-off layer in Alcator C-Mod,” Journal of Nuclear Materials, vol. 337–339, no. 1-3 SPEC. ISS., pp. 281–285, 2005.