Evidence and modeling of turbulence bifurcation in L-mode confinement transitions on Alcator C-Mod

Analysis and modeling of rotation reversal hysteresis experiments show that a single turbulent bifurcation is responsible for the Linear to Saturated Ohmic Conﬁnement (LOC/SOC) transition and concomitant intrinsic rotation reversal on Alcator C-Mod. Plasmas on either side of the reversal exhibit diﬀerent toroidal rotation proﬁles and therefore diﬀerent turbulence characteristics despite proﬁles of density and temperature which are indistinguishable within measurement uncertainty. Elements of this bifurcation are also shown to persist for auxiliary heated L-modes. The deactivation of subdominant (in linear growth rate and contribution to heat transport) ion temperature gradient (ITG) and trapped electron mode (TEM) instabilities is identiﬁed as the only possible change in turbulence within a reduced quasilinear transport model across the reversal which is consistent with the measured proﬁles and inferred heat and particle ﬂuxes. Experimental constraints on a possible change from strong to weak turbulence, outside the description of the quasilinear model, are also discussed. These results indicate an explanation for the LOC/SOC transition that provides a mechanism for the hysteresis through the dynamics of subdominant modes and changes in their relative populations, and does not involve a change in the most linearly unstable ion-scale drift-wave instability.


I. INTRODUCTION
Confined plasmas are active, turbulent nonlinear systems which demonstrate a complex dependency between external actuation and the plasma response. Of prime importance to magnetically confined fusion plasmas is the energy confinement time, the ratio of the plasma stored energy (the plasma response) to the applied heating power (the external actuation), which demonstrates a dependency on plasma parameters that cannot be described with a single scaling law 1,2 . Thus, to understand the behavior of fusion plasmas, it is necessary to understand the dynamics which govern the plasma response, many of which are internal to the plasma, and how they could lead to different regimes of plasma behavior.
This work focuses on two transitions observed universally in tokamak plasmas: the Linear to Saturated Ohmic Confinement (LOC/SOC) transition and the intrinsic rotation reversal.
The LOC/SOC transition refers to a break in the slope of the scaling of energy confinement with respect to density in Ohmic L-mode plasmas. There has been extensive experimental work characterizing the confinement transition and its link to L-mode density scaling 3 . The intrinsic rotation reversal refers to a spontaneous reorganization of the plasma rotation in the absence of external momentum input. Typically the rotation profiles transition from mostly flat co-current rotation profiles to hollow counter-current rotation profiles. In response to ramps of plasma density, these two transitions are found to occur at the same critical density, the reason for which is unknown. This critical density is seen to correlate with a critical collisionality ν * = (ν ee /ǫ)/ω be , where ω be is the bounce time of electron banana orbits caused by trapping in the magnetic well, and ν ee /ǫ corresponds to the rate of electron detrapping due to collisions. There is also a wealth of other phenomenology associated with the transitions, including density profile peaking, non-diffusive cold pulse propagation, fluctuation spectrum changes, and impurity density asymmetries [4][5][6] . The concurrence of these two transitions suggests a link between the heat, particle, and momentum transport channels in tokamak core plasmas, reminiscent of how heat-flux driven turbulent Reynolds stresses are thought to trigger the L-H transition, and lead to shear layer formation. While there has been success in modeling the LOC/SOC transition in silico 7-9 , a selfcontained physical picture of how changes in turbulence lead to the concurrence of these transitions and explain its observed properties remains elusive. Since drift wave turbulence is responsible for most of the heat, particle, and momentum transport in tokamaks, it has 3 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE Here, we provide evidence the two transitions are linked to a single nonlinear bifurcation of the plasma state in Alcator C-Mod, advancing the notion that LOC and SOC are representative of different states of nonlinearly saturated turbulence. We also show how elements of the bifurcation persist in auxiliary heated plasmas, suggesting its importance for understanding L-mode scaling in more generality. Additionally we propose to identify the bifurcation as a partial turbulence population collapse, which involves the deactivation of a We then conclude and discuss open questions and paths forward for further inquiry.

