{'fraction_non_alphanumeric': 0.03258025017406545, 'fraction_numerical': 0.008571222779765192, 'mean_word_length': 4.28981457962916, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Generation of parasitic axial flow by drift wave turbulence with broken symmetry: Theory and experiment

Detailed measurements of intrinsic axial flow generation parallel to the magnetic field in the controlled shear decorrelation experiment linear plasma device with no axial momentum input are presented and compared to theory. The results show a causal link from the density gradient to drift-wave turbulence with broken spectral symmetry and development of the axial mean parallel flow. As the density gradient steepens, the axial and azimuthal Reynolds stresses increase and radially sheared azimuthal and axial mean flows develop. A turbulent axial momentum balance analysis shows that the axial Reynolds stress drives the radially sheared axial mean flow. The turbulent drive (Reynolds power) for the azimuthal flow is an order of magnitude greater than that for axial flow, suggesting that the turbulence fluctuation levels are set by azimuthal flow shear regulation. The direct energy exchange between axial and azimuthal mean flows is shown to be insignificant. Therefore, the axial flow is parasitic to the turbulence-zonal flow system and is driven primarily by the axial turbulent stress generated by that system. The non-diffusive, residual part of the axial Reynolds stress is found to be proportional to the density gradient and is formed due to dynamical asymmetry in the drift-wave turbulence.

sheared azimuthal and axial mean flows develop. A turbulent axial momentum balance analysis shows that the axial Reynolds stress drives the radially sheared axial mean flow.
The turbulent drive (Reynolds power) for the azimuthal flow is an order of magnitude greater than that for axial flow, suggesting that the turbulence fluctuation levels are set by azimuthal flow shear regulation. The direct energy exchange between axial and azimuthal mean flows is shown to be insignificant. Therefore, the axial flow is parasitic to the turbulence-zonal flow system, and is driven primarily by the axial turbulent stress generated by that system.
The non-diffusive, residual part of the axial Reynolds stress is found to be proportional to the density gradient and is formed due to dynamical asymmetry in the drift-wave turbulence. a) R.H. and J.C.L. contributed equally to this paper.

