2.1. Ray tracing and Gaussian beam tracing

This section closely follows notation introduced by Hall-Chen et al. (Reference Hall-Chen, Parra and Hillesheim2022 Reference Hall-Chen, Parra and Hillesheimb ), who rederived Gaussian beam tracing from first principles and presented it in full detail.

The beam-tracing method (or paraxial WKB) is an asymptotic representation of the electric field in anisotropic media for finite size wave packets with characteristic scale in the transverse direction of propagation $\sim W$ . Beam tracing extends ray tracing by including diffraction effects in a set of ODEs. As in ray tracing, the wavevector of the wave $\mathbf {K}$ is ordered inversely proportional to the wavelength $\sim 1/\lambda$ , while the background equilibrium is assumed to vary on a scale $L \gg \lambda$ . In beam tracing, the scale $W$ is intermediate and obeys $L \gg W \gg \lambda$ , with $W \sim (\lambda L)^{\frac {1}{2}}$ . In this manuscript, we denote the wavevector associated with the wave with the capital letter $\mathbf {K}$ . The same notation will apply to its components $K_g$ , $K_x$ , etc., while the wavevector components associated with the turbulence will be in lower case, e.g. $k_x$ and $k_y$ .

In Gaussian beam tracing, the variation of the wave electric field perpendicular to a central ray is assumed to be of Gaussian shape. The central ray propagation obeys the traditional ray-tracing equations for the ray position vector $\mathbf {q}(\tau ) \sim L$ and wavevector $\mathbf {K}(\tau ) \sim 1/\lambda$ , where $\tau$ is a parameter along the central ray. For reference, the ray-tracing equations are

(2.1) \begin{align} \begin{alignedat}{2} & \frac {\text {d}\mathbf {q}}{\text {d}\tau } = \nabla _K H, \\ & \frac {\text {d}\mathbf {K}}{\text {d}\tau } = - \nabla H, \end{alignedat} \end{align}

where $H$ is the cold-plasma dispersion relation and $\mathbf {g} = \nabla _K H$ is the effective group velocity. Note that there are infinite ways to define an $H$ that give rise to the same $\mathbf {q}$ and $\mathbf {K}$ (the different definitions of $H$ only change the dummy parameter $\tau$ ). We use a specific definition that ensures that the wave electric field scales as $\sim 1/g^{\frac {1}{2}}$ , where $\mathbf {g} = \text {d} \mathbf {q}/\text {d}\tau$ (see discussion in § 2.2 of Hall-Chen et al. (Reference Hall-Chen, Parra and Hillesheim2022 Reference Hall-Chen, Parra and Hillesheimb ) for details). Equations (2.1) are routinely solved by codes such as GENRAY (Smirnov et al. Reference Smirnov2009).

In the Gaussian beam-tracing approximation, the beam electric field $\mathbf {E}_b$ is written in amplitude and phase as $\mathbf {E}_b(\mathbf {r}) = \mathbf {A}(\mathbf {r}) \exp [i \psi (\mathbf {r})] +$ c.c., where the phase $\psi$ is expanded to second order about the central ray in the perpendicular direction to $\mathbf {g}$ ,

(2.2) \begin{align} \begin{alignedat}{2} \psi (\mathbf {r}) = s(\tau ) + \mathbf {K}_w \cdot \mathbf {w} + \frac {1}{2} \mathbf {w} \cdot \boldsymbol{\unicode{x1D6F9}\!}_w \cdot \mathbf {w} . \end{alignedat} \end{align}

Here, $s = s(\tau ) = \int ^\tau {\mathbf {K}(\tau ^{\prime}) \cdot \mathbf {g} (\tau ^{\prime})} \ \text {d}\tau ^{\prime} = \int ^\tau {{K}_g g } \ \text {d}\tau ^{\prime}$ is the large phase Eikonal term measuring the variation of the electric field along the central ray. The phase given by $s$ scales as $s \sim L/\lambda \gg 1$ , which follows from $\tau \sim L/\lambda$ , $g \sim \lambda$ and $H \sim 1$ (see (2.1)). Note that the second term in (2.2) is ordered as $\mathbf {K}_w \cdot \mathbf {w} \sim W/\lambda \gg 1$ , while the last term is $(1/2) \mathbf {w} \cdot \boldsymbol{\unicode{x1D6F9}\!}_w \cdot \mathbf {w} \sim 1$ , since the components of $\boldsymbol{\unicode{x1D6F9}\!}_w$ are ordered as $\Psi _{w,ij} \sim 1/W^2$ . The vector $\mathbf {K}_w = (\mathbf {I} - {\hat{\mathbf{g}}} {\hat{\mathbf{g}}} ) \cdot \mathbf {K}$ in the phase (2.2) is the component of $\mathbf {K}$ perpendicular to ${\hat{\mathbf{g}}}$ . The matrix $\boldsymbol{\unicode{x1D6F9}\!}_w$ is a singular $3 \times 3$ matrix and only has non-zero components along the two directions perpendicular to ${\hat{\mathbf{g}}}$ . The matrix $\boldsymbol{\unicode{x1D6F9}\!}_w$ contains information about the phase-front curvature of the Gaussian beam (through the non-zero eigenvalues of its real part = ${K^3}/{(K_g^2 R_{b,i})}$ , where $K = |\mathbf {K}|$ and $K_g = \mathbf {K} \cdot {\hat{\mathbf{g}}}$ ) and beam width $W_i$ (through the non-zero eigenvalues of its imaginary part $ = {2}/{W_i^2}$ ).

Figure 1. Lab frame and beam frame for coordinates with oblique beam incidence (finite $\alpha _0$ ) used throughout this manuscript. Here, $(x_c, y_c)$ denote the coordinates of a point along the central ray and not the cutoff location.

A general point in space is described by

(2.3) \begin{align} \begin{alignedat}{2} \mathbf {r} = \mathbf {q}(\tau ) + \mathbf {w} = \mathbf {q}(\tau ) + X {\hat{\mathbf{X}}}(\tau ) + Y {\hat{\mathbf{Y}}}(\tau ) . \end{alignedat} \end{align}

We define the (orthonormal) beam-frame coordinate system $\{ {\hat{\mathbf{Y}}}, {\hat{\mathbf{g}}}, {\hat{\mathbf{X}}} \}$ as follows:

(2.4) \begin{align} \begin{alignedat}{2} & {\hat{\mathbf{Y}}} = \frac {{\hat{\mathbf{b}}} \times {\hat{\mathbf{g}}}}{| {\hat{\mathbf{b}}} \times {\hat{\mathbf{g}}} |}, \qquad {\hat{\mathbf{g}}} = \frac {\mathbf {g}}{g}, \qquad {\hat{\mathbf{X}}} = \frac {{\hat{\mathbf{Y}}} \times {\hat{\mathbf{g}}}}{| {\hat{\mathbf{Y}}} \times {\hat{\mathbf{g}}} |}, \end{alignedat} \end{align}

where ${\hat{\mathbf{b}}} = \mathbf {B}/B$ is the unit vector along the background magnetic field. This is not to be confused with the Cartesian lab-frame system of coordinates $ \{ {\hat{\mathbf{x}}}, {\hat{\mathbf{y}}}, {\hat{\mathbf{z}}} \}$ (see figure 1 for a comparison between beam frame and lab frame in two dimensions). Using (2.3) and (2.4), a general point in space can be determined from the coordinate $\tau$ along the central ray and the two coordinates $X$ and $Y$ in the directions perpendicular to ${\hat{\mathbf{g}}}$ .

In the beam-tracing formulation which we use in this manuscript, the beam-tracing matrix $\boldsymbol{\unicode{x1D6F9}}$ is a $3 \times 3$ matrix that is convenient to evolve along the central ray. From this matrix, we can obtain $\boldsymbol{\unicode{x1D6F9}\!}_w$ by projection $\boldsymbol{\unicode{x1D6F9}\!}_w = (\mathbf {I} - {\hat{\mathbf{g}}} {\hat{\mathbf{g}}} ) \cdot \boldsymbol{\unicode{x1D6F9}} \cdot (\mathbf {I} - {\hat{\mathbf{g}}} {\hat{\mathbf{g}}} )$ . The beam matrix $\boldsymbol{\unicode{x1D6F9}}$ follows the beam-tracing evolution equations, a set of ordinary differential equations parametrised by $\tau$ , see Hall-Chen et al. (Reference Hall-Chen, Parra and Hillesheim2022 Reference Hall-Chen, Parra and Hillesheimb ) and Appendix A. Given $\boldsymbol{\unicode{x1D6F9}}$ , the electric field $\mathbf {E}_b$ in Gaussian beam tracing is written as

