Electrostatic microturbulence in W7-X: comparison of local gyrokinetic simulations with Doppler reflectometry measurements

The first experimental campaigns of Wendelstein 7-X (W7-X) have shown that turbulence plays a decisive role in the performance of neoclassically optimized stellarators. This stresses the importance of understanding microturbulence from the theoretical and experimental points of view. To this end, this paper addresses a comprehensive characterization of the turbulent fluctuations by means of nonlinear gyrokinetic simulations performed with the code stella in two W7-X scenarios. In the first part of the paper, the amplitude of the density fluctuations is calculated and compared with measurements obtained by Doppler reflectometry (DR) in the OP1 experimental campaigns. It is found that the trend of the fluctuations along the radius is explained by the access of the DR system to different regions of the turbulence wavenumber spectrum. In the second part of the article, frequency spectra of the density fluctuations and the zonal component of the turbulent flow are numerically characterized for comparisons against future experimental analyses. Both quantities feature broad frequency spectra with dominant frequencies of O(1)–O(10)  kHz.


Introduction
Neoclassical transport has traditionally been a major source of concern with regard to energy losses, particle transport and impurity accumulation at the core of stellarators.In generic stellarators, trapped particles experience large radial excursions at low collisionality.In addition, ambipolarity is preserved by the onset of a radial electric field that, in relevant scenarios, leads to a strong inward convection of impurities.As a consequence, turbulent transport, driven by fluctuations with frequency much lower than the gyrofrequency of the plasma species and length scale of the order of the Larmor radius, has been explored in less detail in stellarators than in tokamaks.Nevertheless, the first experimental campaigns of the stellarator Wendelstein 7-X (W7-X) [1,2] -the first large stellarator optimized for neoclassical transport-have demonstrated that turbulence is an essential ingredient to understand how plasmas perform in this device.It has been observed that, in standard plasmas heated via electron-cyclotron-resonance-heating (ECRH), turbulent energy losses exceed by an order of magnitude the neoclassical ones [3].In addition, turbulent diffusion has been shown to underly the practical absence of strong accumulation of impurities, which, in general, exhibit a low peaking factor [4,5].The shape of the density profiles, with no sign of hollowness predicted by neoclassical theory [6], is explained by a robust turbulent pinch that manifests in broad regions of the parameter space [7].In parallel to these observations and partly as a result of them, the interest in modelling gyrokinetic turbulence in stellarators has increased over the recent years.New codes have been developed [8,9,10], verified against each other [11,12,13] and applied to fundamental questions in W7-X and other stellarators [14,15].Nonetheless, comparisons of numerical gyrokinetic simulations against diagnostic measurements of turbulent fluctuations are scarce in stellarators and, in particular, in W7-X [16,17].
In this study, we present a numerical characterization of the turbulent fluctuations in W7-X, comparing them with the existing experimental measurements [18] obtained with the Doppler reflectometer (DR) system [19,18,20,21].This diagnostic is capable of providing the amplitude of density fluctuations (|δn|) and rotation velocity (u ⊥ ) at a specific spatial location in the plasma for a selected perpendicular wavenumber (k DR ⊥ ).In this work we use the gyrokinetic code stella [8] to simulate, in the radial range of measurement of the DR system, two different ECRH discharges by means of flux tube nonlinear simulations, considering both ions and electrons as kinetic species.The amplitude of the density fluctuations computed with stella is compared against measurements by the DR system.In addition, we provide other quantities for future comparisons with the same diagnostic, in particular the zonal component of the turbulent E × B flow.
The rest of the paper is organized as follows.In section 2, the features of the gyrokinetic code stella are briefly summarized.The coordinates used by the code are introduced and the equivalence between the wavenumber of measurement of the DR system and that used in stella is presented.In section 3, the plasma profiles of the simulated ECRH discharges are described, as well as the configuration and parameters used in the simulations.In section 4, the squared amplitude of density fluctuations computed with stella for five different radial positions is compared against the measurements of the DR system analyzed in [18] for the two discharges studied.In section 5, quantities suitable for future comparisons with the DR system are predicted.In particular, in subsection 5.1, the frequency spectra of density fluctuations expected to be measured by the DR system are given.In subsection 5.2, the frequency spectra of zonal flow fluctuations are obtained.Apart from aspects related to DR measurements, in section 6, linear and nonlinear simulations are presented in order to provide a picture of the turbulence in the two scenarios considered throughout this paper.Finally, section 7 contains the summary and conclusions of this work.
2. The code stella, coordinate system and correspondence with k DR ⊥ The gyrokinetic formalism [22] is the theoretical framework to study microturbulence in strongly magnetized plasmas.It is based on the average over the gyrophase, removing the fastest degree of freedom of the system and retaining the effect of the finite Larmor radius size, which is comparable to the wavelength of turbulence.In this work we use the flux tube version of the gyrokinetic code stella.The flux tube domain fully exploits the scale separation between the typical variation lengths of the background magnetic field, plasma profiles and turbulent fluctuations along the magnetic field line with respect to the typical scales of the turbulent fluctuations in the plane perpendicular to the magnetic field.The code solves the gyrokinetic Vlasov and quasineutrality equations [8,11] for an arbitrary number of species using a Fourier spectral treatment in the plane perpendicular to the magnetic field and a mixed implicit-explicit method for the calculation of the parallel streaming and acceleration terms in the gyrokinetic Vlasov equation.This mixed algorithm allows to use larger time steps than explicit methods, being this especially beneficial when kinetic electrons are considered.Usually, standard twist-and-shift boundary conditions [23] are defined in the direction parallel to the magnetic field in flux tube codes.When the global magnetic shear of the device is small, these boundary conditions impose severe restrictions on the minimum wavelength considered in the direction perpendicular to the magnetic field or, equivalently, on the size of the flux tube in that direction.In particular, in W7-X (the device considered in this work), as one moves towards the magnetic axis, the global shear becomes critically low as it is shown in section 3. To avoid these restrictions and have control on the size of the flux tube, the stellarator symmetric twist-and-shift boundary conditions introduced in [24] have been implemented in stella and applied to the analyses presented below.
The coordinates used in the code are denoted by {x, y, ζ, v ∥ , µ}.We assume stellarator configurations with nested magnetic surfaces and define general spatial coordinates {r, α, ζ}, where r ∈ [0, a] is a radial coordinate that labels magnetic surfaces, α ∈ [0, 2π) is an angular coordinate labeling magnetic field lines on each surface and ζ ∈ [ζ min , ζ max ] is a coordinate along magnetic field lines.Here, a is the minor radius of the device and ζ min and ζ max are the ends of the flux tube in the ζ coordinate.In particular, for r we choose where 2πψ t (r) is the toroidal flux enclosed by the surface labeled by r and ψ t,a := ψ t (a).
For the coordinate α, we take where θ and ζ are, respectively, the poloidal and toroidal PEST flux coordinates [25] and ι(r) is the rotational transform.In these coordinates, the magnetic field takes the form B = dψ t dr ∇r × ∇α.
If r 0 and α 0 are the values of r and α that define the flux tube (and, actually, the center of the simulation domain), the spatial coordinates {x, y} in the plane perpendicular to B can be expressed as and Therefore, the modes k x and k y used in the Fourier treatment of the equations can be expressed as where k ⊥ is the wavevector perpendicular to B. Regarding the velocity coordinates, the code employs {v ∥ , µ}, where v ∥ is the component of the velocity parallel to the magnetic field line and µ = m j v 2 ⊥ /2B is the magnetic moment, with v ⊥ the component of the velocity perpendicular to the magnetic field, B = |B| the magnetic field strength and m j the mass of the species j.
On the other hand, the direction of measurement of the DR system is perpendicular both to the radial direction and to the magnetic field [18].Thus, defining the unitary vector pointing along the direction of measurement as e DR = B × ∇x/|B × ∇x|, the perpendicular wavenumber of measurement is expressed as Equation ( 7) can be written in the coordinates used in stella as where we have taken into account that the wavenumbers defined in stella are normalized with the ion thermal gyroradius ρ i , given by with Z i the ion charge number, e the proton charge, B a a reference magnetic field, whose explicit expression can be found in [8,11] and v th,i = 2T i /m i the ion thermal speed, where T i is the ion temperature.Expression (8) implies that, given a certain set of plasma profiles and a certain wavelength of the microwave beam launched by the DR system, the corresponding value of the normalized binormal component of the wavevector defined by the code changes with the spatial location.This is due to the modification from point to point of the geometric quantities |∇x| and B, as well as of T i (r).