I. INTRODUCTION
Plasma flows along the magnetic field play a vital role in the stabilization of MHD instabilities and the development of transport barriers. [1][2][3][4][5][6] In most existing magnetic confinement fusion devices, the parallel flow, or toroidal plasma rotation, is driven directly by external momentum sources, such as neutral beam injection (NBI). However, in large scale devices like ITER, the NBI driven rotation will not be efficient, due to limited neutral beam penetration into high density plasmas.
In order to optimize and improve the confinement regimes in ITER and beyond, it is important to uncover alternative mechanisms that can drive parallel flows.
A phenomenon called intrinsic flow has been identified in magnetically confined tokamak plasmas, [3][4][5][6][7][8] where the plasma rotates toroidally without any input of toroidal momentum. This intrinsic flow can be of the same order of magnitude as that driven by some NBI torques. 2,[5][6][7] Hence, there is strong interest in knowing whether intrinsic flow in future devices is sufficient to affect confinement and MHD stability. Empirical results show that intrinsic torque in H-mode plasmas scales with the plasma stored energy normalized by the plasma current ("Rice scaling"). 2 Further measurements from Alcator C-Mod reveal that the intrinsic torque is proportional to the edge temperature gradient. 8 The production of intrinsic flow can be understood as a process similar to that of a heat engine. 8,9 In this process, temperature gradient, ∇T, excites turbulence, which not only relaxes ∇T but also drives a non-diffusive, residual stress via asymmetry in turbulence spectra k z k θ . 4,10 This residual stress then drives the parallel flow, converting the free energy in ∇T into kinetic energy of macroscopic flow.
As proposed in this heat engine model, the parallel residual stress Π Res r z is the key element that connects radial inhomogeneity to the macroscopic intrinsic flow. It is a component of parallel Reynolds stress, and is not proportional to either flow or flow shear. 4,10 The total parallel Reynolds stress can then be written as 4 The diffusive (− χ z ∂ r V z ) and pinch (V p V z ) terms are strict transport terms which cannot accelerate the plasma from rest. The divergence of this residual stress, −∇·Π Res r z , acts to drive the development of a sheared intrinsic flow via momentum transport. Flows with net momentum can then arise due to exchange of momentum with the surroundings which can give rise to an effective no-slip boundary condition. The residual stress depends on properties of underlying turbulence, and may flip sign when there is a change in the driving radial gradients of the equilibrium profiles.
Evidence for the role of parallel residual stress in driving intrinsic flow has been accumulating.
Probe measurements from the plasma boundary region of TJ-II stellarator confirm the existence of significant turbulent stress which provides a toroidal intrinsic torques. 11 A electrode biasing experiment on J-TEXT achieves a nearly zero toroidal rotation profile, and its results show that the intrinsic torque can be reasonably explained by the measured residual stress. 12 The residual stress profile has also been measured at edge of TEXTOR tokamak by canceling the toroidal rotation using counter-current NBI torque. 13 The observations demonstrate that there is a minimum value for the E r × B flow to trigger the residual stress, and that this stress scales with edge pressure gradient when the E r shear threshold is exceeded. Parallel flow driven by turbulent Reynolds stress has also been observed in a linear device, PANTA. 14,15 Recently, a gyrokinetic simulation predicts that residual stress profile exhibits a dipolar structure and provides the intrinsic torque which is consistent with measured rotation profile in DIII-D. 16 A number of theoretical models based on symmetry breaking in k-space have been proposed to explain the development of the residual stress. 4 In these models, the residual stress is determined by the correlator, k z k θ = k k z k θ φ k 2 / k φ k 2 , which is effectively set by the spatial structure of the k-spectra φ k (r) 2 . Theory suggests that the asymmetry in the k z space can result from the spatial variation of fluctuation intensity profiles, 17 or from the sheared E r × B flow that shifts modes off the resonant surfaces. 10 These mechanisms indicate that the residual stress is related to E r × B flow shear and turbulent intensity gradient, i.e., Π Res r z ∼ V E and Π Res r z ∼ I , respectively. These correlations are consistent with direct measurements from the edge of TEXTOR. 13 Despite these advances, our understanding of the microscopic mechanism is still rather limited.
Until now, there is no direct evidence validating the connection between the requisite symmetry breaking mechanism and the development of residual stress. Moreover, it is also unclear whether the residual stress can efficiently convert the free energy stored in the radial inhomogeneity into kinetic energy of the macroscopic parallel flow.
Due to its turbulence-driven origin, the axial flow must necessarily be coupled to the azimuthal mean flow. The latter is also known as zonal flow and is generated by drift wave turbulence via a modulational instability. 18 A theoretical framework 19 has been proposed to account for the interaction between these two secondary shear flows. However, how to precisely predict what the branching ratio between axial and azimuthal flows remains unknown. Therefore, further studies on how energy is distributed among the turbulence, azimuthal and axial mean flows are of interest.
The dominant branch will have a larger turbulent drive and set the turbulence level through a predator-prey type interaction with turbulent intensity field.
Besides the branching ratio question, the axial and azimuthal flows might also interact with each other directly. For a coupled drift-ion acoustic waves system, a zonal flow can arise from the parallel flow compression due to the effects of acoustic coupling. 20 Specially, when the parallel flow shear is strong enough to trigger parallel shear flow instability (PSFI), the enhanced fluctuating parallel flow compression can act as a source for zonal flow. This mechanism of zonal flow generation differs from conventional models which depend on the potential vorticity (PV) flux, and has not been tested experimentally. On the other hand, the axial flow shear may also be affected directly by its azimuthal counterpart. In the presence of a finite magnetic shear, the E r × B flow shear break parallel symmetry and generate a parallel residual stress Π Res r z , which accelerates the axial flow V z . The effects of azimuthal flows on axial flow generation at zero magnetic shear also remains unclear.
In this study, we discuss axial and azimuthal flow dynamics in CSDX, with a special emphasis on the possible flow interactions discussed above. We begin with a summary of our expectations based upon current theory-based modeling. We then report an experiments in a linear device, the Controlled Shear Decorrelation eXperiment (CSDX). 21,22 We show that the turbulent drive for the axial flow is less than that for the azimuthal flow by an order of magnitude. The turbulence fluctuation level is therefore regulated predominantly by the azimuthal flow shear. The results also show that the axial mean flow is driven by turbulent Reynolds stress. This stress, and particularly the non-diffusive, residual stress, results from a density gradient drive. In agreement with the recently developed dynamical symmetry breaking mechanism, 23 the residual stress emerges from drift wave turbulence with broken k θ − k z spectral symmetry. Note that this dynamical symmetry breaking model is also relevant to zero or weak magnetic shear case, e.g., in devices with straight magnetic fields and in flat-q regime tokamaks. The results presented in this paper validate the theoretical expectations for the link between the residual stress and symmetry breaking in the turbulence k-spectra, as well as the role of residual stress in converting thermodynamic free energy into kinetic energy of macroscopic axial flow.
The rest of the present paper is organized as follows. Section II recapitulates the theoretical background and predictions for turbulence-driven axial and azimuthal shear flows in CSDX. Section III introduces the experimental approach to measurements of mean flows and Reynolds stresses in CSDX. The experimental results and relevant discussions of theory-experiment comparisons are presented in Sections IV to VI, respectively. Section VII summaries the results and findings. In Section VIII, suggestions for future investigations are proposed.