(2.5) \begin{align} \begin{alignedat}{2} & \mathbf {E}_b = &&A_{ant}\exp [{i(\phi _G+\phi _P)}] \bigg [ \frac { \det (\Im [\boldsymbol{\unicode{x1D6F9}\!}_w]) }{ \det (\Im [\boldsymbol{\unicode{x1D6F9}}_{w,ant}]) } \bigg ]^{\frac {1}{4}} \bigg ( {\frac {g_{ant} }{g}} \bigg )^{\frac {1}{2}} \\ & && \times \exp \Big ( is + i \mathbf {K}_w\cdot \mathbf {w} + \frac {i}{2} \mathbf {w} \cdot \boldsymbol{\unicode{x1D6F9}\!}_w \cdot \mathbf {w} \Big ) {\hat{\mathbf{e}}} \ + \ c.c., \end{alignedat} \end{align}

where $\mathbf {E}_b = E_b {\hat{\mathbf{e}}}$ and $\Im$ denotes the imaginary part. The unit vector ${\hat{\mathbf{e}}}$ is the polarisation vector, $\phi _P$ and $\phi _G$ are respectively the polarisation and Gouy phases, which are of limited interest in this manuscript (see Hall-Chen et al. (Reference Hall-Chen, Parra and Hillesheim2022 Reference Hall-Chen, Parra and Hillesheimb ) for their evolution equations). The subscript $(.)_{ant}$ means that quantities are evaluated at the antenna launch location. As previously discussed, the prescription in (2.5) for the electric field shows that $E_b \sim 1/(gW_X W_Y)^{\frac {1}{2}} \propto 1/(KW_X W_Y)^{\frac {1}{2}}$ . This is due to the terms $\det (\Im [\boldsymbol{\unicode{x1D6F9}\!}_w])^{\frac {1}{4}} \propto 1/(W_X W_Y)^{\frac {1}{2}}$ and $1/g^{\frac {1}{2}} \propto 1/K^{\frac {1}{2}}$ . For the purposes of this manuscript, in the 2-D linear layer, only the width $W_Y$ in the direction $\hat {\mathbf {Y}}$ is relevant.

In the next section, we show numerical solutions to the beam-tracing equations (Appendix A) in real tokamak geometry, which exhibit the phenomenon of beam focusing. These will be compared with analytical solutions to beam tracing for the 2-D linear-layer problem. We will see how the analytical solution qualitatively captures the phenomenon of beam focusing.

Figure 2. (a) Trajectory of the central ray of a DBS beam projected on the poloidal $(R,Z)$ plane overlaid by contour lines of the poloidal flux function $\Psi _p$ for JET discharge 97 080 (NBI-heated L-mode). (b) Numerical solution of the two principal widths perpendicular to the central ray propagation, noted here $W_1$ and $W_2$ , using Scotty (Hall-Chen et al. Reference Hall-Chen and Hill2025). The black dots correspond to the turning point in the trajectory (vanishing wavenumber component normal to the flux surface), while the red dots correspond to the plasma exit. (c) Analytic solution of the beam-tracing equations for the 2-D linear layer in slab geometry using experimental parameters corresponding to the case shown in panels (a) and (b): $K_0L\approx 1600, K_0W_0\approx 22$ , $\alpha _0 \approx 10^{\circ }$ , $R_{Y0}=\infty$ (launch at the waist). There is good qualitative agreement between the numerical and analytic solutions.

2.2. Beam focusing in experimentally relevant conditions

In this section, we will show that the beam-focusing phenomenon appears in numerical calculations of beam-tracing modelling of DBS experiments. We will see that the beam width has a tendency to focus as the beam moves towards higher density. A beam propagating in increasing density is routinely encountered in most DBS measurements.

We describe the phenomenon of beam focusing, that is, the decrease or compression of the beam width as the beam propagates near the turning point that characterises a cutoff. This phenomenon was previously observed in numerical simulations (Poli et al. Reference Poli, Pereverzev and Peeters1999, Reference Poli, Pereverzev, Peeters and Bornatici2001 Reference Poli, Pereverzev, Peeters and Bornaticic ; Yu. A. Kravtsov and Berczynski P. Reference Kravtsov and Berczynski2007; Conway et al. Reference Conway, Troester, Schirmer, Angioni, Holzhauer, Jenko, Merz, Poli, Scott and Suttrop2007, Reference Conway, Lechte and Fochi2015, Reference Conway, Lechte, Poli and Maj2019) and in analytic solutions (Maj et al. Reference Maj, Pereverzev and Poli2009, Reference Maj, Balakin and Poli2010; Weber et al. Reference Weber, Maj and Poli2018). To confirm the existence of the beam-focusing phenomenon and its relevance to real-life experimental conditions, in figure 2, we show the result of a numerical calculation with the Scotty code (Hall-Chen et al. Reference Hall-Chen and Hill2025) for the DBS beam propagation with X-mode polarisation in the core of the tokamak JET for L-mode discharge 97 080 in a phase heated only by neutral-beam injection (NBI). Figure 2(a) shows the trajectory of the central ray projected on the 2-D $(R,Z)$ plane overlaid by contour lines of the poloidal flux function $\Psi _p$ . The turning point is at $\Psi _p/\Psi _{p0} \approx 0.7$ , where $\Psi _{p0}$ is the separatrix value. Figure 2(b) shows the perpendicular widths $W_1$ and $W_2$ that correspond to the non-zero eigenvalues of the imaginary part of the beam-tracing matrix $\Im [ \boldsymbol{\unicode{x1D6F9}\!}_w ]$ . We approximate the overmoded waveguide launch by a beam with circular cross-section and initial conditions $W_0 \approx 1$ cm, $R_{Y0} = \infty$ (launch at the waist). The widths $W_1$ and $W_2$ are plotted against the path length $l = \int \text {d}\tau g$ . The beam propagates in vacuum from $l \approx 0$ to $l \approx 0.4$ m, which roughly corresponds to the plasma edge. In vacuum, both principal beam widths follow the same behaviour: expansion following a launch from the waist. In the plasma, both widths initially continue expanding, but one of them (blue) starts contracting and reaching a minimum at $l \approx 0.7$ m (beam focusing). Following the beam focusing, the contracted width starts to expand again. The behaviour observed in figure 2(b) consistently appears in beam-tracing numerical calculations near a turning point, and is confirmed by the beam-tracing code TORBEAM (Poli et al. Reference Poli, Peeters and Pereverzev2001 Reference Poli, Peeters and Pereverzeva , Reference Poli2018) (not shown). This behaviour was also observed in simulations of DBS in other tokamaks (Conway et al. Reference Conway, Troester, Schirmer, Angioni, Holzhauer, Jenko, Merz, Poli, Scott and Suttrop2007, Reference Conway, Lechte and Fochi2015, Reference Conway, Lechte, Poli and Maj2019). This motivates studying the beam-focusing phenomenon using a reduced model that is analytically tractable, the 2-D linear-layer model (Maj et al. Reference Maj, Pereverzev and Poli2009, Reference Maj, Balakin and Poli2010), which we present next.

2.3. Beam focusing in the 2-D linear layer

To understand the beam-focusing phenomenon observed in the Scotty numerical solution in toroidal geometry in figure 2 (b), we first solve the equations for ray tracing (2.1) using a simplified Cartesian-slab geometry in two dimensions (Maj et al. Reference Maj, Pereverzev and Poli2009, Reference Maj, Balakin and Poli2010). Following the ray-tracing solution, we will subsequently solve the beam-tracing equations.

We solve the 2-D linear-layer model in the lab frame, given by the unit vectors $\{ {\hat{\mathbf{x}}}, {\hat{\mathbf{y}}}, {\hat{\mathbf{z}}} \}$ (figure 1). Using (2.4), the lab-frame unit vectors are related to the beam-frame $\{ {\hat{\mathbf{Y}}}, {\hat{\mathbf{g}}}, {\hat{\mathbf{X}}} \}$ in two dimensions by a simple rotation

(2.6) \begin{align} \begin{alignedat}{2} &\hat {\mathbf {Y}} &&= \sin \alpha \ \hat {\mathbf {x}} - \cos \alpha \ \hat {\mathbf {y}} , \\ & \hat {\mathbf {g}} &&= \cos \alpha \ \hat {\mathbf {x}} + \sin \alpha \ \hat {\mathbf {y}} , \\ &\hat {\mathbf {X}} &&= \hat {\mathbf {z}}, \end{alignedat} \end{align}

