Finite T_i Effects on Parallel-Velocity-Gradient Driven Instability

Abstract

This study investigates the influence of finite ion temperature on the linear growth rate of the parallel velocity gradient (PVG)-driven instability. Parallel (toroidal) flows are prevalent in magnetically confined fusion plasmas, where external momentum sources such as neutral beam injection (NBI) can act as a potential driver for PVG instabilities in both conventional and spherical tokamaks. Accurate transport modeling using quasi-linear models such as TGLF or QuaLiKiz requires a fundamental understanding of the characteristics of PVG modes in warm plasmas. The linear stability analysis is conducted to include the finite ion temperature effects. Depending on the radial profiles of ion temperature and parallel velocity, the dominant unstable modes can be categorized as either ion temperature gradient (ITG) or PVG-driven modes. These two instabilities are found to be mutually exclusive in their parameter spaces. Both modes become excited when the compression becomes negative. The ITG and PVG modes are each strengthened by the gradient that drives the other, up to the point where the driving gradient exceeds the instability threshold of the respective mode. In the PVG-dominant regime, increasing the temperature ratio enhances the compression and thereby exerts the stabilizing effects, whereas the temperature gradient contributes to PVG growth until the ITG thresholds is reached.

1.  Introduction

Confinement scaling plays a crucial role in the design of fusion devices. It is now widely recognized that the size of conventional tokamaks can be reduced owing to their advantageous scaling with the toroidal magnetic field [1]. This reduction in size directly implies a potential decrease in construction and operational costs. By contrast, confinement scaling in spherical tokamaks exhibits a stronger dependence on the magnetic field [2], suggesting that spherical configurations may offer a more cost-effective path to fusion energy. Because of fundamental differences in confinement scaling, transport in spherical tokamaks may be governed by turbulence mechanisms beyond the conventional ion temperature gradient (ITG) and trapped electron modes. Indeed, recent gyrokinetic (GK) studies have suggested that kinetic ballooning modes, microtearing modes, and Alfvén eigenmodes may also play a significant role in regulating transport in these systems.

This study focuses on the parallel velocity gradient (PVG)-driven mode. Previous simulations have reported that PVG-driven turbulence can be triggered by neutral beam injection (NBI) in spherical tokamaks [3, 4]. In addition, theoretical analyses have indicated that equilibrium flows in spherical tokamaks can become sufficiently strong to develop significant gradients in parallel velocity, thereby driving PVG instabilities [5]. This mode has been shown to contribute to momentum relaxation, enhanced particle transport, and the formation of secondary flow structures both theoretically [6–9] and experimentally [10–12]. Even when parallel flow shear does not directly drive the instability, it has been found to reduce the nonlinearly computed ion heat fluxes associated with ITG turbulence [13]. Linear gyrokinetic studies have also shown that PVG-driven turbulence can couple not only with drift waves but also with ITG modes, leading to both subcritical excitation and suppression of turbulence [14]. This subcritical behavior can be characterized by the zero-turbulence manifold, where subcritical turbulence cannot be nonlinearly sustained [15]. Finally it is desirable to obtain the transport characteristics of PVG mode using transport simulation in the empirically constrained domain, e.g. using QuaLiKiz [16] and TGLF [17, 18]. However, the basic properties of the PVG mode in the presence of finite ion temperature remain insufficiently understood. The present study is therefore dedicated to investigating the effects of finite $T_i$ on the linear growth of PVG-driven instabilities.

2.  Model

We use a three-field $(\varphi ,p,V_{\parallel })$ fluid model. A linearized system is employed in this study. In this system, we assume an adiabatic electron response, i.e., $\delta n/n=|e|\delta \varphi /T_e$, which reduces the number of fields. The system includes Finite Larmor Radius (FLR) effects. The set of linearized equations for the perturbations, $\delta n$, $\delta \varphi$ and $\delta V_{\parallel }$, is shown below:

  

$$\frac{\mathit{\partial }}{\mathit{\partial }t}[(1-{\rho }_s^2{\nabla }_{⟂}^2)\widehat{\varphi }-\tau {\rho }_s^2{\nabla }_{⟂}^2\widehat{p}]+V_{\star }\frac{\mathit{\partial }\widehat{\varphi }}{\mathit{\partial }y}+{\nabla }_{\parallel }{V}_{\parallel }=0,$$(1)

  