Plasma scenarios and simulation settings
In this section we describe the input for the simulations, including details about the flux tube choice, resolution, magnetic equilibrium and plasma profiles considered throughout the present work.The selected plasma parameters correspond to two gas puff-fueled ECRH scenarios, without NBI injection for heating, from the OP1.2b campaign of W7-X.They belong to the discharges #180920.13 and #180920.17,both using the same heating power of P ECRH ≃ 4.7 MW without electron cyclotron current drive (ECCD) and with the former featuring a lower density than the latter.Hence, in this paper, they will be referred to as 'low density scenario' (#180920.13)and 'high density scenario' (#180920.17).Analyses of the amplitude of the fluctuations measured with the DR system have been presented for these two discharges in [18], which allows a straightforward comparison of the calculations obtained with stella and those already published.The profiles of density (n), ion temperature (T i ) and electron temperature (T e ) of these two scenarios, obtained as fits to the experimental data, are represented in figure 1.From left to right: density (solid line) and its normalized gradient (dashed line), ion temperature (solid line) and its normalized gradient (dashed line), electron temperature (solid line) and its normalized gradient (dashed line).
In this figure, the normalized gradients of these quantities, defined as are also depicted.In order to cover the radial range of measurement of the DR system reported in [18], five radial positions, labeled by r 0 /a = {0.5, 0.6, 0.7, 0.8, 0.9}, have been considered.The values of the plasma profiles and their normalized gradients for the selected radial locations are listed in table 1.Finally, the magnetic geometry used as input for the simulations is provided by the ideal MHD equilibrium code VMEC [26].This configuration is an example of a standard (also known as EIM) W7-X equilibrium in which all planar coil currents are set to zero and the non-planar coil current are set to the same value (see [27] for a detailed description of the different W7-X configurations) with effective minor radius a = 0.514 m, major radius R 0 = 5.513 m, magnetic field on-axis B ax = 2.52 T, and, for normalization purposes, see equation ( 9), B a = 2.28 T.
The global shear, defined as ŝ = − r ι dι dr , and ι profiles corresponding to this configuration are represented in figure 2. During the first W7-X operational phase (OP1), a DR system measured at locations nearby the bean-shaped toroidal plane (ζ = 0) of the plasma.The system launches a microwave beam whose cut-off determines the point of measurement of the perpendicular flow of density fluctuations.Those positions can be obtained applying ray-tracing techniques.In figure 3, the measurement positions for the high density scenario ‡ at the radial locations considered in the gyrokinetic simulations, obtained using the ray-tracing code TRAVIS [28], are represented with stars.The intrinsic restrictions of the flux tube formalism and the boundary conditions used in this study to access those measurement locations should also be noted.The parallel boundary conditions implemented for this work [24] restrict the selected flux tubes to stellarator symmetric ‡ The measurement positions for the low density scenario are very close to those represented in figure 3.
ones, i.e. those fulfilling B(θ, ζ) = B(−θ, −ζ).In W7-X, the suitable flux tubes are then the 'bean' flux tube, which is centered with respect to the position (θ, ζ) = (0, 0) and the 'triangular' flux tube, which is centered at (θ, ζ) = (0, π/5) [11].The intersections of these two flux tubes, assuming they extend one turn in the poloidal direction, with the bean-shaped toroidal plane of W7-X can also be found in figure 3. The bean flux tube position (θ, ζ) = (0, 0) provides the spatial location nearest to the positions of measurements of the DR system.Therefore, the information to compare against the measurements of the DR system will be extracted from that flux tube at (θ, ζ) = (0, 0).In addition, concerning the perpendicular wavenumber of measurement, the relationship between k DR ⊥ and k y ρ i (see expression (8)) is represented as a function of the radial coordinate in figure 4 for a wide range of k DR ⊥ .The k DR ⊥ values selected in the experiment (included in table 1 and reported in [18]) are indicated in figure 4     ⊥ .These plots are the result of evaluating expression ( 8) at (θ, ζ) = (0, 0) considering the T i profiles represented in figure 1 for the low (left) and high (right) density scenarios.The white squares represent the perpendicular wavenumbers measured by the DR system at each radial location.