where we note that $\alpha = \alpha (\tau )$ is a function of $\tau$ . We will ignore the coordinate ${\hat{\mathbf{z}}}$ , which points in the antiparallel direction to the magnetic field. These definitions will be needed in order to subsequently solve the beam-tracing equations. The 2-D linear-layer problem in slab geometry has O-mode polarisation, uniform $\mathbf {B}(\mathbf {r}) = {B}_0 {\hat{\mathbf{b}}} = - B_0 {\hat{\mathbf{z}}}$ and linear density profile $\omega _{pe}^2(x) = \Omega ^2 x/L$ . Here, $x$ is our ’radial’ coordinate, also Cartesian $x$ in figure 1, $\Omega$ is the launch frequency, $\omega _{pe0} = \Omega \cos \alpha _0$ is the electron plasma frequency at the turning point location, $L_n = L \cos ^2\alpha _0$ is the turning point locationFootnote ⁴ and $L$ is the turning point location for zero incident angle ( $\alpha _0=0$ , see figure 1). The ray-tracing equations determine the trajectory of the central ray $\mathbf {q}$ and the wavevector of the central ray $\mathbf {K}$ . In the 2-D linear layer, $\mathbf {q}$ and $\mathbf {K}$ have components

(2.7) \begin{align} \begin{alignedat}{2} & \mathbf {q}(\tau ^{\prime}) &&= \big (x_c(\tau ^{\prime}), y_c (\tau ^{\prime}) \big ) , \\ &\mathbf {K}(\tau ^{\prime}) &&= \big (K_x(\tau ^{\prime}), K_y(\tau ^{\prime})\big ) , \end{alignedat} \end{align}

where $\tau ^{\prime} = \tau /K_0L$ is the normalised parameter along the path and $K_0$ is the wavenumber magnitude at launch. Here, $(x_c, y_c)$ denote the coordinates of an arbitrary point along the central ray and not the cutoff location. In this manuscript, we define the turning point, or cutoff, as the location where $K_x=0$ ( $K_x$ is the Cartesian $x$ -component of the beam wavevector ${\textbf{K}}$ , see (2.9)). We will make use of the dispersion relation $H=0$ , where $H$ takes the form

(2.8) \begin{align} \begin{aligned} H = \frac {K^2}{K_0^2} - 1 + \frac {\omega _{pe}^2}{\Omega ^2} \end{aligned} \end{align}

for the O-mode. The form of $H$ given by (2.8) ensures that the electric field is given by (2.5). Using (2.8), we have $\nabla _K H = \mathbf {g} = 2 \mathbf {K}/K_0^2$ (initial condition $g_{ant} = 2/K_0$ ) and $\nabla H = {\hat{\mathbf{x}}}/L$ . As can be seen from figure 1, the magnetic field is into the page and the angle $\alpha _0$ measures the vertical incidence angle of the central ray at launch. In what follows, $\alpha$ will measure the vertical incidence angle along the central ray. The initial conditions are $(x_{c}, y_{c})_0 = (0, y_0)$ , $(K_{x}, K_{y})_0 = K_0 ( \cos \alpha _0, \sin \alpha _0)$ . The ray-tracing (2.1) can be readily solved in the 2-D linear layer, giving the following components of $\mathbf {q}$ and $\mathbf {K}$ in (2.7):

(2.9) \begin{align} \begin{aligned} & x_c(\tau ^{\prime}) = L \Big ( \cos ^2\alpha _0 - (\cos \alpha _0- {\tau ^{\prime}} )^2 \Big ) = L \big ( \cos ^2\alpha _0 - {K_x^2}/{K_0^2}\big ) , \\ & y_c(\tau ^{\prime}) = y_0 + 2 L \sin \alpha _0 \ \tau ^{\prime} , \\ & K_x(\tau ^{\prime}) = K_0 (\cos \alpha _0 - \tau ^{\prime}) = K_0 (\cos ^2\alpha _0 - x_c/L)^{\frac {1}{2}} , \\ & K_y(\tau ^{\prime}) = K_0\sin \alpha _0. & \end{aligned} \end{align}

Equations (2.9) describe a parabolic trajectory in $(x,y)$ . We will find it useful to parametrise the ray trajectory by the radial wavenumber of the central ray $K_x$ , instead of $\tau$ . We will also parametrise the beam-tracing solution by $K_x$ .

The 2-D linear-layer problem with finite incident angle $\alpha _0$ is also analytically solvable in Gaussian beam tracing, as was previously shown by Maj et al. (Reference Maj, Pereverzev and Poli2009, Reference Maj, Balakin and Poli2010). For details of the derivation, we refer the reader to Appendix A. The beam-tracing equations for the linear-layer problem reduce to ${\text {d} \boldsymbol{\unicode{x1D6F9}}}/{ \text {d} \tau } = -({2}/{K_0^2}) \boldsymbol{\unicode{x1D6F9}}^2$ , which is a nonlinear matrix equation for $\boldsymbol{\unicode{x1D6F9}}$ . This equation belongs to a class of ODEs known as Ricatti equations. The main component of the beam matrix $\boldsymbol{\unicode{x1D6F9}}$ of interest in this manuscript is the $\boldsymbol{\hat {Y}} \boldsymbol{\hat {Y}}$ component $\Psi _{YY}$ , which captures the beam-focusing phenomenon through the perpendicular beam width $W_Y = (2/\Im [\Psi _{YY}])^{\frac {1}{2}}$ . The solution for $\Psi _{YY}$ is

(2.10) \begin{align} \begin{aligned} & \Psi _{YY}(K_x) &&= \frac { - \frac {1}{2} \sin ^2\alpha _0 + \Psi _{yy0}^{\prime} \Big [ \frac {K_x^3}{K_0^3} + 3 \sin ^2\alpha _0 \frac {K_x}{K_0} - \frac {\sin ^2\alpha _0}{\cos \alpha _0}\big(\cos ^2\alpha _0 - \sin ^2\alpha _0 \big) \Big ] }{ \Big [ \sin ^2\alpha _0 + \frac {K_x^2}{K_0^2} \Big ] \Big [ \frac {K_x}{K_0} + 2 \Psi _{yy0}^{\prime} \Big ( \frac {\cos ^2\alpha _0 - \sin ^2\alpha _0 }{\cos \alpha _0} \frac {K_x}{K_0} + \sin ^2\alpha _0 - \frac {K_x^2}{K_0^2} \Big ) \Big ] } \frac {K_0}{L} , \end{aligned} \end{align}

where $\Psi _{yy0}^{\prime} = \Psi _{yy0} ({L}/{K_0})$ is the initial condition for the beam $\Psi _{yy}$ component in the lab-frame direction ${\hat{\mathbf{y}}}$ (vertical direction in figure 1) that contains information about the initial beam width and radius of curvature of the phase front (normalised to $K_0/L$ , see Appendix A for details). Equation (2.10) is the basis for most of the analysis in this manuscript. The reader is referred to Appendix A for details in the calculation of $\Psi _{YY}$ . The analytic solution of $\Psi _{YY}$ contains the beam-focusing phenomenon for the 2-D linear-layer problem and is used to characterise the Doppler backscattering power contributions along the beam trajectory through the filter function $|F_{x\mu }|^2$ (§ 3).