$$\begin{array}{c}\frac{\mathit{\partial }}{\mathit{\partial }t}[\left(1-\frac{10}{3}{\rho }_s^2{\nabla }_{⟂}^2\right)\widehat{p}-\frac{5}{3}{\rho }_s^2{\nabla }_{⟂}^2(1-\tau )\widehat{\varphi }]\hfill \\ +V_{\star p}\frac{\mathit{\partial }\widehat{\varphi }}{\mathit{\partial }y}+\frac{5}{3}{\nabla }_{\parallel }{V}_{\parallel }=0\hfill \end{array},$$(2)

  

$$\frac{\mathit{\partial }{\widehat{V}}_{\parallel }}{\mathit{\partial }t}+V_{\star v}\frac{\mathit{\partial }\widehat{\varphi }}{\mathit{\partial }y}+c_s{\nabla }_{\parallel }(\widehat{\varphi }+\tau \widehat{p})=0,$$(3)

where ${\rho }_s$ is the ion Larmor radius and $c_s$ is the sound velocity. $V_{\star }\equiv c_s{\rho }_s/L_n$, $V_{\star p}\equiv (1+{\eta }_i)c_s{\rho }_s/L_n$ and $V_{\star v}\equiv c_s{\rho }_s/L_v$. $L_n$ and $L_v$ are the scale-lengths for the equilibrium density and velocity. ${\eta }_i\equiv L_n/L_T$ and $L_T$ is the one for equilibrium temperature. $\widehat{\mathrm{\dots }}$ means normalized fluctuating quantities, e.g. $\widehat{n}\equiv \delta n/n_0$, $\widehat{p}\equiv \delta p/p_0$ and ${\widehat{V}}_{\parallel }\equiv \delta V_{\parallel }/V_0$, with $V_0$ taken to be $c_s$. The parameter $\tau \equiv T_{i0}/T_{e0}$ is thermal ratio between the main ion and electron; thus, finite $T_i$ effects enter mainly through ${\eta }_i$ and $\tau$. Without coupling to the pressure equation, the system reduces to a modified Hasegawa-Mima equation, which has been studied extensively in [8]. Hereafter, the correspondence $(r,\theta ,\parallel )\leftrightarrow (x,y,z)$ is assumed.

3.  Local Dispersion Relation

By applying a Fourier transform, the time and spatial derivatives are replaced by, ${\mathit{\partial }}_t\to -i\omega$ and ${\mathit{\partial }}_{\mathit{x}}\to i\mathit{k}$, respectively. Normalized density fluctuation response with adiabatic electron is obtained as:

  

$$\begin{array}{cc}{\widehat{\varphi }}_k=& -\frac{\tau {\rho }_s^2{k}_{⟂}^2}{1+{\rho }_s^2{k}_{⟂}^2-{\omega }_{\ast e}/\omega }{\widehat{p}}_k\hfill \\ & +\frac{k_z{c}_s/\omega }{1+{\rho }_s^2{k}_{⟂}^2-{\omega }_{\ast e}/\omega }{\widehat{V}}_{z,k}\hfill \end{array}.$$(4)

Normalized parallel velocity and pressure perturbation responses are as follows:

  

$${\widehat{V}}_{z,k}=\frac{k_z{c}_s}{\omega }\left(1-\frac{k_y{\rho }_s}{k_z{c}_s}{⟨v_z⟩}^{\prime }\right){\widehat{\varphi }}_k+\frac{k_z{c}_s}{\omega }\tau {\widehat{p}}_k,$$(5)

  

$$\begin{array}{cc}{\widehat{p}}_k=& \frac{1}{1+(10/3)\tau {\rho }_s^2{k}_{⟂}^2}\frac{5}{3}\frac{k_z{c}_s}{\omega }{\widehat{V}}_{z,k}\hfill \\ & -\frac{1}{1+(10/3)\tau {\rho }_s^2{k}_{⟂}^2}\left[\frac{5}{3}(1-\tau ){\rho }_s^2{k}_{⟂}^2-\frac{{\omega }_{\ast p}}{\omega }\right]{\widehat{\varphi }}_k\hfill \end{array}.$$(6)