Radial dependence of density fluctuations: comparison between stella results and DR measurements
The DR system measures the backscattered power (S) of a microwave beam launched into the plasma, which, for low turbulence levels, is proportional to the squared amplitude of the density fluctuations (S ∝ |δn| 2 ) [29].On the other hand, the code stella provides, at each instant and each position along the flux tube, the decomposition in Fourier space of the density fluctuations (δn) on the plane perpendicular to the magnetic field.The density fluctuations at r 0 and α 0 can be expressed as As explained in section 3, the simulated positions nearest to the measurement locations are those of the bean flux tube with ζ = 0. Once | δn kx,ky | 2 (ζ = 0, t) is obtained from the simulations, the quantity to compare against the DR system measurements is postprocessed as follows.First, for each pair of wavenumbers (k x , k y ), the spectrum in k y is obtained by summing over k x as Second, we time-average each (δn) 2 (k y , t) over the saturated nonlinear phase, obtaining the time-averaged squared density fluctuations, ⟨(δn) 2 ⟩ t (k y ).For illustrative purposes, figure 5 shows the time trace of (δn) 2 for the mode k y ρ i = 0.8, which is the one with the largest amplitude at r 0 /a = 0.6 in the high density scenario.In figure 6, the results obtained for ⟨(δn) 2 ⟩ t for the five scanned radial positions and for the two scenarios are represented as a function of k y ρ i .The values of k y that correspond to the wavenumbers of measurement of the DR system, that we denote by k DR y = k y (k DR ⊥ ), are also indicated with vertical dashed lines in these plots.It can be observed that the DR system, depending on the radial location, measures very disparate values of k y ρ i .Whereas at r 0 /a = 0.9 the system accesses nearly scales with k y ρ i ∼ 1, as r 0 decreases it explores different regions of the spectrum, reaching up to k y ρ i ∼ 3.5.Finally, considering the squared density fluctuations for each radial location at the specific value of k y accessed by the diagnostic, the comparison against the backscattered power, S, can be carried out.In figure 7 (left), the numerical results obtained with stella, (δn) 2 DR , are represented as a function of the radial coordinate, while in figure 7 (right), the measurements of S obtained with the DR system (reported in [18]) can be found.It is important to note that the DR measurements represented in figure 7 (right) are expressed in dB, but the reference normalization value is unknown.Hence, S [dB] ∝ 10 log 10 (δn) 2  DR + C, with C an arbitrary constant.For this comparison, C has been chosen to make S (r 0 /a = 0.5) = 10 log 10 (δn) 2 DR (r 0 /a = 0.5) + C or, equivalently, to set the lowest value of the numerical and experimental results at the same level.It can be observed that the numerical results shown in figure 7 (left) and the measurements in figure 7 (right) exhibit a monotonic increase of the squared density fluctuations with r/a.In addition, the difference between their minimum and maximum values is approximately 15 dB in both cases.Finally, the numerical and experimental results corresponding to the low density scenario show lower turbulent fluctuations than those of the high density scenario.These features, in particular the monotonic increase of the density fluctuations towards the edge, result, partially, from the fact that the DR system measures, as the radial position changes, at different locations of the k y spectrum, as we anticipated in figure 6 and related discussion.Indeed, the density fluctuations integrated both in k x and k y , i.e δn = ky ⟨(δn) 2 ⟩ t , loses the monotonic increasing behaviour with r/a, as figure 8 illustrates, for the two discharges analyzed.In that figure, one can observe that the integrated density fluctuations have a maximum at nearly r 0 /a = 0.8 and ranges between 2 and 8 × 10 19 m −3 .This picture aligns remarkably well with what Phase Contrast Imaging (PCI) techniques and GENE [30] simulations have reported for similar plasmas in W7-X [16].On the other hand, returning to the comparison between stella calculations and DR measurements of figure 7, the numerical results do not exhibit the significant reduction that the measurements do for the low density scenario in the range r/a ∼ [0.6, 0.7].A possible source of discrepancy could be that the bean flux tube at ζ = 0 does not correspond to the exact spatial measurement position of the DR system (see figure 3).Nevertheless, as it is explained in Appendix A, when the poloidal deviation from the bean flux tube at ζ = 0 is considered, the numerical results do not change significantly.With respect to the model used, possible extensions, such that accounting for the full flux surface geometry [31] and, in that case, including the radial electric field as well, might be worth addressing in order to attempt a better agreement.With regard to collisions, which have been neglected, the normalized ion collision frequency between r/a = 0.5 − 0.8 (considering ions and electrons with one thermal speed) is aν i /v th,i ≈ 1.3 × 10 −3 − 1 × 10 −2 .This value is significantly smaller than the growth rates of most unstable electrostatic instabilities, typically of a few tenths in units of a/v th,i .Exceptionally, at the position r/a = 0.9, aν i /v th,i ≈ 3.5 × 10 −2 in both scenarios.Electron collisions are higher though (aν e /v th,i ≈ 0.1 − 3, with the lowest values at r/a = 0.5), nevertheless electrons do not seem to play a prominent role except at the innermost studied position (see section 6).The impact of these factors will be addressed in forthcoming works.

