2.1. Ansatz to electric field

Using separation of variables (see Appendix A), the electric field can be written in cylindrical coordinates as

(2.1a)$$\begin{gather} E_{j \gamma} = \exp({\textrm{i}(k_{||}z-\omega t)})\sum_{m={-}\infty}^{+\infty} E_{jm} W_{j \gamma m} \exp({\textrm{i} m \theta}); \quad \gamma = \rho,\theta,z, \end{gather}$$

(2.1b)$$\begin{gather}W_{j\rho m} = \xi_{jx} \textrm{J}_{m}^{\prime}(k_{j \perp}\rho) - \textrm{i}\xi_{jy} \frac{m}{k_{j \perp} \rho} \textrm{J}_{m}(k_{j \perp}\rho), \end{gather}$$

(2.1c)$$\begin{gather}W_{j\theta m} = \textrm{i}\xi_{jx} \frac{m}{k_{j \perp} \rho} \textrm{J}_{m}(k_{j \perp}\rho) + \xi_{jy} \textrm{J}_{m}^{\prime}(k_{j \perp}\rho), \end{gather}$$

(2.1d)$$\begin{gather}W_{jzm} = \textrm{i} \xi_{jz}\textrm{J}_{m}(k_{j \perp}\rho), \end{gather}$$

where $j=0,\ldots,4$ is the wave index, $\boldsymbol{\xi}_j = \{\xi _{jx},\xi _{jy},\xi _{jz}\}$ is the plane-wave polarization of wave $j$, $\textrm {J}_{m}$ is the Bessel function of the first kind and order $m$, and $\textrm {J}_{m}^{\prime }$ is the first derivative of $\textrm {J}_{m}$ with respect to its argument. For the known incident plane wave ($\kern1.5pt j=0$), it is required that

(2.2)\begin{equation} E_{0m} = i^{m-1}, \end{equation}

as $E_{jm}$ for $j>0$ have yet to be determined.

Note that $k_{||}$ is the same for all waves, and is fixed by the incident wave. This is a consequence of Snell's law applied to a medium that is constant along the $z$-direction.

2.2. Boundary conditions

In general, the solutions in (2.1) can have both J and Y terms, where Y is the Bessel function of the second kind. The requirement that $\boldsymbol {E}$ is finite at $\rho = 0$ leads to Y terms being zero for the slow and fast branch inside the filament. For $\rho \rightarrow \infty$, the scattered fields must be radiating away from the filament. For the scattered fast wave, this requires the use of Hankel functions of the first kind, $H^{1}$, instead of J in (2.1). The LH slow wave is backward propagating, which means $\boldsymbol {k}_{\perp }$ and $\boldsymbol {v}_{gr,\perp }$ are anti-parallel. This requires the use of Hankel functions of the second kind, $H^{2}$, for the scattered slow wave.

It should be noted that a backward-propagating incident wave (i.e. the slow wave) has $k_{\textrm {inc},x} = -k_{\perp }$. This flipped sign can most easily be accounted for by substituting $\textrm {J}_{m} \rightarrow \textrm {J}_{-m}$ (Myra & D'Ippolito Reference Myra and D'Ippolito2010).

Lastly, a system of equations must be formulated to determine coefficients $E_{jm}$ for $j=1,\ldots,4$. This is accomplished by imposing the four independent Maxwell boundary conditions at $\rho = a_b$ (see Appendix B). For each poloidal mode number, there are four unknown coefficients and four boundary conditions, which results in a solvable system of equations.

3.1. Modified system of equations

Remember that the ‘flat-top’ ($R=0$) system is solvable because there exist four unknowns $E_1,E_2,E_3,E_4$ and four independent boundary equations for each mode number $m$. A similar system of equations must be derived for the general $R>0$ case. The simplest case ($R=1$) is illustrated in figure 2. In the intermediate layers ($0< r \leq R$), each wave branch is generally a function of both $H_{m}^{1}$ and $H_{m}^{2}$ terms, and are therefore split into these two electric field contributions.

Figure 2. The different contributions of $\boldsymbol {E}$ in the SAS model for a filament with two radial bins ($R=1$). There are two regions of interest: background and filament. The filament is further divided into the central (cntr) and intermediate (mid) regions. Here ‘$S$’ and ‘$F$’ denote slow and fast LH branches. Functions in parentheses denote the type of Bessel function used in (2.1).

With the considerations made above, there are now eight unknown waves for the $R=1$ case. There are also two boundaries, which supply four boundary conditions each. This results in a solvable system for the total electric field everywhere.

For the general case ($R>0$), there are all together $4(R+1)$ equations for each mode number. For more details, see Appendix C. Solving for the total electric field requires inverting $4(R+1) \times 4(R+1)$ matrices a total of $2M+1$ times, where $M$ is the maximum mode number chosen to truncate the series. These matrices are sparse and banded, which results in fast solution times of the order of seconds on a single central processing unit (CPU).

3.2. Poloidal and radial resolution

The electric field will evolve on three possible length scales: $a_b, k_{\text {in} \perp }^{-1}$ or $k_{\text {out} \perp }^{-1}$. Subscripts ‘in’, ‘out’ denote inside, outside the cylinder. Define a characteristic poloidal mode number $\tilde {m} \equiv \text {max}(k_{\text {in} \perp },k_{\text {out} \perp })a_{b}$. If $\tilde {m} \gtrsim 1$, terms with $|m|>\tilde {m}$ will rapidly decay in magnitude. Therefore, for a converged solution, it is necessary that $M \gg \tilde {m}$. If $\tilde {m}\ll 1$, it is necessary that $M \gg 1$.

A rule-of-thumb can also be derived for how many radial bins are required. The source of error stems from discretizing the cylinder's smoothly varying radial profile into homogeneous radial bins. The discontinuity between bins is of the order of $L_{n}^{-1} \Delta r$, where $L_n$ is the characteristic length of the density inhomogeneity. Here, $\Delta r$ is the radial bin width. Assume the cylinder has a monotonically decreasing radial profile with characteristic radial width $a_b$. Then the discontinuity is approximately $L_{n}^{-1} \Delta r \sim ( {n_b}/{n_0 a_b}) ( {a_b}/{R}) = {n_b}/{n_0 R}$. Here ${n_b}/{n_0}$ is the ratio of the peak (centre) cylinder density to the background density, and $R$ is the number of radial bins. To ensure the discontinuities are small requires $R \gg {n_b}/{n_0}$.

4.1. Effective differential scattering width: $\sigma _{\text {eff}}(\theta )$

In the presence of a filament, the resulting scattered field will depend on $a_b$ and $n_{b}$. Therefore, $\sigma (\theta )= \sigma (\theta ; n_b/n_0, a_b)$ for a given $n_0, \boldsymbol {B},N_{||}$ and $\omega$. A joint probability distribution function (p.d.f.), $p(n_b/n_0,a_b)$, is introduced. This is the probability that a filament will have certain parameters $a_b$ and $n_b/n_0$. Taking a weighted average of $\sigma (\theta )$ with this joint-p.d.f. returns the statistically averaged $\sigma (\theta )$ for scattering from a randomly selected filament. This averaged, or ‘effective’ differential scattering width is defined as

(4.5)\begin{equation} \sigma_{\text{eff}}(\theta) = \int_{0}^{\infty} \textrm{d}a_{b} \int_{0}^{\infty}\textrm{d} \left(n_{b}/n_0\right) \sigma_{S}(\theta; n_b/n_0, a_b) p(n_b/n_0,a_b), \end{equation}

where $p(n_b/n_0,a_b)$ is normalized such that $\int _{0}^{\infty } \textrm {d}a_{b} \int _{0}^{\infty }\textrm {d}(n_{b}/n_0)p(n_b/n_0,a_b) = 1$. Again, $\sigma _{F}$ is neglected because the focus is on SOL plasmas.