II. ROLE OF TEM/ITG IN LOC/SOC AND ROTATION REVERSAL
The conventional argument for why LOC/SOC is observed as the density increases is as follows: Ohmic heating primarily deposits power into the electrons. At low collisionality, the ions and electrons are only weakly thermally coupled, and TEMs are particularly virulent, leading to poor confinement. As the density increases, the collisionality serves both to improve the coupling of ions and electrons, and to decrease the growth rate of TEMs through collisional detrapping of electrons. Main ion dilution 1 − (n D /n e ) also tends to decrease and T e /T i approaches 1 with increasing plasma density, leading to increased ITG growth rates.
The TEM/ITG transition picture supposes that once the LOC/SOC transition density is reached, ITGs overtake TEMs as the key contributor to the confinement scaling. Since increasing collisionality does not have as large of an effect on the ITG drive, which tends to pin ion temperature gradients near marginally stability, the confinement scaling saturates.
This argument appears to agree with reduced modeling of LOC/SOC both using analytic estimates of transport 7 and using TGLF 8,9 . The roles of TEM and ITG instabilities are also implicated in producing the observed change in core temperature response to edge cold 5 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444 ies that the QC-TEM feature depended not just on the presence of TEM, but also on the saturated state of the turbulence. Changes in fluctuations have also been seen on Alcator C-Mod. Here we will demonstrate how changes in turbulent fluctuation are seen at experimentally indistinguishable density and temperature profiles, and that circumstantial evidence supports an interpretation of these changes reflecting changes in the underlying turbulence.

III. CHARACTERIZATION OF ROTATION REVERSAL HYSTERESIS ON ALCATOR C-MOD
The experiments presented here were run on Alcator C-Mod, a compact (R = 0.67m, a = 21cm), high-field (B t up to 8.1 T) diverted tokamak with a molybdenum wall 32 . Two methods were used to realize the rotation reversal hysteresis: density control for Ohmic plasmas, and ion cyclotron resonance frequency (ICRF) heating power modulation for auxiliary heated plasmas. The discharges were upper single null, with B t = 5.4T. The three experimental cases performed are shown in Table I. Case I has been reported in previous work 15 . This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. profiles during the reversal evolved on timescales ≈ 4τ e , slower than the modulation period.