II. THEORETICAL PREDICTIONS
In this section, we summarize theoretical predictions concerning the distribution of energy in the ecology of flows and fluctuations in CSDX. In order to investigate the evolution of turbulence and mean profiles in CSDX, we formulated a reduced model that describes the dynamics of the coupled drift-ion acoustic wave plasma. The model is derived from the Hasegawa-Wakatani system with axial flow evolution. 19 It self-consistently describes the variations in the mean profiles of density n, axial and azimuthal flows V z and V θ , as well as fluctuation energy ε = ñ 2 + (∇φ) 2 +ṽ 2 z . The convective derivative due to azimuthal rotation is neglected in the model since k θ v θ /ω k 1 at the location of maximum density gradient in CSDX. However, when ω k → k θ v θ , it could induce a wave-flow resonance, which mainly regulates the perpendicular (i.e., zonal) flows because k θ /k z 1 in systems with collisional drift turbulence. The potential effects of this wave-flow resonance has been discussed elsewhere. 24 The mean field equations are The quantities are normalized as follows: t ≡ t ω ci , v ≡ v /c s , and r ≡ r /ρ s , where ω ci is ion cyclotron frequency, c s is the ion sound speed, and ρ s is the ion Larmor radius at sound speed. The first terms on the RHS of Eqs. (1) to (3) represent the turbulent fluxes of particles and momentum, the terms that contain D c , ν c,⊥ and ν c, represent ion-ion collisional dissipations. In Eqs. (2) and (3), the terms proportional to the ion-neutral collision frequency ν in represent momentum transfer between ions and neutrals, and are significant only in the boundary region. In this study, the Reynolds powers, P Re z = −V z ∂ r ṽ rṽz and P Re θ = −V θ ∂ r ṽ rṽθ , are used to represent the rate of work done by the fluctuations to the mean flows.
In addition to the mean field equations, the evolution of fluctuation intensity ε = ñ 2 +(∇φ) 2 +ṽ 2 z is obtained as The first three terms on the RHS of the previous equation are mean field-fluctuation coupling terms.
They relate variations in ε to the evolution of the mean fields of n, V θ and V z . The energy exchange between fluctuations and mean profiles occurs via the particle flux ñṽ r , and the Reynolds stresses ṽ rṽθ and ṽ rṽz . In the energy equation, the ε 3/2 /l mix term represents energy dissipation by inverse cascade at a rate √ ε/l mix . Dissipated energy is ultimately damped by frictional drag. An energy source term P represents the excitation of drift wave turbulence, which is linear in ε and proportional to γ DW , i.e., P = γ DW ε. This is needed to incorporate turbulence excitation effects.
The flux Γ ε can be traced back to the nonlinear convective terms in the initial Hasegawa-Wakatani system.
Since the density response in CSDX is weakly non-adiabatic, we then calculate turbulent fluxes using quasilinear theory. In the near adiabatic limit, the expression for the particle flux is given by 25 Here D is the particle diffusion coefficient, and is equal to: ν ei and v T he are the electron-ion collision frequency and the electron thermal velocity, respectively.
In addition to the particle flux, an expression for the azimuthal momentum flux is needed. In the near adiabatic limit, and using quasi linear theory, the azimuthal momentum flux is equal to: The first term is the diffusive flux, while the second term is the residual component that accelerates the zonal flow from rest. The pinch term that arises from toroidal effects is neglected for the cylindrical geometry of the experiment. The turbulent viscosity and the residual stress are given as 19 6 In this study, the E r × B flow shearing rate is less than turbulence frequency, i.e., V E ω, so the term 1 ω−kV E x+iγ reduces to |γ| |ω| 2 . The azimuthal residual stress and χ θ thus decouple from azimuthal flow shear.
The axial Reynolds stress is given as 19 The non-diffusive component, i.e, the residual stress Π Res r z , drives the intrinsic axial flow, and is proportional to the correlator k θ k z . We thus write the following expressions for the parallel turbulent diffusivity χ z , and Π Res r z : Note that in order to obtain Π Res r z , we used the expressions for both electron drift frequency ω * e and eigenfrequency ω r = ω * e /(1 + k 2 ⊥ ρ 2 s ) in the adiabatic limit. Here, the axial residual stress and χ z also decouple from V E , since E r × B flow shearing rate is much less than drift wave turbulence frequency in CSDX.