The local dispersion relation of this system originates from Eqs. (4)–(6) and it becomes cubic:

  

$$a{\mathrm{\Omega }}^3+b{\mathrm{\Omega }}^2+c\mathrm{\Omega }+d=0$$(7)

  

$$a=1+\left(1+\frac{10}{3}\tau \right){\rho }_s^2{k}_{⟂}^2+\frac{5}{3}\tau \left(1+\tau \right){\rho }_s^4{k}_{⟂}^4$$

  

$$b=-\left[1-\left({\eta }_i-\frac{7}{3}\right)\tau {\rho }_s^2{k}_{⟂}^2\right]{\rho }_s{k}_y\cdot {\rho }_{\star }$$

  

$$\begin{array}{cc}\hfill c=& -k_z^2{\rho }_s^2[\frac{5}{3}\tau (1+\tau {\rho }_s^2{k}_{⟂}^2)\phantom{\frac{k_y}{k_z}}\hfill \\ & +\left(1+\frac{5}{3}\tau {\rho }_s^2{k}_{⟂}^2\right)\left(1-\frac{k_y}{k_z}\frac{{⟨v_z⟩}^{\prime }}{{\omega }_{ci}}\right)]\hfill \end{array}$$

  

$$d=-[{\eta }_i-\frac{2}{3}]\tau k_y{\rho }_s\cdot k_z^2{\rho }_s^2\cdot {\rho }_{\star },$$

where $\mathrm{\Omega }$ is $\omega$ normalized by ${\omega }_{ci}$. ${\rho }_{\star }\equiv {\rho }_s/L_n$ and ${⟨v_z⟩}^{\prime }\equiv {\mathit{\partial }}_r⟨V_z⟩$ are defined. $⟨\dots ⟩$ means averaged quantity. Dispersion relation for both modes is obtained in Fig. 1. It is worth noting that the ITG and PVG-driven modes are mutually exclusive in terms of destabilization. Especially, ITG mode becomes stabilized with sufficiently strong flow shear. Due to this property, analysis can be separated into two parts: ITG and PVG dominant regions.

Fig. 1.  Contours are growth rates of ITG and PVG modes with $\tau$ = 1.0 and $(k_x,k_y,k_z)=(0.8/{\rho }_s,0.8/{\rho }_s,0.1/L_n)$. Left panel: ITG mode’s growth rate. Right panel: PVG mode’s growth rate.

3.1  ITG mode with flow shear

The slab-ITG (sITG) and PVG modes appear to be decoupled, as shown in Fig. 1. This decoupling allows the subsequent analysis in this section to focus on the effects of parallel flow shear on sITG growth. Notably, sITG growth can be enhanced by parallel flow shear, for example at ${\eta }_i$ ≈ 1.6 in the left panel of Fig. 1. This behavior can be understood by noting that both sITG and PVG modes share a common destabilization mechanism: negative compressibility [8]. In the sITG mode, the negative compression arises from the phase difference between the temperature and density perturbation described by the Boltzmann response: $\mathit{\delta }n/n_0\sim |e|\mathit{\delta }\varphi /T_e$. The effect of parallel sheared flow is introduced through the final term on the left-hand side of Eq. (2), which enhances the growth of $\mathit{\delta }n$ through compression. Before the PVG threshold is reached, the velocity perturbation $\mathit{\delta }{V}_{\parallel }$ helps the temperature perturbation $\mathit{\delta }{T}_i$ in growing by reducing the stabilizing contribution of sound waves. This destabilizing effect is in qualitative agreements with the gyrokinetic results for the Waltz standard case [13].

Within the framework of the present model, this effect is described in Eq. (6). Here, we retain only the terms relevant to the sheared flow ${⟨v_z⟩}^{\prime }$ in the pressure response, so as to highlights its role:

  

