3.1. On computer simulations

We focus on analytics and simulations in the linear stages mostly for validation of our computer model. The computer model is based on the explicit first order in time and second order in space approximation of equations on a rectangular grid. The Poisson brackets are evaluated using the second-order Arakawa (Reference Arakawa1966) method. The inhomogeneous Poisson equation, (2.5), is solved for potential $\varphi$ using preconditioned conjugate gradient algorithm 2 described in Fisicaro et al. (Reference Fisicaro, Genovese, Andreussi, Marzari and Goedecker2016).

The background continuum value provides numerical stability for temperature and density associated effects. Background density is only added to the definition of general vorticity, and background temperature to the denominator of the line-tying term; the rest of the equations remain unchanged with such substitution. More generally, there is a neutral gas in the limiter’s shade that interacts with plasma. Radial flows generated by these gas–plasma interactions are not considered. We keep the background low in our numerical implementation to balance stability and performance. The corresponding numerical value of the background does not exceed 10 % of the peak density and 1 % of the peak electron temperature.

The choice for $m=5$ initial stages of RT and fETG instabilities from figures 2 and 4 is for demonstration purpose only.

3.2. Rayleigh–Taylor instability of rotating plasma

Centrifugal drive in plasma is a consequence of radially inhomogeneous rotation. It can be thought as a RT instability in effective anti-gravitation generated in a rotating reference frame. Interpretation of a centrifugally forced RT instability in liquid can be found in Scase & Hill (Reference Scase and Hill2018). We include FLR effects for the isolated instability configuration in the following system:

(3.6) \begin{align} \dfrac {\text{d}}{\text{d}t}\boldsymbol{\nabla }\boldsymbol{\cdot }\left (n\boldsymbol{\nabla }\varphi \right )&=\left \{ n,\dfrac {\left (\boldsymbol{\nabla }\varphi \right )^{2}}{2}\right \} +U\boldsymbol{\nabla }\boldsymbol{\cdot }\left \{ \boldsymbol{\nabla }\varphi ,n\right \}\! ,\\[-10pt]\nonumber \end{align}

(3.7) \begin{align} \dfrac {\text{d}}{\text{d}t}n&=0. \end{align}

The Lagrangian derivative operator for rigid rotation

(3.8) \begin{equation} d_{t}=-i\omega +\varOmega \partial _{\theta }-imr^{-1}\varphi \partial _{r} \end{equation}

provides relation between perturbations of density and potential $n=\varphi n_{0}^{\prime }/(\varOmega r\mathcal{R})$ from (3.7), where

(3.9) \begin{equation} \mathcal{R}={1-\omega}/{(m\varOmega)} , \end{equation}

is the normalised rotational frequency factor. Linearising the substantial derivative of vorticity, centrifugal and FLR terms results in an ODE in the form (3.4) for the perturbation of the potential

(3.10) \begin{equation} \varphi _{rr}^{\prime \prime }+\varphi _{r}^{\prime }\left [\dfrac {1}{r}+\dfrac {1}{a_{n}}+\widetilde {U}g_{n1}\right ]/\left [1+\widetilde {U}\right ]+\varphi \left [F_{n}-\dfrac {m^{2}}{r^{2}}+\widetilde {U}g_{n2}\right ]/\left [1+\widetilde {U}\right ]=0, \end{equation}

where $a_{n}(r)$ is a density gradient scale and coefficients are

(3.11) \begin{gather} \widetilde {U}\left (r\right )=\dfrac {U}{a_{n}r\mathcal{R}\varOmega },\,\,\,F_{n}\left (r\right )=\dfrac {1}{a_{n}r\mathcal{R}}\left (\dfrac {1}{\mathcal{R}}-2\right )\!,\nonumber \\[9pt] g_{n1}\left (r\right )=\dfrac {1-a_{n}^{\prime }}{a_{n}}-\dfrac {1}{r\mathcal{R}},\,\,\,g_{n2}\left (r\right )=\dfrac {1}{r^{2}\mathcal{R}}-\dfrac {1-a_{n}^{\prime }}{a_{n}r\mathcal{R}}-\dfrac {m^{2}}{r^{2}}. \end{gather}

Following the locality (3.5), the fundamental frequency $\omega _{0}$ for the potential perturbation vanishes, resulting in equation