Π Res
r z contains an expression for k θ k z , which is not easily determined within the scope of this simple, reduced model. To calculate the correlator, we need a spectral model considering the evolution of k θ k z ε , which can be obtained from wave momentum equations. This is beyond the scope of this work. Thus, what we offer here is an empirical approach that relates free energy source, ∇n, to the axial flow shear ∂ r V z . The correlator k θ k z is then expressed in terms of a coefficient that can be used in numerical studies, which is determined as follows. Proceeding in analogy with the treatment of turbulence in pipe flow, 26 the evolution of the fluctuating parallel ion flow is written as where c s denotes the sound speed,ṽ r is the eddy radial velocity,P is the pressure fluctuation, andφ is the potential fluctuation. In a drift wave system with adiabatic electrons like CSDX, one has eφ/T ∼ñ/n 0 andP/P 0 ∼ñ/n 0 as temperature fluctuations are small in this experiment. By introducing the radial mixing length l mix by the familiar relationñ/n 0 ∼ l mix |∇n|/n 0 , the fluctuating parallel flow then can be written as Here L z is the characteristic parallel dimension. The constant σ vT is introduced as a dimensionless scaling betweenṽ z and the density gradient ∇n. Multiplying byṽ r and ensemble averaging, the parallel Reynolds stress then becomes: While the first term represents a diagonal diffusive turbulent viscosity with the remaining part is the residual stress Π Res r z , proportional to ∇n. The coefficient σ vT is written as This coefficient captures the cross phase relation betweenṽ r andṽ z , and calibrates the efficiency of the density gradient in driving the residual stress Π Res r z . σ vT is also a measure of asymmetry in the spectral correlator k θ k z = k k z k θ φ k 2 / k φ k 2 , and encodes information concerning the parallel symmetry breaking that creates the residual parallel stress. An empirical value for σ vT , which can be used in the numerical solution of this model, can be obtained by a least-square fit to the experimental results.
Most of the conventional symmetry breaking mechanisms 4,17 are not applicable to plasmas with weak or zero magnetic shear, since they are usually associated with finite magnetic shears.
To resolve this issue, a dynamical symmetry breaking mechanism has been proposed to explain the development of intrinsic axial flow in absence of magnetic shear. 23 This mechanism does not require a specific magnetic field configuration, and thus it is valid for both finite shear and zero shear regimes. This mechanism is effectively equivalent to the modulational growth of a seed axial flow shear, as in zonal flow generation. In both cases, the initial breaking of symmetry is due to the seed flow.
The dynamical symmetry breaking model 23 was derived from a drift wave system with evolution of axial flow. The axial mean flow introduces a frequency shift to the growth rate of drift wave, i.e., In CSDX, electrons are weakly non-adiabatic, i.e.,ñ = (1 − iδ)φ. The adiabaticity of the electron response is measured by the dimensionless factor α ≡ k 2 where ω * e ≡ k θ ρ s c s /L n is the electron drift frequency. α is directly related to δ, i.e., δ As electrons approach the adiabatic limit, i.e., α → ∞ and δ → 0, drift wave is stabilized yielding γ k → 0. In CSDX, the adiabaticity factor is observed to be α 1, so electrons are weakly non-adiabatic, i.e., δ 1.
A test axial flow shear δV z , i.e., a perturbation to the mean axial flow profile, can break the symmetry of drift wave turbulence through the frequency shift 23 . The real frequency of the drift wave is affected by the test flow shear, and is given as The test flow shear also modifies the drift wave growth rate, which is given as The For a given δV z , the drift wave modes with k θ k z ρ s c s δV z > 0 have a larger frequency shift than the other modes. Thus, these modes grow faster. As a result, a spectral imbalance in the k z − k θ spectra is induced by the test flow shear. Such asymmetry in turbulence spectra can be detected by a joint probability density function of the turbulent velocities in both axial and azimuthal direction. The measurements of spectral imbalance are reported and linked to finite residual stress in this work.
The residual stress set by this dynamical symmetry breaking mechanism provides a negative definite contribution to the total turbulent diffusivity of axial momentum flux, i.e., Π Res r,z = − χ Res z δV z where χ Res z < 0. The negative momentum diffusivity induced by residual stress is Thus, the total Reynolds stress is This process of self-amplification of a test flow shear suggests that intrinsic axial flow can be generated through a modulational instability. When the magnitude of the negative viscosity exceeds the turbulent viscosity driven by drift wave, the total Reynolds stress induces a negative diffusion of axial momentum, thus amplifying the perturbation. In this case, the test shear (i.e., the modulation of mean flow shear profile) becomes unstable. The growth rate of test flow shear is γ q = q 2 r χ Res z − χ z , where q r is the radial mode number of flow shear modulation.
In CSDX, the seed shear is induced by the fact that RF heating is applied to one end of the plasma and, as a result, there is a modest pressure drop along the length of the machine that can drive a seed axial flow, particularly in conditions where the turbulent stress is small (i.e. at lower magnetic fields). Because the power deposition is radially dependent, the pressure drop is inhomogeneous in the radial direction. Thus, the source drives a radially sheared seed axial flow profile, i.e., δV z < 0. The seed shear breaks the spectral imbalance because it sets different growth rates for modes with different k θ k z . Modes that satisfy k θ k z δV z > 0 grow faster than the other modes. With δV z < 0 in CSDX, the saturated spectrum has a larger intensity in the domain where k θ k z < 0 than in the domain where k θ k z > 0, i.e. k z and k θ of dominant fluctuations will eventually become anti-correlated.
The onset threshold of axial flow generation is determined by the balance between residual stress and the turbulent diffusion driven by drift waves. Hence, the ∇n/n 0 threshold can be obtained from χ Res z = χ z . The turbulent viscosity driven by drift wave turbulence is calculated using χ z ∼ l 2 c τ c , where l c is the eddy correlation length and τ c is the eddy correlation time. The critical density gradient is then Using experimentally observed CSDX parameters, we can obtain ∇n crit ∼ 1.5 × 10 20 m −4 , which agrees with the experimental measurements presented below. Here, α = k 2 z v 2 The /ω * e ν ei ∼ 1 is the adiabaticity factor, the perpendicular turbulence scale length is k θ ρ s ∼ 1.5, and the eddy correlation time is τ c ∼ 6 × 10 −5 s.
The density gradient threshold can also be obtained by using the scaling coefficient σ vT of residual stress. The residual stress scales with ∇n as Π Res r,z ∼ σ vT l 2 c c 2 s /(L n L z ). Thus, σ vT is determined by the correlator k θ k z , i.e., σ vT = k θ k z / k 2 θ . Considering the symmetry breaking set by a test flow shear, we can calculate the correlator and thus the coefficient, as Thus, by using the balance between residual stress and turbulent diffusion, i.e., Π Res r,z = χ z δV z , we can also obtain the critical density gradient for onset of axial flow generation, which is the same as is small, such that the leading order of the residual stress is δΠ Res r,z ∼ | χ Res z |δV z . Thus, the axial flow saturates when χ z = | χ Res z |. Ultimately, the flow energy is dissipated by viscous heating and drag dissipation.
In summary, for regimes of moderate azimuthal shear (i.e., |V θ | ω k ), theory predicts that: (1) drift wave fluctuations and azimuthal (i.e., zonal) flows will form a self-regulating system; (2) axial flows will evolve parasitically by Reynolds stress, on the existing drift wave-zonal flow turbulence. Here, the key point is Π Res r,θ Π Res r,z , as k ⊥ k z ; (3) symmetry breaking in the k θ -k z space is required for axial flow generation.
(4) Sheared intrinsic axial flows will be generated when the density gradient exceeds a predicted critical value. Now, we turn to tests of these predictions.