We wish to compare the analytic beam-tracing solution for $W_{Y}$ , computed using (2.10), with the numerical solution in toroidal geometry shown in figures 2(a) and 2(b). To do so, we extract the physical values of $K_0$ , $L$ and the initial width $W_0$ from experimental conditions in the JET discharge 97080: $K_0 \approx 2200$ m $^{-1}$ , $W_0 \approx 1$ cm and $L\approx 0.7$ m. Knowing that, in reality, the beam propagates in a plasma with a varying density gradient, we choose the value of $L$ to be the average value of the density gradient experienced by the beam as it propagates through the plasma, yielding $L \approx 0.7$ m. We use the experimental values of $K_0$ , $L$ and $W_0$ to normalise the initial conditions in the analytic beam-tracing solution, giving $K_0W_0\approx 22$ , $K_0 L \approx 1600$ . The width $W_Y$ that is perpendicular to the central ray propagation and in the $(x,y)$ plane is plotted in blue in figure 2(c), as a function of the beam path length $l$ . The vacuum solution is plotted in orange (absence of plasma),Footnote ⁵ which is also the same solution as the $\hat {\mathbf {z}}$ -component of the width for the linear layer in the $\hat {\mathbf {z}}$ -direction $W_z$ . Despite the limitations of the model (linear density profile, 2-D slab) and the difference in polarisation (X-mode versus O-mode), a qualitative comparison between figures 2(c) and 2(b) clearly shows that the analytic solution for the 2-D linear layer is able to recover the focusing behaviour observed in the numerical solution in toroidal geometry. This confirms that the beam focusing is not a numerical artefact. The beam focusing behaviour in the vicinity of the turning point is physical within the beam-tracing approximationFootnote ⁶ and motivates studying the beam-focusing phenomenon using the analytic solution to the beam-tracing equations in the 2-D linear-layer problem. In what follows, we describe the dependence of the analytic beam-tracing solution for $\Psi _{YY}$ on the initial conditions for the width $W_0$ , the radius of curvature $R_{Y0}$ and the vertical launch angle $\alpha _0$ . The initial condition is at $\tau =0$ , which corresponds to the plasma edge $x=0$ .

Figure 3. Thick coloured curves indicate values of the width $W_Y$ for varying initial incident angle $\alpha _0 = 10^{\circ }, 30^{\circ }, 50^{\circ }$ and $80^{\circ }$ , and different initial $R_{Y0}$ (see legend) and fixed initial $W_0 = 0.40(\lambda L)^{\frac {1}{2}}$ . Corresponding dashed lines of same colour indicate vacuum values of $W_Y$ for same initial conditions as the thick coloured curves. The initial value of $R_{Y0}/L=-2 {\cos \alpha _0}/{\sin ^2\alpha _0}$ (green) corresponds to the particular initial condition employed by Gusakov et al. (Reference Gusakov, Irzak and Popov2014, Reference Gusakov, Irzak, Popov, Khitrov and Teplova2017). The vertical dotted points in grey indicate the location of the turning point, or cutoff (location where $K_x=0$ ).

We start by considering the analytic beam-tracing solution for $\Psi _{YY} = {K}/{R_Y} + i {2}/{W_Y^2}$ (Appendix A). In this beam-tracing solution, the dimensional quantities $K_0 = 2\pi /\lambda$ , $W_0$ and $L$ enter the equations through the beam parameter $W_0/(\lambda L)^{\frac {1}{2}}$ . Therefore, in the rest of the manuscript, only the beam parameter $W_0/(\lambda L)^{\frac {1}{2}}$ will be specified since it is the only parameter needed to determine the normalised solution $\Psi _{YY} ({L}/{K_0})$ . The beam width $W_Y$ is shown in figure 3 as a function of $\tau$ for different initial conditions of the incident angle $\alpha _0 = 10^{\circ }, 30^{\circ }, 50^{\circ }$ and $80^{\circ }$ , and the initial radius of curvature $R_{Y0}$ (coloured curves). In this slab model, we assume that plasma is only present for $x\gt 0$ and the vacuum exists for $x\lt 0$ . Therefore, one needs to calculate the range of $\tau$ for which the central ray trajectory has positive $x_c$ . According to (2.9), for initial conditions at $\tau ^{\prime}=0$ , the plasma exit corresponds to $\tau ^{\prime} = 2\cos \alpha _0$ . The initial width $W_0$ is kept fixed at $W_0 = 0.40(\lambda L)^{\frac {1}{2}}$ . The coloured dashed lines are the vacuum solutions (absence of plasma) with initial conditions corresponding to the respective $W_Y$ of the same colour in the 2-D linear-layer problem. The vacuum solution initially follows the plasma solution until the beam approaches the turning point. The initial value of $R_{Y0}/L=-2 {\cos \alpha _0}/{\sin ^2\alpha _0}$ corresponds to the particular initial condition employed by Gusakov et al. (Reference Gusakov, Irzak and Popov2014, Reference Gusakov, Irzak, Popov, Khitrov and Teplova2017) (see discussion in Appendix A). For that particular case, the initial radius of curvature takes the values $R_{Y0}/L = -65.32, -6.93, -2.19$ and $-0.3581$ respectively for $\alpha _0=10^{\circ }, 30^{\circ }, 50^{\circ }$ and $80^{\circ }$ . Here, $R_{Y0}$ has little effect on the beam focusing around $\tau ^{\prime} \approx 1$ for small angles, but has a non-negligible effect in the vicinity of the initial launch and after the turning point for the outgoing beam. These differences appear for small initial $|R_{Y0}|/L \lesssim 1$ , that is, for strongly focusing or diverging beams (see red and magenta curves in figure 3). The initial radius of curvature $R_{Y0}$ becomes more important in the vicinity of the turning point for increasing incident angles, where the beam focusing location is shifted (see $\alpha _0 = 30^{\circ }$ ). For $\alpha _0 = 50^{\circ }, 80^{\circ }$ , the beam focusing does not take place inside the plasma. The only focusing region corresponds to the beam waist that is captured by the vacuum solution. Note how the initial condition of $R_{Y0}/L = -2 {\cos \alpha _0}/{\sin ^2\alpha _0}$ from Gusakov, Irzak & Popov (Reference Gusakov, Irzak and Popov2014) and Gusakov et al. (Reference Gusakov, Irzak, Popov, Khitrov and Teplova2017) strongly depends on the launch angle $\alpha _0$ : it defines a beam that is closer to a launch at the waist ( $R_{Y0}/L = \infty$ ) for small incident angles, but becomes strongly focusing (smaller, negative $R_{Y0}$ ) for larger $\alpha _0$ .

Figure 4. Similar to figure 3 for varying initial width $W_0 / (\lambda L)^{\frac {1}{2}} = 0.13, 0.40, 1.26$ and fixed $\alpha _0=30^{\circ }$ . For small $W_0$ , the beam is initially strongly focused and experiences a large expansion (due to diffraction) that follows the vacuum solution before focusing slightly past the turning point (for all of $R_{Y0}$ ). For increasing $W_0$ , the initial growth is less severe and the beam focusing depends on the initial $R_{Y0}$ . The beam even focuses twice along the path for particular initial conditions $R_{Y0}/L=-0.5, \alpha _0=30^{\circ }$ , $W_0 / (\lambda L)^{\frac {1}{2}} = 0.4$ (red in panel b), and $1.26(\lambda L)^{\frac {1}{2}}$ (for the same $R_{Y0}/L, \alpha _0$ , red in panel c), where the first focusing region is the waist captured by the vacuum solution. Varying the initial width $W_0$ affects the initial expansion of the beam, from a pronounced initial expansion in panel (a,b) to no initial expansion for some $R_{Y0}/L$ in panel (c). The inset in panel (c) focuses on the initial propagation region, which exhibits converging or diverging beams depending on the initial conditions.

Figure 4 shows the beam width for different initial conditions $W_0/(\lambda L)^{\frac {1}{2}} = 0.13, 0.4, 1.26$ , and the same scan in $R_{Y0}/L$ as in figure 3, all while keeping a fixed initial incident angle $\alpha _0=30^{\circ }$ . For small $W_0=0.13(\lambda L)^{\frac {1}{2}}$ , the beam width $W_Y$ initially follows closely the vacuum solution (dashed line), independent of the initial conditions $W_0$ and $R_{Y0}$ . Small $W_0/(\lambda L)^{\frac {1}{2}}$ corresponds to a highly focused beam. Diffraction produces a strong initial expansion of the beam width $W_Y$ , which is followed by strong focusing slightly past the turning point (for all of $R_{Y0}$ ). For $W_0=0.40(\lambda L)^{\frac {1}{2}}$ , the initial growth is less severe and the beam focusing (value and location) exhibits a noticeable dependence on the initial $R_{Y0}$ . The dependence on the initial $R_{Y0}$ is even more noticeable for $W_0=1.26(\lambda L)^{\frac {1}{2}}$ , especially for the focus location: the beam even focuses twice along the beam for $R_{Y0}/L=-0.5$ (red curve), which is an initially focusing beam at the plasma edge (the first focusing region is the waist captured by the vacuum solution). The double focusing happens for both $W_0 = 0.4 (\lambda L)^{\frac {1}{2}}$ and $W_0 = 1.26 (\lambda L)^{\frac {1}{2}}$ . Comparing the same initial condition $R_{Y0}=\infty$ (launch at waist, black curve) between the different initial $W_0$ , we see that the beam focusing location takes place after the turning point for small initial $W_0 = 0.13, 0.40 (\lambda L)^{\frac {1}{2}}$ and approaches the turning point for large $W_0 = 1.26 (\lambda L)^{\frac {1}{2}}$ (the turning point location is given by the vertical dotted lines). This discussion has important consequences for interpreting the regions in the plasma with predominant contributions to backscattering in DBS measurements, as will be discussed in § 4.