(3.12) \begin{equation} \left [F_{n}\left (r_{0}\right )-\dfrac {m^{2}}{r_{0}^{2}}+\widetilde {U}\left (r_{0}\right )g_{n2}\left (r_{0}\right )\right ]/\left [1+\widetilde {U}\left (r_{0}\right )\right ]=0. \end{equation}

The local spectrum of isolated RT instability is fully defined by the function $F_{n}$ , while $g_{n1,2}$ contribute only in the case of the FLR effect.

3.3. Hydrodynamic rigid rotation

In the absence of FLR effects, the instability reduces to a pure hydrodynamic regime. Early stages of instability onset are illustrated in figure 2. Equation (3.12) becomes a dispersion relation

(3.13) \begin{equation} \left (\mathcal{R}-\alpha \right )^{2}=\alpha ^{2}-\alpha , \end{equation}

where $\alpha =m^{-2}r_{0}/\delta$ .

Figure 2.

The onset of $m=5$ hydrodynamic RT instability.

One of the roots of (3.13) is unstable and has a maximum $\varGamma ^{\star }\equiv \mathcal{I}m(\mathcal{R}_{\alpha =\alpha ^{\star }})$ at

(3.14) \begin{equation} \alpha ^{\star }=\dfrac {1}{2},\,\,\varGamma ^{\star }=\varOmega \left (2\dfrac {\delta }{r_{0}}\right )^{-1/2}\equiv \dfrac {m\varOmega }{2}. \end{equation}

The $\varGamma ^{\star }$ curve passes through all the maxima for a fixed $m\in \mathbb{N}$ (figure 3, dashed black curve). There is no instability for $\delta /r_{0}\leq m^{-2}$ and the equality corresponds to the cutoff $\alpha _{cut}=1$ . As $m\rightarrow \infty$ , an initial perturbation becomes increasingly localised and instability develops at smaller gradient scales down to the limit of an infinitely sharp boundary $\delta \rightarrow 0$ . The global maximum growth rate $\varGamma _{max}$ is the maximum over all possible $m$

(3.15) \begin{equation} \alpha _{max}\rightarrow 0,\,\,\varGamma _{max}=\varOmega \left (\dfrac {\delta }{r_{0}}\right )^{-1/2}, \end{equation}

and always exceeds the maximum growth rate for a fixed azimuthal mode $\varGamma _{max}=\sqrt {2}\varGamma ^{\star }$ . The parameter $\varGamma _{max}$ is an enveloping curve above all $\varGamma _{m\lt \infty }$ (figure 3, solid black curve). The global $m=1$ mode is not excited even if the initial perturbation is “localised” in the width of entire plasma, since its cutoff is on the boundary of computational domain (3.3).

Figure 3.

The RT instability growth rates against $\delta /r_{0}$ . (a) – unstable roots of (3.13) with different $m$ in absence of FLR effect; (b) – unstable roots of the mode $m=7$ for various $U$ ; (c) – numerical increments of the system (3.6, 3.7).

3.4. Stabilisation by FLR

In case of FLR effects the dispersion (3.13) modifies into

(3.16) \begin{equation} \left (\mathcal{R}-\alpha _{\widehat {U}}\right )^{2}=\alpha _{\widehat {U}}^{2}-\alpha _{\widehat {U}}\dfrac {1+\widehat {U}}{1+ {m^{2}}/{2}\widehat {U}}, \end{equation}

where coefficients $\alpha _{\widehat {U}}=\alpha (1+m^{2}\widehat {U}/2)$ and $\widehat {U}=U/(\varOmega r_{0}^{2})$ . Its unstable solution has an additional resonant condition $\omega _{\widehat {U}}\neq m\varOmega (1-\widehat {U}r_{0}/\delta )$ following from the singularity in (3.12) and reaches its maximum at

(3.17) \begin{equation} \alpha _{\widehat {U}}^{\star }=\dfrac {1}{2}\dfrac {1+\widehat {U}}{1+ {m^{2}}/{2}\widehat {U}},\,\,\varGamma _{\widehat {U}}^{\star }=\dfrac {m\varOmega }{2}\dfrac {2\left (1+\widehat {U}\right )}{2+m^{2}\widehat {U}}. \end{equation}