III. EXPERIMENTAL SETUP
In this section, we present the experimental methodology for testing the predictions of model in Section II. The experiments were conducted on the Controlled Shear Decorrelation eXperiment (CSDX), a linear plasma device with an overall length of 2.8 m and a diameter of 0.2 m (Fig. 1).
The working gas was argon at a gas fill pressure of 1.8 mTorr. The argon plasma was produced by a 15 cm diameter 13.56 MHz RF helicon wave source via an m = +1 helical antenna that surrounds a glass bell-jar, and was terminated by insulating end-plates at both ends. The uniform magnetic field is in the axial direction (denoted as the −ẑ direction). In this study 1800W of power was used, and the magnetic field strength was varied from 500 G to 1000 G. A higher magnetic field results in a steepening of the density profile in CSDX. 21,22 Typical plasma parameters are as follows: the peak on-axis electron density of n e ∼ 1 × 10 19 m −3 , the electron temperature of T e ∼ 3 − 5 eV, and the ion temperature of T i ∼ 0.3 − 0.8 eV. More details on this device can be found in previous publications. 21,22,27 A horizontal scanning probe was used to record basic plasma information such as ion saturation currents and floating potentials at port α that is about 1 m downstream from the helicon source. The probe array is a combination of Mach and Langmuir probes and is capable of measuring the axial and radial plasma velocities simultaneously (Fig. 2). The axial velocity, v z , was measured by a Mach give spuriously large axial flows 22 which were later found to be inconsistent with laser-induced fluorescence (LIF). 27 This overestimation of parallel Mach number is found to be related to shadowing effects in Mach probes. 30 In this study, we used small enough tips (D probe ≈ 3 mm) to avoid probe shadowing effects and we verified that the mean flow profile measured by the Mach probe agreed with LIF measurements of the same ion flow taken in the same plasma conditions. 27 The fluctuating E × B velocities are estimated from the floating potential gradients between two adjacent tips (∇φ f ), i.e.,ṽ r = −∇ θφf /B andṽ θ = ∇ rφf /B. The distance between two adjacent floating potential tips is about 3 mm. The sampling rate of the probe data is f s = 500 kHz which gives a Nyquist frequency that is well above the frequency of the observed dominant fluctuations