Numerical characterization of fluctuations in the frequency domain
This section extends the analysis of the fluctuations, performed for the two W7-X discharges described in section 3, to the frequency domain.The aim is to complete the characterization of the fluctuations with features that can be inferred from the time evolution of S and u ⊥ measured by the DR system.In particular, these features correspond to the frequency spectrum of the squared amplitude of density fluctuations, presented in section 5.1, and the frequency spectrum of the zonal component of the E × B flow, discussed in section 5.2.Future analyses, that automate the measurement of highly time-resolved traces of S, will enable a straightforward comparison against the fluctuation frequency spectra provided in section 5.1.On the other hand, the installation in a different toroidal sector of a second DR system -whose data acquisition and analysis will extend during the second W7-X operation phase (OP2)-will allow to identify oscillations in u ⊥ of zonal origin to be verified against the spectra presented in section 5.2.For these future comparisons it will be necessary to subtract the Doppler shift from the experimental measurements of frequency spectra.

Frequency spectra of density fluctuations
In section 4, we have addressed the comparison between the radial dependence of the squared amplitude of the density fluctuations computed with stella and those measured by the DR system.In that case, a time average of (δn) 2 (k y , t) was performed.If, instead, the time dependent squared density fluctuations, (δn) 2 (k y , t), are considered and a Fourier transform in t is taken, we can write leading to the frequency spectrum of the squared density fluctuations.Note that the obtained frequency spectrum differs from the power spectral density of δn.In this section, we assume that once the time evolution of the backscattered power S(t) of the beam launched by the system is known -since S(t) is proportional to the (δn) 2 -the comparison of the frequency spectra of the experimental signal S and the numerically obtained (δn) 2 is straightforward.In figure 9, the amplitude of δn ω is represented as a function of the frequency ω and k y .We define the amplitude of the Fourier harmonics Figure 10 shows δn DR (ω) for the two scenarios under consideration.In general, for all radial positions, the spectrum is broad and extends a few hundreds of kHz.Above frequencies of that order of magnitude, an abrupt drop of the amplitude is observed.
Below ω/2π ∼ 100 kHz, the spectra are rather flat.However, for some radial positionsr 0 /a = {0.5, 0.6, 0.7} in the low density scenario and r 0 /a = {0.6,0.7} in the high density one-peak values of the amplitude are observed for frequencies such that ω/2π ≲ 10 kHz.After a similar Fourier analysis of S(t), the experimental frequency spectra could be compared against those represented in figure 10.