The FLR motion affects modes with $m\gtrsim 2$ and increases with $m$ . Figure 3(b) illustrates the FLR changes to the $m=7$ RT mode, showing its complete stabilisation within the computational domain at $U\approx 0.283$ . The global maximum growth rate $\varGamma _{max}$ is absent as with $m\rightarrow \infty$ the stabilising effect of FLR turns out to be infinite. Asymptotic behaviour of $\varGamma _{max}$ follows from root expansion in (3.16) providing a triggering condition for the instability

(3.18) \begin{equation} U\lt 2\varOmega r_{0}^{2}m^{-1}\left (\dfrac {\delta }{r_{0}}\right )^{1/2}+\mathcal{O}\left (m^{-2}\right)\!. \end{equation}

Estimate (3.18) shows stabilisation of localised modes in rigid rotation for lower $m\sim 1$ with GDT associated parameters $U_{GDT}\approx 5$ and is consistent with numerical modelling. Similarities between pure RT and FLR-stabilised RT increments lead to the interpretation of FLR effect on RT instability given in Appendix A.

3.5. The flute electron-temperature-gradient instability

The fETG instability develops in open traps plasma with axial currents as a coupled effect of electron-temperature gradient and conducting end plates. The GDT axial currents are utilised for vortex transport barrier generation becoming the source for this instability. The fETG drive has been studied in different geometries and regimes. Berk et al. (Reference Berk, Ryutov and Tsidulko1990) considered the isolated fETG in a slab geometry. In later work, Berk, Ryutov & Tsidulko (Reference Berk, Ryutov and Tsidulko1991) derived a dispersion equation in flux coordinates including centrifugal drift and showed regimes with FLR stabilisation. We consider isolated fETG instability with axial loss. Energy dissipation is added to the Lagrangian derivative operator $d_{t}^{(\gamma )}\equiv d_{t}+\gamma H$ , where $\gamma$ is a factor to electron-temperature axial losses

(3.19) \begin{equation} \nu _{\parallel }^{\left (e\right )}=H\left (\dfrac {\rho _{i}}{r_{0}}\right )^{2}\dfrac {\tau _{\parallel i}}{\tau _{\parallel e}}=\gamma H. \end{equation}

This multiple $\gamma \approx 1$ for GDT parameters causes electrons to lose all their energy on a time scale of $\sim 0.3\,\text{ms}$ without external sources of energy. In the linear stages theory and simulations, we consider electron axial losses 10 times those of the warm plasma $\nu _{\parallel }^{(e)}=10\nu _{\parallel GDT}\sim 0.2$ . The fETG-isolated model with substantial derivative operator containing energy dissipation is

(3.20) \begin{align} \left (\dfrac {\text{d}}{\text{d}t}\right )^{\left (\gamma \right )}T_{e}\boldsymbol{\nabla} ^{2}\varphi &=H\left (\varphi -\varLambda T_{e}\right )\!,\\[-10pt]\nonumber \end{align}

(3.21) \begin{align} \left (\dfrac {\text{d}}{\text{d}t}\right )^{\left (\gamma \right )}T_{e}&=0, \end{align}

where we re-assigned the hybrid scalar $n$ to be electron temperature $T_e$ . The following steps are the same as for isolated RT mode. The ODE for the perturbation of potential is as follows:

(3.22) \begin{equation} \varphi _{rr}^{\prime \prime }+\varphi _{r}^{\prime }\dfrac {1}{r}+\varphi \left [F_{T_{e}}-\dfrac {m^{2}}{r^{2}}\right ]=0. \end{equation}

The spectrum is defined by the function

(3.23) \begin{equation} F_{T_{e}}\left (r\right )=-\dfrac {iH\varLambda m}{a_{T_{e}}r\left (m\varOmega \mathcal{R}_{\gamma }\right )^{2}}+\dfrac {iH}{T_{e}^{\left (0\right )}\left (m\varOmega \mathcal{R}_{\gamma }\right )}, \end{equation}

where $a_{T_{e}}(r)$ is the electron-temperature-gradient scale and $\mathcal{R}_\gamma \equiv \mathcal{R}-i\gamma H/(m\varOmega )$ is the frequency factor (3.9) with dissipation. Following the locality, the fundamental frequency $\omega _{0}$ for potential perturbation vanishes, resulting in the equation

(3.24) \begin{equation} F_{T_{e}}\left (r_{0}\right )-\dfrac {m^{2}}{r_{0}^{2}}=0. \end{equation}

3.6. The pure fETG

For isolated instability, the local dispersion (3.24) without dissipation $\gamma =0$ coincides with (23) obtained in Berk et al. (Reference Berk, Ryutov and Tsidulko1991) for a slab geometry