A. Enhanced Shear Flows
In this study, we obtained different equilibrium profiles and fluctuation intensities by changing the magnetic field strength B. As shown in Fig. 3(a), when the B field is raised, the plasma density and its radial gradient increases. During the B scan, the variation in electron temperature is negligible. The axial velocity reverses at edge, and its radial shear increases with increasing B field ( Fig. 3(b)). The axial Reynolds stress, ṽ zṽr (Fig. 3(c)), is estimated using velocity fluctuations in the frequency range of 5 < f < 30 kHz; previous studies have identified these as collisional drift wave fluctuations. 21,22 ṽ zṽr is negligible for r < 3 cm at lower B field, but becomes substantially negative at higher B field (Fig. 3(c)). The Reynolds force, F Re z = −∂ r ṽ zṽr (Fig. 3(d)), increases significantly in the core, and becomes more negative at the edge (3 < r < 6 cm). This negative turbulent force at the edge appears to be matched with the reversed axial mean flow. The parallel Reynolds force is much larger than the force on the ions arising from the parallel electric field. This weak electric field arises from the Boltzmann equilibrium associated with the electron pressure drop along the axial direction (Fig. 3(e)). Thus, the axial shear flow in CSDX reported here is primarily driven by the turbulent Reynolds force. In addition to the evolution of the axial flow, the changes in azimuthal flow have also been measured using a Mach probe during the B scan. As can be seen from Fig. 4(a), the mean azimuthal velocity, V θ , propagates in the electron diamagnetic drift direction (EDD), which is negative in the figure. The magnitude of V θ increases by a factor of two when B is raised from 500 G to 800 G. The azimuthal Reynolds stress, ṽ rṽθ , is also estimated using fluctuations in the frequency range of 5 < f < 30 kHz. ṽ rṽθ is small and flat at lower B, but its magnitude increases when B is increased (Fig. 4(b)). The change in ṽ rṽθ gives rise to substantial turbulent Reynolds force, F Re θ = −∂ r ṽ rṽθ (Fig. 4(c)) at higher B. The generation of sheared azimuthal E × B flow via the turbulent Reynolds stress has been reported in previous studies in CSDX, 27,31,32

B. Axial Force Balance Analysis
To confirm the role of the axial Reynolds force in driving the axial flow, we examine the force balance in axial direction. The azimuthal force balance has been performed in previous studies. 31,33,34 Here, we carry out similar analysis on the axial flow. The axial ion momentum equation is written as where the ion viscosity µ ii = 6 5 ρ 2 i ν ii ∼ 5 − 10 m 2 /s and ion-neutral collision frequency ν in = n gas v ti σ in ∼ 3 − 6 × 10 3 s −1 are estimated from previous studies. 31 µ ii and ν in are likely to have weak spatial variations, i.e., µ ii ∝ nT −1/2 i and ν in ∝ T −1/2 i . Here, we assume the neutral pressure is radially uniform and the neutral temperature is approximated by the ion temperature profile, which has been measured using LIF techniques in previous studies. 27 A no-slip boundary condition is also imposed, justified by strong ion-neutral damping at edge, i.e., V z → 0 at r = 6 cm. Taking the measured profiles of the Reynolds stress and the axial pressure gradient shown in Fig. 3, we can then solve Eq. (17) for V z using a finite difference method. The axial pressure force can also be ignored at higher B field, since it is smaller than turbulence force by a factor of 5. As shown in Fig. 5, the calculated results (curves) are in agreement with the mean axial ion flow profiles measured by the Mach probe (circles). This results confirms that the turbulent stress is responsible for the increased V z and more pronounced flow reversal found at higher magnetic field.