$$\begin{array}{c}\left(1+\frac{10}{3}\tau {\rho }_s^2{k}_{⟂}^2-\frac{5}{3}\tau \frac{k_z^2{c}_s^2}{{\omega }^2}\right){\widehat{p}}_k=\hfill \\ -\left[\frac{5}{3}(1-\tau ){\rho }_s^2{k}_{⟂}^2-\frac{{\omega }_{\ast p}}{\omega }-\frac{5}{3}\frac{k_z^2{c}_s^2}{{\omega }^2}\left(1-\frac{k_y{\rho }_s}{k_z{c}_s}{⟨v_z⟩}^{\prime }\right)\right]{\widehat{\varphi }}_k\hfill \end{array}.$$(8)

The pressure response to the electrostatic potential above shows that sheared flow reduces the compression contribution of sound waves arising from the parallel electric field.

However, the sITG mode becomes stabilized when the parallel flow shear becomes sufficiently strong. This can be understood from Fig. 1 and the remarks above: in the strong-shear regime. This is because the free energy source of the instabilities shifts from the temperature gradient to the flow shear itself. This stabilization effect of sheared flow is consistent with previous GK analysis [14].

Notably, this sITG-PVG regimes need the following assumption, as pointed out in [8, 9]:

  

$$k_y{k}_z{⟨v_z⟩}^{\prime }>0.$$(9)

If Eq. (9) is not satisfied, the parallel flow shear suppresses rather than enhances sITG growth. In this case, the parallel velocity gradient does not act as a free energy source for driving instability. Instead, the parallel flow shear strengthens the competition ion acoustic waves by enhancing the parallel compression, thereby suppressing the sITG mode when Eq. (9) is violated.

3.2  PVG mode with finite T_i

In Fig. 1, the PVG mode appears once the parallel flow shear becomes strong enough to stabilize the ITG mode. As expected, the growth of this mode depends directly on the magnitude of the parallel flow shear. In this Sec. 3.2, the effects of finite ion temperature on the PVG mode are explored from a theoretical perspective.

Both panels of Fig. 2 show the effects of the temperature ratio $\tau$ mainly on PVG growth. First, the temperature gradient ${\eta }_i$ is fixed as 1.6 to isolate the effect of the temperature ratio. The temperature ratio corresponding to Fig. 1 is shown as the orange curves in both panels. Roughly speaking, the temperature ratio plays a role in easing PVG mode growth. Indeed, the left panel illustrates that the PVG mode threshold increases as the temperature ratio increases, while the right panel shows that the PVG growth rate decreases with increasing temperature ratio.

Fig. 2.  PVG growth rate versus parallel velocity gradient (left panel) and ${\rho }_s{k}_y$ (right panel) to show the parameter survey of temperature ratio $\tau \equiv T_{i0}/T_{e0}$. In the left panel, $(k_x,k_y,k_z)=(0.8/{\rho }_s,0.8/{\rho }_s,0.1/L_n)$ are used. In the right panel, $k_x$ and $k_z$ are the same.

Figure 3 shows the effects of the temperature gradient ${\eta }_i$ on PVG growth. Here, the temperature ratio is fixed at unity to isolate the role of the temperature gradient. As shown in Fig. 3, the temperature gradient plays a role in destabilizing the PVG mode. Increasing the temperature gradient enhances PVG growth rate, with the effect being particularly pronounced at ${\rho }_s{k}_y\sim \mathcal{O}(1)$. Interestingly the effects of the temperature gradient differs from that of the temperature ratio. A possible explanation for this distinction is provided below.

Fig. 3.  ${\eta }_i$ scan for growth rate of PVG mode. Thermal ratio $\tau$ is equal to unity. $k_x$ and $k_z$ are the same as other figures. Once the pressure coupled, temperature gradient ${\eta }_i$ helps to enhance PVG growth as the same mechanism discussed in Sec. 3.1.

To identify the mechanism that differentiates the effects of $\tau$ from ${\eta }_i$ effects, we start from the basic destabilization mechanism of the PVG mode. To clarify this mechanism, we temporarily take the limit $\tau \to 0$. Under this assumption, the velocity response reduces to the form reported in previous literatures [7–9].

  

$${\widehat{V}}_{z,k}\approx \frac{k_z{c}_s}{\omega }\left(1-\frac{k_y}{k_z}\frac{{⟨v_z⟩}^{\prime }}{{\omega }_{ci}}\right){\widehat{\varphi }}_k.$$(10)