PLEASE CITE THIS ARTICLE AS
There was no beam injection in these plasmas, so the only particle source is at the edge.
Additionally, on some discharges a perturbed laser blow-off injection of CaF 2 impurity was performed to assess the robustness of the hysteresis to perturbation. Time traces from a representative Case I discharge are shown in Figure 1, where time points with matching line-integrated densities, but differing toroidal rotation, have been marked.
The coexistence of multiple rotation states at a given density can be visualized by plotting the line-average density versus the line-average toroidal velocity of the discharge as a hysteresis curve, as is done in Figure 2. The plotted trajectories overlay closely, showing the robustness of the hysteresis phenomenon to noise and perturbation from LBO. Since the 8 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444
edge rotation is similar for the two rotation states, the change in the core rotation profiles indicates a change in the turbulently generated residual stress, implying that some aspect of the turbulence, yet to be determined, changes between the two states. This shows there exist two distinct states of turbulence, one exhibiting co-current rotation and the other countercurrent rotation, which can stably exist in this range of plasma density. As the density is slowly changed, the co-current rotation states continuously evolve to lower-collisionality LOC states, which corresponds to the 'LOC-like' branch, while the counter-current states continuously evolve to higher-collisionality SOC states, which correspond to the 'SOC-like' branch. The rotation states between the two branches are transient, and do not correspond to steady states of the plasma. The transition densities were found to match those from density ramp experiments, independent of ramp rate. Thus, it makes sense to associate plasmas which exhibit co-current rotation with 'LOC-like' turbulence, and plasmas which exhibit counter-current rotation with 'SOC-like' turbulence. In addition, a hysteresis cycle in core electron temperature versus core line-average toroidal velocity is also plotted for Case III, which demonstrates that the hysteresis in toroidal rotation persists for auxiliary heated L-mode plasmas, where the hysteresis is effected by ICRF heating modulation rather than by density modulation.
Kinetic profiles matched across the reversal for all three cases are shown in Figure 3. The different colors in each plot correspond to different discharges. The discharges all demonstrate clearly separated rotation states at the same density for the ohmic cases, and core electron temperature for the ICRF heated case. These states have been highlighted in red and blue. Note that discharges in cyan and magenta were perturbed using laser blow-off injections. The discharge in the left frame plotted in cyan does not reach the required density to transition from the LOC-like branch to the SOC-like branch, so it does not complete a full hysteresis loop.
shape is used. Note the priors differ from those used for the previously reported profiles, as they were changed to better capture the flattening effect of sawteeth inside the sawtooth mixing radius. For all three cases, the profiles and gradients responsible for determining the linear stability characteristics of the drift-wave turbulence overlay each other very closely.
This confirms previous analysis that drift-wave stability is not responsible for the reversal.
Despite the nearly identical kinetic profiles, differences in fluctuations can be observed in the 1.1 MA Ohmic and 0.8 MA RF-heated cases. Data are presented here from the Phase Contrast Imaging system on C-Mod, which is sensitive to line-integrated density fluctuations with wavevector component in the major radial direction of k R < 30cm −1 , corresponding to k R ρ s 2.5. In these cases, the PCI picks up high-frequency "wing" features in the LOC case, as has been previously reported 10,11,35,36 . This is shown by the PCI spectra in Figure   4, and demonstrates measured fluctuation spectra can change even at matched conditions.
One of the characteristic features observed about these PCI wings is an asymmetry in the ±k R spectrum. The time evolution of a metric capturing this asymmetry is also plotted in  Since PCI measurements are line-integrated, interpreting changes in PCI spectra can be 12 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444 challenging. The PCI sight lines are vertical, and intersect magnetic flux surfaces twice, once towards the top of the tokamak and once towards the bottom. The sign of the projection of k θ on to k R changes between these two intersections, so a ±k R asymmetry is usually interpreted as an up-down fluctuation asymmetry. Note that k R could be either larger or smaller than k θ , depending on the angle of incidence to the flux surface and the radial wavevector k r .
Although it is difficult to precisely localize the radial location of the measured fluctuations, rate is plotted against the E × B shearing rate calculated from force balance for co-and counter-current rotation profiles in Figure 5.
To better understand the impact of the modes on transport and diagnose the turbulent state of the plasma, we adopt a quasilinear transport approximation QLTA. This is the model used in the previous analysis 15 , and will be reintroduced here. In this QLTA, the turbulent fluxes (e.g. the electron heat flux Q e ) are expressed as the sum of a quasilinear mode weight (for electron heat flux, this will be denoted W Qe,k ) times an averaged mode intensity (also called spectral weight, in this case φ2 k ) for each linear eigenmode indexed 14 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444 by wavenumber k. For example, writing this explicitly for the electron heat flux: The resulting anomalous fluxes are plotted in Figure 6. If the spectrum is discretized into N modes, then equation (1) and the corresponding equations for the other fluxes can be viewed as a 3 × N matrix equation, as in equation (2): In the full system, N will be very large or infinite, so a naive application of the flux constraint leaves the mode intensities highly underdetermined. However, for the purpose of broadly determining trends in transport of the plasma, we do not need the detailed shape of the mode 15 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444 spectrum. For example, at electron scales corresponding to k y ρ e 1, electron temperature gradient (ETG) modes will exhaust primarily electron heat flux without exhausting ion heat or particle flux, since ions will behave adiabatically at that scale. Thus the net effect of ETGs on the turbulent fluxes can be well-approximated using only one degree of freedom, instead of using one degree of freedom per mode. This motivates the construction of a reduced transport model where similar modes are lumped together into 'families'. To calculate the total turbulent fluxes, the sum over modes is replaced with a sum over families, using quasilinear weights and spectral weights averaged over the modes in the family. Additionally, here an explicit spectral width ∆k has been included, which captures the number of modes in the original sum that the family replaces. This is represented in matrix form as in the second line of equation (2), where the families are enumerated by Roman numerals I, II, III, ... Note that the product φ2 I ∆k I represents a family-integrated spectral weight, while φ2 I alone is the family-averaged spectral weight. Previously, this reduced family model was applied to the matched profile cases to identify a candidate subdominant mode transition, where the deactivation of an intermediate scale k y ρ s 1 instability was found to be consistent with the observed transport. Here, a refinement of that analysis is presented. Several radial locations were chosen for analysis, starting with the location of maximum ω ′ tor (r), and proceeding radially outward until the points were outside the rotation reversal region. Quasilinear weights calculated from CGYRO are shown in Figure 6, for Case I and Case II. Despite being at different currents, the quasilinear weights look very similar between the two cases, with key features that are present in all analyzed locations: • At ion scales k y ρ s 1, the modes primarily exhaust ion heat flux, although there is a significant electron heat flux component as well. These modes are also thought to be responsible for most of the momentum transport in the plasma, as momentum will primarily be carried by ions. The electron particle transport starts outward at low k, but becomes increasingly inward directed as higher k are approached. For locations within the rotation reversal region, this inward trend is strong enough to reverse the direction of the particle flux. This ITG particle flux trend can be understood through a kinetic picture 47 . Collisions and diffusion cause outward particle fluxes, as is the case at low k. Ion finite larmor radius effects push the phase velocity of higher k modes closer to zero, causing the modes to be resonant with lower energy 16 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444
particles. Since temperature gradients produce an energy dependence in the radial gradient of the background Maxwellian distribution function which is of opposite sign to that produced by a density gradient at low energies, modes which are resonant with lower energy particles produce a stronger inward thermally-driven particle flux. These changes manifest as a change in the phase relation between density and potential fluctuations, leading to different quasilinear weights.
• At intermediate scales k y ρ s ∼ 1, both ITG and a combined TEM/ETG branch are present. Occasionally ubiquitous modes 48 (modes with real frequency ≈ 0) can be observed in linear simulations at this scale. ITG and ubiquitous modes present at this scale tend to have strongly inward particle flux. TEMs which are active at this scale also have inward particle fluxes.
• At electron scales k y ρ s ≫ 1, only a combined TEM/ETG branch is present. These modes are characterized by nearly adiabatic ions due to ion gyroaveraging over the small scales, and hence have weak ion heat and momentum flux, as well as weak particle flux due to the ambipolar nature of transport in axisymmetric systems 49 .
Additionally, it was found that electromagnetic perturbations contributed < 5% to the quasilinear weights, so the transport is dominantly driven by electrostatic turbulence. Note that these weights were not found to be nearly identical between the matched LOC and SOC states in previous work 15 .
Proceeding with the identification of families using the quasilinear weights, here we classify modes into dominant or subdominant families by the following criterion: At the location of maximum ω ′ tor (r), calculate the ratio α ≡ γ E /γ max . γ max is the maximum ion scale linear growth rate, and γ E = r q ∂ω0 ∂r is the shearing rate from ω 0 , the E × B rotation rate calculated using force balance. Using the shearing rates right after the rotation reversal results in α of TEMb k y r b 1 ∇T e driven TEM, where r b is the typical banana orbit width. Also always exhibits inward particle flux, although is particularly strong outside the rotation reversal region. Some metrics of subdominance, such as comparing γ/k 2 ⊥ , find all TEMs to be subdominant and so do not have a TEMb family.
ETG k y r b 1 ETG, exhausts primarily electron heat flux.
grouped into a single family, allowing the identification of six families, shown in Table II, with their main properties summarized. The naming scheme is somewhat arbitrary, as for example there is no discrete delineation between TEM and ETG. Additionally, note that this classification based on the linear growth rate is not the only possible way to group the modes into families, and consideration of other mechanisms including nonlinear mode coupling may suggest different classifications.
With these families identified, the additional constraint that mode intensities must be non-negative distinguishes two qualitatively different solutions consistent with the transport constraints imposed by equation (2), based on whether or not the subdominant modes are active. Since neutral penetration from the edge is low and there is no beam fueling on Alcator C-Mod, the particle flux in the core is nearly zero. The neoclassical particle flux is not large enough to balance a significant turbulent particle flux, so the active modes must have a roughly balanced net particle flux. So, one possible solution is that in the rotation reversal region, only the dominant modes ITGb and TEMb/ETG are active, and the subdominant modes are inactive. The particle flux constraint is satisfied primarily within the modes of ITGb, possibly with a small additional contribution from TEMb. This is possible because the particle flux weight of these modes averages to nearly zero for weights in the rotation reversal region. This leads to a narrower k spectrum of ion-scale turbulence, with electron 18 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444 direction turbulence quenched at ion-scales, corresponding to a "SOC-like" regime. Note that at large enough r/a, an increasingly large outward particle flux is driven by ITG due to increasing collisionality, so it is not possible for the particle balance constraint to be satisfied within ITGb alone. This necessitates significant activity from TEMa/b at these radii. Another possibility is to have both dominant and subdominant modes active. For this other solution, the subdominant ITGa is active, so either or both of ITGc and TEMa/b must have significant activity in order to balance the net particle flux. This leads to a broader k spectrum, with possibly intermixed ion-and electron-direction turbulence at ion-scales, corresponding to a "LOC-like" regime.
Since the "SOC-like" and "LOC-like" regimes continuously connect to the higher collisionality counter-current rotation and lower collisionality co-current rotation branches of the turbulent bifurcation respectively, it is inferred that the rotation physics of ITG dominant or TEM dominant (but possibly ITG active) turbulence continues to be relevant closer to the transition in their respective regimes. The changes suggested by this model are illustrated schematically in Figure 7, and their consequences summarized in Table III. One key implication is that despite heat transport being the principal externally driven means of turbulent free energy release in the plasma, modes which are subdominant in heat transport can play a role in determining the behavior of transport in other channels, such as particle or momentum transport.