In this section, we have seen how the 2-D linear-layer model exhibits focusing of the beam width around the cutoff. This focusing is separate from the focusing caused at the beam waist in vacuum. The beam focusing depends on the initial conditions for the incident angle $\alpha _0$ , normalised beam parameter $W_0/(\lambda L)^{\frac {1}{2}}$ and radius of curvature $R_{Y0}/L$ . Beam focusing tends to be enhanced for small $\alpha _0$ , while it tends to disappear for large $\alpha _0$ . For small $W_0/(\lambda L)^{\frac {1}{2}}$ , strong initial growth precedes strong focusing. For large $W_0/(\lambda L)^{\frac {1}{2}}$ , the beam focuses from the initial condition towards the focusing region, and the initial growth disappears for some initial conditions. These phenomena are not captured by the vacuum solution.

Having gained an understanding of the phenomenon of beam focusing, in what follows, we will use the analytical 2-D model solution to assess the impact of beam focusing on Doppler backscattering measurements.

4.1. The 1-D filter function $|F_{x\mu }|^2$

The filter function $|F_{x\mu }|^2$ in (3.11) provides the $k_x$ -selectivity of the DBS power. This is tied to the localisation along the path: the scattering turbulent $k_x$ can be thought of as a parameter determining the location of scattering along the path $\tau _\mu$ , as $\tau _\mu$ is related to the turbulent scattering $k_x$ via the Bragg condition for backscattering in (3.6), which is used to arrive at (3.13). The radial component $K_{x\mu } = K_0(\cos \alpha _0 - \tau ^{\prime}_\mu )$ from (2.9) is a decreasing function of $\tau ^{\prime}_\mu$ from the launch, becoming negative after the turning point ( $K_{x\mu }=0$ , i.e. $\tau ^{\prime}_\mu = \cos \alpha _0$ ). The vertical component $K_y$ is constant, since the system is homogeneous in $y$ . This means that $K_\mu$ reaches its minimum at the turning point, which should enhance the DBS signal contribution at that location. As we saw in the previous section, the perpendicular beam width $W_{Y\mu }$ has a tendency to focus, which should further enhance the DBS signal power. In this section, we separate the enhancement due to ray and the beam contributions in the filter function.

Figure 5. Filter function $|F_{x\mu }|^2(\tau )/|F_0|^2 = K_0 W_0/K_\mu W_{Y\mu }$ and the corresponding ray component $K_0/K_\mu$ and beam component $W_0/W_{Y\mu }$ for varying values of the incident launch angle $\alpha _0$ , and fixed $R_{Y0} = \infty$ and $W_0 = 0.40(\lambda L)^{\frac {1}{2}}$ . Note that in all cases, the filter function is predominantly affected by the beam term $1/W_{Y\mu }$ which represents the focusing.

Figure 5 shows the ray (blue) and beam (red) components of $|F_{x\mu }|^2$ for varying incident angles $\alpha _0 = 10^{\circ }, 30^{\circ }, 50^{\circ }$ and $ 80^{\circ }$ , and fixed $W_0 = 0.40(\lambda L)^{\frac {1}{2}}$ and $R_{Y0} = \infty$ (launch at the waist). For convenience, we normalise the filter function to its value at $\tau =0$ , $|F_0|^2$ . The filter functions correspond to some of the beam solutions in figure 3. For $\alpha _0=10^{\circ }$ , the filter function $|F_{x\mu }|^2$ is strongly localised at the focus location, which almost matches the location of the turning point, or cutoff (maximum of the ray term in blue). The turning point, or cutoff, is represented in each figure by a vertical grey dotted line. Careful inspection shows that the focus location takes place in the vicinity but after the turning point $\tau ^{\prime}_\mu = \cos \alpha _0$ for $\alpha _0=10^{\circ }$ , consistent with beam focusing taking place after the turning point in figure 3. This shows that the contributions to the DBS power are predominantly from the vicinity (but after) the turning point for these initial conditions.

For larger incident angles, the turning point takes place for smaller $\tau _\mu ^{\prime} = \cos \alpha _0$ , while the beam term peaks at larger $\tau _\mu ^{\prime}$ (beam focus takes place at larger $\tau _\mu ^{\prime}$ ): the ray and beam terms compete to yield a filter function that is less peaked and less localised around the focus location for increasing $\alpha _0$ . The ray term (blue) decreases in amplitude for larger $\alpha _0$ , but less than the beam term (red), which decreases more. The beam term becomes broader in $\tau _\mu ^{\prime}$ for larger angles. These effects make the signal more delocalised along the path. This can be clearly seen in figure 5(c). The localisation becomes broad for $\alpha _0=50^{\circ }$ , where beam focusing becomes less important in favour of contributions from $\tau ^{\prime}_\mu = 0$ . In figure 5, the initial condition is chosen to be at the waist $W_{Y}(\tau =0)=W_0$ . In figure 7, we show the effect of the initial radius of curvature $R_{Y0}$ on the filter function. An initially expanding beam (positive and finite $R_{Y0}/L$ ) is shown in figure 7 not to make a dramatic difference on $|F_{x\mu }|^2$ . For larger incident angles ( $\alpha _0 = 80^{\circ }$ ), the filter peaks around $\tau = 0$ , corresponding to the waist initial condition of figure 5 (figure 7 shows that the behaviour is qualitatively similar for different initial $R_{Y0}$ ). Therefore, for large incident angles, the DBS signal becomes delocalised with predominant contributions before and after the turning point (figures 5 and 7).

Figure 6. Filter function $|F_{x\mu }|^2(\tau )/|F_0|^2 = K_0 W_0/K_\mu W_{Y\mu }$ and the corresponding ray component $K_0/K_\mu$ and beam component $W_0/W_{Y\mu }$ for varying values of the initial width $W_0 = 0.13, 0.40, 1.26(\lambda L)^{\frac {1}{2}}$ , and fixed $R_{Y0} = \infty$ and $\alpha _0 = 30^{\circ }$ . Note that in all cases, the filter function is predominantly affected by the beam term $1/W_{Y\mu }$ which represents the focusing.

Figure 7. Values of the filter function $|F_{x\mu }|^2/|F_0|^2$ for varying incident angles $\alpha _0$ and fixed initial $W_0 = 0.40 (\lambda L)^{\frac {1}{2}}$ as a function of the turbulent, selected $k_x$ component. For small incident angle $\alpha _0=10^{\circ }$ , the filter is strongly peaked near $k_x = 0$ , consistent with the signal being strongly localised near the turning point region. The signal is in fact sensitive to slightly positive $k_x\gt 0$ , which corresponds to a focusing slightly after the turning point. For $\alpha _0=30^{\circ }$ , the peak near $k_x=0$ decreases and shifts to larger $k_x$ values, meaning that the signal starts getting important contributions from finite $k_x$ turbulent fluctuations away from the turning point. For $\alpha _0=50^{\circ }$ , the peak has almost disappeared and the signal receives close-to-uniform contributions in the range of $-1 \lesssim k_x/K_0 \lesssim 2$ , corresponding to a highly delocalised signal along the beam path. For $\alpha _0 = 80^{\circ }$ , the beam expands for most of its path and the filter function approaches the vacuum solution.