(3.25) \begin{equation} \mathcal{R}^{2}+i\mathcal{R}\dfrac {H}{T_{e}^{\left (0\right )}}\left (\dfrac {\delta }{k_{T_{e}}}\right )^{2}+i\dfrac {H\varLambda }{k_{T_{e}}}=0. \end{equation}

We re-assigned $\mathcal{R}=\omega -m\varOmega$ and introduced the azimuthal wave vector $k_{T_{e}}=m\delta /r_{0}$ and $\delta \equiv a_{T_{e}}(r_{0})$ . Figure 4 shows the early stages of fETG instability.

Figure 4.

The onset of $m=5$ fETG instability.

The maximum increment is independent of the azimuthal number and determined by temperature-gradient scale $\delta$ and plasma–wall coupling $H$ . The unstable solution of (3.25) has the following scalings

(3.26) \begin{equation} k_{max}=1.825H^{1/3}\delta ^{4/3}T_{e}^{-2/3}\varLambda ^{-1/3},\,\,\varGamma _{max}=0.384H^{1/3}\delta ^{-2/3}T_{e}^{1/3}\varLambda ^{2/3}. \end{equation}

Numerical solution of the (3.25) and simulations of equations (3.20, 3.21) with $\gamma =0$ are shown in figure 5(a). For some parametric sets $\left \{ H,\delta ,T_{e},\varLambda \right \}$ and azimuthal number $m$ , the value for $k_{max}$ exceeds the borderline $\delta /r_{0}=1$ of our computational domain (3.3). The corresponding increment does not reach its maximum value and continues to grow monotonically towards the borderline. For instance, for the set $\{ \varLambda ,H,T_{e},m\} =1$ the growth of the increment stops at the value ${\sim}0.3$ at the borderline, which is shown in pink in figure 5(a). This is the reason for growing inconsistency in modelling at higher $H$ for different fETG modes showed in figure 5(c). The azimuthal vector value increases with the plasma–wall coupling $k_{max}\sim H^{1/3}$ . The maximum increment of $m=3$ is on boundary $\delta /r_{0}=1$ at $H\sim 20$ , and for $m=5$ at $H\sim 100$ . From this point onward, the dependence of the maximum increment becomes weaker, shown for different $m$ with solid coloured curves in figure 5(c).

Figure 5.

(a,b) The fETG instability increments against $\delta /r_{0}$ , (c) – against wall coupling $H$ . (a) The unstable solution of (3.25) with $\{ \varLambda ,H,T_{e}\} =1$ ; (b) the unstable mode $m=11$ with $\{ \varLambda ,T_{e}\} =1$ and $H=50$ for different values of $\gamma$ ; (c) increments produced in numerical simulations of the system (3.20, 3.21) with $T_{e}=1,\varLambda =5$ .

In the region of low $k$ , the instability increment turns out to be independent of line tying $\varGamma _{k\ll 1}=k\delta ^{2}/(\varLambda T_{e})$ . The growing inconsistency for low $H$ is associated with transverse diffusion of $T_e$ added for numerical stability.

3.7. The fETG with axial loss

Each ion–electron pair escaping the plasma to the endwalls carries away energy, resulting in a decrease of initial perturbation growth. The dispersion relation (3.25) remains unchanged, but the rotation frequency $\mathcal{R}_\gamma$ acquires an imaginary additive resulting in decay of the maximum increment

(3.27) \begin{equation} \varGamma _{max}=0.384H^{1/3}\delta ^{-2/3}T_{e}^{1/3}\varLambda ^{2/3}-\gamma H. \end{equation}

The decrease in maximum increment of the mode $m=11$ with $T_{e}=1$ , $\varLambda =5$ is shown in figure 5(b). Plasma–wall coupling $H$ is no longer decisive for the instability, as any finite $\gamma$ limits the growth rate by the value

(3.28) \begin{equation} H_{crit}\sim 0.238\dfrac {\varLambda \sqrt {T_{e}}}{\delta }\gamma ^{-3/2}. \end{equation}

Damping in the spectrum causes cutoffs on $\delta /r_{0}$ for $\gamma =0.02$ , respectively $\delta /r_{0}\sim 0.3$ on the left and $\delta /r_{0}\sim 2$ on the right. These cutoffs do not affect $\varGamma _{max}$ until the initial perturbation is completely stabilised, or any of them pass outside the computational domain.