Frequency spectra of zonal flow fluctuations
Zonal perturbations of the electrostatic potential φ are constant on flux surfaces and have radial structure.In the flux tube scheme, they correspond to the components of the potential with k y = 0 and finite k x .With a second Doppler reflectometer looking on the same flux surface as the first one but at a different toroidal location, it would be possible to identify zonal oscillations.Specifically, common frequencies found in the time evolution of u ⊥ (t) obtained by each DR could be identified with zonal E × B flow fluctuations projected along the measurement direction of the system, e DR .In this subsection we address the characterization of these frequencies.
The fluctuating E × B velocity can be written as where the fluctuating electrostatic potential can be expressed, for each radial position r 0 , in Fourier series as We consider the modes of the potential with k y = 0, for which expression (16) becomes Considering the projection of δv E along the direction of measurement of the DR diagnostic, e DR , we can compute the zonal flow fluctuations expected to contribute to the total perpendicular flow measured by the DR system, Taking a Fourier transform of δu ZF ⊥ (x, t) in t, allows us to obtain the frequency spectrum, δu ⊥ω , for each x position of our flux tube.
To improve the statistics of our results, we perform an average along x, defined for each frequency of the spectrum as where N x is the number of grid points along the x direction.Figure 11 shows the amplitudes of δu ZF ⊥ω normalized to their maximum value, eliminating the mode ω = 0 (which corresponds to the time-averaged value).For the radial positions r 0 /a = {0.7,0.8, 0.9} dominant frequencies cluster in the range ω/2π ∼ [5,10] kHz.For the innermost positions, r 0 /a = {0.5, 0.6}, dominant frequencies group around the § Here, the field line average is defined as lower boundary of ω with no clear peak in the spectra.In figure 11, we have also depicted with vertical dashed lines the frequency ω 0 above which the amplitude | δu ZF ⊥ω | is always smaller than half the maximum value in each case.It can be observed that, for all radial positions, frequencies ω 0 /2π ≤ 20 kHz dominate the spectrum of the fluctuating zonal E × B flow.