Figure 6 shows the ray (blue) and beam (red) components of $|F_{x\mu }|^2$ for different values of the initial width $W_0 / (\lambda L)^{\frac {1}{2}} = 0.13, 0.40$ and $1.26$ , and fixed incident angle $\alpha _0 = 30^{\circ }$ and $R_{Y0} = \infty$ (launch at the waist). The filter functions in figure 6 correspond to some of the beam solutions in figure 4. In this case, the ray term remains constant. The beam term is narrowest around the initial launch (waist) and around the focusing point for small initial width $W_0=0.13 (\lambda L)^{\frac {1}{2}}$ . This is because a small initial width produces rapid expansion of the beam (see figure 4), which gives negligible contributions to the beam term outside the waist and focus point. This suggests a very localised contribution to the DBS power from the vicinity of the focus point as well as the plasma edge (if enough fluctuations are present). For initial width $W_0=0.40(\lambda L)^{\frac {1}{2}}$ , the location of the peak in $|F_{x\mu }|^2$ is similar to that for $W_0 = 0.13 (\lambda L)^{\frac {1}{2}}$ , but the filter function becomes broader, suggesting that the beam focusing is slower and takes place over a larger region, as shown in figure 4. For $W_0=1.26(\lambda L)^{\frac {1}{2}}$ , the filter function maximum has increased in value (beam focuses more) and has a peak as broad as the filter function for $W_0 = 0.4 (\lambda L)^{\frac {1}{2}}$ . Interestingly, in this case, the filter function $|F_{x\mu }|^2$ does not exhibit an initial decrease after the launch. This is because, in this condition, the beam does not initially expand following the launch, but only contracts from launch to the focus location (see black line in figure 4 c). Moreover, the beam focuses closer to the turning point than for $W_0/(\lambda L)^{\frac {1}{2}} = 0.13, 0.4$ .

Figures 3–7 show that the beam focusing in $W_Y$ tends to take place close to but after the turning point ( $K_{x\mu }=0, \tau ^{\prime}_\mu =\cos \alpha _0$ ) depending on the initial conditions for the incident angle $\alpha _0$ , width $W_0/(\lambda L)^{\frac {1}{2}}$ and especially on the phase front radius of curvature $R_{Y0}$ . This challenges the common belief that the DBS signal is always most sensitive at the turning point and has important consequences for interpreting the DBS scattered power.

It is equally instructive to characterise the filter function in terms of the turbulent scattered $k_x$ . Note that negative $k_x$ corresponds to scattering from the beam in its first pass from launch towards the turning point, and positive $k_x$ to scattering from the beam in its return journey away from the turning point. Figure 7 shows the variation of the filter function $|F_{x\mu }|^2$ for different incident angles $\alpha _0$ , different initial radii of curvature $R_{Y0}$ and $W_0 = 0.40 (\lambda L)^{\frac {1}{2}}$ as a function of the turbulent scattered $k_x$ . For $\alpha _0=10^{\circ }$ , $|F_{x\mu }|^2$ is strongly peaked at $k_x=0$ . This means that the DBS signal is strongly localised near the turning point region and predominantly sensitive to fluctuations with $k_x \approx 0$ . For $\alpha _0=30^{\circ }$ , the filter peak shifts towards larger $k_x$ . For even larger $\alpha _0=50^{\circ }$ , the filter is sensitive almost uniformly to $-1 \lesssim k_x/K_0 \lesssim 2$ , with peaks of the filter function at both positive and negative $k_x$ , and surprisingly a dip in the vicinity of $k_x=0$ . This suggests that the DBS power is sensitive predominantly to specific values of $k_x$ , with a subdominant contribution from $k_x=0$ . For $\alpha _0=80^{\circ }$ , the beam follows closely vacuum propagation and the filter has contributions from small $k_x/K_0$ near $0$ , but this time originating near $\tau =0$ (near the vacuum beam-waist, or edge of the plasma in a real experiment) and decays for larger $k_x$ . The beam expands for most of its path and the filter function approaches the vacuum solution (black dashed line).

Figure 8. Values of the filter function $|F_{x\mu }|^2/|F_0|^2$ for varying initial width $W_0/(\lambda L)^{\frac {1}{2}} = 0.13, 0.40, 1.26$ and fixed $\alpha _0 = 30^{\circ }$ as a function of the turbulent, selected $k_x$ . For small initial $W_0$ , the filter is strongly peaked at two values of $k_x$ : one negative that corrresponds to the entrance to the plasma, and one positive from the beam focusing. For small $W_0$ , the filter is also independent of the initial condition $R_{Y0}$ . For $W_0 = 0.40 (\lambda L)^{\frac {1}{2}}$ , the filter is dominated by a peak of amplitude and location similar to those of the peak due to focusing for small $W_0$ , but broader. For $W_0 = 1.26 (\lambda L)^{\frac {1}{2}}$ , the peak in $k_x$ stays broad and shifts towards $k_x=0$ for launch near the waist (large $|R_{Y0}/L|$ ).

Figure 8 shows the corresponding filter function for varying $W_0 / (\lambda L)^{\frac {1}{2}} = 0.13, 0.40$ and $1.26$ , different radii of curvature and fixed initial angle $\alpha _0=30^{\circ }$ . The behaviour of the filter function for small $W_0$ can be understood from figure 4. For small initial $W_0$ , the filter is strongly peaked at two different values of $k_x$ and is independent of the initial condition $R_{Y0}$ , consistent with the signal being strongly localised at the plasma entrance and near the turning point region. For $W_0 = 0.40 (\lambda L)^{\frac {1}{2}}$ , the filter is dominated by a peak with amplitude and location similar to those of the peak due to focusing for small $W_0$ , but broader, and hence the signal gets enhanced contributions from turbulent fluctuations away from the focusing region, depending on the initial condition $R_{Y0}$ . For $W_0 = 1.26 (\lambda L)^{\frac {1}{2}}$ , the peak in $k_x$ stays broad and shifts towards $k_x=0$ for launch near the waist (large $|R_{Y0}/L|$ ). The small $|R_{Y0}/L|$ cases are different, exhibiting distinctive narrow peaks that even appear twice for $R_{Y0}/L=-0.5$ , consistent with two consecutive focusing regions, as in figure 4.

4.2. Predictions of the measured turbulent spectra

Throughout this manuscript, we have focused on describing the dependence of the filter function $|F_{x\mu }|^2$ on experimentally relevant initial values for $\alpha _0$ , $W_0$ and $R_{Y0}$ . We should note, however, that the intrinsic power spectrum of the turbulence itself has a direct impact on the measurement. The turbulent spectrum is generally a decreasing function of $k_x$ , as demonstrated by direct numerical gyrokinetic turbulence simulations. In conditions where the filter function predominantly selects large $k_x$ (finite to large $\alpha _0$ ), the turbulence spectrum in $k_x$ should be taken into account to understand the dominant $k_x$ contributing to the backscattering signal. In this section, we apply what we have learned about the filter function $|F_{x\mu }|^2$ in the beam-tracing model and use it to understand its effect on the DBS backscattered power from a realistic density-fluctuation spectrum.

We use the electron-density-fluctuation wavenumber power spectrum $ \langle | \delta \hat {n}(k_x, k_y)|^2 \rangle _T$ obtained from nonlinear gyrokinetic simulations resolving strongly developed ETG turbulence in NSTX, which was examined by Ruiz et al. (Reference Ruiz2015, Reference Ruiz, Guttenfelder, White, Howard, Candy, Ren, Smith, Loureiro, Holland and Domier2019, Reference Ruiz, Guttenfelder, White, Howard, Candy, Ren, Smith and Holland2020 Reference Ruiz, Guttenfelder, White, Howard, Candy, Ren, Smith and Hollanda , Reference Ruiz, Guttenfelder, White, Howard, Candy, Ren, Smith and Hollandb ); Ren et al. (Reference Ren2020); Ruiz et al. (Reference Ruiz2022); Guttenfelder et al. (Reference Guttenfelder2022). The turbulent spectrum can be approximated by the following shape:

(4.1) \begin{align} \frac {\big\langle |\delta \hat {n} (k_x, k_y) |^2 \big\rangle _T}{n^2} \sim \frac {A}{1 + \big | \frac {k_x}{w_{k_x} }\big |^{\zeta } + \big | \frac {k_y - k_{y*}}{w_{k_y} } \big |^{\eta } }. \end{align}

This expression is fitted to the specific strongly driven ETG simulations that we have mentioned above (Ruiz et al. Reference Ruiz2022) to find: $\zeta \approx 3.14, \eta \approx 3.19, w_{k_x} \rho _s \approx 0.89, w_{k_y}\rho _s \approx 4.09, k_{y*}\rho _s \approx 6.53, A \approx 1.18 \ 10^{-6}$ (here $\rho _s$ is evaluated at the cutoff of the hypothetical DBS experiment). Ruiz et al. (Reference Ruiz2022) represented the fluctuation spectra in terms of the normal and binormal wavenumber components $k_n$ and $k_b$ perpendicular to the magnetic field, which can be easily calculated from the internal wavenumber components in a gyrokinetic code. The normal and binormal components correspond to the Cartesian $k_x$ and $k_y$ employed throughout this manuscript, as we have explained at the end of § 3. Importantly, from here on, we assume that the turbulent spectrum is uniform in space, that is, the parameters $A, \zeta , \eta , w_{k_x}, w_{k_y}$ and $k_{y*}$ are assumed to be constant throughout the whole plasma volume through which the beam propagates.