In the parallel velocity response, the sheared flow helps reduce the stabilizing effect of the sound wave driven by the parallel electric field. Without a sheared flow, the sound wave carries fluctuations away in the direction parallel to the magnetic field. Naturally, the dispersion relation in Eq. (10) can be reduced to that of the sound wave when no sheared flow is present.

When a finite ion temperature is introduced through the thermal ratio, the parallel pressure term appears as the last term on the left-hand side of Eq. (3). This parallel pressure term drives the ion acoustic wave, which also plays a similar role to the sound wave driven by the parallel electric field. Consequently, it enhances the overall stabilizing effects. The modified velocity response is sketched as below:

  

$$\begin{array}{c}\left(1+\frac{10}{3}\tau {\rho }_s^2{k}_{⟂}^2-\frac{5}{3}\tau \frac{k_z^2{c}_s^2}{{\omega }^2}\right){\widehat{V}}_{z,k}=\hfill \\ \frac{k_z{c}_s}{\omega }\left[1-\frac{k_y}{k_z}\frac{{⟨v_z⟩}^{\prime }}{{\omega }_{ci}}-\tau \left(\frac{5}{3}(1-\tau ){\rho }_s^2{k}_{⟂}^2-\frac{{\omega }_{\ast p}}{\omega }\right)\right]{\widehat{\varphi }}_k\hfill \end{array}.$$(11)

This velocity response shows that the ion acoustic wave appears in the denominator of the coefficient of the electrostatic potential, which suppresses velocity fluctuation in addition to the FLR effects. From this point of view, the dependence on $\tau$ shown in Fig. 2 becomes clear. When $\tau$ ≪ 1, both stabilizing terms become weak, and when $\tau$ increases, the stabilizing effects become strong. This mechanism is straightforward: increasing $\tau$ raises the sound speed, which in turn enhances stabilization of the fluctuations. Then, the PVG threshold increases and the growth rate decreases.

Temperature effects also arise from temperature gradient. The relevant terms appears on the right-hand side of Eq. (11). Among the driving terms, the temperature gradient contributes to distabilization by weakening the sound wave. As shown in Fig. 3, increasing the temperature gradient enhances the PVG growth rate. For example with $\tau$ = 1, Eq. (11) simply expresses that the temperature gradient has the same role as the velocity gradient, i.e., negative compression. Note that the influence of the temperature gradient is somewhat simpler than that of the thermal ratio because the former is regulated by $\tau$, and the dominant unstable mode transitions from the PVG mode to the ITG mode once the temperature gradient becomes sufficiently steep, as shown in Fig. 1.

The characteristic scale of the PVG mode, ${\rho }_s{k}_y\sim \mathcal{O}(1)$, remains consistent though finite temperature effects affect the PVG growth. This implies that PVG mode can still contribute to transport phenomena in heated plasmas. Ongoing transport analysis is currently underway, and will be reported elsewhere.

4.  Summary and Discussion

We investigated the sITG-PVG system as a first step toward understanding confinement scaling in spherical tokamaks. Parallel flow shear is generated both by applying NBI and inherently through equilibrium flow in spherical tokamak [3–5]. Previous GK analyses [14] have reported that parallel flow shear can suppress sITG instability, which we also confirmed in Sec. 3.1. To enable future transport simulations, we performed the present analysis using a reduced fluid model to clarify the parameter dependence of PVG mode, particularly regarding finite $T_i$ effects.

Linear analysis revealed that the sITG and PVG modes become unstable in a mutually exclusive manner, as shown in Fig. 1. This indicates that the dominant growing modes can be distinguished based on underlying plasma profiles. Each mode also exhibits characteristic scales in $k$-space. It is well-established that the sITG mode peaks at around ${\rho }_s{k}_y$ < 1 [19], a trend reproduced by our model as shown in Fig. 4. By contrast, the PVG mode peaks at around ${\rho }_s{k}_y\sim \mathcal{O}(1)$, as shown in the right panel of Figs. 2 and 3. The effects of finite ion temperature originate from the temperature gradient ${\eta }_i$ and temperature ratio $\tau$. While an excessively steep temperature gradient ultimately stabilizes the PVG mode, it enhances the growth rate over a certain range. The temperature ratio contributes quadratically to the stabilizing terms and linearly, to the destabilizing terms, thereby influencing the linear PVG stability in the equations. These dependences are consistent with their underlying destabilization mechanism, i.e. negative compression, as discussed in Secs. 3.1 and 3.2. From this point, increasing $\tau$ strengthens the positive compression (stabilizing) effect, whereas increasing ${\eta }_i$ counteracts it until ${\eta }_i$ reaches sITG instability thresholds. This can aid in identifying the parameter space relevant to the PVG mode in transport simulations.