The choice of joint-p.d.f. for the filament parameters is guided by experimental measurements. In the SOL of C-Mod, mirror Langmuir probe measurements reveal positively skewed p.d.f.s of density fluctuations (Graves et al. Reference Graves, Horacek, Pitts and Hopcraft2005). SOL fluid codes predict that the filament width and density are positively correlated (Decristoforo et al. Reference Decristoforo, Militello, Nicholas, Omotani, Marsden, Walkden and Garcia2020). Therefore, the joint-p.d.f. is reasonably well described by a positively skewed normal p.d.f. for each parameter ($n_b/n_0$ and $a_b$) along with a positive bivariate correlation. The analytic form for this p.d.f., as a function of a two-dimensional (2-D) column vector $\boldsymbol{\zeta} = [ n_b/n_0; a_b ]$ is (Azzalini & Capitanio Reference Azzalini and Capitanio1999)

(4.6a)$$\begin{gather} p(\boldsymbol{\zeta}) = 2 \phi_{2}(\boldsymbol{\zeta} - \boldsymbol{Q},\boldsymbol{\varOmega}) \varPhi(\boldsymbol{\varGamma}^{\text{T}} \boldsymbol{\cdot} (\boldsymbol{\zeta}-\boldsymbol{Q})), \end{gather}$$

(4.6b)$$\begin{gather}\boldsymbol{\varOmega} = \begin{bmatrix} s_{n_b/n_0}^{2} & \eta s_{n_b/n_0} s_{a_b}\\ \eta s_{n_b/n_0} s_{a_b} & s_{a_b}^{2} \end{bmatrix}, \end{gather}$$

(4.6c)$$\begin{gather}\boldsymbol{Q} = \langle \boldsymbol{\zeta} \rangle - \sqrt{\frac{2}{\rm \pi}} \boldsymbol{\delta} \end{gather}$$

(4.6d)$$\begin{gather}\boldsymbol{\delta} = \frac{1}{\sqrt{1 + \boldsymbol{\varGamma}^{\text{T}}\boldsymbol{\cdot} \boldsymbol{\varOmega}\boldsymbol{\cdot} \boldsymbol{\varGamma}}} \boldsymbol{\varOmega} \boldsymbol{\cdot} \boldsymbol{\varGamma}, \end{gather}$$

(4.6e)$$\begin{gather}\boldsymbol{\varGamma} = \begin{bmatrix} \varGamma_{n_b/n_0}\\ \varGamma_{a_b} \end{bmatrix}, \end{gather}$$