Comparisons between linear and nonlinear calculations of φ
The analysis presented in this work has focused, so far, on numerical quantities that can be directly compared against DR measurements.Nevertheless, aspects beyond that comparison, like the nature of the background turbulence, its localization along the flux tube or its linear properties have not been discussed.In the present section we briefly address some of these characteristics, comparing the linear and nonlinear frequency spectra of the electrostatic potential in order to assess to what extent linear frequencies remain during the nonlinear phase.
In linear simulations, the time evolution of each φ kx,ky ζ (t) is assumed to be proportional to exp[(γ − iω r )t], with ω r and γ the linear frequency and growth rate.On the other hand, from the saturated phase obtained from nonlinear simulations, we can calculate the frequency spectrum for each k y mode of ⟨φ⟩, summing over all k x components as This Fourier decomposition allows a direct comparison between the nonlinear frequency spectrum of ⟨φ⟩ and the frequencies obtained from linear simulations.Figure 12 shows, as a function of ω, the amplitude of φ ω normalized to the largest value found (at a certain ω = ω M ) for each k y , i.e.
In addition, this figure also depicts the linear frequencies obtained from k y scans for a fixed k x = 0. Since the modes φ kx,ky are complex, we do not expect a symmetry in | φ ω | with respect to ω = 0 and, in contrast to section 5, we have represented both positive and negative values of ω in figure 12.It can be seen that at high k y , the frequency spectra broadens and a dominant frequency in the nonlinear spectra is lacking.It is interesting though that for the wavenumbers measured by the DR system, changes in the sign of the linear frequency take place.Specifically, for the two analyzed scenarios, the linear frequency is negative for r 0 /a = 0.5 at k y = k DR y , whereas it is positive for all other radii, which points out to a change in the propagation direction of the drift waves driving the instability.On the other hand, for low k y , the linear frequency values are reasonably close to those of the dominant nonlinear ones, which correspond with the darkest regions of the maps in figure 12.This fact can be more clearly observed in figure 13, where the normalized amplitudes |φ ω |(k y ) are represented as a function of ω for k y ρ i = 0.6.In that case, the linear frequency, indicated by a triangle in the figure, is located fairly close to the peak of the nonlinear frequency spectrum for every radial position analyzed.Similar agreements have been reported for the ASDEX Upgrade tokamak in [32], comparing linear and nonlinear simulations carried out with GENE.Finally, we analyze the localization of the fluctuating electrostatic potential along the flux tube for the wavenumber explored by the DR system.In particular, in figure 14 we represent the time-averaged squared amplitude of φ evaluated at y = k DR y , this is To assess the correlation between the structure of nonlinear modes at k DR y and those computed for the most linearly unstable mode for the same k DR y (generally located and obtained at k x = 0), we have depicted the parallel structure of the electrostatic potential from these linear simulations in figure 15.It is observed that both linear and nonlinear results exhibit similar structures, with slight differences possibly due to the contribution of every k x mode considered in the nonlinear analysis and, of course, nonlinear effects.In particular, (φ) 2  DR is strongly localized around the center of the flux tube, ζ = 0, for every radial position except at r 0 /a = 0.5.That localization in ζ overlaps with regions of bad curvature in this configuration of W7-X.This suggests that ion-temperature-gradient (ITG) modes [33] play a leading role at r 0 /a > 0.5.In addition, we can conclude that the linear modes at k y = k DR y propagate in the ion diamagnetic direction, which in the present work correspond with positive values of ω (see white dots in figure 12), as it is expected for ITG modes.On the other hand, the fluctuating electrostatic potential for the position r 0 /a = 0.5 is localized at regions of magnetic field wells (see the insets of figure 14, where the structure of this fluctuating electrostatic potential is represented together with the magnetic field strength).This fact, and the change in the sign of the linear frequency, points out to a prominent role of trapped electrons on the turbulence at r 0 /a = 0.5.In summary, looking at k y = k DR y , where the standard linear stability analysis finds a change in the sign of the frequency, a very different localization of the turbulent electrostatic potential is found nonlinearly.