Figure 9. (a) Values of the normalised beam width $W_Y/W_0$ along the path $\tau$ for varying values of the incident angle $\alpha _0$ , $W_0 \approx 1.26(\lambda L)^{\frac {1}{2}}$ and $R_{Y0}/L = 0.5$ . The curves are plotted for the values of $\tau$ in which the beam is traversing the plasma slab, $x\gt 0$ . We assume vacuum for $x\lt 0$ and the beam is launched at $x=0$ . (b) Corresponding filter function $|F_{x\mu }|^2/|F_0|^2 = {K_0 W_0}/{ K_\mu W_{Y\mu } }$ plotted as a function of the scattered turbulent $k_x$ within the plasma slab, that is, for $|k_x/K_0 \cos \alpha _0| \lt 2$ .

DBS experiments routinely vary the angle of incidence $\alpha _0$ to select different wavenumbers from the turbulence spectrum. In what follows, we show how the beam width $W_Y$ , the filter function $|F_{x\mu }|^2$ and ultimately the measured turbulence spectrum are impacted by varying $\alpha _0$ .

Figure 9(a) shows the normalised beam width $W_Y/W_0$ from the beam-tracing equations as a function of $\tau$ for a range of incident angles $\alpha _0 = (10^{\circ }, 30^{\circ }, 50^{\circ }, 80^{\circ })$ , $W_0 = 1.26(\lambda L)^{\frac {1}{2}}$ and $R_{Y0}/L = 0.5$ . The beams in figure 9 have the same initial frequency $\Omega$ and reach a turning point location at $x=L_n = L\cos ^2\alpha _0$ (see figure 1). Figure 9(a) exhibits beam focusing for $\alpha _0 = 10^{\circ }, 30^{\circ }$ within the plasma while the beam focusing takes place outside the plasma for $\alpha _0 = 50^{\circ }$ . Both the value of the minimum beam width and the focusing location increase with $\alpha _0$ in this particular case.

Figure 9(b) shows the DBS filter function $|F_{x\mu }|^2$ corresponding to figure 9(a) plotted as a function of the turbulent scattered $k_x$ normalised to $K_0\cos \alpha _0$ . The scattered $k_x$ at every $\tau$ location is computed using the Bragg condition for backscattering written in Cartesian coordinates, $k_x = - 2 K_{x\mu } = - 2 K_0 (\cos \alpha _0 - \tau ^{\prime}_\mu )$ (see § 3 and Appendix B). The turning point takes place where $K_{x\mu }=0$ , i.e. $\tau ^{\prime}_\mu =\cos \alpha _0$ . Note how the filter function peaks closer to $k_x=0$ for decreasing $\alpha _0$ , corresponding to the beam focusing location approaching the turning point for small angles. The intensity of the enhancement decreases for increasing $\alpha _0$ . The dominant $k_x$ contributing to scattering is always positive $k_x\gt 0$ in this situation, which means that the signal is predominantly originating from plasma locations after the turning point. Note how this depends on the initial condition: figure 4(c) shows how the beam focusing can take place arbitrarily close to the turning point for the same $\alpha _0=30^{\circ }$ and $W_0 = 1.26 (\lambda L)^{\frac {1}{2}}$ values but different $R_{Y0}$ . Increasing values of $\alpha _0$ will move the filter function localisation of the DBS signal further and further away from the turning point, reaching the plasma exit for $\alpha _0 \approx 50^{\circ }$ in these conditions (see red line in figures 9 a and 9 b).

Figure 10. Density-fluctuation power $\langle |\delta \hat {n}|^2 \rangle _T$ and the product of the density-fluctuation power and the filter function $\langle |\delta \hat {n}|^2 \rangle _T |F_{x\mu }|^2/|F_0|^2$ as a function of the normalised scattered turbulent radial wavenumber $k_x\rho _s$ for different values of the scattered turbulence $k_{y0}\rho _s = -2K_{0}\rho _s \sin \alpha _0$ . The different $\alpha _0 = (10^{\circ }, 30^{\circ }, 50^{\circ }, 80^{\circ })$ correspond to scattered wavenumbers $k_{y0}\rho _s = -(1.50, 4.92, 10.15, 48.29)$ . Solid lines correspond to $k_x \gt 0$ , while dashed lines correspond to $k_x \lt 0$ . The red curve shows that the beam focusing has a strong effect in the measurement $k_x$ selection for $\alpha _0 \lesssim 30^{\circ }{-}40^{\circ }$ , while the effect becomes negligible for $\alpha _0 \gtrsim 40^{\circ }$ .

Figure 11. (a) Integrated filter function $\int \text {d}k_x\rho _s |F_{x\mu }|^2/|F_0|^2$ for the scattered $k_x$ within the plasma ( $|k_x/K_0| \lt 2\cos \alpha _0$ ). Note how the integrated weight $\int \text {d}k_x\rho _s |F_{x\mu }|^2/|F_0|^2$ is not constant, meaning that the filter function can affect the measured DBS $k_y$ spectrum. (b) Turbulent spectrum $\langle |\delta \hat {n}|^2(k_x=0, k_y) \rangle _T$ (black), $k_x$ -integrated spectrum $\int \text {d}k_x\rho _s \langle |\delta \hat {n}(k_x, k_y)|^2 \rangle _T$ (magenta) and synthetic DBS spectrum $\int \text {d}k_x\rho _s ( |F_{x\mu }|^2 /|F_0|^2 ) \langle | \delta \hat {n} |^2 \rangle _T (k_x, k_{y})$ (blue) as a function of $k_y$ of the turbulence. The turbulence spectrum has been fitted to a gyrokinetic simulation (Ruiz et al. Reference Ruiz2022). Note how the predicted DBS measurement is not able to capture the spectral peak of the turbulent spectrum (injection scale), but it is able to quantitatively reproduce the spectral exponent of the $k_x$ -integrated spectrum. The red dots indicate the specific angles $\alpha _0 = 10^{\circ }, 30^{\circ }, 50^{\circ }, 80^{\circ }$ given for reference.

Next, we quantify the combined effect of the filter function and its dependence on $k_x$ , in conjunction with a realistic, power-law turbulence spectrum. In figures 10 and 11, we normalise $k_x$ , $k_y$ and the initial wavenumber magnitude $K_0$ to the local ion sound gyroradius at the cutoff $\rho _s = c_s/\Omega _D$ , where $c_s = (T_e/m_D)^{1/2}$ is the local ion sound speed, $m_D$ the deuterium mass and $\Omega _D = eB/m_D$ is the deuterium gyro-frequency. Using the relation $\Omega = \omega _{pe0}/\cos \alpha _0$ between $\Omega$ and the local value of the electron plasma frequency at the cutoff $\omega _{pe0}$ (see § 2.3), we have $K_0\rho _s = (\beta _e m_D/2m_e)^{1/2}/\cos \alpha _0 \approx 4.26/\cos \alpha _0$ , where $\beta _e \approx 1\, \%$ is the electron beta using the local magnetic field, electron density and temperature of the NSTX experiment analysed by Ruiz et al. (Reference Ruiz2015, Reference Ruiz, Guttenfelder, White, Howard, Candy, Ren, Smith, Loureiro, Holland and Domier2019, Reference Ruiz, Guttenfelder, White, Howard, Candy, Ren, Smith and Holland2020 Reference Ruiz, Guttenfelder, White, Howard, Candy, Ren, Smith and Hollanda , Reference Ruiz, Guttenfelder, White, Howard, Candy, Ren, Smith and Hollandb ); Ren et al. (Reference Ren2020); Ruiz et al. (Reference Ruiz2022); Guttenfelder et al. (Reference Guttenfelder2022).