4.1. Inhomogeneous plasma inside vortex transport barrier

We start by formulating the governing equations for the radial transport simulations. First, we include current dissipation as the only loss term neglecting direct axial losses for vorticity and density. Second, we avoid radial transport associated with fETG, reducing the model to density distributed only

(4.1) \begin{align} \dfrac {\partial \varpi }{\partial t}+\left \{ \varphi ,\varpi \right \} &=\left \{ n,\dfrac {\left (\boldsymbol{\nabla }\varphi \right )^{2}}{2}\right \} +H\left (\varphi -\varphi _{w}\right )+\nu _{\perp }\boldsymbol{\nabla} ^{2}\varpi ,\\[-10pt]\nonumber \end{align}

(4.2) \begin{align} \dfrac {\partial n}{\partial t}+\left \{ \varphi ,n\right \} &=\nu _{\perp }^{n}\boldsymbol{\nabla} ^{2}n+Q_{n}. \end{align}

The density source $Q_{n}\equiv -\nu _{\perp }^{n}\boldsymbol{\nabla} ^{2}n_{st}$ is defined from the stationary density distribution $n_{st}$ using the profile (3.1).

We use the density radial transport coefficient $D_{n}$ to describe the quasi-stationary rotation. In our model, the flow (2.2) is divergence free and radial transport cannot be evaluated directly from electric potential. The $D_{n}$ can be calculated using the particle flux density $\boldsymbol{\mathcal{F}}_{n}$ in the continuity equation with source

(4.3) \begin{equation} \dfrac {\partial n}{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\mathcal{F}}_{n}=Q_{n} \end{equation}

and we are interested in its radial component $\mathcal{F}_{n}$ . The flux density is determined by integrating (4.3) over the area enclosed in a circle of radius $r$

(4.4) \begin{equation} \underset {S\lt \pi r^{2}}{\int }\text{d}S\boldsymbol{\nabla }\boldsymbol{\mathcal{F}}_{n}=2\pi r\mathcal{F}_{n}\left (r\right ), \end{equation}

where $\mathcal{F}_{n}$ is radial particle flux through the boundary of the circle. The source and time derivative terms are expressed from the density convection (4.2) and integrated over the same circle. Expressing the particle flux gives two contributions

(4.5) \begin{equation} \mathcal{F}_{n}=\dfrac {1}{2\pi r}\int \text{d}S\left \{ \varphi ,n\right \} -\dfrac {\nu _{\perp }^{n}}{2\pi r}\int \text{d}S\boldsymbol{\nabla} ^{2}n\equiv \mathcal{F}_{V}+\mathcal{F}_{D}, \end{equation}

the vortex $\mathcal{F}_{V}$ and diffusive $\mathcal{F}_{D}$ . On the other hand, the particle flux is proportional to classical transport coefficient $D_{n}$ and density gradient

(4.6) \begin{equation} {\mathcal{F}_{n}}=-D_{n}{\nabla} n. \end{equation}

To interpret the numerical experiments, we perform time and angle averaging in the transport coefficients. Convection forms density filaments with local extrema, where the transport is undefined. For the same reason, radial components of the particle flux and the density gradient are averaged separately

(4.7) \begin{equation} D_{n}=-\dfrac {\left \langle \mathcal{F}_{n}\right \rangle _{t,\theta }}{\left \langle \nabla n\right \rangle _{t,\theta }}\equiv D_{V}+D_{D}, \end{equation}

where, similarly, $D_V$ and $D_D$ are the vortex and classical transport coefficients. Radial transport is not described by effective diffusion coefficient in general. Naulin (Reference Naulin2007) points out that the flux-gradient relation (4.6) is valid for passive scalar convection only. The turbulence is fed from the energy stored in gradients of distributed quantities. In our case there is an active energy transfer between RT instability and density convection. A turbulent flow reaching quasi-equilibrium contains coherent fluctuations that are constantly born and dissipated within the flow width. This provides an average eddy spectrum, so the radial transport coefficient can be evaluated as a time and angle average over quasi-stationary stage.