A. Turbulent Flow Drive Scales with Density Gradient
The magnetic field scan yields a clear rise in ∇n, which is much larger than ∇T e and has been identified in previous work as the primary free energy source driving the fluctuations. 21,22 This change presents us an opportunity to determine the link between ∇n, the turbulent drive, and the macroscopic intrinsic flow. In this study, we did a shot-by-shot B field scan, and used the Reynolds power, P Re z = − V z ∂ r ṽ zṽr , to represent the rate of work performed by the turbulent fluctuations on the mean axial flow. The axial shear flow and the Reynolds power are plotted as a function of ∇n (Fig. 6). The magnitude of axial flow shearing rate, |V z | = |∂ r V z |, increases sharply when the density gradient exceeds a critical value, ∇n e 1.6 × 10 20 m −4 ( Fig. 6(a)). This critical density gradient is in agreement with the theoretical prediction shown in Eq. (15). Concurrently, the Reynolds power also increases substantially when this threshold is exceeded (Fig. 6(b)). Here, we used volumed-averaged Reynolds power, P av z = ∫ − V z ∂ r ṽ zṽr r dr/ ∫ r dr where 1 < r < 5 cm.
These observations show that the axial shear flow and its Reynolds power increase consistently as ∇n increases, indicating that the turbulence acts as a converter, transferring the free energy to the intrinsic flow. These results are consistent with the heat engine model. 9 Here, the free energy due to ∇n is converted into kinetic energy of macroscopic parallel flow.
The azimuthal flow and its turbulent drive are also driven by the density gradient. Similar to the analysis of the axial flow case, we use the azimuthal Reynolds power, P Re θ = − V θ ∂ r ṽ rṽθ , to represent the nonlinear kinetic energy transfer into the mean azimuthal flow. We then plot the axial flow shear and azimuthal Reynolds power as a function of the density gradient. As shown in Fig. 6(c), there is a clear threshold effect in the density gradient, which is the same as the axial flow case. After the threshold, the azimuthal flow shear, V θ = |∂ r V θ − V θ /r |, and the azimuthal Reynolds power, P Re θ , increase with the density gradient ∇n (Fig. 6(d)). The similar trends of V θ and P Re θ suggest that the underlying turbulence also converts the free energy from the density gradient into kinetic energy of azimuthal mean flow.
The results above show that both the axial and azimuthal mean flows are turbulence-driven in CSDX. However, the nonlinear kinetic energy transfer to the two secondary shear flows are not equally distributed. The axial Reynolds power is smaller than the azimuthal one by an order of magnitude, i.e., P Re z P Re θ , since k z k ⊥ for turbulent fluctuations in CSDX. Therefore, we  conclude that the azimuthal shear flow sets the turbulent fluctuation level through predator-prey type interaction, while the axial flow evolves in this intensity field. The disparate magnitudes of nonlinear energy transfer also suggest that there is no significant direct energy exchange between axial and azimuthal shear flows. The axial flow is then parasitic to the turbulence-zonal flow system, and is driven by the turbulent Reynolds stress, especially the non-diffusive, residual stress.
The weak axial to azimuthal flow coupling allows us then to simplify the 4-field model in Section II to a 2-field predator-prey model.

B. Residual Stress Driven by Density Gradient
As discussed in Section II, it is the residual stress that converts the thermodynamic free energy to the kinetic energy of the axial mean flow. 4,23 The residual stress can be synthesized from the measured total Reynolds stress (Fig. 3(c)) and the diffusive stress inferred from experimental measurements, 32 i.e., Π Res r z = ṽ rṽz + χ z ∂ r V z with the diffusivity χ z = ṽ 2 r τ c expressed in terms of the measured eddy radial velocityṽ r and eddy correlation time τ c . Here, the pinch term (V p V z ) is ignored, since it arises from toroidal effects and thus is not significant in a linear device. As shown in Fig. 7, the magnitude of the synthesized residual stress increases as the B field, as well as ∇n, is increased. The magnitude of the residual stress, Π Res r z , is then plotted against the normalized density gradient in Fig. 8. At smaller density gradient, the magnitude of residual stress, Π Res r z , is small, and is almost independent of the normalized density gradient. At larger ∇n, Π Res r z increases in proportion to the normalized density gradient, with a slope σ vT ≈ 0.10. Here, Π Res r z is volumeaveraged in the range of 1 < r < 5 cm. This finding is consistent with the hypothesis that the residual stress is driven by the density gradient. Also, a finite σ vT ≈ 0.1 indicates the existence of a k θ − k z symmetry breaking mechanism at higher ∇n.

TURBULENCE SPECTRA
The development of residual stress is also proposed to be correlated with symmetry breaking in k-space, 4 i.e., k z k θ = k k z k θ φ k 2 / k φ k 2 0. The symmetry breaking can be assessed by investigating the joint probability density function (PDF) of radial and axial velocity fluctuations, P (ṽ r ,ṽ z ). Note that in CSDX we haveṽ z ∼ ∇ P ∼ k zφ andṽ r ∼ k θφ , due to the adiabatic electron response and negligible temperature fluctuations. By normalizing the velocity fluctuations using their standard deviations, P (ṽ r ,ṽ z ) can represent the correlator k z k θ . As shown in Fig. 9, the anisotropy of P (ṽ r ,ṽ z ) grows with increasing B field strength and ∇n. The critical density gradient occurs at B ≈ 650 G, and P (ṽ r ,ṽ z ) starts to tilt ( Fig. 9(b)) at slightly higher B and ∇n. At higher ∇n, P (ṽ r ,ṽ z ) is strongly elongated along the diagonal, suggesting large asymmetry in k z k θ . As proposed by the dynamical symmetry breaking model, 23 the mean axial flow shear modifies the drift wave growth rate, by introducing a frequency shift proportional to k z k θ V z . In our experiments, the seed axial flow shear is negative, V z < 0, because V z (r) is initially driven by the axial pressure drop and hence decreases from the core to the edge. As a result, the modes with k z k θ < 0 grow faster than modes with k z k θ > 0, and eventually become dominant. This in turn induces a spectral imbalance, with predominance of the spectral intensity in quadrants II and IV of the k θ − k z plane, as shown in the right panel of Fig. 10. The predicted spectral imbalance, k θ k z < 0, is consistent with the tilted contour of P (ṽ r ,ṽ z ), as shown in left panel of Fig. 10.
Since larger residual stress occurs at higher ∇n, we can therefore infer that this symmetry breaking is related to a finite residual stress.