Figure 10 shows the radial-wavenumber dependence of the density-fluctuation power spectrum (blue lines) that would be sampled as the microwave beam propagates through the plasma. The product of the filter function $|F_{x\mu }|^2$ and the density-fluctuation power is shown in red. These quantities are plotted as a function of the turbulent $k_x$ normalised by $\rho _s$ using $k_x \rho _s = (k_x/K_0)(K_0\rho _s)$ . The solid lines correspond to $k_x$ scattered after the turning point ( $k_x \gt 0$ in figure 9 b), while the dashed lines correspond to $k_x$ scattered before the turning point ( $k_x\lt 0$ in figure 9 b). As expected, the filter function has an important effect on the DBS signal power for smaller incidence angles $\alpha _0 = 10^{\circ }, 30^{\circ }$ , while it becomes unimportant for larger $\alpha _0$ , at which point, the spectral falloff of the density spectrum becomes the dominant effect determining the $k_x$ selection in the DBS measurement. Note, for example, that the peaking of the filter function near the plasma exit for $\alpha _0=50^{\circ }$ becomes negligible because it occurs at large $k_x$ . Importantly, figure 10 contains both information about the radial localisation of the DBS power (if plotted as a function of $\tau$ ) as well as the $k_x$ -selectivity of the power spectrum in the DBS measurement (plotted as a function of $k_x$ ) through the combined effect of the filter function and the turbulent spectrum. We stress again that, in this discussion, we are neglecting the spatial dependence of the turbulence intensity and spectrum, $\langle |\delta \hat {n}|^2 \rangle _T$ , which is assumed not to vary along the path.

Ultimately, one of the main purposes of DBS diagnostics as a turbulence measurement is to obtain the density-fluctuation power spectrum in the binormal wavenumber of the turbulence, known as $k_\perp$ in the DBS jargon. In the context of this manuscript, $k_\perp$ is $k_y$ . To understand the impact of the filter function in such a measurement, in figure 11(a), we plot the integrated filter function $|F_{x\mu }|^2$ over the relevant $k_x$ which would be sampled as the beam propagates through the plasma as a function of the selected wavenumber $k_{y} = 2K_0 \sin \alpha _0$ . Since $|F_{x\mu }|^2$ multiplies the density-fluctuation spectrum $\langle |\delta \hat {n}|^2 \rangle _T$ , a constant $\int \text {d}k_x\rho _s |F_{x\mu }|^2/|F_0|^2$ as a function of $k_{y}$ would suggest that the DBS measurement is able to reproduce the true shape of the density-fluctuation wavenumber power spectrum $\int \text {d}k_x\rho _s \langle |\delta \hat {n}|^2(k_x, k_{y}) \rangle _T$ . Figure 11(a) shows that the $k_{y}$ dependence of $\int \text {d}k_x\rho _s |F_{x\mu }|^2/|F_0|^2$ varies by a factor $\sim 3$ for $k_{y}\rho _s \sim 0.4{-}40$ ( $\alpha _0 \sim 3^{\circ }{-}80^{\circ }$ ). The red dots correspond to the values of the integral of the filter function for $\alpha _0 = 10^{\circ }, 30^{\circ }, 50^{\circ }$ and $80^{\circ }$ . For small $k_{y}$ (small $\alpha _0$ ), the filter function is a decreasing function of $k_{y}$ . This is due to the fact that for small $\alpha _0$ , the beam width focuses and produces an overall enhancement upon integration over the sampled $k_x$ . In fact, $W_Y \sim \alpha _0$ for small $\alpha _0$ (Belrhali et al. Reference Belrhali, Ruiz Ruiz, Parra and Hall-Chen2025), that is, $|F_{x\mu }|^2 \sim 1/\alpha _0^2$ ( $K_\mu \sim \alpha _0$ ). Since beam tracing assumes $W_Y \gg \lambda$ , the decrease of $W_Y$ at small $\alpha _0$ means that, for small enough angles, beam tracing is bound to break down: the prediction of a signal enhancement for small angles might overestimate the true enhancement in reality, as discussed in detail by Maj et al. (Reference Maj, Pereverzev and Poli2009, Reference Maj, Balakin and Poli2010). The precise angles for which this might happen and corrections to the beam-tracing model for small angles will be the subject of future publications.

For $k_{y}\rho _s \approx 4{-}5$ ( $\alpha _0 \approx 30^{\circ }$ ), the integrated filter function plateaus before abruptly decreasing for $\alpha _0 \approx 50^{\circ }$ . A similar plateau behaviour is observed in detailed measurements of the scattered power in DIII-D, reported in a recent publication by Pratt et al. (Reference Pratt2023). In our case, this is due to the fact that for such large angles, beam focusing has started to take place close to the plasma exit (see figure 9 b). This eliminates the signal enhancement due to the beam focusing inside the plasma (the $k_x$ integration is only performed inside the plasma), therefore decreasing the value of the integrated filter function. For larger $k_{y}$ ( $\alpha _0 \gtrsim 60^{\circ }$ ), the integrated filter increases. This last increase is due to the fact that the launch condition at the plasma edge starts to become important. For such large incidence angles, the beam is practically glancing and can be treated as a high-k backscattering system (Rhodes et al. Reference Rhodes, Peebles, Nguyen, VanZeeland, deGrassie, Doyle, Wang and Zeng2006; Hillesheim et al. Reference Hillesheim, Crocker, Peebles, Meyer, Meakins, Field, Dunai, Carr and Hawkes2015). Figure 11(a) shows that the range of sampled radial wavenumbers in the density-fluctuation spectrum can have an important effect on the measured DBS power spectrum.

Figure 11(b) shows the turbulence binormal-wavenumber power spectrum $\langle |\delta \hat {n}|^2(k_x=0, k_y) \rangle _T$ (black) corresponding to $k_x=0$ . This follows the traditional interpretation that dominant contributions to the DBS signal originate from the cutoff. The magenta line shows the $k_x$ -integrated wavenumber spectrum $\int \text {d}k_x\rho _s \langle |\delta \hat {n}|^2 (k_x, k_y)\rangle _T$ . This could be interpreted as a line-integrated measurement of the density-fluctuation spectrum. The DBS system can be interpreted as a line-integrated measurement for large angles $\alpha _0$ , since for large $\alpha _0$ , $|F_{x\mu }|^2$ varies slowly and is not strongly peaked (see figure 7 d). The extent to which it can be interpreted as a line integral rather than a localised measurement is quantified by the filter function $|F_{x\mu }|^2$ , which we discuss next.

The blue dots in figure 11(b) show expected synthetic scattered power spectrum in a DBS measurement, given by the quantity $\int \text {d} k_x\rho _s \left ( |F_{x\mu }|^2/|F_0|^2 \right ) \langle |\delta \hat {n}|^2 \rangle _T (k_x, k_{y})$ in this beam-tracing model (see (3.9)). For each incidence angle $\alpha _0$ (each $k_{y}$ ), the synthetic DBS scattered power includes the effect of the filter function and the radial wavenumber dependence of the density-fluctuation spectrum. As expected from the non-monotonic nature of $\int \text {d}k_x\rho _s |F_{x\mu }|^2 / |F_0|^2$ from figure 11(a), the synthetic scattered power does not reproduce the true density-fluctuation wavenumber spectrum for $k_x=0$ nor the $k_x$ -integrated spectrum. The peak in the density-fluctuation power spectrum ( $k_y\rho _s \approx 6{-}8$ ), which is the driving scale of the turbulence, cannot be identified in the synthetic power spectrum. The peak should be visible for $\alpha _0 \approx 30^{\circ }{-}40^{\circ }$ , where the beam model is expected to be quantitatively accurate (Maj et al. Reference Maj, Pereverzev and Poli2009, Reference Maj, Balakin and Poli2010). The enhancement of the synthetic scattered power with respect to the true turbulence spectrum for low $k_{y}$ (small incidence angles $\alpha _0 \lesssim 20^{\circ }$ ) is due to the enhancement of the signal by beam focusing. For larger $k_{y}\rho _s \gtrsim 10$ ( $\alpha _0 \gtrsim 40^{\circ }$ ), the synthetic spectrum exhibits a power-law spectral decay $\propto k_{y}^{-3.05}$ , which is different from that of the true turbulent spectrum for $k_x=0$ ( $\propto k_{y}^{-4.16}$ ). The spectral decay is shown to quantitatively agree with the $k_x$ -integrated spectrum, for which $\propto k_{y}^{-2.84}$ . This suggests that DBS scattered power measurements could accurately capture the spectral exponent of the $k_x$ -integrated spectrum, in contrast to the traditional belief that DBS measurements only select the $k_x=0$ component (Hillesheim et al. Reference Hillesheim, Holland, Schmitz, Kubota, Rhodes and Carter2012; Holland et al. Reference Holland2012). For this particular NSTX-inspired case, the DBS scattered power measurements can quantitatively reproduce the true spectral exponent of the $k_x$ -integrated, density-fluctuation power spectrum.