We normalise the transport coefficient $D_{n}$ to the density diffusion coefficient $\nu _{\perp }^{n}$ to distinguish the nature of the flow. For pure collisional transport with $D_{n}/\nu _{\perp }^{n}\approx 1$ and $D_{V}\approx 0$ a flow reaches a stationary state with axially symmetric density and vorticity distributions. A quasi-stationary state corresponds to a transport different from pure collisional that remains unchanged over a long time period. The experimental base of average radial transport coefficients (4.7) is collected over all possible combinations of model parameters presented in table 1. The projection of the biased limiter in these simulations coincides with plasma radius $r_{0}$ and the initial flow width is set $\delta =0.2$ . The total number of unique quasi-stationary scenarios is 600, including 60 simulations with uniform density.

Table 1.

Model parameters in numerical experiments.

Figure 6.

(a) Radial transport coefficient $D_{n}/\nu _{\perp }^{n}$ in different scenarios. (b) Linear stage of $m=7$ KH instability with positive inner (red dashed) and negative outer (blue dashed) vortex rings symmetrical about maximum speed location (black dashed).

Three different types of $D_{n}$ observed in numerical simulations are shown in figure 6(a). The blue line is stationary plasma rotation with an axially symmetric distribution with purely collisional transport. The other two are different regimes of anomalous transport that we refer to as single hump (green) and two hump (red).

From the linear stage perspective the two-hump regime has to be more relevant. The KH instability of a sheared flow generates a Kármán (Reference von Kármán1911) vortex street with two lines of alternating vortices corresponding to each inflection point on the potential shown in figure 6(b). Vortex shedding preserves the number of vortices and net vorticity making each vortex ring contain the same number of individual vortices, as in the absence of sheared flow their total amount is zero. The inner ring is smaller in radius, its vortices are packed denser and we can expect the associated transport to be greater. Such behaviour is common in simulations with a two-hump regime, however, around 10 % of instances have reversed outer-ring-dominated radial transport. Aside from that, there is a single-hump quasi-equilibrium that is not consistent with the sheared rotation.

The observed quasi-stationary regimes are determined by their spatial parameters – humps heights, widths and positions, etc., in total 8 unique values per sample. The model variables are a part of parameter space and can be correlated with transport establishment. Another measurable value is angle average of the quasi-stationary flow width $\xi _{d}$ related to normalised potential $\xi _{d}=1/( 2|\boldsymbol{\nabla }\varphi |_{max})$ . The value of $\xi _{d}$ determines the flow broadening as the flow tends to its equilibrium state $\varphi _{w}$ . The angle and time mean value of $\xi _{d}$ is calculated for both distributed and uniform density models.

The parameter space of two-hump simulations has at least 15 explicit entries, with three linear classifiers to break down the dependencies and classify parameters determining regimes of anomalous transport. The flow regime is set to be the target variable and the parameter space is left as features. The Pearson correlation matrix on two-hump simulations showed the following dependencies:

1. (i) the flow width $\xi _{d}$ decreases with line-tying $H$ ;

2. (ii) the inner-ring width acts like an independent variable with $\left \langle R^{2}\right \rangle \approx 0.09$ , while the outer-ring width is defined by $H$ and correlated with flow width $\xi _{d}$ ;

3. (iii) the fall between the rings is positioned by collisional viscosity $\nu _{\perp }$ with $R^{2}=0.85$ ;

4. (iv) radial locations of inner- and outer-ring peaks shift towards limiter projection as H increases. Correlation between $\xi _{d}$ and the peaks’ interval is $R^{2}\approx 0.78$ .

Most of them result from the stationary potential profile characterised by its flow width $\xi _{d}=(4\nu _{\perp }/H)^{1/4}$ . Among less clear dependencies is a possible explanation for outer-ring-dominated transport simulations.

1. (i) Density diffusion $\nu _{\perp }^{n}$ provides numerical stability for inner-ring-dominated transport. However, this correlation breaks in all outer-ring-dominated instances. Density diffusion coefficient in these simulations determines the density gradient scale in the same manner as collisional viscosity determines flow width $a_{n}=(4\nu _{\perp }^{n}/H)^{1/4}$ . Visuals showed active RT instability in these numerical experiments. This can be a consequence of reaching a certain local density gradient scale sufficient for RT generation outside the limiter. The flow width is established by line tying, which is ${\sim}30$ times higher beyond the limiter. The KH instability of this region generates eddies ${\sim}2.3$ times smaller than inside that transfer energy from potential to density gradients much faster producing RT instability associated transport outside the limiter.