Summary and conclusions
In the present work we have used the electrostatic flux tube version of the gyrokinetic code stella with the aim of calculating important turbulent quantities that can be directly compared with Doppler reflectometry measurements in the W7-X stellarator.Five radial positions of two ECRH discharges corresponding to the first operation phase (OP1) of W7-X have been considered throughout the paper.
In the first place, nonlinear simulations have been carried out to compute the amplitude of the density fluctuations at the spatial region and perpendicular wavenumbers explored by the DR system.Numerical results and DR measurements of the squared density fluctuations cover a range of about 15 dB and feature in both cases a monotonic increase with the radial coordinate.This is due to the different regions of the wavenumber spectrum accessed by the diagnostic at each radial location.Discrepancies between the numerical and experimental results in the different fluctuation levels between the two scenarios have also been observed.To expand the numerical characterization of the turbulence towards other quantities measurable by DR, the frequency spectra of the density fluctuations and of the zonal flow fluctuations have been provided.With regard to the density fluctuations, in general, the largest amplitudes are found for frequencies around a few tens of kHz.For higher frequencies, a slight decrease follows up to approximate ω/2π ≃ 500 kHz, where an abrupt fall of the amplitude is found.Concerning the frequency spectra of the zonal flow fluctuations, they are dominated by frequencies with ω/2π ≤ 20 kHz and find their peak values for ω/2π ≲ 10 kHz at most radial positions.Future highly time-resolved characterization of the backscattered power and measurements of u ⊥ fluctuations at distant locations over the same flux surface will allow to validate these numerical findings.Finally, the turbulent electrostatic potential has been provided and compared against linear simulations, showing, for low k y values, that the dominant frequencies of the nonlinear spectrum cluster nearly around the frequency of the most unstable mode linearly.In addition, these results have also shown that, at r 0 /a = 0.5, trapped electrons seem to play a leading role on the background turbulence.

DR
The numerical study presented in this work has been performed considering the position (θ, ζ) = (0, 0), with respect to which the bean flux tube is centered.However, as it is shown in figure 3, that position does not correspond to the exact location where the beam launched by the DR system reflects and, consequently, where the measurement is performed.In order to assess the impact of this deviation in the results presented in section 4, we have obtained poloidal profiles of (δn) 2  DR , that are represented in figure A1.For this, (δn) 2  DR has been interpolated considering its value at the points of the plane ζ = 0 crossed by the bean and triangular flux tubes (see figure 3).The maximum poloidal deviation of the DR measurement points with respect to θ = 0, which is found at the innermost radii, is also indicated as a shadowed region in figure A1.From these results, deviations > 20% between (δn) 2 DR evaluated at the exact measurement location and that at (θ, ζ) = (0, 0) are not expected.Such differences are negligible compared to the quantitative discrepancy between the simulations and the experimental measurements shown in figure 7. Hence, considering the position (θ, ζ) = (0, 0) is a good approximation to the real location where the DR system measures.

