Quantitative comparisons of electron-scale turbulence measurements in NSTX via synthetic diagnostics for high-k scattering

Plasma turbulence gives rise to anomalous transport of particles and heat in magnetic confinement fusion devices [1], which is detrimental to confinement. The complex, kinetic nature of the turbulence has led to the development of sophisticated gyrokinetic models implemented in state-of-the-art numerical simulations [2] (non-linear gyrokinetic simulation) to study the turbulence and consequent turbulence-driven transport. The gyrokinetic model requires extensive validation in today's fusion experiments before achieving a predictive capability for future fusion devices such as ITER [3], FNSF [4] and beyond. Confidence in the predictions from gyrokinetic simulations can be gained via a thorough validation process, which should include detailed comparisons of turbulence characteristics in addition to the traditional comparisons of turbulent fluxes [5, 6]. In this article we make direct comparisons of density fluctuation spectra between experimental turbulence measurements by high-k scattering and non-linear gyrokinetic simulations, which are part of an extensive validation study of electron thermal transport in NSTX [7].

Turbulence in the tokamak core can live at ion gyro radius scales ($k_\perp \rho_s \lt 1$) and at electron gyro radius scales (${k_ \bot }{\rho _s} > 1$), where $k_\perp$ is the perpendicular wavenumber of the turbulence and $\rho_s = c_s/\Omega_{c,i}$ is the ion sound gyro radius evaluated at the ion sound speed $c_s = \sqrt{T_e/m_i}$, and $\Omega_{c,i} = eB/m_i$ is the ion gyro-frequency. Experimentally, coherent scattering techniques implemented through Doppler backscattering (DBS) [8–16] and high-k scattering [17–33] can probe turbulent density fluctuations at electron-scales $({k_ \bot }{\rho _s} > 1)$. These can be of special importance in spherical tokamak plasmas where electron thermal transport dominates heat transport losses in H-mode scenarios [34–39], and electron-scale turbulence can be the main heat loss mechanism in some operating regimes [40, 41]. Previous work to validate non-linear electron-scale gyrokinetic simulations against experimental turbulence measurements has focused on establishing qualitative turbulence comparisons, with more recent efforts expanding to quantitative comparisons through synthetic diagnostic development [42–46]. However, the challenge remains to achieve direct quantitative comparisons of the high-k frequency and wavenumber turbulent fluctuation spectrum between experiment and simulation.

Coherent scattering diagnostics are sensitive to a specific wavenumber $\vec{k}_+$ of the turbulent fluctuations. As a result, previous synthetic diagnostic work has been formulated in wavenumber space (k-space), via the selection of $\vec{k}_+$ by use of a filter in wavenumber space [42–46]. However, the scattering signal for coherent scattering is fundamentally calculated from first principles via integration of the electron density fluctuation amplitude in real space [47, 48]. In this article we build on past work to show how the wavenumber formulation can be naturally derived from real space. We propose two equivalent formulations, in real space and in k-space, for the computation of the scattering signal from coherent scattering in realistic field-aligned coordinates. The quantitative agreement shown between the real space and k-space based synthetic diagnostics provides improved confidence on the validity of the computed synthetic spectra shown in this work.

One of the difficulties in developing synthetic diagnostics for coherent scattering is the different wavenumber definitions employed in experiments and gyrokinetic codes. Experiments generally use cylindrical or Cartesian coordinates for the components of the measured wavenumber $\vec{k}_+$, which are provided by ray-tracing or equivalent calculations. Gyrokinetic codes operate in field-aligned coordinates and use internal wavenumber definitions. An 'apples-to-apples' comparison between experimental turbulence measurements and simulated turbulence requires a mapping of the measured wavenumber to field-aligned coordinates, implemented as part of this work. The wavenumber mapping is also an important step in the development of synthetic diagnostics in wavenumber space. This complexity is absent in the real space formulation, but is necessary to understand the measurement wavenumber range of the diagnostic in the density fluctuation spectra. This motivates the implementation of two equivalent formulations of the synthetic diagnostic, in real space and in k-space.

The rest of this article proceeds as follows. In section 2 we discuss some theoretical considerations of coherent scattering measurements of turbulence fluctuations. In section 3 we outline the implementation of synthetic diagnostics for coherent scattering systems in realistic field-aligned tokamak geometry and introduce the wavenumber mapping. In this section we highlight the importance of geometric effects, such as the normalizing B-field, the effect of plasma elongation and Shafranov shift, which can strongly affect the interpretation of the measured wavenumber components from scattering measurements (by up to factors of ∼5 for the present NSTX case). In section 4 we apply the synthetic diagnostic to compute numerically generated synthetic spectra for the high-k scattering system in NSTX [26–28], using realistic electron-scale gyrokinetic simulations based on a modest-β NSTX NBI heated H-mode plasma. Finally, we compare synthetic high-k frequency and wavenumber spectra with experimental spectrum measurements. The main outcome of this work is the successful validation of electron-scale gyrokinetic simulations in the core-gradient region of a modest-β NSTX H-mode plasma via direct comparison with measured high-k density fluctuation spectra.

Coherent scattering from turbulence fluctuations inherently takes place in a confined region known as the scattering volume V_s, which is generally delimited by the size of the electromagnetic wave beam input in the plasma and by the magnetic field geometry. This leads one to interpret the scattering process as the integration of fluctuations in real space within the scattering volume. However, scattering measurements are usually interpreted in wavenumber space, based on the measurement of a specific turbulence wavenumber $\vec{k}_+$, which is determined by the launching and receiving geometries of the electromagnetic wave beam. This leads one to interpret the scattering process as a selection of a specific wavenumber $\vec{k}_+$ from the density fluctuations. In this section we show how the scattering process can be interpreted equally in real space as well as in k-space.

In the coherent scattering process, plasma electrons are exposed to an external source of electromagnetic radiation (e.g. a laser or microwave source, considered here to be a beam of radius a₀). Accelerated by the incoming electric field, electrons radiate electromagnetic energy in the form of a radiation field. The expression for the scattered power per unit frequency and unit solid angle $\frac{{\mathrm{d}}^2P_s}{{\mathrm{d}}\omega {\mathrm{d}}\Omega}$ is related to the frequency signal of density fluctuations $\delta{n}_{\textrm{u}} (\vec{k}_+, \omega)$ by the textbook formula (appendix B)

$$\begin{equation}\frac{{{{\rm{d}}^2}{P_s}}}{{{\rm{d}}\omega {{d\Omega }}}} = \frac{{{P_0}}}{{{A_i}}}r_0^2|\hat s \times (\hat s \times \hat e){|^2}\frac{1}{{2\pi T}}\left\langle {|\delta {n_{\rm{u}}}({{\vec k}_ + },\omega ){|^2}} \right\rangle \end{equation} \tag{ 1 }$$

where the subscript 'u' indicates that the density fluctuation signal has been properly filtered in the scattering process by a filter U in real space. In expression (1) dΩ is the solid angle, P₀ is the incident beam power in watts, A_i is the incident beam area $A_i = \pi a_0^2$, $r_0 = \frac{e^2}{m_ec^2}$ is the classical electron radius, T is the collection time, $\hat{s}$ is the direction of scattering, $\hat{e}$ is the direction of the scattered electric field and $\left\langle . \right\rangle$ denotes an ensemble average. The incident radiation oscillating at a frequency ω_i with wave-vector $\vec{k}_i$ can be related to the scattered frequency ω_s and wavenumber $\vec{k}_s$ by the scattering matching conditions $\omega = \omega_s - \omega_i$ and $\vec{k}_+ = \vec{k}_s - \vec{k}_i$ (Bragg condition), where ω and $\vec{k}_+$ are the matching turbulence frequency and wavenumber (the contributions from the matching wave-vector $\vec{k}_- = \vec{k}_s+\vec{k}_i$ are negligible and are ignored in this work).

In the case of DBS or Doppler reflectometry [8–10], ray-tracing or beam-tracing methods break down near the cut-off encountered by the incident electromagnetic wave, and a full-wave treatment might be necessary to accurately model the propagation of the electromagnetic wave in the plasma. Despite this, much of the work presented here could still prove useful to help interpret density fluctuation measurements using DBS in some operating regimes (e.g. the linear response regime). In the case of high-k scattering, the incident frequency is typically much higher than any other frequencies in the plasma (${\omega _i},{\omega _s} > > {\omega _{pe}},{\omega _{ce}}, \ldots$). The electromagnetic wave propagates above any cut-off and resonance in the plasma, validating the use of the ray-tracing methods employed in this work.

The quantity $\delta{n}_{\textrm{u}} (\vec{k}_+, \omega)$ in equation (1) has been Fourier decomposed from $\delta{n}_{\textrm{u}} (\vec{k}_+, t)$, which is the synthetic time signal of electron density fluctuations for the selected turbulence wavenumber $\vec{k}_+$. $\delta{n}_{\textrm{u}} (\vec{k}_+, t)$ can be formally computed in real space as well as in wavenumber space:

$$\begin{equation} \begin{aligned} \delta{n}_{\textrm{u}} (\vec{k}_+, t) &= \int_{V_s} {\mathrm{d}}^3\vec{r'}\ \delta{n}(\vec{r'}, t) U(\vec{r'}) {e^{-i\vec{k}_+\cdot\vec{r'}} }\\ &= \frac{1}{(2\pi)^3} \int{{\mathrm{d}}^3\vec{k}\ \delta{n}(\vec{k}, t)} W(\vec{k}-\vec{k}_+). \end{aligned} \end{equation} \tag{ 2 }$$

The real space filter U determines the shape of the scattering volume V_s from the incident and scattered beam profiles and the magnetic field geometry. The quantity $\delta n(\vec{r'}, t)$ is the real electron density fluctuation field, and $\delta n(\vec{k}, t)$ is the raw electron density fluctuation spectrum, computed from $\delta n(\vec{r'}, t)$ by Fourier analysis. W is the scattering filter in wavenumber space corresponding to the weights for each wavenumber $\vec{k}$, and is directly related to the Fourier transform of the scattering volume shape U (appendix B).

Equation (2) states the equivalence between the computation of the synthetic signal of fluctuations in real space and in wavenumber space. In wavenumber space, the synthetic signal is a sum over all turbulence wavenumber contributions around the detected wavenumber $\vec{k}_+$, where a filter $W(\vec{k}-\vec{k}_+)$ is applied to the wavenumber spectrum of fluctuations $\delta n(\vec{k}, t)$. W peaks around the measurement wave-vector $\vec{k}_+$, and down-selects a range of wavenumbers neighboring $\vec{k}_+$ within the range $\Delta \vec{k}^3 \sim 1/V_s$ (V_s is the scattering volume extent, m³). In real space, the synthetic signal can be interpreted as the Fourier component $\vec{k}_+$ of the real quantity $\delta{n}(\vec{r'}, t) U(\vec{r'})$. As a result of U in the real integration, the scattering signal has not only contributions from one lone $\vec{k}_+$(obtained for U = 1), but also from an array of wavenumbers around $\vec{k}_+$ in the range $\Delta \vec{k}^3 \sim 1/V_s$ ($\Delta \vec{k}^3$ is the spectral width of the filter $W(\vec{k}-\vec{k}_+)$). Full information about the detected turbulence wavenumber $\vec{k}_+$ and the spectral width $\Delta \vec{k}^3$ is preserved in the computation of $\delta{n}_{\textrm{u}} (\vec{k}_+, t)$ according to both formulations, and motivates their implementation for realistic tokamak scattering experiments.

In this section we implement synthetic diagnostics for coherent scattering turbulence measurements in the toroidal geometry characteristic of tokamak scattering experiments, both in real space and in wavenumber space. The expression of the synthetic signal δn_u in axisymmetric, toroidal geometry is provided in section 3.1. In section 3.2 we give a brief outline of the derivation and formulation of the synthetic diagnostic signal. Only a succinct derivation is presented, highlighting the most important points. The reader is referred to appendix E for additional details. In section 3.3 we give a specific example corresponding to the experimentally relevant case of scattering at the outboard mid-plane and in the 2D approximation. Three main geometric effects will prove to be crucial for accurate 'apples-to-apples' comparisons between the experiment and simulation of the measured wavenumber $\vec{k}_+$: the normalizing magnetic field entering the definition of the sound gyroradius ρ_s, the Shafranov shift Δ and the flux-surface elongation κ. Not taking into account these effects could lead to systematic errors in the interpretation of the measured wavenumber components, up to a factor of 5 in the present NSTX case. These might be particularly important in the high-β, strongly shaped geometries characteristic of high performance tokamak scenarios, and particularly in spherical tokamaks.

The choices of the magnetic geometry parametrization, the field-aligned wavenumber coordinates, the scattering volume shape, etc, might all differ depending on specific experiments and subsequent modelling tools. The goal of this section is not to be general, but to provide guidelines that one might follow for developing synthetic diagnostics based on coherent scattering from turbulent fluctuations.

3.1. Formulation in real-space versus k-space

Before proceeding to formulate a synthetic diagnostic, one needs information about the scattering location $\vec{r}_0$, the measurement wavenumber $\vec{k}_+$, as well as the scattering volume extent V_s. This information can generally be provided by ray-tracing or beam-tracing calculations. Full wave calculations might be needed close to cut-offs and/or resonances, however, these are not relevant in the context of high-k scattering and are omitted in this work.

We start by writing the synthetic signal of density fluctuations $\delta{n}_{\textrm{u}} (\vec{k}_+, t)$ in real space cylindrical coordinates (R, Z, ϕ):

$$\begin{equation} \begin{aligned} &\delta{n}_{\textrm{u}} (\vec{k}_+, t) = \int_{V_s} {\mathrm{d}}^3\vec{r'}\ \delta{n}(\vec{r'}, t) U(\vec{r'}) {e^{- i\vec{k}_+\cdot\vec{r'}} }\\ &=\int{{\mathrm{d}}\varphi } {R{\mathrm{d}}R{\mathrm{d}}Z \delta{n}(R, Z, \varphi, t) U(R - R_0, Z - Z_0, \varphi - \varphi_0) e^{- i\vec{k}_+\cdot\vec{r} } } \end{aligned} \end{equation} \tag{ 3 }$$