The single-hump regimes have lesser parameter space and the only dependence is found – the flow width decreases both with line tying and collisional viscosity. Spectral analysis showed the absence of mode structure beyond the limiter and global $m=1$ inside. This regime is established when outer plasma lies in the stable region of parameter space.

Cross-validation with Lasso and backward elimination methods confirmed the Pearson results in all simulations, but one stood out. Initial state defined by $r_{0}$ being almost an independent variable with mean $\langle R^{2}\rangle \approx 0.12$ showed an impact on geometric properties of the flow inside the limiter. In our case this is natural, as we used $r_{0}$ -characterised source (3.1) for a stationary density profile. Using an exponential model source with gradient scales dependent on $r_{0}$ followed the tendency, however, numerical experiments with a constant radial flux source free from $r_{0}$ -dependency invalidated it.

The highest correlation with $R^{2}=0.99$ is observed between the flow widths of distributed $\xi _{d}$ and homogeneous $\xi _{d0}$ simulations. The linear regression for changes in flow width is illustrated in figure 7. All quasi-stationary points represented by circles belong to the trend with slope 0.89, meaning that the width of the distributed flow shrinks in comparison with the homogeneous flow. Each point has the same model parameters and only differs by featuring inhomogeneous density.

Figure 7.

Flow width in distributed and homogeneous simulations. Circles represent instances that reached quasi-equilibrium, squares – instances with unsettled anomalous transport, colour denotes line tying. A group of crosses below the main trend is a model with exponentially distributed source.

Initial states determine the equilibrium via the source model. Using an exponential source model verified the linear relation with a different slope as shown in figure 7 with yellow circles. The constant radial flux source model free of $r_{0}$ dependency also followed the linear trend. The linear shrink turns out to be universal independently from the source, in general

(4.8) \begin{equation} \xi _{d}=A\left (Q_{n}\right )\xi _{d0}, \end{equation}

where $A\lt 1$ depends on the equilibrium model.

4.2. The gap effect in vortex confinement

In these simulations we juxtapose two models – original (1.2)–(1.3) and hybrid (2.7)–(2.6). Transversal diffusion coefficients are set $\nu _{\perp }^{n}=\nu _{\perp }=0.002$ and axial losses $\nu _{\parallel }=A_{E}\nu _{\parallel }^{n}=0.02$ with $A_{E}=8$ . The comparison is built around two different biasing schemes in GDT experiments – 1) with a positive edge potential applied to both the last endwall plate and limiter, and 2) to the limiter only. In the case of positive biasing of the endwall and limiter relative to the vacuum chamber, there is a spatial gap between their projections in the plasma. Experimentally, instability is not observed in the first case with a gap, but in the gapless case with limiter-only biasing additional transport occurs. This spatial difference in biasing schemes results in different flow regimes outside and inside the biased potential radius that might interact producing scrape-off layer (SOL) instability. Visual representation of the gap is shown in figure 8.

Figure 8.

Sketch for visual differences in biasing schemes. Here, $\varphi _w$ is the biased potential (positive edge in simulations), $H$ line tying that is proportional to axial losses.

The line tying contains the net ion current to the walls that is much higher beyond the limiter. Current paths of this region are impossible to trace, so the change in $H$ can be modelled as an effective increase outside the limiter by a mean constant factor of 10–20, which is fairly consistent with experimental data.

A measure for confinement quality in hybrid model is particle confinement time

(4.9) \begin{equation} \tau _{n}\left (t\right )=\dfrac {\left \langle n\right \rangle }{\left \langle Q_{n}\right \rangle -\left \langle \dot {n}\right \rangle }, \end{equation}

where averaging is over the cross-section. The initial distributions of hybrid scalar and potential are set to satisfy purely axial losses. Comparison of different quasi-stationary regimes is given in figure 9. In the gap-based simulations the vortex confinement appears to be stable even when the gradients of plasma density and electron temperature are included. The additional polarising effects associated with the distributed ion density and electron temperature are insignificant, making the relative decrease in particle confinement time around 3 %.

Figure 9.

Particle confinement time in units of $\tau _{_{GDT}}\approx 30\mu$ s for distributed and homogeneous models in two biasing schemes. Blue lines are homogeneous model, red lines are distributed model. Dotted lines are simulations with biased limiter, plane lines are with biased endwall.