Figure 1 :
Figure 1: Radial profiles for the low (light blue) and high (dark blue) density scenarios.From left to right: density (solid line) and its normalized gradient (dashed line), ion temperature (solid line) and its normalized gradient (dashed line), electron temperature (solid line) and its normalized gradient (dashed line).

Figure 2 :
Figure 2: Radial profiles of ι (left) and ŝ (right) of the standard (or EIM) configuration considered throughout this paper.

Figure 3 :
Figure 3: Points belonging to the bean (blue dots) and triangular (red squares) flux tubes of W7-X that lay on the plane ζ = 0.The measurement positions for the high density scenario, obtained with TRAVIS for the radial locations considered in the gyrokinetic simulations, are represented as stars.

Figure 4 :
Figure 4: Relation between k y ρ i and k DR ⊥ as a function of the radial position for a wide range of k DR⊥ .These plots are the result of evaluating expression (8) at (θ, ζ) = (0, 0) considering the T i profiles represented in figure1for the low (left) and high (right) density scenarios.The white squares represent the perpendicular wavenumbers measured by the DR system at each radial location.

Figure 5 :
Figure 5: Time evolution of (δn) 2 (k y ρ i = 0.8) at r 0 /a = 0.6 in the high density scenario.The range considered for the time-average is represented by the shadowed area and the mean value of this quantity is indicated with the horizontal dashed line.

Figure 6 :
Figure 6: k y spectrum of the time-averaged squared density fluctuations, ⟨(δn) 2 ⟩ t , obtained with stella for the low (left) and high (right) density scenarios.The dashed vertical lines indicate the wavenumber measured by the DR system, k DR ⊥ , at each radial position.The values of k DR ⊥ are normalized following the convention employed in [18], where ρ s := √ T i m i /eB.

Figure 7 :
Figure 7: Squared amplitude of the density fluctuations as a function of the radial coordinate for the low (light blue circles) and high (dark blue triangles) density scenarios.Left: Numerical results of (δn) 2 DR obtained with stella.The error bars represent the standard deviation from the mean value evaluated over the saturated nonlinear phase.Right: backscattered power measured by the DR system.

Figure 8 :
Figure 8: Amplitude of the density fluctuations integrated over k x and k y computed with stella as a function of the radial coordinate for the low (light blue circles) and high (dark blue triangles) density scenarios.

Figure 9 :
Figure 9: Frequency spectra of (δn) 2 for the low (top) and high (bottom) density scenarios computed with stella as function of k y .

Figure 10 :
Figure 10: Density spectra of δn DR computed with stella for the low (top) and high (bottom) density scenarios for the case k y = k DR y .

Figure 12 :
Figure 12: Nonlinear frequency spectrum of kx ⟨ φ kx,ky ⟩ ζ for the low (top) and high (bottom) density scenarios computed with stella as a function of k y .The amplitudes | φ ω | are normalized to their maximum value at each k y .The white dots are the results of linear k y scans assuming k x = 0.

Figure 13 :
Figure 13: Nonlinear results of |φ ω | for the low (top) and high (bottom) density scenarios, evaluated at k y ρ i = 0.6.The linear result for the mode (k x ρ i , k y ρ i ) = (0, 0.6) is represented as a black triangle in every situation.

Figure 14 :
Figure 14: Normalized (φ) 2 DR as a function of ζ, computed with stella for the low (left) and high (right) density scenarios.The insets show the results corresponding to r 0 /a = 0.5, together with the magnetic field strength as a gray solid line.

Figure 15 :
Figure 15: Normalized (φ) 2 DR for the linearly unstable mode with k x = 0 as a function of ζ, computed with stella for the low (left) and high (right) density scenarios.

Figure A1 :
Figure A1: (δn) 2 DR as a function of the poloidal coordinate θ for a fixed ζ = 0 for the low (left) and high (right) density scenarios.Computed as an interpolation of the (δn) 2 DR calculated for the bean (circles) and triangular (squares) flux tubes.The shadowed region shows approximately the maximum deviation of the DR measurement position from θ = 0.

Table 1 :
Local parameters used in the simulations at each radial position for the low and high density scenarios.From left to right: density, normalized density gradient, ion temperature, normalized ion temperature gradient, electron temperature, normalized electron temperature gradient and perpendicular wavenumber measured by the DR system.