where R is the major radial direction, Z is the vertical direction and ϕ is the toroidal direction. The filter in real space U is centered around the scattering location $\vec{r}_0 = (R_0, Z_0, \varphi_0)$, and the product $\vec{k}_+\cdot\vec{r}$ needs to be written in cylindrical coordinates (appendix E). In the context of magnetized plasma turbulence, fluctuations perpendicular to the magnetic field are expressed in wavenumber space components that depend on the field-aligned coordinates, and are routinely employed by gyrokinetic codes. The formulation of the synthetic diagnostic in k-space field-aligned coordinates is more cumbersome than in real space, but yields a direct map of the measured wavenumber in the density fluctuation wavenumber power spectrum. This is useful for correctly interpreting the measurement range of current coherent scattering measurements, as well as for the projection of future measurements. Figure 1 schematically shows the synthetic diagnostic procedure in real space (figure 1(a)) versus k-space (figure 1(b)) in the 2D approximation for a realistic NSTX H-mode discharge.

Zoom In Zoom Out Reset image size

To compute the synthetic signal of density fluctuations δn_u we follow the expansion of fields as implemented in the gyrokinetic codes GYRO/CGYRO (appendix E, [49–51]). In the field-aligned coordinates (r, θ, ϕ), the electron density field δn(r, θ, ϕ, t) is expanded as a function of the toroidal and radial mode number components (n, p) as shown by equation (4) (note GYRO internally computes δn_n(r, θ, t) in real space while CGYRO computes δn_np(θ, t) spectrally). The real density field δn(r, θ, ϕ, t) can be substituted into equation (3) to compute the synthetic signal of density fluctuations $\delta{n}_{\textrm{u}} (\vec{k}_+, t)$, leading to

$$\begin{equation} \begin{aligned} &\delta{n}_{\textrm{u}} (\vec{k}_+, t) \approx \sum_{n, p} {U_{np} \delta n_{np}(\theta_0, t) e^{in\overline{\omega}_0t}} \\ \text{where } &\delta n_{np} \text{is expanded in GYRO/CGYRO as} \\ &\delta n(r,\theta, \varphi, t) = \sum_{n,p}{\delta n_{np}(\theta,t)e^{-in\alpha}e^{i2\pi pr/L_r}e^{in\overline{\omega}_0t}}. \end{aligned} \end{equation} \tag{ 4 }$$

Here $\delta n_{np}(\theta_0, t)$ are the (n, p) components of the real electron density field at the poloidal location θ₀ of scattering, $\overline{\omega}_0$ is the toroidal rotation frequency of the background plasma (producing the Doppler shift) and α is the field line label (appendix E and [49, 50]). U_np is the scattering matrix, defined as

$$\begin{equation} \begin{aligned} &U_{np} = \int {{\mathrm{d}}^3\vec{r} U(\vec{r})e^{-in\alpha}e^{i2\pi p r/L_r} e^{{-i\vec{k}_+}\cdot\vec{r} } } \end{aligned} \end{equation} \tag{ 5 }$$

where L_r is the radial box-size of the simulation (used to define p). The scattering matrix is a filter in (n, p), and is the representation of the k-space filter $W(\vec{k} - \vec{k}_+)$ from equation (2) when expressed in (n, p) mode numbers. The scattering matrix U_np peaks around specific toroidal and radial mode numbers, which can be calculated and mapped from a turbulence wavenumber $\vec{k}_+$ via a mapping in wavenumber space (equation (8)). Similar to the relation between U and W from equation (2), U_np can be interpreted as a 'Fourier-like' transform of the scattering volume shape U, and can be analytically computed for simple cases such as at the outboard midplane and a separable U(R, Z, ϕ), as we will show in the next section. An example of a specific shape of U_np is shown in figure G4 (appendix G). Equations (3) and (4) are the equivalent formulations of the synthetic signal in real versus k-space in toroidal geometry. Specifics about the computation of the synthetic signal in toroidal, field-aligned geometry are shown in the next section.

3.2. Computation of the synthetic signal $\delta n_\textrm{u}$

In this section we give a brief outline of the procedure to compute the synthetic signal of density fluctuations δn_u. We restrain ourselves to Gaussian scattering volume shapes U for simplicity in section 3.2.1. In section 3.2.2 we introduce the wavenumber mapping from Cartesian coordinates to field-aligned coordinates which are routinely used by gyrokinetic codes. In section 3.2.3 we give analytical expressions for the filters to be applied in k-space field-aligned coordinates in the full 3D formulation and for arbitrary poloidal locations of scattering.

3.2.1. Scattering volume shape

The specific shape of the scattering volume U entering in the computation of the scattering matrix U_np can vary depending on the specific scattering geometry of particular scattering experiments. Here we assume the scattering volume envelope is separable in filter functions $\Psi_R, \Psi_Z, \Psi_\varphi$, and is characterized by a Gaussian shape centered around $(R_0, Z_0, \varphi_0 = 0)$

$$\begin{equation} \begin{aligned} & U(R, Z, \varphi) && = \Psi_R(R-R_0)\Psi_Z(Z-Z_0)\Psi_\varphi(\varphi) \\ & \text{with} \\ & \Psi_R(R-R_0) && = \exp[-(R-R_0)^2/\Delta R^2] \\ & \Psi_Z(Z-Z_0) && = \exp[-(Z-Z_0)^2/\Delta Z^2] \\ & \Psi_\varphi(\varphi) && = \exp[-\varphi^2/\Delta \varphi^2] && 3{\mathrm{D}}\\ & && = \delta(\varphi/\Delta \varphi) && 2{\mathrm{D}}\\ \end{aligned} \end{equation} \tag{ 6 }$$

where we have made the distinction between a 3D and a 2D implementation in the toroidal filter $\Psi_\varphi$ (δ(ϕ/Δϕ) is the Dirac delta function). R = R(r, θ) and Z = Z(r, θ) are specified by the flux surface parametrization. ΔR, ΔZ and Δϕ are the dimensions of the scattering volume shape U along the major radius, vertical and toroidal directions, respectively. In the 2D approximation we neglect any toroidal variation and the fluctuations will be filtered at a fixed toroidal slice. Although we will show the full 3D formulation of the synthetic diagnostic, in the practical example shown in section 4 we will restrict ourselves to 2D as we will justify. More details about the 2D versus 3D approximation can be found in appendices E and F.

At the outboard midplane (θ₀ ≈ 0), the radial and vertical filters Ψ_R and Ψ_Z can be expressed as

$$\begin{equation} \begin{aligned} \exp[- (R {-} R_0)^2/\Delta R^2] \approx \exp[ - (r-r_0)^2/\Delta r^2] \equiv \Psi_r(r-r_0) \\ \exp[- (Z {-} Z_0)^2/\Delta Z^2] \approx \exp[ - (\theta-\theta_0)^2/\Delta \theta^2] \equiv \Psi_\theta(\theta-\theta_0) \\[6pt] \end{aligned} \end{equation} \tag{ 7 }$$

where $\Delta r = |\nabla r|_0 \Delta R$ and Δθ = ΔZ/(r₀κ). $|\nabla r|_0$ is a local gradient related to Shafranov shift via $|\nabla r|_0 \approx 1/(1+\Delta)$ ([49]) and κ is the flux-surface elongation. This allows U to be written as $U(R, Z, \varphi) \approx \Psi_r(r-r_0) \Psi_\theta(\theta-\theta_0) \Psi_\varphi(\varphi)$. For a local flux-tube simulation the radial filter Ψ_r can take the value of 1, with an equivalent radial extent of the scattering volume $\Delta R \approx L_r/2|\nabla r|_0$. R₀ and Z₀ are directly provided by ray-tracing calculations, and r₀ and θ₀ can be computed from R₀ and Z₀ by use of the flux-surface shape parametrization R(r, θ), Z(r, θ). Using the particular shape of the filters in real space $\Psi_r, \Psi_\theta, \Psi_\varphi$ from equation (7), one can compute the corresponding filters in wavenumber space and the scattering matrix U_np (section 3.2.3).

3.2.2. Wavenumber mapping

Before proceeding to compute the equivalent filters in wavenumber space from those in real space, one needs to map the measured wavenumber components from those provided in experiments (typically Cartesian coordinates) to the field-aligned geometry definitions employed by gyrokinetic codes. As shown in greater detail in appendix D, the wavenumber components in Cartesian coordinates $(k_x, k_y, k_z)_+$ are mapped to toroidal and radial mode number components $(n^\theta_+, n^\varphi_+, p_+)$ via the wavenumber mapping