19
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Although not systematically ruled out, the robustness of the hysteresis to perturbation from LBO injection also suggests that the crossing of an undetected linear stability boundary is not responsible for the transition. This is due to the large effect of LBO on the linear stability of the plasma 17,18 . Additionally, while profile shear has been identified as being important to intrinsic rotation 29,51 , these experiments suggest that changing profile shear is not responsible for the transition itself. However, the experiments do not rule out profile shear being responsible for setting the final steady-state shape of the rotation profile.

PLEASE CITE THIS ARTICLE AS
Local mean-field mechanisms for bistability are not entirely ruled out by the profiles, as change of the mean rotation profile could feedback on the plasma turbulence. The E × B flow shear is not large enough to entirely quench the ion-scale turbulence, but the shearing rate is a significant fraction of the maximum linear growth rate at ion-scales, as shown in 20 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444 One such change is suggested by Figure 8, where the perpendicular group velocity in the flux surface v gr = ∂ω(ky) ∂ky is compared to the phase velocity in the same direction v ph = ω(ky) ky . Visually, v ph ≈ v gr for an intermediate scale k y ρ s from 2-5, implying waves in that range are weakly dispersive. The existence of this range can be understood from the fact that the ratio of the thermal banana orbit width to the ion sound gyroradius ∆r b /ρ s ≈ 0.062 is small, so there exists a range of scales where k y ρ s is large enough that ions are largely adiabatic while k y ∆r b is small enough that finite banana orbit width effects have yet to become important. Thus, TEMs at ion-scales are not expected to have this weak dispersion.
Since the wave group velocity matches the phase velocity over a wide range in k, wave packets 21 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444 would travel with resonant particle trajectories in the absence of nonlinear scattering. As has been pointed out in the past 59,60 , since trapped electron bounce centers precess toroidally and do not decorrelate from the wave due to parallel streaming, they can remain at a fixed phase with the wave. This is contrasted with electrons and ions in passing orbits, which rapidly decorrelate with the wave phase due to parallel streaming. This long field-particle interaction time could lead to a large deflection of the trapped electron bounce center from its unperturbed trajectory. Such large deflections would lead to the breakdown of quasilinear theory, which is based on 'small' kicks due to the response of particles to the perturbing fields.
This breakdown is signaled quantitatively by the Kubo number K ≡ τ ac /τ S crossing unity.
Here τ ac is the autocorrelation time of the fields (in this case, electrostatic potential) as seen by resonant particles on unperturbed trajectories, and τ S is the 'bounce' time This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5144444 satisfies the requirements for quasilinear theory to hold. In contrast to this situation, for TEMs in the absence of collisions, the decorrelation will be due to the relative difference of the bounce center toroidal precession and the group velocity, and can be estimated as to the maximum ion scale linear growth rate, and show that ion-direction turbulence remains dominant over the entire plasma radius shortly before the rotation reversal. Next a plot of the maximum ion scale growth rate against the E × B shearing rate calculated from force balance is shown, where red and blue correspond to LOC and SOC respectively. Note the shearing rate reaches a significant fraction of the maximum ion scale growth rate. Finally, the linear spectrum for a single radius is shown.

33
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Note the exact ratios between families are underconstrained, so these are only illustrative solutions.

PLEASE CITE THIS ARTICLE AS
In SOC the subdominant families ITGa/c and TEMa are quenched in a way which respects particle flux constraint. (Bottom) An example ion heat flux spectrum which is consistent with the above mode families. Note the SOC spectrum in blue is narrower than the LOC spectrum in red.

35
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Case I (top) and II (bottom) both at r/a = 0.6. The ITG branch (magenta) and TEM/ETG branch (green) are shown in different colors. Note that v gr follows v ph more closely for the TEM/ETG branch than it does for the ITG branch. The text discusses this more quantitatively.

36
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.