The gap-based biasing is more challenging to use, as end plates tend to discharge. Positive potential edge is then created only at the limiter, making the biasing gapless while eradicating the discharge in the experiment. As there is no longer a spatial gap between limiter and endwall projections, different equilibrium is reached on the outside. For the original model, this results in coupling of the $m=1$ mode inside the limiter and $m=2$ outside oscillating at $\omega \sim 65$ kHz.

In case of hybrid model the interaction between inside and outside is more complex. Visuals show that plasma becomes unstable after reaching the axially symmetric phase. These are sawtooth-like oscillations with $\omega \sim 23$ kHz having breakdowns at $t\approx 0.01,0.012,0.014$ s. Large-scale waves with comparable frequency have been observed in GDT experiments by Prikhodko et al. (Reference Prikhodko, Bagryansky, Beklemishev, Kolesnikov, Kotelnikov, Maximov, Pushkareva, Soldatkina, Tsidulko and Zaytsev2011). It is important to note that in the experiment, these oscillations are observed at a relatively high mirror ratio.

Reaching critical gradient scale during the axially symmetric phase for either of fETG or RT instability might be the cause of vortex confinement degradation, as we have seen in simulations of quasi-stationary transport. One could distinguish different flute modes in linear stages that are impossible in developed edge turbulence. However, RT mode should be effectively stabilised by FLR effects, while fETG has a faster growth rate outside the limiter as $H$ is effectively higher. Visuals show moving waves along the plasma boundary carrying small filaments slower than induced rotation, implying a frequency delay, which is a strong indicator of fETG onset. There can also be a contribution from KH instability (Beklemishev Reference Beklemishev2024), as the effective increase in $H$ lowers the inner layer width determining the KH-generated vortex size in the SOL. One could speculate about the interplay between KH and fETG in the SOL that Lotov, Ryutov & Weiland (Reference Lotov, Ryutov and Weiland1994) show to be regulated by the parameter

(4.10) \begin{equation} R= \frac{l}{\rho _i}\left (\frac{\varepsilon l}{\varLambda L}\right )^{1/3}, \end{equation}

where $\varepsilon$ is the fraction of impinging ions absorbed by the conducting wall and $l$ is the electron-temperature-gradient scale. Evaluating (4.10) for GDT parameters with $\varepsilon =0.01$ we get $R\sim 10^{-5}$ , meaning that the edge turbulence might be dominated by KH instability, rather than by fETG.

The influence of biasing scheme can be studied further by introducing the spatial gap $\varDelta =r_{lim}-r_{w}$ as a parameter and inspecting how the flow degrades with gap varying into the gapless. The limiter radius $r_{lim}$ is fixed, while the biased potential radius $r_{w}$ varies. The mean particle confinement times retrieved from simulations in four different source models depending on the relative gap size are shown in figure 10. Arrows point to the corresponding values of energy confinement time evaluated from the time evolution of the plasma density measured by the Thomson scattering system in Bagryansky et al. (Reference Bagryansky, Lizunov, Zuev, Kolesnikov and Solomachin2003). Experiments before the year 2011 with controllable potential applied to the limiters and endwalls showed a plasma decay time comparable to pure axial gas-dynamic loss. The biasing scheme changed to limiter only after around the year 2011, which led to decrease in confinement time by ${\sim}10$ % that is consistent with our modelling.

Figure 10.

Mean particle confinement times against gap size. Colour denotes different equilibrium source models.

The gap dependency turns out to have an optimal plateau around two decreasing branches. The low- $\varDelta$ is the earlier observed unstable regime. The high- $\varDelta$ unstable branch correspond to a thin plasma limit, where density gradients inside the biased plates become sufficient to produce additional convection. The plateau is laying around $0.2\lesssim \varDelta /r_{w}\lesssim 0.4$ that can be recalculated into the safety zone, where additional instability is not generated

(4.11) \begin{equation} \dfrac {5}{7}\lesssim \dfrac {r_{w}}{r_{lim}}\lesssim \dfrac {5}{6}, \end{equation}

and the optimal relation is proposed as $r_{w}=0.8r_{lim}$ . The optimal constant of 0.8 is only valid for GDT. It is specific to any open trap that utilises externally controlled potentials and has to be recalculated in a suitable model. More complex biasing schemes can include a gradual increase in potential using sectioned endwalls, e.g. as is done in Ivanov et al. (Reference Ivanov2023). In that case only the last biased section is enough to satisfy safety conditions as (4.11), regardless of the number of biased steps.