where $\phi _{2}(\boldsymbol{\zeta},\boldsymbol{\varOmega})$ is a 2-D normal p.d.f. with zero mean and correlation matrix $\boldsymbol{\varOmega}$. Here $\varPhi (x)$ is the 1-D normal cumulative distribution function (CDF) for scalar input $x$, $\boldsymbol{\varGamma}$ is a 2-D column vector of skewness factors and $\langle \boldsymbol{\zeta} \rangle \equiv [ \langle n_b/n_0 \rangle ; \langle a_b\rangle ]$ is the 2-D column vector of mean $n_b/n_0$ and $a_b$. Similarly, $\boldsymbol{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 parametrize $p(\boldsymbol{\zeta})$ as a function of $\langle \boldsymbol{\zeta} \rangle, \boldsymbol{\varGamma}, \boldsymbol{s}$ and $\eta$.

Filament mean width $\langle a_b \rangle$ is well bounded by gas puff imaging (GPI) measurements as well as theory/simulation (Zweben et al. Reference Zweben, Stotler, Terry, LaBombard, Greenwald, Muterspaugh, Pitcher, Group, Hallatschek and Maqueda2002; Krasheninnikov, D'Ippolito & Myra Reference Krasheninnikov, D'Ippolito and Myra2008; Keramidas Charidakos et al. Reference Keramidas Charidakos, Myra, Ku, Churchill, Hager, Chang and Parker2020). Langmuir probe and GPI measurements provide a rough lower bound on the filament mean relative density $\langle n_b/n_0 \rangle$, though this value will vary significantly at different radial locations in the SOL (Zweben et al. Reference Zweben, Stotler, Terry, LaBombard, Greenwald, Muterspaugh, Pitcher, Group, Hallatschek and Maqueda2002; Terry et al. Reference Terry, Zweben, Hallatschek, LaBombard, Maqueda, Bai, Boswell, Greenwald, Kopon and Nevins2003).

5.1. Solution to RTE using a Markov chain

Equation (5.2) is an integro-differential equation, and cannot, in general, be solved analytically. One numerical method is to discretize the wave spectrum into photons/rays, and stochastically evolve their trajectories, as per the standard Monte Carlo technique. This method is rather slow in providing a converged wave spectrum near $\theta = \pm {\rm \pi}/2$, where tally counts are usually low. Given how simple the slab geometry is, a more elegant absorbing Markov chain method can be employed. This Markov chain (MC) method is deterministic, so it avoids the low tally count problem. It is commonly used to solve for a reflected and transmitted wave spectrum through a turbid slab (e.g. solar rays interacting with Earth's atmosphere). The present study closely follows the formalism by (Esposito & House Reference Esposito and House1978) and Xu et al. (Reference Xu, Davis, West and Esposito2011).

Power is incident on the turbulent slab from the left, directed along the $x$-direction ($\theta _{0}=0$). It is convenient to define a $\zeta \equiv |\cos {\theta }|$ such that $\zeta _{0}=1$. It is simple to calculate the fraction of transmitted power that does not scatter in the slab. This is the ‘ballistic’ fraction:

(5.3)\begin{equation} P_{\textrm{ball}} = \exp\left({-\frac{L_x}{\zeta_{0}}\varSigma_{\text{eff}}}\right).\end{equation}

It is also straightforward to calculate the transmitted and reflected fraction that only scatter once in the slab:

(5.4)\begin{align} P_{T,1}(\theta) &= \hat{\sigma}(\theta-\theta_{0}) \exp\left({-\frac{L_{x}}{\zeta}\varSigma_{\text{eff}}}\right)\nonumber\\ &\quad \times \begin{cases} 0 & \cos{\theta} < 0,\\ \dfrac{\zeta \zeta_{0}}{\zeta-\zeta_{0}}\left[1- \exp\left({-L_{x}\varSigma_{\text{eff}}\left(\dfrac{1}{\zeta_0}-\dfrac{1}{\zeta}\right)}\right)\right] & \theta \neq \theta_{0},\\ \varSigma_{\text{eff}}L_{x} & \theta = \theta_0, \end{cases} \end{align}

(5.5)\begin{equation} P_{R,1}(\theta) = \hat{\sigma}(\theta-\theta_{0}) \times \begin{cases} 0 & \cos{\theta} > 0,\\ \dfrac{\zeta \zeta_{0}}{\zeta+\zeta_{0}}\left[1-\exp\left({-}L_{x}\varSigma_{\text{eff}}\left(\dfrac{1}{\zeta_0}+\dfrac{1}{\zeta}\right)\right)\right] & \text{otherwise}. \end{cases}\end{equation}

The above are the first-order scattering terms. To compute the higher-order terms (fraction of power undergoing $> 1$ scattering events), it is necessary to use the MC method. The slab is discretized into $n=1,\ldots,{N}$ segments of width $\Delta x=L_{x}/{N}$. The angular spectrum is also discretized into $m=1,\ldots,{M}$ segments of width $\Delta \theta =2{\rm \pi} /{M}$. Next, the $\text {NM}\times \text {NM}$ ‘transition’ matrix $\boldsymbol{T}(x_n,\theta _m;x_{n^{\prime }}\theta _{m^{\prime }})$ is generated, which accounts for the probability of a photon in segment $n$ and directed along $\theta _m$ to scatter in segment $n'$ into $\theta _{m^{\prime }}$. The ${N}\times {M}$ ‘source’ matrix $\boldsymbol{\varPi} (x_n,\theta _m)$ is defined as the probability distribution of photons in segment $n$ directed along $\theta _{m}$ right after the first scattering event. Lastly, the $\text {NM}\times {M}$ ‘absorption’ matrix $\boldsymbol{R}_{T/R}(x_n,\theta _m;\theta _{m^{\prime }})$ is the probability of a photon to escape the slab via transmission/reflection following its final scattering event in segment $n$ from $\theta _m$ to $\theta _{m^{\prime }}$. The form for these matrices are as follows. The ‘transition’ matrix can be broken into four components:

(5.6)\begin{equation} \boldsymbol{T}(x_n,\theta_m;x_{n^{\prime}}\theta_{m^{\prime}}) = p_{\text{esc}}(x_n,\theta_m)\,p_{\text{trvl}}(x_n, \theta_m, x_{n'})\, p_{\text{sct}}(x_{n'},\theta_m)\,\hat{\sigma}(\theta_{m'}-\theta_{m}), \end{equation}

where

(5.7a)$$\begin{gather} p_{\text{esc}}(x_n,\theta_m) = \frac{\zeta_m}{\varSigma_{\text{eff}}\Delta x} \left(1-\exp\left({-\frac{\Delta x}{\zeta_m}\varSigma_{\text{eff}}}\right)\right), \end{gather}$$

(5.7b)$$\begin{gather}p_{\text{trvl}}(x_n,\theta_m,x_{n'}) = \begin{cases} \exp\left({-\dfrac{x_{n'}-x_n}{\zeta_m}\varSigma_{\text{eff}}}\right) & \dfrac{x_{n'}-x_n}{\cos{\theta_m}} \geq 0,\\ 0 & \dfrac{x_{n'}-x_n}{\cos{\theta_m}} < 0,\end{cases} \end{gather}$$

(5.7c)$$\begin{gather}p_{\text{sct}}(x_{n'},\theta_m) = 1 - \exp\left({-\frac{\Delta x}{\zeta_m}\varSigma_{\text{eff}}}\right). \end{gather}$$

The escape probability, $p_{\text {esc}}(x_n,\theta _m)$, is the probability for a photon to travel through segment $n$ without scattering. The travel probability, $p_{\text {trvl}}(x_n,\theta _m,x_{n'})$, is the probability of travelling between segments $n$ and $n'$ without scattering. Note that $p_{\text {trvl}}$ is set to zero in cases where the photon in segment $n$ with $\theta _m$ is oriented such that it is travelling away from $n'$. The scatter probability, $p_{\text {sct}}(x_{n'},\theta _m)$, is the probability of scattering within segment $n'$. Lastly, $\hat {\sigma }(\theta _{m'}-\theta _{m})$ is the probability of the photon rotating from $\theta _m$ to $\theta _{m'}$ given that it undergoes a scattering event. The ‘source’ matrix is

(5.8)\begin{gather} \hspace{-18pc} \boldsymbol{\varPi} (x_n,\theta_m) = \hat{\sigma}(\theta_m-\theta_{0})C^{{-}1}\nonumber\\ \hspace{2.5pc}\times \begin{cases} \dfrac{\zeta_m}{\zeta_m-\zeta_0} \exp\left({-\left(\dfrac{x_n}{\zeta_0}+\dfrac{\Delta x}{\zeta_m}\right)\varSigma_{\text{eff}}}\right) & {}\\ \quad\times\left(1-\exp\left({-\varSigma_{\text{eff}}\Delta x \left(\dfrac{1}{\zeta_{0}}- \dfrac{1}{\zeta_m}\right)}\right)\right) & \cos{\theta_m} > 0,\\ \dfrac{\zeta_0 \zeta_m}{\zeta_0 + \zeta_m} \exp\left({-\dfrac{x_n}{\zeta_0}\varSigma_{\text{eff}}}\right) & {}\\ \quad\times\left(1 - \exp\left({-\varSigma_{\text{eff}}\Delta x \left(\dfrac{1}{\zeta_0}+\dfrac{1}{\zeta_m}\right)} \right)\right) & \cos{\theta_m} < 0,\\ \dfrac{\varSigma_{\text{eff}} \Delta x }{\zeta_0} \exp\left({-\dfrac{x_n + \Delta x}{\zeta_0}\varSigma_{\text{eff}}}\right) & \theta_m = \theta_0, \end{cases} \end{gather}

(5.9)\begin{gather} C = \frac{\zeta_m}{\varSigma_{\text{eff}}\Delta x} \left(1-\exp\left({-\frac{\Delta x}{\zeta_m}\varSigma_{\text{eff}}}\right)\right),\end{gather}

where the coefficient $C$ is required to properly volume-average the source over segment $n$. The transmission and reflection ‘absorption’ matrices are

(5.10)\begin{align} \boldsymbol{R}_{T}(x_n,\theta_m;\theta_{m'}) &= \hat{\sigma}(\theta_{m'}-\theta_{m}) \nonumber\\ &\quad \times \begin{cases} 0 & \cos{\theta_{m'}} < 0,\\ \dfrac{\zeta_{m'}}{\zeta_{m'}-\zeta_m} \exp\left({-\varSigma_{\text{eff}}\left(\dfrac{L_x}{\zeta_{m'}} - \dfrac{x_n}{\zeta_{m}}\right)}\right) & {}\\ \quad\times\left(\exp\left({-\varSigma_{\text{eff}}x_n \left(\dfrac{1}{\zeta_{m'}} - \dfrac{1}{\zeta_{m}}\right)}\right) \right. & {}\\ \quad \left.- \exp\left({-\varSigma_{\text{eff}}L_x \left(\dfrac{L_x}{\zeta_{m'}} - \dfrac{x_n}{\zeta_{m}}\right)}\right)\right) & \cos{\theta_{m}} > 0,\\ \dfrac{\zeta_{m'}}{\zeta_m+\zeta_{m'}}\exp\left({-\varSigma_{\text{eff}} \left(\dfrac{L_x}{\zeta_{m'}}+\dfrac{x_n}{\zeta_m}\right)}\right) & {}\\ \quad \times\left(\exp\left({\varSigma_{\text{eff}}x_n\left(\dfrac{1}{\zeta_m}+ \dfrac{1}{\zeta_{m'}}\right)}\right)\right. & {}\\ \quad - \left.\exp\left({\varSigma_{\text{eff}}L_x\left(\dfrac{1}{\zeta_m}-\dfrac{1}{\zeta_{m'}}\right)}\right)\right) & \cos{\theta_{m}} < 0,\\ \dfrac{\varSigma_{\text{eff}}}{\zeta_m} \exp\left({-\varSigma_{\text{eff}}\left(\dfrac{L_x}{\zeta_{m'}}-\dfrac{x_n}{\zeta_m}\right)}\right) (L_x - x_n) & \theta_m = \theta_{m'}, \end{cases} \end{align}

(5.11)\begin{align} \boldsymbol{R}_{R}(x_n,\theta_m;\theta_{m'}) &= \hat{\sigma}(\theta_{m'}-\theta_{m}) \nonumber\\ &\quad \times \begin{cases} 0 & \cos{\theta_{m'}} > 0,\\ \dfrac{\zeta_{m'}}{\zeta_{m}-\zeta_{m'}} \exp\left({-\varSigma_{\text{eff}}\dfrac{x_n}{\zeta_{m}}}\right) & {}\\ \quad \times\left(1 - \exp\left({-\varSigma_{\text{eff}}x_n \left(\dfrac{1}{\zeta_{m'}} - \dfrac{1}{\zeta_{m}}\right)}\right)\right) & \cos{\theta_{m}} < 0,\\ \dfrac{\zeta_{m'}}{\zeta_m+\zeta_{m'}} \exp\left({\varSigma_{\text{eff}}\dfrac{x_n}{\zeta_{m'}}}\right) & {}\\ \quad \times \left(\exp\left({-\varSigma_{\text{eff}}x_n\left(\dfrac{1}{\zeta_m}+\dfrac{1}{\zeta_{m'}}\right)}\right) - 1\right) & \cos{\theta_{m}} > 0,\\ \dfrac{\varSigma_{\text{eff}}}{\zeta_m} \exp\left({-\dfrac{\varSigma_{\text{eff}}}{\zeta_{m'}}x_n}\right)x_n & \theta_m = \theta_{m'}. \end{cases} \end{align}

Using these matrices, one can calculate the higher-order transmitted/reflected wave spectrum terms:

(5.12a)$$\begin{gather} \boldsymbol{P}_{T/R,2}(\theta_m) = \boldsymbol{\varPi} \boldsymbol{\cdot} \boldsymbol{I} \boldsymbol{\cdot} \boldsymbol{R}_{T/R}, \end{gather}$$

(5.12b)$$\begin{gather} \boldsymbol{P}_{T/R,l}(\theta_m) = \boldsymbol{\varPi} \boldsymbol{\cdot} \boldsymbol{T}^{l-2} \boldsymbol{\cdot} \boldsymbol{R}_{T/R} \quad \text{for } l \geq 3, \end{gather}$$

where $l$ is the order of the scattering term and $\boldsymbol{I}$ is the $\text {NM}\times \text {NM}$ identity matrix. In summing all scattering terms, the total transmitted and reflected wave spectrum is

(5.13a)$$\begin{gather} \boldsymbol{P}_{T}(\theta_m) - \frac{P_{\textrm{ball}}}{\Delta \theta}\delta_{\theta_m,\theta_0} = \boldsymbol{P}_{T,1}(\theta_m) + \sum_{l=1}^{\infty} \boldsymbol{P}_{T,l}, \end{gather}$$

(5.13b)$$\begin{gather} \boldsymbol{P}_{R}(\theta_m) = \boldsymbol{P}_{R,1}(\theta_m) +\sum_{l=1}^{\infty} \boldsymbol{P}_{R,l}, \end{gather}$$

and $\delta _{i,j}$ is the Kronecker delta. Furthermore, (5.13a) can be rewritten as

(5.14)\begin{equation} \boldsymbol{P}_{T}(\theta_m) - \frac{P_{\textrm{ball}}}{\Delta \theta}\delta_{\theta_m,\theta_0} - \boldsymbol{P}_{T,1}(\theta_m)= \boldsymbol{\varPi} \boldsymbol{\cdot} \left( \boldsymbol{I} + \sum_{l=1}^{\infty} \boldsymbol{T}^{l} \right) \boldsymbol{\cdot} \boldsymbol{R}_{T} = \boldsymbol{\varPi} \boldsymbol{\cdot} \left( \boldsymbol{I} - \boldsymbol{T}\right)^{{-}1}\boldsymbol{\cdot} \boldsymbol{R}_{T}.\end{equation}

A similar form applies to (5.13b). The second relation in (5.14) produces the solution following a matrix inversion. In practice, it is often faster to evaluate the first relation and truncate the series at a finite $l$ when the solution is sufficiently converged (Yang et al. Reference Yang, Xiao, Li and Ma2018).

In deriving $p_{\text {esc}}$ and $p_{\text {sct}}$, the possibility of multiple scattering events in segment $n$ is neglected. This is a reasonable assumption as long as $\Delta x \equiv {L_x}/{{N}} \ll {\zeta }/{\varSigma _{\text {eff}}}$. The population of photons with $\zeta \approx 0$ is the largest source of error for any finite $N$. Nevertheless, in practice, $\boldsymbol{P}_{T/R}$ is found to converge as long as $\varSigma _{\text {eff}} \Delta x \ll 1$. A criterion for the angular resolution is not as straightforward. It depends on the smoothness of $\hat {\sigma }(\theta )$. Naturally, a fine resolution is needed to accurately resolve sharp peaks in $\hat {\sigma }(\theta )$.

6.1. Scattered field for single filament

Consider an incident slow wave with a prescribed frequency $f$ and parallel refractive index $N_{||}\equiv {c k_{||}}/{\omega }$. Also assume a filament with Gaussian radial profile such that

(6.1)\begin{equation} n(\rho) -n_0= n_0 \left( \frac{n_b}{n_0} -1 \right) \exp\left({-\left(\frac{2\sqrt{\ln(2)}\rho}{a_b}\right)^{2}}\right), \end{equation}

where $n_{b}/n_0$ is the relative density at the filament's peak ($\rho = 0$) and $a_b$ is re-defined as the full-width at half-maximum of the filament. A case with $n_0=1\times 10^{19}\ \text {m}^{-3}$, ${B}=4\ \text {T}$, $f=4.6\ \text {GHz}$, $N_{||} =2$, $n_b/n_0 = 4.8$ and $a_b=1\ \text {cm}$ is simulated using the SAS model. Simulation resolution is $R=22$ and $M=100$. Figure 3(a–c) shows the $(x,y,z)$ components of the time-averaged Poynting flux $\boldsymbol {S}$ exterior to the filament. The $\boldsymbol {S}$ is calculated using the relation $\boldsymbol {S} = ({1}/{2 \mu _0 \omega })\textrm {Im}\,{\boldsymbol {E}^{*} \times (\boldsymbol {\nabla } \times \boldsymbol {E})}$. The normalized field $\boldsymbol {P} \equiv {\boldsymbol {S}}/{|\boldsymbol {S}_{\textrm {inc}}|}$ is introduced, where $\boldsymbol {S}_{\textrm {inc}}$ is the Poynting flux of the incident wave. Figure 3 reveals a shadowing effect downstream of the filament. The striations in the field indicate strong back and side scattering of the incident wave. Note that $\boldsymbol {P}$ has been numerically computed from the interpolated $\boldsymbol {E}$-field on a grid in the $(x,y)$-plane. This method of plotting $\boldsymbol {P}$ is susceptible to large errors inside the filament where the gradients of $\boldsymbol {E}$ are large. Therefore, $\boldsymbol {P}$ inside the filament is not plotted.

Figure 3. Gaussian filament: normalized Poynting flux computed using the SAS model (a–c) and PETRA-M (d–f), where $n_0=1\times 10^{19} \ \text {m}^{-3}$, ${B}=4 \ \text {T}$, $f=4.6\ \text {GHz}$, $N_{||} =2$, $n_b/n_0 = 4.8$, and $a_b=1\ \text {cm}$.

Figure 3(d–f) shows this case repeated in PETRA-M, and result in an excellent agreement with the SAS model. This 2-D simulation is done in a circular domain. The incident wave is excited to the left of the filament using an external current source term. To approximate an infinite background plasma, a perfectly matched layer (PML) is modelled at the perimeter of the circular simulation domain.

6.2. Reflection coefficient for turbulent slab

Next, the MC step of the SAS-MC model is compared with a turbulent slab modelled in PETRA-M. This will reveal whether the far-field limit, a critical approximation in the MC model, is valid for the treatment of LH scattering in the SOL. In theory, the far-field approximation should break down as filaments are packed closer together (Mishchenko Reference Mishchenko2014).

For the MC model, $\sigma _{\text {eff}}(\theta )$ must be calculated, which first requires prescribing a joint-p.d.f. of filaments. Figure 4 plots an example joint-p.d.f. A bivariate skewed-normal distribution is assumed for $a_b$ and $n_b/n_0$ (see (4.6)). In this case, $\langle a_b \rangle = 0.48\ \textrm {cm}$ and $\langle n_b/n_0 \rangle = 2.6$. These values are bounded by experimental SOL measurements (Zweben et al. Reference Zweben, Stotler, Terry, LaBombard, Greenwald, Muterspaugh, Pitcher, Group, Hallatschek and Maqueda2002; Terry et al. Reference Terry, Zweben, Hallatschek, LaBombard, Maqueda, Bai, Boswell, Greenwald, Kopon and Nevins2003). (Assuming SOL turbulence is predominantly filamentary, the approximation $\langle n_b/n_0 \rangle \approx 1 + ({n_{\text {RMS}}}/{n})f_{p}^{-1}$ is made, where $f_{p}$ is the packing fraction and is assumed to be $0.2$.) Filament relative density and size skewness are prescribed via the skewness parameters: $\boldsymbol{\varGamma} = [7,10]$. Filament relative density and size standard deviation are prescribed as $\boldsymbol{s} = [1.35, 0.1\ \text {cm}]$. The bivariate correlation coefficient $\eta =0.9$.

Figure 4. Example joint-p.d.f. of filament parameters $n_b/n_0$ and $a_b$. A bivariate skewed-normal distribution is assumed (see (4.6)), where $\langle n_b/n_0 \rangle = 2.6$ and $\langle a_b \rangle = 0.48\,$cm. Filament density and size skewness $\boldsymbol{\varGamma} = [7,10]$. Filament density and standard deviation $\boldsymbol{s}=[1.35, 0.1\ \text {cm}]$. Correlation coefficient $\eta =0.9$. Dashed black line denotes $n_b/n_0 =1$.

Figure 5 shows the simulation set-up for a slab turbulent geometry in PETRA-M. A slow wave travelling in the $+x$-direction interacts with a turbulent layer populated with Gaussian filaments. The filaments are randomly generated in the slab using a Monte Carlo approach (Sierchio et al. Reference Sierchio, Cziegler, Terry, White and Zweben2016; Biswas et al. Reference Biswas, Baek, Bonoli, Shiraiwa, Wallace and White2020). Each filament is generated with a randomly picked $n_b/n_0$ and $a_b$ with a probability that satisfies the prescribed joint-p.d.f. $p(n_b/n_0, a_b)$. In this way, the turbulent slab used in the SAS-MC model and in PETRA-M are made statistically equivalent. The incident slow wave is excited with an external current density source function upstream of the turbulence. Top and bottom boundaries are periodic. The turbulence is also periodic in the $y$-direction. To minimize the effect of the periodic geometry on the wave, the $y$ length of the solution domain is much larger than $\langle a_b \rangle$ or $k_{\textrm {inc}\perp }$. To mimic an infinite domain in the $\pm x$-direction, the left and right boundaries can be modelled as either a perfectly matched layer (PML) or an absorbing boundary condition (ABC). The PML, while more computationally efficient, does not work when both the slow and fast wave can propagate in the background plasma.

Figure 5. Set-up for PETRA-M simulation with turbulent slab. Here, ‘PML’ denotes perfectly matched layer and ‘PEC’ denotes perfect electric conductor.

In PETRA-M, the reflection coefficient $F_{\text {ref}}=1-P_{x}/P_{x,0}$ is calculated and directly compared with the SAS-MC value. Here, $P_{x}$ is the $x$-component of the Poynting flux and $P_{x,0}$ is the ‘nominal’ value when no turbulence is present. In the SAS-MC model, $F_{\text {ref}}=\int P_{{R}}(\theta )\,\textrm {d}\theta$.

Table 1 compares $F_{\text {ref}}$ computed in the PETRA-M and SAS-MC models. Both models follow the same general trend. Cases 1–4 and 8–11 reveal that $F_{\text {ref}}$ increases with $f_{p}$. Cases 3,5,6,8 reveal that $F_{\text {ref}}$ increases with $\langle n_b/n_0 \rangle$ and decreases with $\langle a_b \rangle$. This is consistent with previous scattering theories (Ott Reference Ott1979; Andrews & Perkins Reference Andrews and Perkins1983). Another way to analyse the trends in $F_{\text {ref}}$ between models is by inspecting $\varSigma _{\text {eff}}L_{x}$ calculated in the SAS-MC model. This is the attenuation factor for the ballistic power (see (5.3)). As $\varSigma _{\text {eff}}L_{x}$ increases, so should $F_{\text {ref}}$. Indeed, this is true for both models.

Table 1. Comparison between SAS-MC model and PETRA-M. Slow wave launched at $4.6\ \text {GHz}$ and $N_{||}=2$ with ${B}=4\ \text {T}$. Filament joint-p.d.f. parameters are same as in figure 4 unless otherwise noted. The $\varSigma _{\text {eff}}L_{x}$ is calculated using the SAS-MC model. For each case, results from multiple iterations (with different turbulence realizations) are averaged until $F_{\text {ref}}$ is statistically converged.

In general, the SAS-MC model over-predicts $F_{\text {ref}}$, such that the absolute error $\Delta F_{\text {ref}} \equiv F_{\text {ref,SAS-MC}} - F_{\text {ref,PETRA-M}} \geq 0$. This error increases with $f_p$. Again, this arises from the far-field approximation breaking down. While far-field validity is dependent on $f_{p}$, the aggregate error depends on $\varSigma _{eff} L_{x}$. For example, the SOL width is varied between cases 2 and 7, while the turbulence is kept statistically identical. The $L_{x}=5\ \textrm {cm}$ case results in $\Delta F_{\text {ref}} = 0.00$. For the $L_{x}=15\ \textrm {cm}$ case, $\Delta F_{\text {ref}} = +0.13$.

SOL measurements indicate $f_p \approx 0.05\text {--}0.25$ (Agostini et al. Reference Agostini, Zweben, Cavazzana, Scarin, Serianni, Maqueda and Stotler2007; Carralero et al. Reference Carralero, Artene, Bernert, Birkenmeier, Faitsch, Manz, de Marne, Stroth, Wischmeier and Wolfrum2018; Zweben et al. Reference Zweben, Terry, LaBombard, Agostini, Greenwald, Grulke, Hughes, D'Ippolito, Krasheninnikov and Myra2011). SOL widths are also ${<}5\ \textrm {cm}$ in present-day devices. As a result, cases 2 and 9 are most representative of a C-Mod SOL, depending on whether the background density is evaluated at the far-SOL or the separatrix, respectively. At low-density (case 2), the two models agree well ($\Delta F_{\text {ref}} = 0.00$). At high density (case 9), the SAS-MC model over-predicts $F_{\text {ref}}$, such that $\Delta F_{\text {ref}} = +0.08$. A possible reason for the disagreement at high density may be because $\varSigma _{\text {eff}}L_{x}$ is greater (compared with similar cases at low density).

6.3. Comments on computational cost

Generally, the semi-analytic scattering method has three key advantages compared with finite-element Maxwell solvers: (1) the large (in fact infinite) background plasma region does not need to be meshed; (2) it exactly solves scattering problems, because the infinite exterior domain does not need to be artificially truncated; (3) analysing the scattered wave spectrum is straight-forward, because the solution is already deconvolved into the constituent poloidal mode numbers for each branch.

Points 1 and 2 result in the SAS-MC model being considerably less expensive than slab turbulence simulations in PETRA-M. Using the SAS technique, computing a single differential scattering width $\sigma (\theta ;n_b/n_0, a_b )$ takes ${\sim }10\ \textrm {s}$ on a single CPU, and considerably less time if parallelized between poloidal mode numbers. Computing $\sigma _{\text {eff}}(\theta )$ may require sampling a few hundred combinations of $(n_b/n_0, a_b)$, depending on the filament joint-p.d.f. Fortunately, each sampled $\sigma (\theta ;n_b/n_0, a_b )$ needs to be computed only once. Any number of $\sigma _{\text {eff}}(\theta )$ can then be generated from the sampled differential scattering widths.

The most expensive process in the MC routine is generating the transition matrix T. This takes ${\sim }20\ \textrm {s}$ on a single CPU, depending on poloidal and radial bin resolution.

Using PETRA-M, each $n_0 = 2.25 \times 10^{19}\ \text {m}^{-3}$, $L_{x} = 5\ \textrm {cm}$ slab case required ${\sim }25$ CPU-hours and ${\sim }300\ \textrm {GB}$ of RAM on the MIT Engaging computing cluster. The size of PETRA-M simulations is primarily limited by available RAM. In comparison, all computations for the SAS-MC model were conducted on a PC with 8 GB of available RAM.

6.4. Caveats to the SAS-MC model

The SAS-MC model offers higher physics fidelity than ray tracing, while being computationally less expensive than numeric full-wave solvers. This is possible owing to a number of assumptions made in the model that makes it less universally applicable than numeric full-wave solvers.

The SAS model assumes a homogeneous background plasma with a cylindrical scattering object (the filament) that is poloidally and azimuthally symmetric. This allows an efficient solution scheme using separation of variables. In reality, filaments usually develop a shock front as they convect outward (Krasheninnikov et al. Reference Krasheninnikov, D'Ippolito and Myra2008), and the resulting crescent-like filament shape can lead to significantly modified scattering behaviour, at least for ion cyclotron waves (Tierens, Zhang & Manz Reference Tierens, Zhang and Manz2020a). Furthermore, the filament in the SAS model is assumed to be aligned with the magnetic field, so that $\boldsymbol {\nabla }_{||} ( {n}/{n_0} )=0$. This is likely a reasonable assumption because filaments introduce a $\boldsymbol {\nabla }_{||} ( {n}/{n_0} )$ that is much smaller than $k_{||}$ of the LH wave (Grulke et al. Reference Grulke, Terry, Cziegler, LaBombard and Garcia2014). As a result, the effect of $k_{||}$ broadening arising from a typical SOL filament is small (Madi et al. Reference Madi, Peysson, Decker and Kabalan2015).

The MC model introduces additional assumptions. The SOL is treated as a slab, which means the effect of toroidal geometry is neglected. This is a reasonable assumption for the treatment of first-pass scattering in front of the antenna. In addition, the background plasma and turbulence parameters are constant within the slab, when in reality, they are sensitive to the radial coordinate in a tokamak. The MC model also assumes the reflected wave spectrum is lost, when in reality, a fraction of this power may once again reflect at a cutoff and re-enter the core plasma. Lastly, the MC model assumes the filaments are far enough apart so that the RTE is valid. This assumption is increasingly poor as $f_p$ rises.

To quantify the inaccuracies introduced by the MC model, it may be worthwhile to model scattering along the full extent of a ray trajectory in a realistic tokamak geometry. This can be done by employing $\sigma (\theta ; n_0, B, N_{||}, \langle n_b/n_0 \rangle, \langle a_b \rangle )$ as a scattering probability in a Monte Carlo ray tracing simulation, similar to what has been done for the $k$-scattering model (Bonoli & Ott Reference Bonoli and Ott1982; Bertelli et al. Reference Bertelli, Wallace, Bonoli, Harvey, Smirnov, Baek, Parker, Phillips, Valeo and Wilson2013). Alternatively, three-dimensional (3-D) PETRA-M simulations of LH launching in a turbulent SOL, with realistic geometry, would address all the caveats mentioned. These tests are outside the scope of this paper.

7.1. Parametric scan of SOL density and filament parameters

Figure 7(a,c,e) plots $\sigma _{S}$ as a function of background density $n_0$, relative filament density $n_b/n_0$ and filament width $a_b$. The incident wave is at 4.6 GHz with $N_{||}=2$, typical for LH launching in C-Mod. In accordance with the low-field-side SOL in C-Mod, ${B}=4\ \text {T}$. As expected, $\sigma _{S}$ increases as $n_b/n_0$ deviates from unity. Notably, scattering resonances can be seen, as indicated by bands of higher $\sigma _{S}$. These arise from standing wave resonances excited within the filament, which result in stronger coupling to the scattered waves. Expressions for these resonances can be analytically derived for the ‘flat-top’ filament case (Myra & D'Ippolito Reference Myra and D'Ippolito2010; Tierens et al. Reference Tierens, Zhang and Myra2020b), and are related to radial and poloidal harmonics in cylindrical geometry. To the best of the authors’ knowledge, these analytic calculations are intractable for Gaussian (and more general) filament cross-sections.

Figure 7. Plots of (a) scattering width $\sigma _{S}$ and (b) asymmetric scattering metric $\alpha$ for an incident slow wave at 4.6 GHz, $N_{||}=2$, ${B}=4\ \text {T}$ and $a_b =0.30\ \text {cm}$. (c,d) Plot of $\sigma _{S}$ and $\alpha$, respectively, for $a_b=0.85\ \text {cm}$. (e,f) Plot of $\sigma _{S}$ and $\alpha$, respectively, for $a_b=1.30\ \text {cm}$. The white region in the upper-right corner of each subplot has no data plotted.

7.2. Asymmetric scattering

The quantity $\alpha$ is introduced as a metric for asymmetric scattering:

(7.1)\begin{equation} \alpha = \int_{0}^{\rm \pi} \hat{\sigma}_{S}(\theta) \,\textrm{d} \theta - \tfrac{1}{2}. \end{equation}

Here $\alpha =-0.5,+0.5$ denote that power is only scattered downwards ($-{\rm \pi} < \theta < 0$), upwards ($0 < \theta < {\rm \pi}$). If $\alpha =0$, then an equal fraction of power is scattered downwards and upwards. Figure 7(b,d,f) reveals that $\alpha$ is not guaranteed to be zero. This means, in general, $\sigma (\theta )$ is not an even function and the scattered power is not equally distributed downwards or upwards. This effect is most noticeable at high background densities ($n_0 \gtrsim 2\times 10^{19}\ \text {m}^{-3}$). For positive density modifications ($n_b/n_0 > 1$), power is predominantly scattered upwards. The reverse is true for negative density modifications ($n_b/n_0 < 1$). In the typical tokamak SOL, filaments are predominantly denser than the background plasma (Graves et al. Reference Graves, Horacek, Pitts and Hopcraft2005). As a result, the incident wave is preferentially scattered upwards. The strength of this asymmetry depends on the statistical properties of the filaments.

Asymmetric scattering is possible in an anisotropic medium. Specifically, the dielectric dyadic tensor, $\varepsilon$, has off-diagonal components $\varepsilon _{12} = -\varepsilon _{21} \neq 0$, which permit asymmetric scattering (Wu Reference Wu1994), so $\varepsilon _{12}=-\textrm {i} \varepsilon _{xy}$. In the LH limit, $\varepsilon _{xy}\approx \omega _{pe}^{2}/{\omega \varOmega _{ce}}$, where $\varOmega _{ce} \equiv eB/m_e$ is the electron cyclotron frequency. It is clear that the sign of $\varepsilon _{xy}$ is dependent on the sign of $B$ (that is, whether the magnetic field is oriented co-parallel or counter-parallel with $\hat {e}_{z}$). Correspondingly, when the direction of $\boldsymbol {B}$ is flipped, $\boldsymbol{\varepsilon} \rightarrow \boldsymbol{\varepsilon}^{T}$ and $\sigma (\theta ) \rightarrow \sigma (-\theta )$.

Notably, this asymmetric scattering effect is not accounted for in previous treatments for LH wave scatter. It is easy to see why this is the case for models that assume drift-wave-like turbulence. These models assume that density fluctuations are equally likely to be above or below the background density. Therefore, this asymmetric scatter effect is statistically cancelled out.

There also exist LH scattering models that assume coherent turbulent structures that can, on average, be denser than the background (Hizanidis et al. Reference Hizanidis, Ram, Kominis and Tsironis2010; Biswas et al. Reference Biswas, Baek, Bonoli, Shiraiwa, Wallace and White2020). These models also do not account for asymmetric scattering, because they make the ray tracing approximation.

The reason why ray tracing cannot model asymmetric scattering is subtle. It is related to the breakdown of the ray tracing approximation. Consider the ray tracing equations for an LH ray initially propagating with $\boldsymbol {N}_{\perp }$ aligned along the $x$-direction and background $\boldsymbol {B}$ aligned along the $z$-direction. The ray tracing equations involve partial derivatives of $\text {det}(\boldsymbol {D})$, where $\boldsymbol {D}$ is now the dielectric tensor and

(7.2)\begin{equation} \boldsymbol{D} \boldsymbol{\cdot} \tilde{\boldsymbol{E}}= \begin{bmatrix} \varepsilon_{{\perp}} - N_{||}^{2} & -\textrm{i}\varepsilon_{xy} & N_{{\perp}}N_{||} \\ i\varepsilon_{xy} & \varepsilon_{{\perp}} - N^{2} & 0 \\ N_{{\perp}}N_{||} & 0 & \varepsilon_{||} - N_{{\perp}}^{2} \end{bmatrix} \boldsymbol{\cdot} \begin{bmatrix} \widetilde{E_x} \\ \widetilde{E_y} \\ \widetilde{E_z} \end{bmatrix} = 0,\end{equation}

where $\tilde {\boldsymbol {E}}(\boldsymbol {r})$ is the slowly varying part of the electric field. Here, det($\boldsymbol {D}$) has terms that are quadratic in $\varepsilon _{xy}$, but no linear $\varepsilon _{xy}$ terms. As a result, information about the sign of $B$ along $\hat {e}_{z}$ is lost. Compare this to the full EM wave equation, with no ray tracing approximation, written in a form similar to that of (7.2).

(7.3)\begin{equation} \begin{bmatrix} \varepsilon_{{\perp}} - \dfrac{c^{2}}{\omega^{2}}(F_{yy} + F_{zz}) & \dfrac{c^{2}}{\omega^{2}}F_{yx} - \textrm{i}\varepsilon_{xy} & \dfrac{c^{2}}{\omega^{2}}F_{zx}\\ \dfrac{c^{2}}{\omega^{2}}F_{xy} + \textrm{i}\varepsilon_{xy} & \varepsilon_{{\perp}} - \dfrac{c^{2}}{\omega^{2}}(F_{zz} + F_{xx}) & \dfrac{c^{2}}{\omega^{2}}F_{zy}\\ \dfrac{c^{2}}{\omega^{2}}F_{xz} & \dfrac{c^{2}}{\omega^{2}}F_{yz} & \varepsilon_{||} - \dfrac{c^{2}}{\omega^{2}}(F_{xx} + F_{yy}) \end{bmatrix} \boldsymbol{\cdot} \begin{bmatrix} \widetilde{E_x} \\ \widetilde{E_y} \\ \widetilde{E_z} \end{bmatrix} = 0,\end{equation}

where

(7.4a)$$\begin{gather} F_j = k_j - \textrm{i} \frac{\partial}{\partial j} , \end{gather}$$

(7.4b)$$\begin{gather}F_{jl} = F_j(F_l). \end{gather}$$

Equation (7.3) accounts for $\boldsymbol {\nabla } \boldsymbol {k}$ and $\boldsymbol {\nabla } \tilde {\boldsymbol {E}}$ terms, which are usually neglected in ray tracing because the plasma is assumed to be sufficiently homogeneous, such that $k_{\perp }L_n \gg 1$, where $L_n\equiv |{\boldsymbol {\nabla } n}/{n_0}|^{-1}$ is the characteristic length of the density inhomogeneity. This heuristic validity criterion is actually too lax for magnetized plasma. A perturbation analysis of (7.3) reveals that these higher-order gradient terms can be comparable to the zeroth-order terms (7.2) even if $k_{\perp }L_n \gg 1$. Following some algebra, it is found that the two leading higher-order terms are linear in $\varepsilon _{xy}$ and quadratic in $N_{||}$, such that the actual validity criterion for ray tracing is $(|\varepsilon _{xy}| + N_{||}^{2}) {1}/{k_{\perp } L_{n}} \ll 1$. (More details about this perturbation analysis can be found in Appendix A of Biswas et al. (Reference Biswas, Baek, Bonoli, Shiraiwa, Wallace and White2020), although in that derivation, $\boldsymbol {\nabla }\boldsymbol {k}$ terms were erroneously neglected, which led to the dropping of the $N_{||}^{2}$ term in the final ray tracing validity criterion.) At initial launching of the LH wave $N_{||}^{2} , |\varepsilon _{xy}| \sim 1$, but they can both grow to be much larger as the ray continues to propagate. Specifically, $|\varepsilon _{xy}|\propto n_e$, and so it rapidly increases as the ray propagates into the plasma. The restriction that $\varepsilon _{xy}$ places on LH ray tracing has been commented on before (Ott Reference Ott1979). The following discussion is the first time it has been linked to asymmetric scattering in the context of LH waves.

In deriving the new ray tracing criterion for LH waves in a magnetized plasma, it was revealed that one of the leading higher-order gradient terms neglected in ray tracing is linear in $\varepsilon _{xy}$. This is precisely the term with information about the orientation of $\boldsymbol {B}$. In accordance, as $|\varepsilon _{xy}|{1}/{k_{\perp } L_n}$ grows and becomes comparable to unity, asymmetric scattering also becomes important. This is shown numerically by simulating the scatter of LH waves from four increasingly dense Gaussian filaments. Figure 8 plots the validity regime of ray tracing in the presence of a Gaussian filament as a function of $n_b/n_0$ and $a_b$. The incident wave is launched at 4.6 GHz with $N_{||}=2$ and ${B} = 4\ \text {T}$. It is assumed that $L_{n}^{-1} \approx ({n_b}/{n_0}-1)/a_b$. The black line denotes the validity limit for $n_0 = 1 \times 10^{19}\ \text {m}^{-3}$. To the right of this line, $|\varepsilon _{xy}|{1}/{k_{\perp } L_n} > 1$ and to the left $|\varepsilon _{xy}|{1}/{k_{\perp } L_n} < 1$. (For simplicity, the $N_{||}^{2}$ term is ignored.) The four starred points denote filaments with $a_b = 1$ cm and $n_b/n_0 = [1.24, 1.6, 2.44, 4.6]$ (plotted left to right). Alternatively, these filaments satisfy $n_b = [0.1, 0.25, 0.6, 1.5]\times n_{b,\text {max}}$, where $n_{b,\text {max}}$ satisfies $|\varepsilon _{xy}|{1}/{k_{\perp } L_n} = 1$. Qualitatively, the green point signifies a filament that is validly treated with ray tracing because $|\varepsilon _{xy}|{1}/{k_{\perp } L_n} \approx {n_b}/{n_{b,\text {max}}} = 0.1 < 1$. The yellow points are marginally valid, because ${n_b}/{n_{b,\text {max}}} < 1$ but also $\mathscr {O}(1)$. The red point, for which ${n_b}/{n_{b,\text {max}}} > 1$, certainly cannot be treated using ray tracing.

Figure 8. Validity of ray tracing for Gaussian filaments as a function of relative density ($n_b/n_0$) and radial width ($a_b$). Black line denotes validity limit at $n_0 = 1 \times 10^{19}\ \text {m}^{-3}$ (see § 7.2 for more details). Stars denote filaments that are used in scattering studies in figures 9 and 10. Numbers in colour denote the ratio $n_b/n_{b,\text {max}}$.

Figure 9 plots the ray trajectories of LH rays incident from the left and interacting with a filament. As $n_b/n_0$ increases, rays are more strongly refracted, which results in a shadowing effect downstream of the filament. Notably, these ray trajectories are always perfectly symmetric with respect to $y=0$.

Figure 9. Ray tracing simulations of LH waves scattering from a filament. Here, $f=4.6\ \textrm {GHz}$, $N_{||}=2$, ${B}=4\ \text {T}$ and $n_0 = 1 \times 10^{19}\ \text {m}^{-3}$. Each subplot corresponds to a star on figure 8. Blue lines denote ray trajectories.

In contrast, figure 10 plots $\sigma (\theta )$ calculated using the SAS method for the same four cases simulated in figure 9. The SAS method (like all full-wave treatments) implicitly accounts for all terms in (7.3). For $n_b/n_0=1.24$ and 1.6, $\sigma (\theta )$ is symmetric about $\theta = 0$. The profiles peak at $\theta =0$, signifying predominantly forward scatter. At $n_b/n_0=2.44$, $\sigma (\theta )$ is slightly asymmetric because the side lobe at $-40^{\circ }$ is larger than the one at $+40^{\circ }$. At $n_b/n_0=4.6$, $\sigma (\theta )$ is clearly asymmetric. Notably, the largest lobe is centred at $+5^{\circ }$. While a direct quantitative comparison between figures 9 and 10 is not possible, it is clear that as $|\varepsilon _{xy}|{1}/{k_{\perp } L}$ grows (filaments get denser), ray tracing becomes less accurate because the increasingly important asymmetric scattering effect is ignored. It is important to note that filaments with $n_b/n_0 \gtrsim 2$ are common in the SOL (Zweben et al. Reference Zweben, Stotler, Terry, LaBombard, Greenwald, Muterspaugh, Pitcher, Group, Hallatschek and Maqueda2002), which signifies ray tracing is inadequate for the treatment of LH wave scattering in realistic SOL turbulence.

Figure 10. Polar plots of differential scattering width $\sigma (\theta )$ calculated using the SAS method. Simulation parameters are the same as in figure 9.

7.3. Modified wave spectrum in front of LH antenna

The incident wave parameters, background plasma parameters and joint-p.d.f. of filament parameters determine $\sigma _{\text {eff}}(\theta )$. The MC method is used to compute the modified wave spectrum after the incident wave interacts with a slab layer of thickness $L_{x}$ and packing fraction $f_p$. Figure 11 plots the modified wave spectra resulting from the $\sigma _{\text {eff}}(\theta )$ shown in figure 6. Here, $L_{x} = 2.5\ \textrm {cm}$, which is the typical gap between the LH antenna and the separatrix in C-Mod. Green, black and red lines denote $f_p=[0.1,0.25,0.5]$. Note that the ballistic power fraction is not plotted (if it were, it would be a Dirac delta plotted at $\theta = 0$). In the low-density ($n_0 = 0.55 \times 10^{19}\ \text {m}^{-3}$) case, the modified wave spectrum is smoothly broadened in $\theta$-space, with a peak centred at $\theta = 0$. Increase in $f_p$ leads to a decrease in ballistic power and increase in reflected power. This is expected, because $\varSigma _{\text {eff}}$, the inverse mean-free-path to scatter, is linearly proportional to $f_p$. In the high-density ($n_0 = 4.8 \times 10^{19}\ \text {m}^{-3}$) case, the modified wave spectrum is significantly asymmetric, with net power scattered in the $+\theta$-direction. Naturally, this is the result of $\sigma _{\text {eff}}(\theta )$ in the high-density case being very asymmetric. Again, ballistic power decreases and reflected power increases with $f_p$. In comparing the low- and high-density cases, it is found that the high-density case results in significantly greater reflected power. This arises from $\sigma _{\text {eff}}$ being larger in the high-density cases, as can be seen from figure 6.

Figure 11. Modified LH wave spectrum after interacting with turbulent slab. Here, $|\theta | < {\rm \pi}/2$ denotes transmitted power and $|\theta | > {\rm \pi}/2$ denotes reflected power. Ballistic power is not plotted and $f_p =$ 0.1 (green), 0.25 (black), 0.5 (red). Left and right plots assume low and high background density, respectively. The $F_{\text {bal}}$, $F_{{T}}$ and $F_{{R}}$ denote fractional ballistic, transmitted and reflected power, respectively. The $\sigma _{\text {eff}}(\theta )$ used is shown in figure 6.

7.4. SAS-MC compared with ray tracing

A comparison study between the SAS-MC model and ray tracing model is conducted. Rays are launched in a slab geometry and are incident normal to a slab composed of randomly generated filaments (see § 6.2). A ray terminates when it leaves the slab (either reflected backward or transmitted forward), at which point the angle between the ray's perpendicular group velocity ($\boldsymbol {v}_{gr_{\perp }}$) and $\hat {e}_x$ is tallied. This poloidal angle is the direction that the ray continues to propagate and radiate power away from the slab. It is therefore equivalent to $\theta$ in the SAS-MC model. Following multiple ray launches, a histogram of these tallies is constructed. This histogram, once properly normalized, is equivalent to a modified wave spectrum that can be compared with the wave spectrum computed with the SAS-MC model.

Figure 12 plots modified wave spectra computed using the SAS-MC model and the ray tracing model for statistically identical turbulent slabs. Three different cases are run. All cases assume LH rays incident at 4.6 GHz and $N_{||}=2$. Here, ${B}=4\ \text {T}$ and $L_x=2.5\ \textrm {cm}$. In the first case, $n_0 = 1 \times 10^{19}\ \text {m}^{-3}$, and the joint-p.d.f. in figure 4 is assumed. The SAS-MC model and ray tracing model show good agreement in the wave spectra. Both predict low reflected power fractions. The wave spectrum computed with the SAS-MC model is fairly symmetric. This means that asymmetric scatter was weak, and therefore the ray tracing approximation was valid. Thus, the good agreement between the two models. In the next case, $\langle a_b \rangle$ is halved to 0.5 cm. The SAS-MC model results in an asymmetric wave spectrum, such that the peak is shifted to $+0.2$ rad. The onset of significant asymmetric scattering is caused by the decrease in $L_{n}$. At the same time, the ray tracing approximation begins to break down. As a result, the ray tracing results (which are symmetric) begin to deviate from the SAS-MC results. Notably, the ray tracing model severely under-predicts the fraction of reflected power compared with the SAS-MC model. The last case increases the background density to $n_0 = 4.8 \times 10^{19}\ \text {m}^{-3}$. Now, $\varepsilon _{xy}$ is large enough that asymmetric scattering is quite strong and the ray tracing approximation is surely invalid. As a result, the two models result in very different wave spectra. Ray tracing predicts a scattered wave spectrum with a large central peak at $\pm 0.2\ \textrm {rad}$. In contrast, the SAS-MC model predicts a smaller central peak slightly shifted in the $+\theta$-direction. The tail to the right of this peak is significantly larger than the one on the left. Lastly, the SAS-MC model results in $\sim$50 % more power being reflected than the ray tracing model.

Figure 12. Comparison between the SAS-MC and ray tracing slab models. Blue bars denote a histogram of ray $\boldsymbol {v}_{gr_{\perp }}$ angles after leaving slab. Red line denotes a modified wave spectrum $P_{sct} (\theta )$ calculated using the SAS-MC method. Here, $L_x = 2.5\ \textrm {cm}$ and $f_p = 0.25$. Filament joint-p.d.f. parameters are the same as in figure 4 unless otherwise noted in the subplot title. The $F_{\text {ref}}$ denotes fractional power reflected in the ray tracing (blue) and SAS-MC (red) models. Ballistic power and unscattered rays are not plotted.