Fig. 4.  ITG growth rate in k-space. Relevant values are the same as in Fig. 2 except for the parallel velocity gradient. ${⟨v_z⟩}^{\prime }/{\omega }_{ci}$ = 0.0 is used here.

Notably, this analysis is restricted to the linear regime and therefore omits many important effects, making quantitative predictions difficult. However, this study provides a foundation for future studies on PVG-driven transport in tokamak. In particular, it identifies the characteristic scale of the PVG mode—peaking at ion-scales in $k$-space—clarifies finite $T_i$ effects by separating the contributions from $\tau$ and ${\eta }_i$. In future devices, parallel flow shear may play a role in the impure plasmas [20, 21], in the excitaion of breather wave discussed in [22] and in the coupled dynamics with drift-Alfvén turbulence in finite β plasmas [23, 24]. As the next step in this research, a transport analysis based on gyro-Landau fluid theory is currently underway.

Acknowledgement

This work was supported by JST SPRING, Grant Number JPMJSP2136, the Grants-in-Aid for Scientific Research of JSPS of Japan (Nos. 21H01066, 23K20838) and the joint research project in RIAM, Kyushu University.

References

• [1]   A.J.  Creely et al., J. Plasma Phys. 86, 865860502 (2020).
• [2]   S.M.  Kaye et al., Plasma Phys. Control. Fusion 63, 123001 (2021).
• [3]   W.X.  Wang et al., Phys. Plasmas 22, 102509 (2015).
• [4]   I.T.  Chapman et al., Nucl. Fusion 52, 042005 (2012).
• [5]   Y.  Kosuga, J. Phys. Soc. Jpn. 93, 064501 (2024).
• [6]   S.L.  Newton et al., Plasma Phys. Control. Fusion 52, 125001 (2010).
• [7]   I.  Oyama et al., Plasma Fusion Res. 19, 1403019 (2024).
• [8]   Y.  Kosuga et al., Phys. Plasmas 24, 032304 (2017).
• [9]   Y.  Kosuga et al., Plasma Fusion Res. 10, 3401024 (2015).
• [10]   T.  Kaneko et al., Phys. Rev. Lett. 90, 124001 (2003).
• [11]   S.  Inagaki et al., Sci. Rep. 6, 22189 (2016).
• [12]   T.  Kobayashi et al., Phys. Plasmas 23, 102311 (2016).
• [13]   S.  Mazzi et al., Nucl. Fusion 62, 096024 (2022).
• [14]   A.  Schekochihin et al., Plasma Phys. Control. Fusion 54, 055011 (2012).
• [15]   E.G.  Highcock et al., Phys. Rev. Lett. 109, 265001 (2012).
• [16]   J.  Citrin et al., Plasma Phys. Control. Fusion 59, 124005 (2017).
• [17]   G.M.  Steabler et al., Phys. Plasmas 14, 055909 (2007).
• [18]   G.M.  Steabler et al., Phys. Plasmas 17, 122309 (2010).
• [19]   A.M.  Dimits et al., Phys. Plasmas 7, 969 (2000).
• [20]   J.  Bourgeois et al., Open Plasma Science 1, 1 (2024).
• [21]   G.  Cuerva-Lazaro et al., Plasma Phys. Control. Fusion 67, 075019 (2025).
• [22]   Y.  Kosuga et al., Plasma Phys. Control. Fusion 66, 075018 (2024).
• [23]   V.V.  Mikhailenko et al., Phys. Plasmas 23, 092301 (2016).
• [24]   V.V.  Mikhailenko et al., Phys. Plasmas 27, 112103 (2020).