$$\begin{equation} \left\{ \begin{alignedat}{3} & \frac{2\pi}{L_r}p_+ - \frac{\partial \alpha}{\partial r}\bigg|_0 n^{\theta}_+ && = \frac{\partial R}{\partial r}\bigg|_0 k_{x+} + \frac{\partial Z}{\partial r}\bigg|_0 k_{z+} \\ & -\frac{1}{r_0} \frac{\partial \alpha}{\partial \theta}\bigg|_0 n^{\theta}_{+} && = \frac{1}{r_0}\frac{\partial R}{\partial \theta}\bigg|_0 k_{x+} + \frac{1}{r_0}\frac{\partial Z}{\partial \theta}\bigg|_0 k_{z+} \\ & -\frac{1}{R_0}n^{\varphi}_+ && = k_{y+} \\ \end{alignedat} \right. \end{equation} \tag{ 8 }$$

involving r and θ derivatives of the flux surface coordinates R = R(r, θ), Z = Z(r, θ) and the field line label α = α(r, θ) (appendix E and [49, 50]). The Cartesian coordinates of $\vec{k}_+ = (k_x, k_y, k_z)_+$ are defined as: k_x is along the major radius direction (for ϕ₀ = 0), k_y is along the toroidal direction and k_z is along the vertical direction (figure D2 and appendix E). Given a particular measured wavenumber $\vec{k}_+$ provided by ray-tracing or beam-tracing calculations, equation (8) states what are the corresponding mode numbers $n^\theta_+, n^\varphi_+$ and p₊ in a gyrokinetic simulation following the magnetic field line. The mapping given by equation (8) is local, denoted by a subscript 0 indicating local values at $(R_0, Z_0, \varphi_0 = 0)$. A simple explanation of some of these terms entering the wavenumber mapping (8) can be found in the simplified case of 2D and outboard midplane, discussed in section 3.3 and in appendix G.

The mapped mode numbers $n^\theta_+, n^\varphi_+, p_+$ are generally non-integers. We make the distinction between a toroidal mode number $n_+^{\theta}$ associated to the vertical component k_z+, and a different toroidal mode number $n_+^{\varphi}$ associated to the toroidal component k_y+. $n_+^{\theta}$ and $n_+^{\varphi}$ are in principle independent of each other, however, it will be shown in appendix E how the condition for successful scattering $\vec{k}\cdot\vec{B} \approx 0$ restricts them to have a similar value $n^\theta_+ \approx n^\varphi_+$. For reference, in s − α geometry, $n_+^\theta$ would reduce to $n_+^\theta q_0/r_0 = k_{z+}$ at the midplane (θ₀ = 0) while $2\pi p_+/L_r = k_{x+}$, and $n_+^\varphi$ would be ignorable in the 2D approximation. In section 4.2.3 and appendix G examples are given of the application of the wavenumber mapping.

3.2.3. Scattering matrix U_np

Under the assumption of separable scattering volume shape (equation (6)) and at the outboard midplane θ₀ ≈ 0, the scattering matrix U_np can be decomposed as a product of toroidal and radial mode number filters $\Phi_n, \Theta_n, \Pi_p$ as follows (appendix E):

$$\begin{align} & U_{np} = e^{-in\alpha_0}e^{i2\pi p r_0/L_r} \Phi_n(n - n^\varphi_+) \Theta_n(n - n^\theta_+) \Pi_p(p - p_+) \notag\\ &\text{with}\quad \notag\\ &\Phi_n(n - n^\varphi_+) \approx R_0 \bigg( \frac{4\pi }{(\Delta n^\varphi )^2} \bigg)^{1/2} e^{-(n-n^\varphi_+)^2/(\Delta n^\varphi )^2}\qquad \text{3D} \notag\\ & \qquad\qquad\enspace\ \approx R_0\Delta \varphi \qquad\qquad\qquad\qquad\qquad\qquad\quad\ \text{2D} \notag\\ &\Theta_n(n-n^\theta_+) \approx \sqrt{\pi}(\Delta Z) e^{-(n- n^\theta_+)^2/(\Delta n^\theta)^2} \notag\\ &\Pi_p(p-p_{+}) \approx\sqrt{\pi}(\Delta R) e^{-(p- p_{+})^2/\Delta p^2}\qquad \text{Global simulation} \notag\\ & \qquad\qquad\enspace\ \approx \frac{L_r}{|\nabla r|_0} {\mathrm{sinc}}[{\pi(p- p_{+})}]\qquad\ \ \text{Local simulation} . \end{align} \tag{ 9 }$$

The toroidal and radial mode number filters Φ_n, Θ_n and Π_p are, respectively, centered around the mapped mode numbers $(n^\varphi_+, n^\theta_+, p_+)$, which can be calculated using equation (8). The toroidal mode number filter Φ_n takes different expressions in the 2D-approximation versus in the full 3D treatment. The 2D approximation relies on a fixed toroidal slice and U has no toroidal dependence, which translates to an infinitely thin toroidal extent $\Delta \varphi \rightarrow 0$, or equivalently $\Psi_\varphi = \delta(\varphi/\Delta \varphi)$ (equation (6)). The resulting toroidal mode number filter Φ_n is simply constant = R₀Δϕ.

Equation (9) also shows a different radial mode number filter Π_p in a local versus global simulation. Local, flux-tube simulations are characterized by constant background profile gradients along the full radial domain, justifying the radial filter to take the value Ψ_r = 1 (resulting in the sinc function in equation (9)). However a global simulation retains radial profile variation, and the radial filter in real space Ψ_r has to take the shape dictated from experiment. In this article the gyrokinetic code GYRO is run in local, flux-tube mode.

In equation (9), the toroidal, poloidal and radial mode number resolutions take the following values (appendix E):

$$\begin{equation} \begin{alignedat}{3} &(\Delta n^\varphi )^2 && = 4/\Delta \varphi^2 - 2iR_0k_{x+} \\ &\Delta n^\theta && = \frac{2}{\Delta \theta \frac{\partial \alpha}{\partial \theta}|_0} \\ &\Delta p && = \frac{L_r}{\pi|\nabla r|_0\Delta R} \\ \end{alignedat}\label{} \end{equation} \tag{ 10 }$$

where ΔR is the radial extent of the scattering volume in a global simulation, but $\Delta R \approx L_r/2|\nabla r|_0$ for a local simulation. The resolution associated with the toroidal filter Φ_n is complex in nature and depends on the toroidal extent of the scattering volume Δϕ and the x component of the sampled wavenumber $R_0k_{x+}$. The product $(\Delta \varphi)^2R_0k_{x+}$ will determine the importance of 3D effects in the computation of the synthetic signal, as discussed in appendix E.

3.3. 2D and outboard midplane approximation

In this subsection we build on the intuition behind the formulas presented in the previous section in the particular example of scattering at the outboard midplane and in the 2D approximation (neglecting the toroidal variation of the scattering volume). This is motivated by the fact that most scattering turbulence measurements take place at the outboard midplane, in which traditional ballooning drift wave instabilities and consequent microturbulence fluctuations tend to exhibit the highest amplitude. It is useful here to introduce the radial and poloidal wavenumber components of the turbulence $(k_r, k_\theta)$, since their normalizations by ρ_s are the physically meaningful quantities characterizing microturbulence fluctuations. In GYRO [49, 50] and CGYRO [51] these are related to the toroidal and radial mode numbers (n, p) by $k_r = 2\pi p/L_r$ and $k_\theta = nq_0/r_0$, where q₀ is the local safety factor. At the outboard midplane and in the 2D approximation, the corresponding $(k_r\rho_s, k_\theta\rho_s)_+$ values mapped from a $(k_x, k_z)_+$ couple can be simplified from equation (8) to take the following form:

$$\begin{equation} \left\{ \begin{aligned} (k_{r+}\rho_s)^\textrm{sim} &\approx \frac{k_{x+}}{|\nabla r|_0} \rho_s^\textrm{sim} \\ (k_{\theta+}\rho_s )^\textrm{sim} &\approx - \frac{\kappa q}{\frac{\partial \alpha}{\partial \theta}}\bigg|_0 k_{z+} \rho_s^\textrm{sim} \\ \end{aligned} \right. \end{equation} \tag{ 11 }$$

where we employed the Miller flux surface parametrization [56] and κ is the flux-surface elongation. The wavenumber mapping in equation (11) highlights three main geometric effects affecting the mapping: the effect of the normalizing magnetic field entering the definition of ρ_s, the effect of Shafranov shift Δ affecting the radial wavenumber component k_r through $|\nabla r|_0$ and the effect of flux-surface elongation κ affecting the poloidal wavenumber $k_\theta$. In unshifted s − α geometry we have α ≈ ϕ − qθ, Δ = 0 and κ = 1, resulting in $(k_r, k_\theta)_+ = (k_x, k_z)_+$ as expected. However, realistic flux-surface geometries and off-midplane locations can significantly modify the mapping with respect to the s − α midplane approximation. Appendix G presents additional intuition behind these effects.

Within the 2D approximation and at the outboard midplane the scattering matrix U_np can be expressed as a product of separate filter functions $W_{k_r}$ and $W_{k_\theta}$, written now in terms of $(k_r, k_\theta)$

$$\begin{equation} \begin{alignedat}{3} &U_{np} = (R_0 \Delta \varphi) && e^{-in\alpha_0}e^{ir_0k_r} W_{k_r}(k_r-k_{r+}) W_{k_\theta}(k_\theta - k_{\theta+}) \\ &\text{with} \\ &W_{k_\theta}(k_\theta-k_{\theta+}) &&{\approx}\,\sqrt{\pi}(\Delta Z) e^{-(k_\theta- k_{\theta+})^2/\Delta k_\theta^2} \\ &W_{k_r}(k_r-k_{r+}) &&{\approx}\, \sqrt{\pi}(\Delta R) e^{-(k_r- k_{r+})^2/\Delta k_r^2}\quad \text{Global simulation} \\ & &&{\approx}\, \frac{L_r}{|\nabla r|_0} {\mathrm{sinc}}[{(k_r{-} k_{r+})L_r/2}]\quad \!\!\!\!\text{Local simulation}. \\ \end{alignedat} \end{equation} \tag{ 12 }$$

where we recall $(k_r, k_\theta) = (2\pi p/L_r, nq/r)$. We have expressed the toroidal and radial mode number filters $\Theta_n(n-n_+^\theta)$ and $\Pi_p(p-p_+)$ in equation (9) as poloidal and radial wavenumber filters $W_{k_\theta}$ and $W_{k_r}$. Figure 4 shows two examples of radial and poloidal wavenumber filters corresponding to realistic geometry from the high-k scattering diagnostic in NSTX.

It is useful to say a few words about the extent of the scattering volume U and how it might affect the measured wavenumbers. Assume a scattering measurement sensitive to a scattering vector with components $(k_r, k_\theta)_+$ and having a scattering volume with a characteristic length along the major radius ΔR and vertical dimension ΔZ. In the outboard midplane approximation, this will result in a wavenumber resolution $\Delta k_r, \Delta k_\theta$ given by

$$\begin{equation} \left\{ \begin{aligned} \Delta k_{r} \rho_s &\approx \frac{2}{\Delta R |\nabla r|_0} \rho_s \\ \Delta k_{\theta} \rho_s &\approx \frac{2\kappa q}{\Delta Z \frac{\partial \alpha}{\partial \theta}\big|_0 } \rho_s \\ \end{aligned} \right. \end{equation} \tag{ 13 }$$

which corresponds to the resolution of the filters in equation (12). The resolutions are inversely proportional to the scattering volume dimensions, namely $\Delta k_r \propto 1/\Delta R$ and $\Delta k_\theta \propto 1/\Delta Z$. Equation (13) indicates that a wide scattering volume extent will result in spectrally localized measurements in k-space. On the other hand, a narrow scattering volume extent will result in a spatially localized measurement, having contributions from a wide array of wavenumbers. This feature is reminiscent of Heisenberg's uncertainty principle in quantum mechanics.

In this section we show the implementation of a synthetic diagnostic for high-k scattering in NSTX [26]. We introduce the diagnostic in section 4.1 and present the numerical resolution details in section 4.2. Simulation spectra outputs from GYRO simulations are shown in section 4.3, and synthetically generated spectra are shown in section 4.4. Comparisons between experimental and simulated spectra are shown in section 4.5.

4.1. High-k diagnostic in NSTX

A high-k scattering diagnostic designed for the measurement of electron density fluctuations on the electron gyro radius scale ($k_{\perp} \rho_e \lesssim 0.6$) was designed, built and operated in NSTX [26]. This high-k scattering system used a 280 GHz microwave beam source of 15 mW, propagating close to the midplane in a tangential geometry with respect to the flux surfaces, as can be seen in figure 2. In this geometry, the measured wave vectors are primarily radial k_x, with a smaller vertical component k_z satisfying $k_{z}/k_{x} \approx$ 0.2−0.3. The scattering system consisted of five collection channels that simultaneously measure five different wave numbers in the range $5 \lesssim k_{\perp} \lesssim 30$ cm⁻¹. Heterodyne receivers installed on each channel allowed us to determine the direction of propagation of the observed fluctuations. The wavenumber resolution of the observed electron density fluctuations is Δk ≈ ±0.7 cm⁻¹ and the radial resolution is ΔR ≈ ±3 cm. The near mid-plane trajectory of the probe beam and the k-response are computed using a ray-tracing code. Figure 2 shows the trajectory of four channels of the high-k scattering system for NSTX shot 141 767, which has been extensively analyzed in [7, 52, 53]. The scattering system is sensitive to fluctuations taking place at R ≈ 135 cm (r/a ∼ 0.7). For reference, the major and minor radii of NSTX are, respectively, $R_\text{maj} = 0.85$ m, minor radius a = 0.68 m. Channels 1, 2 and 3 measure $k_{x} \rho_s \sim$ 8–13 and $k_{z} \rho_s \sim$ 1.5–2.5, which in physical units correspond to k_x ∼ 11−19 cm⁻¹ and k_z ∼ 2.4−3.5 cm⁻¹ (ρ_s is computed using local values of electron temperature T_e and magnetic field from LRDFIT equilibrium reconstruction). Additional details can be found in table 2. The electron and ion gyro-radii typically have values ρ_e ≈ 0.1 mm and $\rho_i \approx \rho_s \approx 0.7$ cm in these NSTX plasmas. Note the $(k_x, k_z)$ definitions employed in this manuscript are identical to those defined in [7], but do not correspond to the $(k_r, k_b)$ definitions employed in [52].

Zoom In Zoom Out Reset image size

4.2. Nonlinear gyrokinetic simulation set-up

In our attempt to establish quantitative comparisons of electron-scale turbulence we present two types of non-linear gyrokinetic simulations: standard electron-scale gyrokinetic simulation featuring a box size characteristic for the resolution of electron-scale modes $(L_r, L_\theta) = (4.5, 4) \rho_s$ and 'big-box' electron-scale simulation with an increased simulation domain $(L_r, L_\theta) = (20, 20.6) \rho_s$. The increased simulation domain results in a finer wavenumber grid resolution, which proves necessary to resolve the experimental wavenumbers from the high-k system.

4.2.1. Physics parameters and numerical resolution

The physics parameters employed in both simulation types are taken from NSTX H-mode plasma shot 141 767. Standard electron-scale and 'big-box' electron-scale simulations model three gyrokinetic species ($e-, D^{+}, C^{+6}$). Simulations are performed in the local, flux-tube limit at the scattering location r/a  ∼ 0.7, including electron collisions (ν_ei ∼ 1 c_s/a, but not ion collisions), background flow and flow shear (M ∼ 0.2− .3, γ_{E ∼} 0.1−0.2 c$_s/a, \gamma_p \sim 1$ c_s/a) and fully electromagnetic fluctuations $(\delta \phi, \delta A_{||}, \delta B_{||})$. Linear background profiles were simulated employing non-periodic boundary conditions in the radial direction with typical buffer widths $\Delta_b \sim 1/2.5 \rho_s$, respectively, for standard electron-scale and 'big-box' electron-scale simulation. Parallel resolution employed 14 poloidal grid points (× 2 signs of parallel velocity), 12 energies and 12 pitch-angles (6 passing + 6 trapped). This choice of numerical grids was made according to previous convergence and accuracy tests for the GYRO code simulating microinstabilities in the core of NSTX [40] and was also tested for convergence for the present conditions [7].

4.2.2. Radial and poloidal wavenumber resolution

The radial and poloidal wavenumber resolution of the non-linear simulations is of crucial importance in order to accurately resolve the measured wavenumbers by the high-k system. Standard electron-scale simulation resolves only electron-scale turbulence wavenumbers $k_r\rho_s \ \epsilon \ [1, 50]$ and $k_\theta\rho_s \ \epsilon \ [1.5, 65]$. 'Big-box' electron-scale simulation resolves electron-scale turbulence, similarly as electron-scale simulation, but includes modes characteristic of low-k instabilities, typically $k_r\rho_s \ \epsilon \ [0.3, 40]$ and $k_\theta\rho_s \ \epsilon \ [0.3,$ 65− 85] depending on the plasma condition. However 'big-box' electron-scale simulation does not correctly resolve the full spectrum of ion-scale turbulence (which would require a poloidal wavenumber grid spacing dk_θs ∼ 0.05−0.1). In addition, simulations are only run for electron time scales (T ∼ 20−30 a/c_s, when ions have not had time to reach a fully saturated state). Consequently, 'big-box' electron-scale simulation should not be considered as a multiscale simulation, such as the ones documented in [57–61] by Howard et al and Maeyama et al. All electron-scale simulations presented are converged in radial and poloidal box-sizes L_r and $L_\theta$, radial resolution dr, poloidal wavenumber resolution $\textrm{max}(k_\theta\rho_s)$, as well as in simulation run-time T [7]. Additional numerical resolution details can be found in table 1.

Table 1. Numerical resolution parameters typical of a standard and a 'big-box' electron-scale simulation: dr is the radial resolution ($\rho_s, \rho_e$ are, respectively, the ion and electron sound gyro radius using electron temperature T_e), $L_r [\rho_s]$ is the radial box size, n_r is the number of radial modes, max$(k_r\rho_s)$ is the maximum radial wavenumber resolved, $L_\theta[\rho_s]$ is the poloidal box size, d$k_\theta\rho_s$ is the poloidal wavenumber resolution, max$(k_\theta\rho_s)$ is the maximum poloidal wavenumber resolved, n_n is the number of toroidal modes, T is the simulation run time, dt is the simulation time step. Both simulation models only resolve electron-scale modes.

	dr	$L_r (\rho_s)$	nr	max$({k_r\rho_s})$	$L_\theta (\rho_s)$	d$k_\theta \rho_s$	max$({k_\theta \rho_s})$	nn	T(a/cs)	dt(a/cs)
std. e-scale	2ρe	4.5	≈ 200	50.5	4	1.5	65	42	30	${10^{ - 3}} - {5.10^{ - 4}}$
'big-box' e-scale	2.5ρe	20	≈ 500	30	20.6	0.3	65	≈ 200	30	${10^{ - 3}} - {5.10^{ - 4}}$

Table 2. Comparison between the measured wavenumber components via high-k scattering from channels 1, 2 and 3 in Cartesian coordinates $(k_x, k_z)$, and the corresponding values mapped to the field-aligned coordinate definitions $(k_r, k_\theta)$ in GYRO/CGYRO [49–51]. Equation (8) was used to perform the mapping.

	Cartesian		Field-aligned
	$\mathbf{{\mathit{k_x}} (cm^{-1})}$	$\mathbf{{\mathit{k_z}} (cm^{-1})}$	$\mathbf{{\mathit{k_r\rho_s}}^{unit}}$	$\mathbf{{\mathit{k_\theta\rho_s}}^{unit}}$
Ch 1	–18.84	–3.47	–2.57	–5.34
Ch 2	–14.87	–2.82	–2.05	–4.21
Ch 3	–11.00	–2.37	–1.54	–3.45

Figure 3 displays the radial and poloidal wavenumber simulation grid from a typical electron-scale simulation (left) and a 'big-box' electron-scale simulation (right), along with mapped wavenumbers detected by channels 1, 2 and 3 of the high-k scattering system. The black dots denote the dominant wavenumber detected by each diagnostic channel $\vec{k}_+$, and the ellipses surrounding them are the wavenumber resolution, denoting the 1/e amplitude of the effective wavenumber filter (scattering matrix). Ideally, one would want to simulate several radial and poloidal wavenumbers inside each $(k_r, k_\theta)$ ellipse to accurately replicate the experimental fluctuation measurement. However, due to a coarse wavenumber grid spacing, standard electron-scale simulation can at most resolve one radial and two poloidal wavenumbers inside the effective measurement range delimited by the elliptical shape of the wavenumber filter. This poor resolution due to the diagnostic requirements and numerical resolution requirements will result in inaccurate synthetic frequency spectra computed from electron-scale simulation (figure 8). By decreasing the wavenumber grid spacing (dk_r, d$k_\theta$), a 'big-box' electron-scale simulation can effectively filter a handful of poloidal wavenumbers inside the measurement range from each channel, yielding it adequate for attempting quantitative turbulence spectra comparisons. This results in computationally intensive simulations, typically running on 10−20 thousand parallel CPU cores taking ∼1− 2 M CPU hours to completion on leadership high-performance supercomputers such as NERSC's Edison.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Values for the total scattered power P_tot [a.u.], spectral peak <ω>  [c_s/a] and spectral width $\sigma_\omega\ [c_s/a]$ corresponding to synthetic frequency spectra from figure 8. Similar values of the total power P_tot and spectral peak <ω> are obtained between the two simulation models (standard versus 'big-box' electron-scale simulation). The spectral width $\sigma_\omega$ is wider for 'big-box' electron-scale simulation, in closer agreement to the experimental value (section 4.5).

	Ptot [a.u.]	<ω> [cs/a]	$\sigma_\omega\ [c_{s}/a]$
Std. e- scale	1.82 10−11	−24.18	3.04
'Big-box' e- scale	2.17 10−11	−23.29	5.13

4.2.3. Measured wavenumbers

Ray-tracing provides the measured wavenumbers in the components $(k_x, k_y, k_z)_+$. These are mapped to $(k_r, k_\theta)$ in figure 3 by use of the full 3D mapping (equation (8)). Since the measurement is local and close to the outboard midplane (poloidal location $\theta_0 \sim -4^o$), in what follows we present alternative calculations of the wavenumber mapping using the 2D and outboard midplane approximation (equation (11)), as well as a 'naive' mapping (equivalent to unshifted, circular flux surface geometry). These are presented for channel 1, which is sensitive to $\vec{k}_+ = (-1\ 883.5, -347.2, 138.2)$ m⁻¹.

The 2D and outboard midplane approximation requires computing the geometric factors $|\nabla r|_0, \kappa, q, \frac{\partial \alpha}{\partial \theta}|_0$ and the GYRO normalizing $\rho_s^{\textrm{unit}}$. We use the flux-surface geometry of NSTX H-mode plasma 141 767, and find $|\nabla r|_0 \approx 1.43$, κ ≈ 2.11, q ≈ −3.79, $\frac{\partial \alpha}{\partial \theta}|_0 \approx 1.33$ and $\rho_s^{\text{unit}} \approx 0.2$ cm. Note here how $\rho_s^{\text{unit}}$ is ∼ 3× smaller than the experimental, local value of ρ_s∼0.7 cm (due to the normalizing field $B_{\textrm{unit}}$ in GYRO, which does not correspond to the local value). Using the 2D and outboard midplane approximation (equation (11)) we find that the wavenumber components from channel 1 map to $(k_r\rho_s^\textrm{unit}, k_\theta\rho_s^\textrm{unit}) \approx (-2.59, -4.11)$. For comparison, the full 3D mapping of equation (8) gives $(k_r\rho_s^\textrm{unit}, k_\theta\rho_s^\textrm{unit}) \approx (-2.57, -5.34)$, which corresponds to the plotted values for channel 1 in figure 3. The k_r component is very well reproduced, but an error of ≈ 20%−25% is produced in the $k_\theta$ component. This discrepancy emphasizes the importance of using the full mapping for realistic tokamak geometry, even for rather small off-midplane poloidal angles $\theta_0 \sim -4^o$. This is partly due to the larger k_x + component of the high-k scattering system and the high flux surface shaping of this spherical tokamak plasma (which make the $\frac{1}{r_0}\frac{\partial R}{\partial \theta}\big|_0k_{x+}$ term non-negligible in the second line of equation (8)). The result of the mapping applied to channels 1, 2 and 3 is given in table 2.

Using the local ρ_s value as in experimental measurements gives the normalized $(k_x, k_z)$ components from channel 1 $(k_x\rho_s, k_z\rho_s)_+ \approx (-13, -2.4)$. If one were to make the 'naive' mapping $(k_x, k_z) = (k_r, k_\theta)$ and ignore the different ρ_s definitions, a significant systematic error of factor of ∼5× would be performed in the interpreted k_r, and ∼ 2× in the interpreted $k_\theta$. As a result, we learn that the high-k scattering system is sensitive to fluctuations with a lower k_r and larger $k_\theta$ than predicted by the 'naive' mapping, bringing it closer to the streamer peak of fluctuations (figure 6). This makes the NSTX high-k scattering system more transport relevant than previously thought, and emphasizes the importance of performing the wavenumber mapping (equation (8) or (11), depending on conditions) in order to correctly interpret the measurement range of the high-k diagnostic. As an illustration, figure 6 shows the mapped wavenumbers from channels 1, 2 and 3 using the full 3D mapping in black dots (equation (8)), and in white dots using the 'naive' mapping $(k_x, k_z) = (k_r, k_\theta)$. These are superimposed on the 2D density fluctuation wavenumber spectrum $S(k_r, k_\theta)$.

4.2.4. Filters in wavenumber space and real space

For the same plasma discharge condition, figure 4 shows the shape of the radial and poloidal wavenumber filters $W_{k_r}$ and $W_{k_\theta}$ in the θ₀ ≈ 0 approximation (equation (12)), using a simulation grid from a standard electron-scale simulation (red) and from a 'big-box' electron-scale simulation (blue). The Gaussian shape comes from the Gaussian shape of the scattering volume U (equations (6)). Notice the lack of resolution when using a standard simulation domain (red), particularly in k_r, and the improved resolution when using a bigger simulation domain (blue), due to the finer $(k_r, k_\theta)$ grid resolution. The 'virtual' dashed line shows the theoretical Gaussian expression of the filter. Note we chose ΔZ = 3 cm (the experimental value) and ΔR = 1 cm (reduced from the experimental ΔR = 3 cm due to the reduced simulation domain, even for the increased box size). The reduced radial extent of the scattering volume ΔR in a local simulation only scales the fluctuation amplitude by a constant value of irrelevance in the current work.

Figure 5 shows a snapshot of the raw 2D electron density fluctuation field δn(R, Z, ϕ₀ = 0) for a standard electron-scale simulation in 5(a) and for a 'big-box' electron-scale simulation in 5(b). The elongated streamer structures are tilted by the strong E × B shear flow. (c) and (d) show the corresponding electron density fluctuation field filtered by a 2D filter in real space $U(R-R_0, Z-Z_0)$ and selecting the $\vec{k}_+$ wave-vector, respectively for a standard (c) and 'big-box' electron-scale simulation (d). Volume integration of $\textrm{Re}[\delta n(R, Z, \varphi_0 = 0)U(R-R_0, Z-Z_0)e^{-i \vec{k}_+\cdot \vec{r}}]$ directly yields the synthetic signal of fluctuations $\delta n_\textrm{u}$ as in equation (3) (Re denotes the real part). Figures (c) and (d) give a spatial representation in real-space of the detected structures by the high-k system, and illustrate the effect of the filtering by U and wavenumber selection via the complex exponential $e^{-i \vec{k}_+\cdot \vec{r}}$. Since simulations are run in the local approximation, profile parameters are constant within the radial domain and the radial filter is chosen to be constant Ψ_r = 1. The poloidal filter shape is Gaussian in θ and mapped to (R, Z), having maximum amplitude at the thick black line passing through Z₀ ≈ −0.06 cm. The additional black dashed lines denote the 1/e, 1/e² and 1/e³ amplitude of the filter in the poloidal direction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3. Electron-scale simulation spectra

In this subsection we show electron-scale simulation spectra to gain insight into the measurement range of the high-k diagnostic. Spectral differences between a standard electron-scale and a 'big-box' electron-scale simulation are also discussed.

Figure 6 shows the GYRO 2D electron density fluctuation power spectrum from a standard electron-scale simulation in (a) and a 'big-box' electron-scale simulation in (b), proportional to the spectral density $S(k_r, k_\theta)$. $(k_r, k_\theta)$ are the internal field-aligned definitions in GYRO. We define $S(k_r, k_\theta) = \frac{\langle |\delta n_{np}|^2 \rangle_{\theta, T} }{(\text{d} k_r\rho_s) (\text{d} k_\theta \rho_s)}$, where $\langle . \rangle _{\theta, T}$ denotes the θ and time averages, and d$k_r\rho_s$, d$k_\theta\rho_s$ are the simulation radial and poloidal wavenumber grid resolutions. The black dots surrounded by ellipses correspond to the wavenumber measurement range from channels 1, 2 and 3 of the high-k diagnostic, the same as figure 3 shows. The spectrum is not symmetric in $k_\theta$ and is tilted due to the high E × B flow shear. The highest spectral power given by streamers is characterized by finite $({k_r} > 0,{k_\theta } < 0)$ and $({k_r}\left\langle {0,{k_\theta }} \right\rangle 0)$, consistent with the tilt of streamers in real space (figures 5(a) and (b)). Figure 6 is only intended to give a qualitative idea of the measurement wave numbers in the simulated fluctuation spectrum, since it is the amplitude δn_np that should be filtered in the scattering process (preserving phase information), but not the spectral density $S \propto |\delta n_{np}|^2$. This can be clearly seen from equations (2), (3) and (4).

Figure 7(a) and (b) shows the k_r and $k_\theta$ electron density fluctuation power spectrum S(k_r) and $S(k_\theta)$ from a standard electron-scale simulation (red) and from a 'big-box' electron-scale simulation (blue). We define

$$\begin{equation} \left\{ \begin{aligned} S(k_r) = \text{d} k_\theta \rho_s \sum_{k_\theta\lt 0}S(k_r, k_\theta) \\ S(k_\theta) = \text{d} k_r \rho_s \sum_{k_r}S(k_r, k_\theta) \\ \end{aligned}. \right. \end{equation} \tag{ 14 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The choice of $k_\theta \lt 0$ is made due to the symmetry property of the density spectrum $\delta n_{np} = \delta n_{-n-p}^*$ (where ^* indicates the complex conjugate [49]). Vertical black lines indicate the measurement wavenumbers by the high-k system. Inspection of S(k_r) in figure 7(a) emphasizes that the measurement is not aligned in k_r with the highest amplitude streamer fluctuations (for $k_\theta \lt 0$, streamers have positive k_r > 0 while the measurement is made for $k_r\lt 0$). This suggests that it may have been possible to detect streamer fluctuations in the present experiment, had the measurement been designed for $k_\theta \lt 0$ and k_r > 0. Since k_r changes sign, the k_r > 0 and $k_r \lt 0$ branches are plotted, exhibiting close to an order of magnitude difference in spectral power with respect to the streamer branch (k_r > 0) near $k_r\rho_s\sim 1-2$. The k_r spectra between standard and 'big-box' electron-scale simulation exhibit quantitative agreement, from the low-k wavenumber peak of the spectrum to the spectral slope at higher k_r.

Figure 7(b) shows the $k_\theta$ density fluctuation power spectrum from a standard electron-scale simulation (red) and from a 'big-box' electron-scale simulation (blue). Due to the logarithmic scale and the symmetry property in δn_np (translating to $S(k_\theta) = S(-k_\theta)$), the $k_\theta \lt 0$ branch only is shown here and $k_\theta$ should be interpreted as having a negative sign (similarly to the negative k_r branch in (a)). Figure 7(b) shows that the 'big-box' electron-scale simulation exhibits a quantitatively similar spectrum to a standard electron-scale simulation, showing similar wavenumber peaking and spectral slopes. However, the predicted spectral power is about ∼20% lower for the 'big-box' electron-scale simulation. This difference, however, lies within the simulation standard deviation of the total turbulent power. To summarize, figures 6 and 7 show how the 'big-box' electron-scale simulation spectra is quantitatively similar to that of a standard electron-scale simulation, providing ultimate confidence that the resolved electron temperature gradient (ETG) physics are very similar between the two simulation models.

4.4. Synthetic spectra and Doppler shift

In this section we deploy the 2D synthetic diagnostic for high-k scattering in k-space to show some of the spectral features of the synthetic spectra. The equivalence between the real space and the k-space implementation of the synthetic diagnostic is essentially identical, as shown in appendix C.

Figure 8 shows the spectral density $S(\vec{k}_+, \omega)$ predicted from a standard electron-scale simulation (red) and 'big-box' electron-scale simulation (blue), exhibiting the same plasma physics parameters but different wavenumber grid-resolution. The synthetic spectra exhibit qualitatively similar features such as frequency response and power levels (within ∼20% agreement), which are quantified in table 3. However, the spectra exhibit at least two appreciable differences. First, the spectrum in red exhibits a 'double-peak' structure in frequency, whereas the blue spectrum shows only one hump. Second, the spectrum in blue is wider than the spectrum in red. The particular values of the total scattered power P_tot, the spectral peak <ω> and the spectral width $\sigma_\omega$ are computed by analysis of the turbulence frequency spectrum $S(\vec{k}_+, \omega)$ in figure 8 (throughout this manuscript we denote 'spectral peak' as the frequency value <f > or <ω> at the peak of the spectral power). Table 3 shows the specific numeric values of P_tot, <ω> and $\sigma_\omega$.

As discussed in figures 3 and 4(b), the reduced wavenumber grid-resolution from a standard electron-scale simulation only allows it to sample a maximum of two simulation toroidal mode numbers contributing to the high-k signal. The two peaks in the red curve of figure 8 correspond to the two dominant toroidal mode numbers within the measurement $k_\theta$ range. Each mode has its own propagation frequency and is additionally Doppler-shifted by a different amount ($\omega_{Dop} = \vec{k}_+\cdot\vec{v} \sim n_+ \overline{\omega}_0$, where n₊ is the sampled mode number and $\overline{\omega}_0$ is the plasma toroidal rotation frequency). This results in a separation of spectral peaks in the frequency spectrum due to the two dominant modes contributing to the synthetic signal when using a small simulation domain (red spectrum). This phenomenon is not present when using an increased simulation domain (blue curve) due to the increased number of sampled modes. Using a bigger simulation domain increases the number of toroidal modes sampled within the measurement $k_\theta$ range, which tends to 'fill-in' the frequency spectrum and yield a single frequency feature in figure 8. This last point also contributes to a widening of the spectrum from a value of $\sigma_\omega \sim 3 a/c_s$ (red) to a value of ~ 5a/c_s (blue, table 3), yielding an improved agreement with the experimentally detected spectra (section 4.5).

In the present conditions, the toroidal rotation level (Mach number M ∼ 0.2) added to the relatively high poloidal wavenumbers sampled ($k_{\theta+}\rho_s \sim 3$–5) contribute to a Doppler shift frequency ω_Dop that can exceed the plasma-frame frequency of fluctuations by factors of ∼10× or more (recall $\omega_{Dop} \sim k_{\theta+} \omega_0 r_0/q_0$). The frequency spectrum of the high-k system is dominated by Doppler shift. To illustrate this, figure 9 shows the frequency power spectrum from a 'big-box' electron-scale simulation in which no Doppler-shift was applied to the fluctuations (gray), and with the experimental Doppler shift value applied (blue). The different frequency response between the two spectra highlights the important effect of Doppler shift in these conditions, shifting the peak of spectral power from $\langle\omega\rangle_{M = 0} \sim 1.54\ c_s/a$ with no Doppler shift applied (ω > 0 is in the electron diamagnetic drift direction) to $\langle\omega\rangle_{M = M^{\textrm{exp}}} \sim -23.29\ c_s/a$ when Doppler shift is applied (Doppler shift shifts the frequency spectrum to the ion diamagnetic drift direction). Other frequency spectra quantities such as the total power P_tot and spectral width $\sigma_\omega$ are quantitatively similar (table 4). A difference of ∼30% can be observed in the total spectral power P_tot, which might be related to numerical errors, but also to a pure effect of Doppler shift.

Table 4. Values for the total scattered power P_tot [a.u.], spectral peak <ω> [c_s/a] and spectral width $\sigma_\omega [c_s/a]$ corresponding to synthetic frequency spectra from figure 9, showing the impact of Doppler shift on the frequency spectrum (M is the Mach number). Similar values of the total power P_tot and spectral peak $\sigma_\omega$ are obtained between the two simulation models, but the spectral peak is completely different: fluctuations propagate in the electron diamagnetic direction in the plasma frame (ω > 0 for GYRO), but are shifted strongly towards the ion diamagnetic drift direction by Doppler shift.

	Ptot [a.u.]	<ω> [cs/a]	$\sigma_\omega\ [c_{s}/a]$
Plasma frame: M = 0	1.62 10−11	1.54	4.89
Lab frame: $M = M^{\textrm{exp}}$	2.17 10−11	−23.29	5.13

Zoom In Zoom Out Reset image size

One can gain further insight into the effect of Doppler shift by studying the 2D frequency spectra plots $(\omega, k_\theta)$ and (ω, k_r). Figures 10(a) and (b) show the synthetic frequency power spectrum of fluctuations $S_{k_{r+}}(k_\theta, \omega)$ computed from a 'big-box' electron-scale simulation, where spectra have been filtered in k_r around the radial component k_r + (from channel 1), but not in $k_\theta$. The black vertical band shows the measurement range in $k_\theta$. No Doppler shift is applied in (a), and the experimental Doppler shift value is applied in (b). Figures 10(a) and (b) show that the effect of Doppler shift is primarily a shift in frequency ω proportional to $\omega \propto k_\theta \overline{\omega}_0$ as expected. A smaller effect is a widening of the spectrum $\sigma_\omega$, which is negligible for the present conditions but becomes more important for higher toroidal rotation values and higher $k_\theta$. The spectrum also shows how the positive ${k_\theta } > 0$ part of the spectrum exhibits higher spectral power than the negative $k_\theta\lt 0$ counterpart, consistent with the spectrum shape from figures 6 and 7.

Zoom In Zoom Out Reset image size

Figures 10(c) and (d) show the synthetic frequency power spectrum of fluctuations $S_{k_{\theta+}}(k_r, \omega)$ computed from a 'big-box' electron-scale simulation, where spectra have been filtered in $k_\theta$ around the poloidal component $k_{\theta+}$ (channel 1), but not in k_r. The black vertical band shows the measurement range in k_r and the white dashed line denotes the k_r = 0 line. Differently to the $(\omega, k_\theta)$ spectra, the (ω, k_r)-spectra show that Doppler shift essentially produces a similar frequency shift for all k_r (corresponding to the same $k_\theta$ since it has been previously filtered in $k_\theta$), as expected for the current conditions (this may not be the case in far off-midplane scattering locations). As observed for the corresponding $(\omega, k_\theta)$ plots, one can notice the asymmetry in $S_{k_{\theta+}}(k_r, \omega)$ for positive versus negative k_r: a higher spectral power is observed for positive k_r > 0 fluctuations (to the right of the vertical white dashed line). This is once more consistent with the $(k_r, k_\theta)$ spectra from figures 6 and 7.

4.5. Comparisons with high-k scattering fluctuation measurements from NSTX

In this subsection we compare synthetically generated frequency and wavenumber spectra with experimentally detected high-k spectra. The present experimental conditions correspond to a highly unstable ETG regime from NSTX H-mode plasma discharge 141 767, which has been extensively analyzed in [7, 52, 53]. In this regime, ETG was shown to be the only microturbulence process able to predict experimentally relevant levels of electron heat flux, and ion-scale turbulence was shown to be fully suppressed by strong E × B shear flow.

Synthetic spectra are obtained from two GYRO 'big-box' electron-scale simulations (details in appendix A). The first one uses the nominal experimental profile values as input (blue), and predicts ∼30% of the experimental electron heat flux value. In the second one (red), the values of the normalized electron density gradient a/L_ne, safety factor q and magnetic shear $\hat{s}$ are scaled from the experimental values in a sensitivity scan to maximize the ETG drive, which was able to reproduce the experimental electron heat flux value within uncertainty. The values are: a/L_ne = 0.502 4 is the ETG stabilizing mechanism, scaled by 1σ experimental uncertainty, q = 3.410 3 (−10%) and $\hat{s} = 2.165\,6$ (+20%). The uncertainty σ(a/L_ne) was computed from uncertainty in the background electron density profile followed by a Monte Carlo analysis approach. More details on this discharge condition and the corresponding turbulent transport fluxes can be found in appendix A and [7].

Figures 11(a), (b) and (c) show the frequency power spectrum of the high-k scattering system from channels 1, 2 and 3, respectively, plotted in real frequencies f (MHz). In black are the experimental frequency spectra, and in blue the corresponding synthetic spectra from a 'big-box' electron-scale simulation using nominal inputs. Two main differences can be observed: (i) experimental spectra exhibit a high spectral peak at zero frequency, corresponding to spurious reflections of the input microwave beam in the plasma; and (ii) experiment exhibits an increased background noise level with respect to simulation. In fact, experimental spectra are affected by electronic noise as well as by additional electromagnetic emission from the plasma that does not directly correspond to a scattering process. For these reasons, experiment and simulation comparisons are done in a prescribed frequency band, delimited by the black vertical dashed lines in figure 11. Experimental spectra are not absolutely calibrated, and are consequently rescaled by a constant. The particular value of the scaling constant is chosen in order to minimize the total integrated power from the different channels with respect to the synthetic frequency spectra (a least-squares minimization of the power difference). All channels are scaled by the same constant, preserving the fluctuation level ratio and the k-spectrum shape.

Zoom In Zoom Out Reset image size

Qualitatively, experimental and synthetic spectra in figure 11 exhibit similar frequency response. The spectral peak <f> shows quite good agreement with experiment, lying within ∼10% for all channels. The spectral width σ_f lies within ∼25% of the experimental value (for channel 3, the reduced Doppler shift from the experimental spectrum with respect to the f = 0 peak does not allow a reliable measure of the spectral width σ_f). Table 5 shows the particular values of <f > and width σ_f (kHz). With respect to the spectral power, channels 1 and 2 in figure 11 exhibit reasonable agreement, however, the synthetic power level from channel 3 overpredicts the experimental power level by over an order of magnitude. The disagreement in the spectral power is also clearly seen in the wavenumber spectrum of fluctuations in figure 13, and suggests this simulation is likely missing some necessary physics ingredient that is able to reproduce the shape of the wavenumber spectrum.

Table 5. Summary of the main frequency spectrum characteristics corresponding to figures 11 and 12. The spectral peak <f> and the spectral width σ_f (kHz) from the different channels are compared, corresponding to experiment and simulation. Two simulations are used, one using the nominal experimental profile values as input, and also a simulation using scaled gradients (within 1σ(∇n_e), −10% q and $+20\% \ \hat{s}$). d_P denotes a validation metric of distance between experiment and simulation total fluctuation power P_tot, as suggested by [62, 63]. The σ_f = 196 (kHz) value from channel 3 is given in parentheses due to the unreliability of measurement from the reduced Doppler shift (figures 11 and 12).

	dP	ch 1		ch 2		ch 3
/(kHz)		<f >	σf	<f >	σf	<f >	σf
Exp.	—	–895	214	–810	199	–613	(196)
Sim. (nominal)	5.08	–843	186	–735	157	–664	127
Sim. (scan)	1.53	–878	196	–746	147	–672	126

Figure 12 shows similar experiment and synthetic frequency spectra comparisons, where the synthetic spectrum this time is computed from a simulation with scanned inputs in a/L_ne, safety factor q and magnetic shear $\hat{s}$, chosen to maximize the ETG drive. A different scaling constant than in figure 11 is applied to the experimental frequency spectra in order to minimize the difference spectral power between experiment and simulation (however, all channels are scaled by the same constant). A similar agreement within ∼10% and ∼25% is observed in the spectral peak <f> and width σ_f, respectively (table 5). Recalling that <f > is completely dominated by Doppler shift, agreement in <f> does not imply agreement with the intrinsic frequency of fluctuations in the plasma frame. With respect to σ_f, as shown in table 4, σ_f appears to be insensitive to Doppler shift in the present conditions. This means that the measured spectral width is essentially the same as the intrinsic plasma frame value, which is reproduced by the synthetic spectra within 25% agreement. Interestingly, the spectral peak and width from both simulations exhibit very similar values, suggesting the frequency spectra characteristics are insensitive to the turbulence drive.

Zoom In Zoom Out Reset image size

Figure 13 shows the wavenumber spectrum comparisons between experiment and simulation using the nominal experimental parameters (blue, corresponding to figure 11) and using the scanned input values (red, corresponding to figure 12). The wavenumber spectrum is calculated from the frequency spectrum by integration of the frequency spectra within the prescribed frequency band for each of the three available channels. In figure 13 only, it is the synthetic spectra that are scaled by a constant, not the experimental spectra. A different normalization constant minimizing the distance with respect to experimental spectra is applied to the blue and red synthetic spectra to yield a comparison of the shape of the wavenumber spectrum, but not the total fluctuation power.

Simulation using scaled gradients displays improved agreement in the shape of the k-spectrum. This can be quantified via a validation metric $d_P$ between experiment and simulation. We define the validation distance $d_P$ of the total fluctuation power P_tot as a metric ${d}_{P} = \sqrt{\sum_i({P_i^{\textrm{exp}}-P_i^{\textrm{syn}}})^2/((\Delta P_i^{\textrm{exp}})^2 + (\Delta P_i^{\textrm{syn}})^2)}$, as suggested in [62, 63]. Here the sum ∑_i is over the different channels (i = {1, 2, 3}), $P_i^{\textrm{exp}}$ and $P_i^{\textrm{syn}}$ are the total fluctuation power from experiment and simulation, respectively, and Δ indicates the standard deviation of the time series used to compute P_i. In the definition of $d_P$, $P_i^{\textrm{syn}}$ has been scaled by a constant in order to minimize the difference with respect to the experimental spectrum. $d_P$ is thus the minimum distance between experiment and simulation, where the minimization is carried out by the scaling constant (recall all channels are scaled by the same constant).

Using this metric, we find a value of d_P = 5.08 for the simulation spectra using the experimental profile values. The simulation using scaled inputs yields an improved distance of d_P = 1.53. For reference, a value of $d_P \approx 5$ is roughly equivalent to a difference between experiment and simulation of 7σ, while a value of $d_P \approx 1.5$ roughly equates to a difference of 2σ. This substantial improvement in the shape of the wavenumber spectrum shows the strong sensitivity of the k-spectra to input parameters, and highlights how reasonable agreement in the shape of the spectrum can be obtained within small variations in the simulation input drive terms ($a/L_{ne}, q, \hat{s}$). Reference [7] contains additional information about these and additional simulations performed for this NSTX discharge.

Zoom In Zoom Out Reset image size

We have presented a formulation of two synthetic diagnostics for coherent scattering of microwaves and applied it to the particular case of the high-k scattering diagnostic in NSTX [26]. This has yielded direct comparisons between experiment and simulation of high-k frequency and wavenumber turbulence spectra. The principles outlined in sections 2 and 3 remain quite general and applicable to additional fluctuation diagnostics such as DBS, reflectometry, and even cross-polarization scattering (CPS) measurements. Although we have formulated a full 3D synthetic diagnostic in sections 2 and 3, we have only deployed a 2D synthetic diagnostic in section 4. Appendix E includes additional insight on the differences expected in a 3D formulation.

This work has built on previous synthetic diagnostic efforts of high-k scattering [42] and DBS [45], which were based on the standard interpretation of scattering in k-space. We have shown the equivalence of the formulation in k-space to a formulation in real space. Agreement between the two formulations is also achieved numerically in realistic flux-surface geometry (appendix C), providing improved confidence in the synthetic calculations presented. Additional insight into the measurement can be gained via the implementation in k-space, which provides precious information about the measurement wave-vector $\vec{k}_+$ mapped to the field aligned $(k_r, k_\theta)$ components. In particular, the k-space implementation shows that a 'big-box' electron-scale simulation is more suitable for attempting quantitative comparisons with experimental turbulence measurements, as figure 3 shows. A detailed understanding of purely numerical artifacts in the frequency spectra (figure 8) was only possible with the wavenumber formulation of the synthetic diagnostic. This information would not have been available if only the real space formulation was implemented. We highlight that full understanding of the specific scattering measurement and corresponding numerically generated synthetic spectra can only be gained with the combined implementations in real space and in k-space.

Part of this work has shown how the measurement wavenumbers from the high-k scattering system are much closer to the peak of the fluctuation spectrum (streamers) than was previously thought. The mapping between cylindrical and field-aligned wavenumber components (equation (8)) has highlighted three main geometric effects affecting the interpretation of the measured wavenumber: the effect of the normalizing magnetic field entering the definition of ρ_s, the effect of Shafranov shift Δ affecting the radial wavenumber component k_r through compression of flux surfaces at the outboard midplane, and the effect of flux-surface elongation κ affecting the poloidal wavenumber $k_\theta$ through flux surface stretching poloidally. These three effects combine to yield systematic errors in the interpretation of the measured wavenumber components of up to factors of 5 in the present conditions. These systematic errors are amplified in conditions of strong Shafranov shift and highly shaped geometries in spherical tokamaks.

We find that the high-k scattering system in NSTX is sensitive to –3, 3–6) . Although the high-k scattering system was initially designed to be sensitive to high-k_r and low-$k_\theta$ fluctuations, we find that the corresponding $(k_r\rho_s)_+$ values are not as high-k_r as originally thought (as predicted by the 'naive' mapping), while the $(k_\theta\rho_s)_+$ component is higher (the high-k diagnostic is sensitive to high-k_x in the lab frame, but this does not necessarily map to high-k_r in the field-line frame). The smaller $k_r\rho_s$ and larger $k_\theta\rho_s$ values with respect to a 'naive' mapping indicate the measurement k by the high-k scattering system is more transport relevant than previously thought. Figure 6 shows how close the high-k measurement is to the streamer peak of the fluctuation spectrum. This hints at the possibility that the high-k scattering system may in fact already have been sensitive to streamer fluctuations for other experimental conditions and high-k scattering geometry. A careful analysis of additional NSTX plasma discharges to the one presented here would be needed in order to confirm this speculation. Additionally, a newly designed high-k scattering diagnostic is planned to be installed in NSTX-U. Projected to be sensitive to smaller k_r and higher $k_\theta$, it is also expected that this diagnostic will be able to detect density fluctuations from streamers.

This work has highlighted the profound impact of Doppler shift on the measured high-k signal, making the frequency spectrum of fluctuations completely dominated by Doppler shift. In particular, figure 9 and table 4 show how the spectral peak of fluctuations shifts from a value of < ω> ∼ 1.5 c_s/a in the plasma frame, to a value of ∼−23 c_s/a in the lab frame. This suggests the spectral peak of the high-k measurement is completely opaque to the intrinsic plasma frame value, making it nearly useless as a quantitative discriminator on the turbulence model. In fact we have observed very close agreement in the spectral peak for the different turbulence models tested (standard electron-scale simulation, 'big-box' electron-scale simulation, simulations with scaled inputs, etc), as long as a correct Doppler shift value is applied. This can also be observed from tables 3 and 5.

With respect to the spectral width $\sigma_\omega$, it is shown to be less affected by Doppler shift for the present case, however, one should be careful to interpret the measured width as an intrinsic turbulence value. As we have seen in figure 8 and table 3, a simulation with an increased simulation domain can modify the predicted spectral width from a value of $\sigma_\omega \sim 3$ to a value of ∼5 c_s/a, for the same simulation physics parameters. Even with this improvement, a 'big-box' electron-scale simulation tends to underpredict the measured spectral width of fluctuations by ∼25% (figures 11, 12 and table 5). We interpret the improved agreement achieved by the 'big-box' electron-scale simulation to be due to the increased numerical resolution in k-space as suggested by figure 3. This also suggests that simulations having even higher wavenumber resolution and including additional physics, i.e. as multiscale simulation, could in fact provide a closer match to the measured spectral width of fluctuations.

This discussion aims to make the reader aware of the difficulty in interpreting the measurement frequency spectra characteristics from high-k scattering. Fortunately, we find that all the previous effects polluting the high-k frequency spectrum appear to have a less profound impact on the total fluctuation power P_tot, implemented in the wavenumber spectrum. Integration of the frequency spectrum along a prescribed frequency band seems to 'erase' artificial numerical artifacts present in the frequency spectrum. The relative fluctuation power also provides a precious constraint on the simulations through the characterization of the shape of the fluctuation wavenumber spectrum. This allows the possibility of model selection and discrimination, as we have seen from figure 13 and table 5 (even in the absence of absolute diagnostic calibration). Even though the comparisons between experiment and simulation presented remain at a preliminary stage, they beg for additional quantitative comparisons of turbulence wavenumber spectra in future high-k scattering experiments. Improved diagnostic capability (only three channels are available in this work) as well as absolute diagnostic calibration would provide highly valuable validation constraints.

There is no doubt that the synthetic high-k turbulence predictions presented here suffer from uncertainties and inaccuracies emanating from the approximations made in the synthetic model. One important approximation lies in the 2D implementation, which is discussed in appendices E and F. An additional approximation in the synthetic diagnostic described here is based on a constant k: the same turbulence wavenumber is sampled within the whole simulation domain. In fact, the measurement wavenumber provided by ray-tracing calculations is only representative of the central ray of the input microwave beam. However, in reality a slightly different wavenumber is sampled by the diagnostic at different radial, poloidal and toroidal locations within the scattering volume. To assess the impact of this constant-k approximation, additional ray-tracing calculations for non-central rays were carried out within the scattering volume, showing that the measured k can vary at most by ∼20%. This would have a small impact on the synthetic frequency and k-spectrum characteristics when compared to other factors such as the simulation wavenumber resolution or Doppler shift. However, it is possible that taking into account the spatial variation of $\vec{k}_+ = \vec{k}_+(\vec{r})$ within the scattering volume could recover the underpredictions observed in the frequency spectral width in table 5 and figures 11 and 12. The spatial variation of $\vec{k}_+$ could be computed through ray-tracing of non-central rays, beam-tracing, or even full-wave simulation, and is left for future work.

The work presented here is part of the broader framework of validation of turbulent transport models. The ultimate goal of the validation effort is to provide confidence that current models are able to explain the transport processes observed in present experiments, in order to be able to predict the plasma profiles and ultimately fusion performance of future fusion reactors. Although seemingly far from this ultimate goal, establishing quantitative turbulence comparisons via synthetic diagnostics is a necessary step for the validation of the fundamental first principles based simulations. From these, reduced models can be optimized and be subsequently used for profile prediction and performance prediction of future devices. Great strides have been made so far to predict the plasma density and temperature profiles of the spherical tokamak, but work is still in its early stage. This work serves as a stepping stone towards validating first principles based electron-scale simulation in the core-gradient region of modest-β NSTX NBI-heated H-modes. We have found with reasonable confidence that electron-scale simulation is in fact able to reproduce the detected frequency and wavenumber spectra by the high-k scattering system in a highly unstable ETG regime. Improved confidence in the current transport models will only be possible by placing additional constraints via additional diagnostic measurements and further testing of the models in higher-β, lower collisionality and higher performance plasmas. Combining high-k measurements to low-k and intermediate-k fluctuation diagnostics (DBS, reflectometry, BES), magnetic field fluctuation measurements (CPS), etc, would provide invaluable information to constrain our turbulence models at all relevant spatial and temporal scales characteristic of microturbulence fluctuations. More importantly, they are imperative to test and validate our current models in the wake of future fusion generating devices such as ITER and beyond.

The authors would like to thank all co-authors for fruitful discussions and invaluable advice. This work is supported by the US Department of Energy Contract DE-AC02-09CH11466, and partially by the Engineering and Physical Sciences Research Council (EPSRC) [EP/R034737/1]. Computer simulations were carried out at the National Energy Research Scientific Computing Center, supported by the Office of Science of the US DOE under Contract No. DE-AC02-05CH11231, and at 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. The digital data for this paper can be found in http://arks.princeton.edu/ark:/88435/dsp018p58pg29j.

We define here all the relevant quantities necessary to obtain the real-space and k-space computation of the synthetic signal $\delta{n}_{\textrm{u}}(\vec{k}_+,t)$ in section 2, from knowledge of the density perturbation $\delta{n}(\vec{r'},t')$ and the filter in real space $U({\vec{r'}})$:

$$\begin{equation} \begin{array}{l} \delta n(\vec{k},t) = \int{{\mathrm{d}}^3\vec{r} \ \delta n(\vec{r},t)e^{-i\vec{k}\cdot \vec{r} }} \qquad \text{Fourier transform in space of } \delta{n}(\vec{r},t). \\[6pt] \delta n(\vec{k},\omega) = \int{{\mathrm{d}}^3\vec{r}\int{{\mathrm{d}}t \ \delta n(\vec{r},t)e^{i(\omega t - \vec{k}\cdot \vec{r} )}}} \qquad \text{Fourier transform in space and time of } \delta{n}(\vec{r},t). \\[6pt] W(\vec{k}-\vec{k}_{+}) = \int_{V_s}{{\mathrm{d}}^3\vec{r} \ U(\vec{r}) e^{i(\vec{k}-\vec{k}_{+})\cdot \vec{r}}} \qquad \text{Filter in wavenumber space, or }k{\textrm{-filter of fluctuations.}} \\[6pt] \delta{n}_{\textrm{u}}(\vec{k}_+,t) = \frac{1}{(2\pi)^3} \int{{\mathrm{d}}^3\vec{k}\ \delta n(\vec{k},t)} W(\vec{k}-\vec{k}_+) \qquad \text{Synthetic time signal of electron density fluctuations.} \\[6pt] \delta{n}_{\textrm{u}}(\vec{k}_+, \omega) = \frac{1}{(2\pi)^3} \int{{\text{d}}^3\vec{k}\ \delta n(\vec{k},\omega)} W(\vec{k}-\vec{k}_+) \qquad \text{Synthetic frequency signal of electron density fluctuations.} \\[6pt] \end{array} \end{equation} \tag{ B1 }$$

From the computation of the synthetic signal $\delta{n}_{\textrm{u}}(\vec{k}_+,t)$, here is a summary of the formulas needed to achieve the final expression for the scattered power P_s as a function of the spectral density $S(\vec{k}_+, \omega)$

$$\begin{equation} \begin{array}{l} \delta{n}_{\textrm{u}}(\vec{k}_+,t) = \int{{\mathrm{d}}^3\vec{r} \ U(\vec{r}) \delta{n}(\vec{r},t) e^{-i\vec{k}_+\cdot \vec{r} } }\\ \qquad\qquad = \frac{1}{(2\pi)^3}\int{{\mathrm{d}}^3\vec{k}\ W(\vec{k}-\vec{k}_+)\delta{n}(\vec{k},t) } \\[6pt] \delta{n}_{\textrm{u}}(\vec{k}_+, \omega) = \int{{\mathrm{d}}t\ \delta{n}_{\textrm{u}}(\vec{k}_+,t)e^{i\omega t}} \\[6pt] S(\vec{k}_+,\omega) = \frac{1}{TV_s} \Big< \frac{|\delta n_{\textrm{u}}(\vec{k}_+,\omega)|^2}{n_0} \Big> \\[6pt] S(\vec{k}_+) = \frac{1}{2\pi}\int{{\mathrm{d}}\omega S(\vec{k_+},\omega)} \\[6pt] \frac{{\mathrm{d}}^2P_s}{{\mathrm{d}}\omega {\mathrm{d}}\Omega} = \frac{P_0}{A_i} r_0^2|\hat{s} \times ( \hat{s} \times \hat{e} )|^2 \frac{n_0 V_s}{2\pi}S(\vec{k}_+,\omega) \\[6pt] \frac{{\mathrm{d}}P_s}{{\mathrm{d}}\Omega} = \frac{P_0}{A_i} r_0^2|\hat{s} \times ( \hat{s} \times \hat{e} )|^2 n_0 V_s\ S(\vec{k}_+) \\[6pt] \end{array} \end{equation} \tag{ B2 }$$

In this appendix we compare the output of the synthetic spectra predicted from the real-space formulation versus the k-space formulation. Figure C1 shows the synthetic frequency power spectrum of fluctuations $S(\vec{k}_+, \omega)$ corresponding to a filtered wavenumber $\vec{k}_+$ from channel 1 of the high-k scattering system. Simulations use a big numerical domain $(L_r, L_\theta) = (20, 20.6) \rho_s$. The formulation based on real-space (dashed line) is compared to the formulation based on k-space filtering (continuous line). Quantitative agreement is obtained between the two synthetic diagnostic formulations, achieving 15% agreement in the total power P_tot, and a further improved agreement in the spectral peak <ω> and spectral width $\sigma_\omega$. This agreement is not coincidental but is generally observed, and validates the implementation of both synthetic diagnostic methods in the context of realistic flux-surface geometries.

Zoom In Zoom Out Reset image size

In this section we present a derivation of the wavenumber mapping to field-aligned geometry. The wavenumber components $(k_x, k_y, k_z)_+$ mapped to field-aligned geometry allow for a direct interpretation of the measurement range of the high-k diagnostic. Figure D2 shows the reference geometry definitions in real space and in k-space that are used in this article.

Zoom In Zoom Out Reset image size

In real-space cylindrical coordinates (R, Z, ϕ), we define the wavenumber components $(k_x, k_z, k_y)$ by

$$\begin{equation} \left \{ \begin{array}{l} k_x = -i \frac{\partial }{\partial R} \qquad \text{along the major-radial direction} \\[6pt] k_z = -i \frac{\partial}{\partial Z} \qquad \text{along the vertical direction} \\[6pt] k_y = -i \frac{1}{R} \frac{\partial}{\partial \varphi} \qquad \text{along the toroidal direction} \end{array} \right. \end{equation} \tag{ D3 }$$

Since we assume axisymmetry, the reference frame is chosen at a fixed toroidal angle ϕ = 0 as figure D2 indicates. Associated to the field-aligned coordinates (r, θ, ϕ), the wavenumber components (k_r, $k_\theta$, $k_\varphi$) are internally defined in terms of n, p and $\frac{\partial }{\partial \varphi}$

$$\begin{equation} \left \{ \begin{array}{l} k_r = \frac{2\pi p}{L} \\[6pt] k_\theta = \frac{n q}{r} \\[6pt] k_\varphi = -i\frac{1}{R} \frac{\partial}{\partial \varphi} \end{array}. \right. \end{equation} \tag{ D4 }$$

Recall the expansion of a generic fluctuating field in GYRO/CGYRO (ignoring Doppler shift) is $f (r, \theta, \varphi, t) = \sum_{n,p} f_{np}(\theta, t) e^{-in\alpha} e^{i\frac{2\pi p}{L} r }$, and where α = ϕ + ν(r, θ) [49–51]. The radial wavenumber component k_r is based on a Fourier decomposition of f in the radial coordinate (note $k_r \neq -i \frac{\partial}{\partial r}$). The poloidal wavenumber $k_\theta$ is a flux surface quantity, independent of θ. The toroidal wavenumber $k_\varphi$ has the same definition as in cylindrical coordinates, here $k_\varphi = - n/R$. Note how $k_\varphi = k_y$, also consistent with figure D2(b). With these definitions, the wavenumber mapping from Cartesian coordinates $(k_x, k_y, k_z)$ to field-aligned coordinates $(k_r, k_\theta, k_\varphi)$ is computed for a fixed toroidal mode number n and a fixed radial mode number p (recall equation (4)). Defining $z_{np}(r,\theta, \varphi, t) = f_{np}(\theta,t) e^{-in(\varphi + \nu(r,\theta))} e^{i\frac{2\pi p}{L} r }$, we have

$$\begin{align} \frac{\partial }{\partial r} z_{np} &= i \bigg( \frac{2\pi p}{L} - n \frac{\partial \nu}{\partial r} \bigg) z_{np} \notag\\ \frac{1}{r}\frac{\partial}{\partial \theta} z_{np} &= \frac{1}{r} \bigg( \frac{\partial \log(f_{np}) }{\partial \theta} - i n \frac{\partial \nu}{\partial \theta} \bigg) z_{np} &= -i \frac{n}{r} \frac{\partial \nu}{\partial \theta} z_{np} \notag\\ \frac{1}{R}\frac{\partial}{\partial \varphi} z_{np} &= -i \frac{n}{R} z_{np} . \end{align} \tag{ D5 }$$

The last equality shows the definition of $k_\varphi = k_y$. Note how $\partial f_{np}/\partial \theta = 0$ since $\partial/\partial \theta$ only applies to the rapidly oscillating part of the fields. Making use of the chain rule yields

$$\begin{equation} \left \{ \begin{array}{l} \frac{\partial}{\partial r} = \frac{\partial R}{\partial r} \frac{\partial}{\partial R} + \frac{\partial Z}{\partial r} \frac{\partial}{\partial Z} + \frac{\partial \varphi}{\partial r}\frac{\partial}{\partial \varphi} = \frac{\partial R}{\partial r} i k_x + \frac{\partial Z}{\partial r} i k_z \\[6pt] \frac{\partial}{\partial \theta} = \frac{\partial R}{\partial \theta} \frac{\partial}{\partial R} + \frac{\partial Z}{\partial \theta} \frac{\partial}{\partial Z} + \frac{\partial \varphi}{\partial \theta}\frac{\partial}{\partial \varphi} = \frac{\partial R}{\partial \theta} i k_x + \frac{\partial Z}{\partial \theta} i k_z \end{array}. \right. \end{equation} \tag{ D6 }$$

Combining D5 and D6 yields

$$\begin{equation} \left \{ \begin{array}{l} \frac{\partial}{\partial r} z_{np} = i \bigg(k_r - n \frac{\partial \nu}{\partial r} \bigg) z_{np} = i \bigg( \frac{\partial R}{\partial r} k_x + \frac{\partial Z}{\partial r} k_z \bigg) z_{np} \\ \frac{1}{r}\frac{\partial}{\partial \theta} z_{np} = -i \frac{n}{r} \frac{\partial \nu}{\partial \theta} z_{np} = i \bigg( \frac{1}{r}\frac{\partial R}{\partial \theta} k_x + \frac{1}{r}\frac{\partial Z}{\partial \theta} k_z \bigg) z_{np} \\ \frac{1}{R}\frac{\partial}{\partial \varphi} z_{np} = -i \frac{n}{R} z_{np} = i ( k_y ) z_{np} \end{array}. \right.\label{} \end{equation}$$

This relation being true for any arbitrary function z_np, we arrive to the wavenumber mapping relations, relating a measured wavenumber $\vec{k}_+$ with components $(k_x, k_y, k_z)_+$ (in the form typically issued by ray-tracing codes or equivalent) to the field-aligned components $(k_r, k_\theta, k_\varphi)_+$

$$\begin{equation} \left\{ \begin{alignedat}{2} &k_{r+} - \frac{r}{q}\frac{\partial \nu}{\partial r} k_{\theta+} && = \frac{\partial R}{\partial r} k_{x+} + \frac{\partial Z}{\partial r} k_{z+} \\ &- \frac{r}{q}\frac{\partial \nu}{\partial \theta} k_{\theta+} && = \frac{\partial R}{\partial \theta} k_{x+} + \frac{\partial Z}{\partial \theta} k_{z+} \\ &k_{\varphi+} && = k_{y+} \end{alignedat}. \right. \end{equation} \tag{ D8 }$$

Equation (D8) is over constrained, and one needs to define two separate mapped toroidal mode numbers $n^\theta_+$ and $n^\varphi_+$, in addition to p₊ (which are generally not integers). Using α = ϕ + ν(r, θ) we recover formula (8) used in the main text:

$$\begin{equation} \left\{ \begin{alignedat}{3} & \frac{2\pi}{L_r} p_+ - \frac{\partial \alpha}{\partial r}\bigg|_0 n^{\theta}_+ && = \frac{\partial R}{\partial r}\bigg|_0 k_{x+} + \frac{\partial Z}{\partial r}\bigg|_0 k_{z+} \\ & -\frac{1}{r_0} \frac{\partial \alpha}{\partial \theta}\bigg|_0 n^{\theta}_{+} && = \frac{1}{r_0}\frac{\partial R}{\partial \theta}\bigg|_0 k_{x+} + \frac{1}{r_0}\frac{\partial Z}{\partial \theta}\bigg|_0 k_{z+} \\ & -\frac{n^{\varphi}_+}{R_0} && = k_{y+} \\ \end{alignedat}. \right. \end{equation} \tag{ D9 }$$

Using these relations, the wavenumber components $(k_x, k_y, k_z)_+$ typically given by experiments can be mapped to the field-aligned wavenumber definitions $(k_r, k_\theta, k_\varphi)_+$ employed by simulation codes using (D8), or equivalently to the mode number components $(n^\theta_+, n^\varphi_+, p_+)$ using (D9). In equation (D9) $n^\theta_+$ denotes a toroidal mode number mapped from a vertical component k_z + (not quite a poloidal component), while $n^\varphi_+$ is a toroidal mode number mapped from a toroidal component k_y +. These, in principle, need not be equal to each other. However, as we will see in the following appendices, the condition for successful scattering $\vec{k}\cdot\vec{B} \approx 0$ is equivalent to $n^\theta_+ \approx n^\varphi_+$.

In this appendix we show a more detailed derivation of the computation of the synthetic signal of density fluctuations $\delta{n}_{\textrm{u}} (\vec{k}_+, t)$ than the one presented in section 3. We start by rewriting the scattering signal in cylindrical coordinates (R, Z, ϕ), equation (3), repeated once more for reference:

$$\begin{equation} \begin{array}{l} \delta{n}_{\textrm{u}} (\vec{k}_+, t) = \int_{V_s} {\mathrm{d}}^3\vec{r'}\ \delta{n}(\vec{r'}, t) U(\vec{r'}) {e^{- i\vec{k}_+\cdot\vec{r'}} } = \\ \int{{\mathrm{d}}\varphi } {R{\mathrm{d}}R{\mathrm{d}}Z\ \delta{n}\ (R, Z, \varphi, t) U(R - R_0, Z - Z_0, \varphi - \varphi_0) e^{- i\vec{k}_+\cdot \vec{r} } } \end{array} \end{equation} \tag{ E10 }$$

where $\vec{k}_+$ is the measurement wavenumber in a particular scattering experiment, assuming constant in time and space in this work. The scattering volume shape U is centered around $(R_0, Z_0, \varphi_0)$, and has the spatial extent ∼(ΔR, ΔZ, Δϕ). We place ourselves at the outboard midplane and assume Z₀ ≈ 0. As in appendix D and in the main text, we choose a coordinate system in which the toroidal angle of scattering is set to ϕ₀ = 0 (figure D2). The quantities ΔR, ΔZ, RΔϕ are much smaller than minor radius a but much greater than the corresponding wavelengths associated to $\vec{k}_+$.

The product $\vec{k}_+\cdot\vec{r}$ written in cylindrical coordinates is

$$\begin{equation} \vec{k}_+\cdot\vec{r} = k_{x+} R \cos\varphi + k_{y+}R \sin\varphi + k_{Z+}Z. \end{equation} \tag{ E11 }$$

Expanding the density field as $\delta n(r,\theta, \varphi, t) = \sum_n{\delta n_{n}(r,\theta,t)} e^{-in(\varphi + \nu)} = \sum_{n,p}{\delta n_{np}(\theta,t)e^{-in(\varphi + \nu)}e^{i2\pi pr/L_r}}$ (ignore here Doppler shift $\overline{\omega}_0$) and plugging into equation (E10), we recover equations (4) and (5):

$$\begin{equation} \begin{aligned} &\delta{n}_{\textrm{u}} (\vec{k}_+, t) \approx \sum_{n, p} {U_{np} \delta n_{np}(\theta_0, t) }\\ &\text{where }U_{np} \text{ is defined by} \\ &U_{np} = \int {{\mathrm{d}}^3\vec{r} U(\vec{r})e^{-in\alpha}e^{i2\pi p r/L_r} e^{{-i\vec{k}_+}\cdot\vec{r} } }. \end{aligned} \label{} \end{equation}$$

We have assumed a slowly varying $\delta n_{np}(\theta_0,t)$ in θ. At the outboard midplane it is reasonable to assume a scattering volume envelope to be separable in (R, Z, ϕ), such as $U(R, Z, \varphi) \approx \Psi_R(R-R_0)\Psi_Z(Z-Z_0)\Psi_\varphi(\varphi)$, and assume these functions are well approximated by Gaussians (equation (6)). Now the scattering signal $\delta n_\textrm{u}$ can be written as separate integrals over R, Z, and ϕ as

$$\begin{equation} \begin{alignedat}{2} &\delta n_{\textrm{u}} = \int{{\mathrm{d}}R{\mathrm{d}}Z e^{-i Zk_{z+}}\Psi_R\Psi_Z \sum_n{e^{-i n\nu} \delta n_n(r,\theta,t) J_{n^\varphi}} } \\ &\text{where} \\ &J_{n^\varphi} = R\int{{\mathrm{d}}\varphi \Psi_\varphi e^{-i n\varphi} e^{-i(k_{x+}R\cos\varphi + k_{y+}R\sin\varphi )} } \\ \end{alignedat} \end{equation} \tag{ E13 }$$

where we used expression (E11) for $\vec{k}_+\cdot \vec{r}$. Since we assume a localized measurement in ϕ, the sin and cos terms in the exponent can be expanded about ϕ ≈ 0, leading to

$$\begin{equation} k_{x+} R \cos\varphi + k_{y+} R \sin\varphi \approx k_{x+} R(1-\varphi^2/2) + k_{y+}R\varphi. \end{equation} \tag{ E14 }$$

Using this expansion, the integral $J_{n^\varphi}$ can then be written as

$$\begin{equation} \begin{alignedat}{2} &J_{n^\varphi} &&\approx e^{-iRk_{x+} } \Phi_n(n{-}n^\varphi_+) \\ &\text{where} \\ &\Phi_n(n {-} n^\varphi_+) &&\approx R_0\int{{\mathrm{d}}\varphi \Psi_\varphi e^{i(\varphi^2/2)Rk_{x+}} e^{-i\varphi(n{-}n^\varphi_+)} } . \end{alignedat} \end{equation} \tag{ E15 }$$

We assumed R ≈ R₀ in Φ_n due to the slow spatial dependence on (r, θ), and we made use of the wavenumber mapping of equation (D9). The expression for Φ_n is nothing but a Gaussian integral, leading to

$$\begin{equation} \begin{alignedat}{2} &\Phi_n(n - n^\varphi_+) &&\approx R_0 \bigg( \frac{4\pi }{(\Delta n_\varphi)^2} \bigg)^{1/2} e^{-(n {-} n^\varphi_+)^2/(\Delta n_\varphi)^2} \\ &\text{with} \\ & (\Delta n_\varphi)^2 && = 4/\Delta \varphi^2 - 2iR_0k_{x+}. \end{alignedat} \end{equation} \tag{ E16 }$$

We recover a Gaussian shape for Φ_n. The toroidal mode number resolution $\Delta n_\varphi$ is complex in nature, and depends on the toroidal extent of the scattering volume Δϕ and the x component of the sampled wavenumber k_x+. The combination $\Delta \varphi^2R_0 k_{x+}$ is dependent on the specific scattering experiment and geometry and would have to be analyzed case by case. Now the scattering signal $\delta n_\textrm{u}$ can be written as

$$\begin{equation} \begin{aligned} \delta n_{\textrm{u}} &= \int{{\mathrm{d}}R{\mathrm{d}}Z e^{-i Zk_{z+}}e^{-i Rk_{x+}}\Psi_R\Psi_Z}\\ &{\sum_n{e^{-i n\nu} \delta n_n(r,\theta,t) \Phi_n(n-n^\varphi_+) } } . \end{aligned} \end{equation} \tag{ E17 }$$

The complex exponential part in $\delta n_\textrm{u}$ can also be expanded about $(r,\theta) \approx (r_0, \theta_0)$, leading to

$$\begin{equation} \begin{aligned} & Rk_{x+} + Zk_{z+} + n\nu \approx [\text{phase}]\\ & - r\bigg[\frac{2\pi }{L_r}(p - p_+) - (n - n_+^\theta)\frac{\partial \nu}{\partial r}\big|_0 \bigg] + \theta \bigg[(n - n_+^\theta)\frac{\partial \nu}{\partial \theta}\big|_0 \bigg]\\ & + n\bigg[\nu_0 -r_0 \frac{\partial \nu}{\partial r}\big|_0 - \theta_0 \frac{\partial \nu}{\partial \theta}\big|_0 \bigg]. \end{aligned} \end{equation} \tag{ E18 }$$

We can also expand (R, Z) about $(r_0, \theta_0)$ in the expressions for $\Psi_R, \Psi_Z$ (equation (6)), leading to

$$\begin{equation} \begin{aligned} U(R, Z, \varphi) &\approx \Psi_r(r-r_0)\Psi_\theta(\theta-\theta_0)\Psi_\varphi(\varphi) \\ \text{with}& \\ \Psi_r(r-r_0) &= \exp[-(r-r_0)^2/\Delta r^2] \\ \Psi_\theta(\theta-\theta_0) &= \exp[-(\theta-\theta_0)^2/\Delta \theta^2] \\ \Psi_\varphi(\varphi) &= \exp[-\varphi^2/\Delta \varphi^2] \\ \end{aligned} \end{equation} \tag{ E19 }$$

where $\Delta r = |\nabla r|_0 \Delta R$ and Δθ = ΔZ/(r₀κ). Next we can write the Jacobian J_r from $(R, Z, \varphi) \rightarrow (r, \theta, \varphi)$ approximately as $J_r(r,\theta) \approx R_0r_0\kappa/|\nabla r|_0$. Putting equations (E16), (E17), (E18) and (E19) together, we arrive to the final expression of the scattering matrix U_np

$$\begin{align} U_{np} & = e^{-in\alpha_0}e^{i2\pi p r_0/L_r} \Phi_n (n - n^\varphi_+) \Theta_n(n - n^\theta_+)\notag\\ \text{with}& \quad\Pi_p(p - p_+) \notag\\ \Phi_n(n - n^\varphi_+) &\approx R_0 \bigg( \frac{4\pi }{(\Delta n^\varphi )^2} \bigg)^{1/2} e^{-(n-n^\varphi_+)^2/(\Delta n^\varphi )^2} \qquad\text{3D} \notag\\ & \approx R_0\Delta \varphi \qquad\qquad\text{2D} \notag\\ \Theta_n(n-n^\theta_+) &\approx \sqrt{\pi}(\Delta Z) e^{-(n- n^\theta_+)^2/(\Delta n^\theta)^2} \notag\\ \Pi_{p}(p-p_{+}) &\approx\sqrt{\pi}(\Delta R) e^{-(p- p_{+})^2/\Delta p^2}\ \text{Global simulation} \notag\\ &\approx \frac{L_r}{|\nabla r|_0} {\mathrm{sinc}}[{\pi(p- p_{+})}]\ \text{Local simulation} . \end{align} \tag{ E20 }$$

Equation (E20) gives the general expression of the scattering matrix U_np in 3D and 2D, at the outboard midplane and assuming axisymmetry. The scattering signal can be computed similarly as before via $\delta{n}_{\textrm{u}} (\vec{k}_+, t) = \sum_{n,p} {U_{np} \delta n_{np}(\theta_0, t) }$.

The different radial mode number filter Π_p in a local versus global simulation stems from the fact that local simulation has the same profiles along the full radial domain. One can choose a radial filter in real space Ψ_r = 1, resulting in the sinc function. The toroidal, poloidal and radial mode number resolutions are

$$\begin{equation} \begin{aligned} (\Delta n^\varphi )^2 &= 4/\Delta \varphi^2 - 2iR_0k_{x+} \\ \Delta n^\theta &= \frac{2}{\Delta \theta \frac{\partial \alpha}{\partial \theta}|_0} \\ \Delta p &= \frac{L_r}{\pi|\nabla r|_0\Delta R}. \\[-3pt] \end{aligned} \end{equation} \tag{ E21 }$$

With respect to the toroidal filter, in the 2D approximation $\Delta \varphi \rightarrow 0$, or equivalently $\Psi_\varphi = \delta(\varphi/\Delta \varphi)$, is a delta function, and the toroidal mode number filter Φ_n is simply constant. Taking into account the toroidal variation of the scattering volume leads to a combination of two toroidal mode number filters Φ_n and Θ_n. The toroidal filter in real space $\Psi_\varphi$ gives rise to a toroidal mode number filter Φ_n about $n^{\varphi}_+$. The poloidal filter in real space (Ψ_Z or $\Psi_\theta$) gives rise to a toroidal mode number filter Θ_n about $n^{\theta}_+$. However, recall how $n^{\varphi}_+$ and $n^{\theta}_+$ are a priori independent, since they are separately computed from $(k_x, k_z)_+$ and k_y+, respectively. The question remains whether $n^{\theta}_+$ and $n^{\varphi}_+$ have similar values in actual scattering experiments. Assuming simple s−α geometry at the outboard midplane, equation (D9) simplifies to

$$\begin{equation} \left\{ \begin{alignedat}{1} &n^\theta_+ = \frac{\frac{\partial R}{\partial r} k_{x+} + \frac{\partial Z}{\partial r}k_{z+} }{-\frac{\partial \nu}{\partial \theta}} \bigg|_0 \rightarrow \frac{r_0}{q_0}k_{z+} \\ &n^\varphi_+ = -R_0k_{y+} \end{alignedat}. \right. \end{equation} \tag{ E22 }$$

Equation (E20) suggests the condition $n^{\varphi}_+ \approx n^\theta_+$ is needed for achieving a finite amplitude scattering signal. Using equation (E22), the condition $n^{\varphi}_+ \approx n^\theta_+$ translates to the condition $\vec{k}_+ \cdot \vec{B} \approx 0$, i.e. fluctuations are aligned with $\vec{B}$. The condition $\vec{k}_+ \cdot \vec{B} \approx 0$ is a necessary requirement for scattering experiments in magnetized plasmas, since fluctuations have maximal amplitude when aligned along the magnetic field. This information is directly encoded in the magnetic geometry along the field lines via the definition of ν and ultimately in the scattering matrix U_np. Successful, dedicated scattering experiments should be designed to satisfy the scattering condition $\vec{k}_+ \cdot \vec{B} \approx 0$, which is recovered here in the context of the computation of the scattering signal.

In this appendix we question whether toroidal effects are important for the high-k scattering system in NSTX. The scattering spectra presented in this article have been computed in the 2D approximation, i.e. setting $\Phi_n \approx R_0\Delta \varphi$. In fact, the synthetic diagnostic signal is scaled by R₀Δϕ. In this appendix we relax that restriction to understand if toroidal, 3D effects will have important consequences on the scattering signal.

As is suggested by equation (E20), the most important dependence on toroidal mode number n from the scattering matrix U_np comes from the filters Φ_n and Θ_n. We wish to understand how Φ_n and Θ_n compare to each other for different values of the toroidal scattering volume length R₀Δϕ. This is an unknown quantity so far, although good estimates could be found following the procedures outlined in [54, 55]. The 2D approximation can be recovered for $R_0\Delta \varphi \rightarrow 0$.

Figure F3 summarizes the preliminary analysis performed to assess the effect of the toroidal scattering length R₀Δϕ on the filters. These tests were performed for scattering wave-vector components corresponding to channel 1 of the scattering system and using realistic NSTX geometry from H-mode plasma 141 767. The vertical extent of the scattering volume is taken to be ΔZ = 0.03 m, consistent with the high-k scattering system.

Zoom In Zoom Out Reset image size

Figure F3(a) shows the toroidal mode number filter Φ_n for varying values of the toroidal scattering length R₀Δϕ [0.001−0.6] m. Small values of $R_0\Delta \varphi \lesssim 0.05$ m correspond to highly toroidally localized measurements. The toroidal mode number resolution $(\Delta n^\varphi)^2$ (equation (E21)) is dominated by the 4/Δϕ² contribution, and decreases with increasing R₀Δϕ from R₀Δϕ = 0.001 m to R₀Δϕ = 0.05 m. In this situation Φ_n is a real quantity. As the toroidal scattering length R₀Δϕ further increases for $R_0\Delta\varphi \gtrsim 0.05$ m, we are in the opposite situation and $(\Delta n^\varphi)^2 \approx -2iR_0k_{x+}$ is now a complex quantity. In this situation, an increasing toroidal scattering length for values larger than 0.05 m has the opposite effect of widening the resulting toroidal mode number width $\Delta n^\varphi$, increasing from R₀Δϕ = 0.05 m to R₀Δϕ = 0.6 m. Note how the narrowest width in this scan was found for a value of R₀Δϕ = 0.05 m. We now wish to compare the extent of these filters Φ_n to the toroidal mode number filter Θ_n.

Figure F3(b) shows an actual comparison of the toroidal mode number filters Φ_n (black) and Θ_n (red), defined in equation (E20), and corresponding to the narrowest toroidal width found in (a) for R₀Δϕ = 0.05 m. The toroidal mode number filter Φ_n is wider than the filter Θ_n, even in this 'worst-case' scenario in which $\Delta n^\varphi$ was narrowest. This suggests than one can assume Φ_n to be constant in the regime of variation of Θ_n. The overall product $\Phi_n \times \Theta_n$ is in orange. Comparing $\Phi_n \times \Theta_n$ (orange) and Θ_n (black), the main impact of Φ_n is simply a scaling factor. However, the overall width of the filter seems not to be much affected by the inclusion of Φ_n.

Figure F3(c) shows the overall product $\Phi_n \times \Theta_n$ for all the corresponding values R₀Δϕ performed in (a). As can be seen, varying R₀Δϕ has negligible impact on the overall toroidal mode number filter product $\Phi_n \times \Theta_n$. This can be understood from figures F3(a) and (b), since the filter Φ_n is wider than Θ_n for essentially all values of R₀Δϕ. Compared to Θ_n, Φ_n is well approximated by a constant and results in an excellent approximation $\Phi_n \times \Theta_n \approx \textrm{const.} \times \Theta_n$. This recovers the 2D formulation of the synthetic diagnostic.

Although only a qualitative assessment, we have not numerically implemented a 3D synthetic diagnostic and cannot compare actual frequency power spectra from 2D and 3D synthetic diagnostics. However, this preliminary assessment on the wavenumber filters suggests that the effect of the toroidal extent of the scattering volume is expected to be negligible for the high-k scattering system in NSTX, and the 2D approximation holds to a high degree of accuracy. A more detailed analysis comparing the actual frequency spectra from 2D and 3D synthetic diagnostics, as well as additional assessments to the one presented here for other coherent scattering experiments, might be the object of a future publication.

In this appendix we expand on the consequences of having off-midplane scattering locations on the wavenumber mapping. We also present intuitive pictures behind the effect the normalizing magnetic field, Shafranov shift and elongation affecting the wavenumber mapping $(k_x, k_y, k_z)_+ \rightarrow (k_r, k_\theta)_+$.

Figure G4 shows an example mapping of equation (8), corresponding to a pair $(k_x, k_z)_+ = (-500, 1)$ m⁻¹ when the mapping is computed at different poloidal locations along the flux surface, and neglecting the influence of the toroidal component (2D approximation). A circular scattering volume cross-section in U at a fixed toroidal slice gives rise to a circularly shaped filter in $(k_x, k_z)$ as shown in G4(a). This is characteristic of the high-k scattering diagnostic in NSTX [26]. Figure G4(b) shows how a circularly shaped filter in $(k_x, k_z)$ is mapped to an elliptical shape in $(k_r, k_\theta)$, depending on the flux-surface location where scattering takes place. The ellipse is elongated, and becomes slanted for off-midplane locations. Even a small k_z + = 1 m⁻¹ value maps to finite $k_{\theta+}\rho_s$, which can itself be positive or negative depending on the location of scattering along the flux surface. This means that the same measurement of $(k_x, k_z)_+$ at different poloidal locations will be sampling different turbulent wavenumbers when expressed in $(k_r\rho_s, k_\theta\rho_s)$, which are the properly normalized wavenumber components characterizing microturbulence fluctuations. This will result into different scattering amplitudes at each poloidal location. The particular case of outboard midplane location is discussed in subsection 3.3. The flux-surface geometry is taken from NSTX H-mode plasma discharge 141 767.

Zoom In Zoom Out Reset image size

Within the 2D approximation and at the outboard midplane (equation (11)), one can also easily build an intuitive picture behind the effect of the normalizing magnetic field, Shafranov shift and elongation on the wavenumber mapping. We recall the 2D approximation of the wavenumber mapping here once more

$$\begin{equation} \left\{ \begin{aligned} (k_{r+}\rho_s)^\textrm{sim} &\approx \frac{k_{x+}}{|\nabla r|_0} \rho_s^\textrm{sim} \\ (k_{\theta+}\rho_s )^\textrm{sim} &\approx - \frac{\kappa q}{\frac{\partial \nu}{\partial \theta}}\bigg|_0 k_{z+} \rho_s^\textrm{sim} \\ \end{aligned}. \right. \end{equation} \tag{ G23 }$$

Experiments tend to use local values of the electron temperature T_e and magnetic field B entering the normalizing ρ_s, while gyrokinetic codes tend to employ internal definitions. In expression (G23), the subscript $(.)^{\textrm{sim}}$ means the values of the mapped wavenumbers have been properly normalized by ρ_s, using T_e and B consistent with the normalizations in gyrokinetic codes. In this work we used the GYRO $B_\textrm{unit}$ normalizing magnetic field [50] in $\rho_s^{\text{sim}}$. Not using the same normalizing magnetic field in experiments as in gyrokinetic codes will lead to a systematic error in the interpreted measured wavenumber components $(k_r, k_\theta)$ by a scattering experiment.

Shafranov shift Δ primarily affects k_r due to the compression of flux-surfaces at the outboard midplane through $|\nabla r{|_0} \approx 1/(1 + \Delta) > 1$ (Δ is negative at the outboard midplane, but positive at the inboard midplane). For a given radial perturbation of wavelength $\lambda_x = 2\pi/k_x$ at the outboard midplane given by experiment, the strong radial compression due to Shafranov shift will mean a larger λ_r in the field-line frame than in the absence of compression (where $|\nabla r| = 1, \Delta = 0$). Note how the same perturbation of wavelength λ_x in the lab frame fits in the radial domain in the absence of Shafranov shift in figure G5(a), but not with strong Shafranov shift. This will translate to a smaller k_r component with respect to the absence of compression, resulting in $k_r \approx k_x/|\nabla r|_0$ in equation (G23). The opposite will take place at the inboard midplane.

Zoom In Zoom Out Reset image size

With respect to elongation, it primarily affects $k_\theta$ due to flux-surface 'stretching' poloidally. For a given vertical perturbation of wavelength $\lambda_z = 2\pi/k_z$ given by experiment, a strongly elongated plasma will have a smaller $\lambda_\theta$ in the field-line frame than in the absence of elongation (κ = 1). This will translate to larger $k_\theta$ component when compared to the absence of elongation, resulting in $k_\theta \propto \kappa k_z$ in equation (G23). Note how the same perturbation of wavelength λ_z in the lab frame has a smaller poloidal extent (smaller $\lambda_\theta$) in the presence of elongation in figure G5(b) than in the absence of elongation. The factor $q/(\frac{\partial \nu}{\partial \theta})$ in equation (G23) appears due to the definition of $k_\theta$ as a flux function, and not a local quantity in θ.