VII. CONCLUSIONS
In this work, we study axial and azimuthal flow dynamics in drift wave turbulence in CSDX.
We focus on possible interactions between azimuthal and axial flows. The principal results of this study are: • Turbulent azimuthal Reynolds stresses ṽ rṽθ drive zonal flows which regulate the turbulence.
• Turbulent axial Reynolds stresses ṽ rṽz drive axial flows-akin to intrinsic rotation. However, the azimuthal Reynolds power is much larger than the axial Reynolds power, i.e. P Re θ P Re z , so one may regard the axial flow evolution as parasitic to the drift wave-zonal flow system.
• Spectral symmetry breaking was observed and measured-i.e., k θ k z 0. The observed broken symmetry is consistent with that required for axial flow generation. The symmetry breaking is dynamical, and is not produced by magnetic field geometry.
• Azimuthal and axial flows as well as the symmetry breaking scale with ∇n, consistent with the scenario of the engine model of the system.
• Experimental results support the predictions of the reduced model discussed in this paper.
We emphasize that conclusions pertinent to azimuthal-axial flow coupling are limited to magnetic field in range from 500 G to 1000 G. In this range, V θ ω k and L −1 , which are fundamental to the system dynamics observed and modeled here.

VIII. AXIAL-AZIMUTHAL FLOW INTERACTION-A FUTURE DIRECTION
A plausible physical picture of the system of flows and turbulence discussed in this paper is summarized in Fig. 11. In this study, the axial Reynolds power is smaller than the azimuthal one by an order of magnitude. Thus, the azimuthal flow-turbulence interaction is the primary branch in the turbulence-flow system. The axial mean flow is then parasitic to such system, and is driven by the residual stress. The azimuthal flow shearing rate is much less than the drift wave frequency, so the residual stress decouples from the effect of azimuthal flow (dashed line in Fig. 11). This axial residual stress results from a dynamical symmetry breaking mechanism, i.e., driven by drift wave turbulence with broken symmetry in k-space. This spectral imbalance in k z k θ is induced by the seed axial flow shear, which is in turn amplified by the axial residual stress. These observations are consistent with the causal link proposed by the heat engine model, i.e., a pathway from symmetry breaking to the development of residual stress and the onset of axial mean flow. Although the axial-azimuthal flow coupling appears to be weak in this study, it needs not always be so. There are at least two ways to enhance the interaction between axial and azimuthal flows in CSDX. The proposed mechanisms are illustrated in Fig. 12. One way is to increase the power of the plasma source, such that ∇n drives stronger drift wave turbulence and thus leads to enhanced zonal flows via the Reynolds force. When the zonal flow shear is comparable to drift wave frequency, it will regulate the axial flow production and dissipation by entering explicitly-and reducing-the axial residual stress and turbulent diffusivity. The enhanced zonal flow shear will then increase the axial flow shear by reducing the cross-field momentum transport, i.e., thus forming a transport barrier.
The other way to enhance the coupling between axial and azimuthal flows is to increase the parallel momentum source. The enhanced axial flow can increase the zonal flow production via the acoustic coupling. 20 The parallel flow compression can be converted to zonal flow by coupling with potential vorticity (PV) fluctuations. This coupling, i.e., q∇ ṽ , breaks PV conservation, and thus forms a source for zonal flow. This conversion occurs when parallel flow compression is significant, especially near the PSFI threshold. With increased axial and azimuthal flow shears, a transport barrier can be formed by increasing the axial momentum source. CSDX will be equipped with an axial gas-puff system that provides an axial momentum source. The axial flow then can also be driven by a strong axial momentum source, and thus V z would be adjustable within a wide range. In our current experiments, the peak value of the axial Mach number is about 0.2, which is well below the PSFI threshold. The upgraded system will present us an opportunity to investigate the role of PSFI in parallel flow saturation as well as axial-azimuthal flow coupling.
In conclusion, we remark that CSDX offers an excellent venue to study the detailed physics of transport barrier formation with turbulent-driven transverse and parallel shear flows at zero magnetic shear. In tokamaks, it has been observed that coexistence of large toroidal rotation and low magnetic shear, i.e., flat-q regime, leads to enhanced confinement states, and profile "destiffening". 35 This regime is under intensive study in the magnetic fusion energy community, and it is worthwhile to note that basic experiments can produce substantial insights into the relevant physics.