Gyrokinetic simulation of short wavelength ion temperature gradient instabilities in the ADITYA-U tokamak

In this work, linear and nonlinear collisionless electrostatic simulation studies of the standard and short wavelength ion temperature gradient mode (SWITG) for experimental profiles and parameters of ADITYA-U tokamak are performed using the linear global eigenvalue gyrokinetic code GLOGYSTO and the nonlinear global gyrokinetic particle-in-cell code ORB5. All simulations are carried out with non-adiabatic ions and adiabatic electrons. The ADITYA-U tokamak which has recently been upgraded to divertor configuration, is small in size and well suited for investigation of micro-instabilities in the presence of density and temperature gradients. Due to steep density and temperature gradients, simulation shows that the SWITG mode naturally exists along with the standard ion temperature gradient (ITG) mode in ADITYA-U. In this work, the experimental shot# 33536 of the ADITYA-U tokamak is considered as a reference. There is good agreement in the growth rate and the real frequency values between GLOGYSTO and ORB5 with variations of less than 25% . Two maxima of growth rate versus mode number are obtained, the first around kθρs≃0.4 is the standard ITG, the second around kθρs≃1.2 is the SWITG. Additionally, using linear stability analysis, it is observed that the SWITGs are suppressed for low values of R0/LT i.e. only the standard ITG mode remains unstable. For the ADITYA-U tokamak, nonlinear global simulations with ORB5 are also carried out. Nonlinearly, SWITG dominating case results are compared with the conventional ITG case, where SWITG is relatively suppressed. The nonlinear contribution of the SWITG mode to the total thermal ion heat transport is found to be minimal due to an increased zonal flow shearing effect on the SWITG mode suppression, even though it may be linearly more unstable than the conventional long wavelength ( kθρs<1 ) ITG mode.


Introduction
One of the crucial concerns in current fusion research is the mitigation of anomalous transport to enable improved plasma confinement. In modern tokamaks [1], high confinement mode plasmas are frequently created. The creation of pedestals and transport barriers in this scenario, which minimize the transport of particles and heat from the system, is a significant aspect. Strong pressure gradients throughout the pedestal's plasma make it a hotspot for a variety of instabilities of different flavours. Similar to this, the plasma also displays sharp gradients across internal transport barriers, which gives free energy to small-scale instabilities [2]. In tokamak experiments, these barriers are produced in a self-consistent way due to nonlinear interactions between multi-scale structures. The nature of the conventional instabilities may vary when strong gradients are present. Instabilities at low frequency as compared to the ion gyro-motion frequency and scale length comparable to the ion Larmor radius instabilities [3][4][5][6][7][8] are thought to be the cause of confinement's degradation, which results in anomalous transport of energy and particles. The density and temperature inhomogeneities that are present in a magnetically confined plasma are the source of free energy for these modes. For example, it is observed that even at wave lengths k θ ρ s > 1.0, the ion temperature gradient (ITG) mode, which is driven by the temperature gradient of ions, becomes unstable, when the background gradients (density and temperature) are extremely sharp [3,4]. These background gradients tend to drive short scale ITG modes (SWITGs) unstable [3][4][5][6]. Similarly, trapped electron modes can also exhibit a shorter wavelength branch when significant gradients are present [9,10]. Using gyrokinetic simulations, it has been shown that this short wavelength branch of micro-instabilities is vital for experimental parameters [11].
Therefore, it is crucial to comprehend these modes' nonlinear characteristics and their contribution to the anomalous transport of energy and particles. The contribution of the shorter-scale ITG modes to the overall heat flux is weaker than that of the standard wavelength ITG branch, as determined by a nonlinear flux tube simulation of SWITG modes with adiabatic electrons [5]. In this case, local flux tube calculations might not be applicable because the region across the steep gradients might be relatively small. This is due to the fact that the local computations are built on the assumption that the equilibrium quantities and fine-scale fluctuations have disparate scale-lengths and are hence separable. However, in the regions with narrow gradients, such an assumption breaks down. This therefore necessitates the use of global computations that address small scale fluctuations and large scale equilibrium variations on an equal basis. Therefore, it is crucial to use a global gyrokinetic solver to examine how this mode acts nonlinearly and determine whether it contributes significantly to the net ion transport in the core of the system. The self-consistent dynamics of SWITGs driven by extremely sharp background gradients, their linear and nonlinear evolution and their saturation after the onset of zonal flows are addressed in the present work, possibly for the first time using a global, gyrokinetic, electrostatic solver that includes adiabatic electrons and non-adiabatic ions. The scope of the paper is not a direct comparison with the experiment but is merely a first step, focusing on strong gradient cases. More realistic modelsin particular, non-adiabatic electron response, will be studied in future works. In order to do this, we perform a systematic linear and nonlinear analysis of the mode for ADITYA-U [12] using ORB5 [13,14] and GLOGYSTO [15,16]. However, nonlinearly, we observe that the SWITG mode only makes a very minor contribution to the total amount of ion thermal transport while having a growth rate that is linearly comparable to that of the conventional ITG mode.
The present manuscript is arranged as follows. Section 2 briefly describes the simulation models used, namely ORB5 and GLOGYSTO. In section 3, ADITYA-U profiles, plasma parameters used in the present simulations and the results and discussion are described. Finally, conclusions are drawn in section 4.

Simulation model
In the following, we briefly describe the global nonlinear gyrokinetic particle-in-cell (PIC) code ORB5 [13,14] and global linear eigenvalue gyrokinetic code GLOGYSTO [15,16]. Note that the linear SWITG mode for ADITYA-U is studied using ORB5 and GLOGYSTO while the nonlinear SWITG study is carried out using the ORB5 code in the present work.

ORB5 Gyrokinetic model
For an axisymmetric toroidal plasma, the Vlasov-Poisson system is solved using the ORB5 code in the gyrokinetic limit. Labels for magnetic surfaces are either the poloidal flux ψ or the radial coordinate s = √ ψ ψ edge . The magnetic field is expressed as B = I(ψ)∇φ + ∇ψ × ∇φ, Here, I(ψ) is the poloidal current flux function, and φ is the toroidal angle. Circular concentric magnetic surfaces (ad hoc equilibrium) and real MHD equilibria are two different types of magnetic equilibria that are implemented in ORB5. For the real MHD equilibrium, the Grad-Shafranov equation with a fixed plasma boundary is solved using the CHEASE code [17] and it is coupled with the ORB5 code. In ORB5, a straight-field-line coordinate system is used with the poloidal coordinate θ ⋆ defined by where q(s) is the safety factor and θ is the geometrical poloidal angle. The 5D phase space (R, v ∥ , µ) is sampled in ORB5 using markers or pseudo-particles for charges. In the absence of collisions, the magnetic moment (µ) of a marker is constant and the equations are independent of the gyrophase. The guiding center position is represented by R, v ∥ and v ⊥ , respectively, stand for the parallel and perpendicular components of the velocity with respect to the magnetic field. The dynamics is determined by using a gyrokinetic framework for ions, where the time evolution of the particle distribution function is given by: where the gyro-position center's in real space is X, v ∥ , its velocity parallel to the equilibrium magnetic field B = Bb, C accounts for the effects of collisions and S accounts for the sources [18,19]. Noting that this work is restricted to electrostatic perturbations and collisionless (C = 0) regime despite the fact that ORB5 can account for both electromagnetic fluctuations [20] and collisions [18,19]. In ORB5, a PIC scheme with control variates is used to reduce statistical sampling noise. It consists in formally splitting the full distribution f into a time-independent part f 0 and a time-dependent part δf.
Only the δf part is discretized using numerical particles, called markers, while the f 0 terms are treated analytically. As indicated earlier, S stands for the sources that can be added to regulate numerical noise and/or maintain temperature and density profiles; S = S k + S h , where S h is a heating source term (considered zero in the present simulation) and S k is a Krook operator. In ORB5 [21], a Krook operator is employed with corrections for energy, momentum, and zonal flow (ZF) conservation; however, in the present simulation, a Krook operator is used without any corrections. Therefore, all of the following equations are given in this limit. For a species with mass m and charge q, the equations of motion are given bẏ where v ∇B = (µ/ (mΩB)) B × ∇B represents the grad-B drift velocity, v E×B = ( 1/B 2 ) B × ∇⟨ϕ ⟩ G represents the E × B drift velocity, and v c = (v 2 ∥ /Ω) (∇ × b) ⊥ denotes the curvature drift velocity. Here, ⟨ϕ ⟩ G represents the gyroaverged electrostatic potential. Finally, the effective magnetic field (B ⋆ ) can be written as The system of equations (2)-(4) is closed by the gyrokinetic Poisson equation which assuming a quasi-neutrality, reads as where δn α is the perturbed density and the sum extends over all α plasma species. Gyro-density and polarisation density are the two components that are produced when δn α is evaluated as the zeroth order moment of the corresponding distribution function; the latter is caused by the difference between the coordinates of the guiding centre and the gyro-center, i.e., a result of finite amplitude field fluctuations. As a result, equation (6) is transformed into a linear integral equation for the electrostatic potential. The ion polarization term that appears in the quasi-neutrality equation has been linearized and can be expressed in various ways in ORB5: (a) a longwavelength approximation [13] (b) a Pade approximation [14,22,23] and (c) a solver valid at arbitrary order in k ⊥ ρ [24]. In equation (7), ⟨ϕ ⟩ FS represents a flux-surface-averaged electric potential defined as´ϕ is the Jacobian of the straight-field-line coordinate system (s, θ ⋆ , φ). The solver valid at arbitrary order is used in the current work because it is our goal to explore both conventional ITG and SWITG modes (0.0 ⩽ k θ ρ s ⩽ 1.4). Four parameters-the ion mass (m i ), the ion charge (q i = Z i e, Z i the atomic number and e the electric charge), the electron temperature at a specified reference position s 0 , T e (s 0 ), and the magnetic field on axis (B 0 ) are used to normalise all the quantities in the code. These parameters are used to determine all other normalized quantities. The time units are given as the inverse of the ion-cyclotron frequency, Ω ci = q i B 0 /(m i ), the velocity units are normalized through the sound velocity of the ions (c s = √ eT e (s 0 )/m i , the temperature T e (s 0 ) is in eV), the length units through the ion sound Larmor radius (ρ s = c s /Ω ci ), and the densities are normalized by their averaged value in space.

GLOGYSTO gyrokinetic model
With the use of the Nyquist method, the GLOGYSTO code [15,16], a global spectral code, determines the eigen mode structure as well as the real frequency and growth rates of unstable modes for a given equilibrium. Ion and electron species are both taken as fully gyrokinetic. However, in the present simulation, electrons are considered adiabatic. The sum of the adiabatic and non-adiabatic components can be used to express the perturbed electrostatic density for a species α as follows: The charge and temperature for the species α are represented by q α and T α , respectively, while N accounts for the equilibrium density. The perturbed electrostatic potential is represented byφ. The diamagnetic drift frequency is given by and v thα is the thermal velocity of species α. The full finite Larmor radius effect is taken into account in the Bessel func- The safety factor is defined as q = m n , m and n are poloidal and toroidal wave numbers, respectively, B p is the poloidal magnetic field, and k θ is the poloidal wave vector. The local Maxwellian f Mα for the species α of mass m α is given by N(ψ) and T(ψ) are the density and temperature at a given magnetic surface. The term U α represents the guiding center propagator for the passing particles. Using the quasi-neutrality condition, it provides the following equation ∑ αñ α (r; ω) = 0.
This results in an eigenvalue problem that can be solved in Fourier space, in the case of electrostatic fluctuations only, where ω andφ are the eigenvalues and eigenvectors, respectively. For fully gyrokinetic ions and electrons with only passing particles, the following equation is given: HereM α k , k ′ is the convolution matrix [6] in Fourier space and it has been given in the appendix B. The radial and poloidal wave numbers κ and m, respectively, are represented by the wave vector k = (κ, m) when the toroidal mode number n for an axisymmetric system is fixed. Thus, k = (κ, m) and k ′ = (κ ′ , m ′ ). Here, m ′ = nq(s 0 ) ± δm, where q(s 0 ) denotes the q value at s = s 0 and δm determines the range of poloidal mode numbers.

Results and discussion
In this section, we will discuss our results for both nonlinear runs with ORB5 and linear runs with GLOGYSTO and ORB5.

Linear gyrokinetic simulations with ORB5 and GLOGYSTO
The ADITYA-U tokamak [12,25,26] which has recently been upgraded to divertor configuration [26], is small in size and well suited for the investigation of micro-instabilities in the presence of density and temperature gradients. The gradient in the plasma profile can drive several temperature and density gradient instabilities, such as ITG [6,[27][28][29][30], trapped electron mode (TEM) [31,32], universal drift modes [33,34] etc which are electrostatic in nature. Similar profile gradients may also produce electromagnetic instabilities such as kinetic ballooning mode (KBM) [35][36][37][38][39][40] and micro tearing mode [41][42][43][44][45] if the plasma β is high [46,47]. In the present study, the core plasma beta β = 2µ 0 n 0e T 0e /B 2 0 , where T 0e , n 0e and B 0 are the on-axis electron temperature, density and the magnetic field strength with the values of n 0e = 2.3 × 10 19 m −3 , T 0e = 250 eV and B 0 = 1.0 T is approximately 0.1%. Furthermore, the threshold β for the KBM is higher than the β value for the current discharge of ADITYA-U for hydrogen plasma [48]. At such a low β value, the microturbulence responsible for the heat and particle transport in a tokamak is expected to be electrostatic [39,47,49].
Hence, in this work, we perform global linear and nonlinear electrostatic simulation studies of the conventional (ITGs) mode and SWITGs, which are both unstable, resulting in a multi-scale gyrokinetic transport, for profiles and parameters relevant to the ADITYA-U tokamak. Conventionally, in ADITYA-U discharges, it is observed that the electron temperature is almost three times greater than the ion temperature [25]. In the presence of external auxiliary heating mechanisms, it is expected that the ion temperature would become comparable to the electron temperature, i.e, T i = T e . Figure 1 shows the plasma profiles and the corresponding normalized gradient (R 0 /L g ) used in the simulations, where L g is the profile gradient length scale of any quantity g(r), given by L g = − (d ln g/dr) −1 , r is the local minor radius. We first benchmark the ORB5 code with the familiar Cyclone DIII-D base case (CBC) [50] for ITG modes with adiabatic electrons. This is presented in the appendix A. After this exercise, we carry out linear simulations using GLOGYSTO and ORB5 codes for parameters and the density profile of shot# 33 536 of the ADITYA-U tokamak [12,25,26]. However, as indicated earlier, we consider temperature profiles with T i = T e . The simulations are performed on spatial grids N s = 448, N θ⋆ = 512, N φ = 256, (s, θ ⋆ , φ) representing the radial, poloidal, and toroidal directions. The number of markers is N p = 15 × 10 8 , the time step is ∆t = 10Ω −1 ci and ρ * = ρ s /a = 0.00365. The ITG instability is measured at s 0 = 0.6, where the temperature and density gradients peak. For the linear simulation 25 particles per cell are used [51]. As described earlier, in ORB5, a field-aligned Fourier filter [13,14] is invoked using θ ⋆ and φ coordinates. In line with the gyrokinetic order, i.e. k ∥ /k ⊥ ≪ 1, the effective number of markers per Fourier mode in the filter and per radial mode then becomes Np Ns(2×∆m) ∼ 3 × 10 5 (∆m is the width of the field-aligned filter), for the parameter used. Each linear simulation corresponds to a single toroidal mode number. Accordingly, k θ ρ s is computed by using the relationship k θ ρ s = nq r0 ρ * a/s with s = 0.6a and q s0 is the value of the safety factor at s 0 . The safety factor, shear, temperature, density, temperature and density gradient scale lengths and η = L n /L T profiles for ADITYA-U circular plasma are depicted in figure 1. Here, we use the circular concentric model [13, 52] to describe the equilibrium even though it is not a real MHD equilibrium. This model is very useful for numerical studies since it offers straightforward formulations for geometrical and magnetic field parameters. Flux-tube and global gyrokinetic studies using this model have found linear growth rates and a heat diffusivity fairly similar to those determined using an ideal MHD equilibrium [52,53]. This equilibrium with circular cross-section magnetic configurations is defined in radial(r), poloidal(θ), and toroidal(φ) coordinate system, which is related to the cylindrical coordinate-system (R, φ, Z) by the relationships R = R 0 + r cos θ = R 0 (1 + ϵ cos θ) and Z = r sin θ, here ϵ = r/R 0 is the inverse aspect ratio and R 0 the major radius in the cylindrical co-ordinate system, but keeping all geometrical terms consistent, which was shown in order to get quantitatively correct results [52]. The magnetic flux is chosen as circular concentric surfaces defined by dψ dr = B0r q(r) , with a specifiedq(r) profile given in table 1 The required plasma parameters and profiles for the simulation are taken from [25,26]. The profiles and parameters that are considered in the simulation of ADITYA-U are given in table 1 for which R 0 /L n = 5.52 and R 0 /L T = 26.8. The real frequencies and growth rates are calculated using GLOGYSTO and ORB5 for ADITYA-Us for different toroidal mode numbers and are shown in figure 2.
As seen in figure 2 for R 0 /L T = 26.8 case, the growth rate has two peaks rather than the single peak often seen in linear conventional ITG modes, which is around k θ ρ s ≃ 0.4 (shown by the orange dashed vertical line). The second peak appears around k θ ρ s ≃ 1.2 (shown by the blue dashed vertical line) and is characteristic of the SWITG mode [3-6, 54, 55]. For k θ ρ s ⩽ 1, the real frequency increases monotonically with k θ ρ s , but then stays almost constant for 1 ⩽ k θ ρ s ⩽ 2.
The following dispersion relation for SWITG [56] for adiabatic electrons is obtained by the quasi-neutrality equation in the context of a local gyrokinetic formulation in the limit here τ = T e /T i , I 0 is the modified Bessel function of order zero, ω * i = − (v thi /L n ) (k θ ρ s ) is the ion diamagnetic drift frequency, and L n is the density scale length. It can be seen from the expression of the dispersion relation, equation (12) that the mode frequency ω scales as k θ ρ s for small k 2 θ ρ 2 s Table 1. Parameters and equilibrium profiles.
Parameters: Equilibrium Profiles: • B-field : B 0 = 1.0 Tesla, mi me = 1836 • N-profile and T-profile: shearŝ is positive and at s = s 0 ,ŝ = 1.  and for larger k 2 θ ρ 2 s (k θ ρ s ≫ 1) scales almost as constant as The toroidal magnetic drift term ω di of the ions resonates with the mode frequency ω and results in the double hump behaviour in the growth rate in toroidal geometry. It is crucial to remember that the finite Larmor radius effects also have an impact on the SWITG mode [56,57].
From figure 2, we also see that there is, in general, a good agreement in growth rate and real frequency values between ORB5 and GLOGYSTO with quantitative differences that are less than 25%. Figure 3 shows the mode structure obtained from ORB5 for toroidal mode numbers n = 35 (first peak) and n = 100 (second peak), respectively. With GLOGYSTO, we obtain the shifted peaks of the growth rates compared to ORB5, which are located at toroidal mode numbers n = 40 (first peak) and n = 105 (second peak), respectively, as shown in figure 2. The mode structures obtained using GLOGYSTO for toroidal mode numbers n = 40 (first peak) and n = 105 (second peak), respectively, are also shown in figure 4. For comparison with a case without SWITG, we have looked at R 0 /L T = 13.1 keeping R 0 /L n fixed. Both ORB5 and GLOGYSTO are run for R 0 /L T = 13.1 case to make sure that SWITG is actually suppressed, and it does indeed suggest that for R 0 /L T = 13.1, SWITG contribution would be small as shown in figure 2. The growth rate peaks at n = 30 (conventional ITG peak) as obtained from ORB5. The mode structures of R 0 /L T = 13.1 case for toroidal mode number n = 30 (conventional ITG peak) using ORB5 is shown in the left panel of figure 5. Also in this case, we find a shifted peak of growth rate with GLOGYSTO, which is located at n = 35 (conventional ITG peak). The right panel of figure 5 depicts mode structure using GLOGYSTO. The mode structures obtained from ORB5 and GLOGYSTO exhibit good agreement. Additionally, for R 0 /L T = 13.1 case, we observe from figure 2 that there is generally reasonable agreement between ORB5 and GLOGYSTO in terms of growth rate and real frequency values, with quantitative deviations that are less than 10%.
Before discussing our nonlinear simulations, it is crucial that the zonal flows respond correctly in nonlinear runs because of ZFs' significant contribution to the evolution of plasma turbulence and their ability to effectively suppress plasma turbulence and reduce turbulent transport due to their low frequency characteristic [58]. Therefore, the Rosenbluth-Hinton test for ADITYA-U using ORB5 is performed and shown in the appendix C, which was analytically detailed in [59], is the common benchmark [13,48,60,61] that addresses the physics of the zonal flow/geo-acoustic mode.
Following the linear simulation, in the below section we present nonlinear global simulations with ORB5 for R 0 /L T = 26.8 and R 0 /L T = 13.1 keeping R 0 /L n = 5.52 fixed.

Nonlinear gyrokinetic simulations with ORB5
We will outline our findings for gradient-driven nonlinear runs in this section. The following numerical parameters are chosen the number of markers is N p = 15 × 10 8 , which corresponds to ∼ 2360 particles per Fourier modes and per radial interval using the expression N p /(N s (2∆m + 1)(n max + 1)), where ∆m is the width of the field-aligned filter and n max is the maximum value of the toroidal mode number. This leads to a signal-to-noise ratio sufficiently high (above 80) for all nonlinear simulations shown later in this paper. The 3D grid resolution for the field solver is N s × N θ⋆ × N φ = 448 × 512 × 256. Note that we have 6 × 10 7 grid points, resulting in 25 particle/cell on average and a 32-point averaging technique is employed for gyro-averaging. The arbitrary wavelength field solver is used for the quasi-neutrality equation and the time step size is ∆t = 10Ω −1 ci . We conduct a convergence test for the grid and time steps, and the selected number of grid points and time steps are sufficient for nonlinear simulation. Toroidal   Γ is the gyrocenter particle flux, Q kin is the kinetic energy flux, Q pot is the potential energy flux, and q H is the heat flux.
Heat diffusivity χ is defined as The gyro-Bohm heat diffusivity is defined as χ GB = ρ 2 s0 c s0 /a, with a is the minor radius, ρ s0 the sound Larmor radius of the ion, and c s0 the sound speed of the ion where the temperature is taken at the reference surface s 0 and the magnetic field on the axis. Using the plasma parameters taken in this work, χ GB = 0.29 m 2 s −1 .
The E × B ZF shearing rate is given by [62].
3.2.1. SWITG dominant case: R 0 /L T = 26.8. All of the simulations reported in this paper have a signal-to-noise ratio that is high enough, S/N > 80, as shown in figure 6(a). Further information on the S/N estimation can be found in [63]. With a decay rate γ K = 2.0 × 10 −4 Ω ci , the Krook operator, as described in [21], is used to control noise over long time scales without corrections for energy, momentum and ZF conservation [64]. The time evolution of the corresponding average particle weight squared, ⟨w 2 ⟩, is calculated with a different number of markers, 15 × 10 8 and 45 × 10 8 , respectively, as shown in figure 6(b). The average particle weight squared ⟨w 2 ⟩ in both cases is almost identical and does not depend on whether 15 × 10 8 or 45 × 10 8 markers are used in these simulations. The particle weight is found to saturate during the simulations and remains at a low level (⟨w 2 ⟩ < 0.12) throughout the simulations. These convergence tests clearly show that the noise-induced transport in our simulations is insignificant compared to the turbulence driven transport. Our simulations are therefore 'gradient-driven', however small variations are permitted, as will be seen in the results presented below. The Krook source term acts as a source of heat. As a consequence, while the temperature profile is allowed to relax due to the turbulent heat flux, it is maintained close to the given initial profile. In the long timescales, the time-averaged temperature gradients reach a quasi steadystate. In figure 7(a), the temperature profile relaxation, computed using equation (19) is shown. Since the adiabatic electron model is used, there is no particle transport. In figure 7(b) initial and final R 0 /L T profiles are shown.
The plasma volume is split into N b radial bins, with each having a volume of V j , j ∈ [1, N b ], in order to compute the profiles. The evolving plasma temperature and density are computed as follows for species α , with f α = f 0α + δf α , using the expressions [65]: where v ∥ is the parallel velocity along the magnetic field and v ⊥ is the perpendicular velocity across the magnetic field. f 0α is the equilibrium distribution function and δf α is the perturbed distribution function. From figure 8(a), we can see an initial exponential growth phase of the electrostatic field energy as a function of time, followed by a nonlinear saturation phase. The electrostatic field energy E field , defined as where ⟨ ⟩ G indicate a gyro-averaged quantity and n 0 is the equilibrium density. From figure 8(a), it can be seen that after the initial linear phase the mode amplitude tends to saturate from time t ∼ 20.0 × 10 3 Ω −1 ci . The normalized field energy is given as E field /(⟨n⟩VT) and the normalized heat flux as Q/(⟨n⟩c s T e (s 0 ), where V is the plasma volume, n is the density and ⟨n⟩ its averaged value in space and T e = T e (s 0 ) is the electron temperature at s 0 = 0.6. Figure 8(b) displays the volumeaveraged heat flux Q, defined by where V is the volume of the torus. It is clear from the figure that the heat flux initially increases exponentially in the linear phase and peaks around t ∼ 1.0 × 10 4 Ω −1 ci . However, as time progresses the simulation enters the nonlinear phase around t = 2.5 × 10 4 Ω −1 ci where the zonal-flow sets in. Due to the interaction between turbulence and zonal flow, the overall heat flux reduces and tends towards a steady state. The timeaveraged heat flux in the steady state is around 0.04 in the normalized unit. As mentioned in appendix C, when conducting nonlinear simulations, it is crucial to accurately describe the structure of the zonal flow component i.e., .(n = m = 0), as it is crucial for the nonlinear saturation in the ITG regime. Due to its ability to shear radially coherent turbulent formations, zonal flows have an impact on microturbulence. This effect is governed by the shearing rate ω E×B , which is determined by equation (18). The shearing rate is proportional to the first radial derivative of the electric field or the second radial derivative of the electrostatic potential. The structure of turbulent eddies in the nonlinear phase can be seen in figure 8(c) taken at t = 2 × 10 5 Ω −1 ci , where the ϕ − ⟨ϕ ⟩ FS (where ⟨ϕ ⟩ FS is the flux surface averaged potential) is plotted in a poloidal cross section. From the snapshot, it is clear how zonal flow shearing affects the potential. The zonal flow tears the global structures and regulates the turbulence.
In figure 9(a), we plot the simulated heat flux as a function of the radial coordinate and time. Figure 9(b) shows the spatio-temporal evolution of the ion heat diffusivity χ i . From figure 9(b), we can see that a number of bursts (localized features in space and time) as well as avalanches frequently occur following an initial turbulence overshoot commencing near s ∈ [0.55 − 0.75], which is in the region of maximum R 0 /L T . The averaged heat diffusivity over the time interval tΩ ci = [1.0e5 − 2.0e5], and the radial interval s = [0.5 − 0.7] is χ i /χ GB = 6.04. In physical units, this is χ i ≃ 1.7 m 2 s −1 . Figure 10 shows the spatio-temporal behaviour of the electrostatic field energy, equation (24), separately for its zonal (top panel) and non-zonal (bottom panel) components. The mode intensity peaks at time t ∼ 1.0 × 10 4 Ω −1 ci and around s = 0.65. The turbulence exists over a wide radial domain approximately from s = 0.5 to s = 0.8. The zonal flow caused by the turbulence reaches its peak, as expected, in the region of significant ITG and SWITG activity, as can be seen in the upper panel of figure 10. As a result, a saturation phase of turbulence is reached. The radial profiles of R 0 /L T at initial (red dashed line) and final (magenta dashed line) times are shown in figure 10. The equilibrium R 0 /L n profile is depicted by the white dashed line and η = L n /L T at initial (orange dashed line) and final (green dashed line) are displayed in figure 10. From figure 10, we can also see that the turbulent energy and ZF energy have their maximum amplitude pushed a bit outwards from the position of maximum log gradient. It looks much more active in the interval s ∈ [0.6, 0.8]. This can be attributed to temperature profile relaxation.  Here, we consider the same density profile as before, R 0 /L n = 5.52, but a lower temperature gradient profile with R 0 /L T = 13.1 instead of 26.8. The time evolution of the electrostatic field energy is shown in figure 11(a). It is clear that the initial exponential growth phase of the electrostatic field energy as a function of time is followed by a nonlinear saturation phase. After the initial linear phase, the mode amplitude tends to saturate from time t ∼ 25.0 × 10 3 Ω −1 ci . Compared to figure 8(a) the mode amplitude increases slowly due to the reduced strength of the gradients. Figure 11(b) displays the volume-averaged heat flux. The heat flux peaks around t ∼ 1.5 × 10 4 Ω −1 ci and then enters the nonlinear steady state phase around t = 2.5 × 10 4 Ω −1 ci . The time averaged heat flux in the steady state is around 0.02 in the normalized unit. The steady state heat flux is quite small we also see that the turbulence is more global (wide spread in radial direction) as compared to R 0 /L T = 26.8 case.
In figure 12(a), a spatio-temporal plot of the heat flux is depicted. The spatio-temporal evolution of the ion heat diffusivity χ i is depicted in figure 12 The upper and bottom panels of figure 13 show the spatiotemporal behaviour of the ZF and turbulence energy. The mode intensity peaks at time t ∼ 1.5 × 10 4 Ω −1 ci and around s = 0.65. The turbulence exists and spreads over a wide radial domain approximately from s = 0.5 to s = 0.8. As could be expected, the zonal flow caused by the turbulence reaches its peak in the region of the strong ITG activity, as can be seen in the upper panel of figure 13. As a result, the turbulence enters into a saturation regime. It can also be seen from the figure that ZF and turbulence energy are lower for R 0 /L T = 13.1 than the one for R 0 /L T = 26.   also equally unstable. As we also see from figure 14 that the zonal flow shear tears the global structures to regulate the turbulence for both R 0 /L T = 26.8 and R 0 /L T = 13.1 cases. Therefore, it is essential to investigate the role of the zonal flow shearing rate [58,66,67] in both cases.
In order to achieve this, in both simulations (R 0 /L T = 13.1 and R 0 /L T = 26.8), we have determined the zonal flow shearing rate. In both situations, the shearing rate and temporal evolution are depicted in figure 15. The time and radially averaged shearing rate ω tot E×B = ⟨⟨|ω E×B |⟩ s ⟩ t is 0.0033Ω ci for , respectively. The quantity ω tot E×B is a measure of the total absolute value of the ZF shearing rate. In the case of R 0 /L T = 26.8, the zonal flow shearing rate for the SWITG mode is approximately twice as high as in the case of R 0 /L T = 13.1. The primary mechanism for saturating the SWITG mode turbulence is zonal flows, as shown by the fact that the shearing rate is higher in both cases than the linear growth rate. It has also been found that ZF shearing is a generally accepted mechanism for turbulence suppression [68]. Maximum linear growth rates of the SWITG mode are 0.00214Ω ci and 0.00065Ω ci , respectively, for R 0 /L T = 26.  conventional ITG mode is taken into account and the SWITG modes are artificially suppressed from the simulations. The aim is to determine the role of SWITG modes on turbulent heat transport.
We can see from figures 16(a) and (b) that the contribution of the SWITG mode to the normalized heat flux is very small when compared to the contribution of the standard ITG mode (the region 0.0 ⩽ k θ ρ s ⩽ 0.8, which mostly contributes to the transport) in both k θ ρ s scenarios because we obtain roughly the same heat fluxes 0.024 for R 0 /L T = 13.1 and 0.039 for R 0 /L T = 26.8, shown in green and magenta dashed lines, respectively, whether SWITG mode is taken into account or not. One can calculate the difference between the higher k θ tail, k θ ρ s ⩾ 1.0, which corresponds to the SWITG mode and the lower k θ component, k θ ρ s ⩽ 1.0, which is related to the conventional ITG, in terms of their relative contributions to the net ion heat transport.
where Q frac, SWITG is the relative contributions to the net ion heat transport. From equation (26), we obtain, Q frac, SWITG = 0.032 for R 0 /L T = 13.1 and Q frac, SWITG = 0.035 for R 0 /L T = 26.8. Thus, one finds that despite the SWITG mode for R 0 /L T = 26.8 having a linear growth rate that is more than thrice that of the SWITG mode for R 0 /L T = 13.1 (see figure 2), the SWITG branch of the k θ spectrum's net contribution to the total ion heat flux is approximately 3.5 % in both cases. As a result, the lower k θ components of the fluctuation related to the conventional ITG mode dominate in determining the thermal ion heat transport and the SWITG mode only makes a very small contribution to the overall ion thermal transport. Mixing length estimation of heat diffusivity is also carried out. The estimation of the heat diffusion coefficient for conventional and short-wavelength ITG from mixing length arguments (quasi-linear theory) is given as follows: χ ∼ γL , where k 2 ⊥ = k 2 r + k 2 θ and γ L is the linear growth rate; here we have used k ⊥ ≈ k θ , which is a valid approximation in ITG relevant scales in the absence of non-adiabatic electrons [69].
In figure 9(b), the total heat diffusivity calculated from nonlinear simulations including mode numbers up to the  short-wavelength range (ITG+SWITG) is found to be χ NL ITG+SWITG = 1.7 m 2 s −1 . Total heat diffusivity including mode numbers up to the conventional ITG range (n ⩽ 70) is χ NL ITG = 1.64 m 2 s −1 . Hence, the ratio of the total heat diffusivity coming from the ITG+SWITG branch to the conventional ITG branch is χ NL ITG+SWITG /χ NL ITG = 1.04. The current heat transport coefficient is estimated from the mixing length calculation for all the toroidal mode numbers using the linear calculations as shown in figure 17. For the conventional ITG mode (denoted by the orange dashed vertical line at k θ ρ s ≃ 0.4 in figure 17), the transport coefficient is found to be ≃ 0.92 m 2 s −1 whereas the same for the SWITG mode (denoted by the blue dashed vertical line at k θ ρ s ≃ 1.2 in figure 17) is ≃ 0.12 m 2 s −1 . Hence, from quasi-linear theory, the ratio of the diffusivity coming from the conventional ITG branch to that coming from the SWITG branch is about 7.7, whereas it is 1.04 from the nonlinear simulations. This shows that the mixing length theory is unable to explain the small contribution from SWITG and overestimates it. This necessitates the requirement of nonlinear simulation. From the present linear and nonlinear simulations of both cases, the following conclusions emerge: 1. In the linear simulations, ORB5 and GLOGYSTO exhibit good agreement in growth rate and real frequency values, with quantitative differences that are less than 25%. 2. We find two peaks of the linear growth rate versus mode number as opposed to the single peak that is often seen in the linear analysis of the conventional ITG modes. The second peak, which is characteristic of the SWITG mode, appears at k θ ρ s ≃ 1.2. 3. Using linear analysis, it is also observed that the SWITGs are suppressed at low values of R 0 /L T , implying that only the conventional ITG mode remains unstable. 4. The nonlinear SWITG mode for R 0 /L T = 26.8 exhibits a higher zonal flow shearing rate than the one for R 0 /L T = 13.1. 5. The nonlinear contribution from the SWITG modes to the net ion heat flux is very small in both cases, despite the fact that the growth rate of the SWITG mode is significantly higher for R 0 /L T = 26.8 than for R 0 /L T = 13.1 in the linear phase. This implies that the bulk part of the ion heat fluxes comes from the standard ITG mode.
The significance of secondary instabilities [70,71] can be used to clarify the fifth point. The secondary instability caused by the E × B shear scales as k 2 θ Φ L [72], where Φ L is the amplitude of the linear SWITG mode. It should be noted that the values of k θ 's where the linear growth rates of the conventional ITG and SWITG modes peak are approximately 0.4 and 1.2, respectively (see figure 2). The corresponding growth rates also have similar magnitudes. Due to the fact that secondary instabilities' nonlinear growth rates (∝ k 2 θ ) exceed their linear growth rates, hence SWITG modes are sub-dominant nonlinearly. The strength of the secondary instability is weaker in the conventional ITG component of the k θ spectrum than in the SWITG part because the conventional ITG part's k θ is smaller. Since SWITG has a far stronger suppression level than ITG in the nonlinear regime, conventional ITG contributes the majority of the total ion heat flux. Therefore, the SWITG mode's contribution to the total ion heat flux for R 0 /L T = 26.8 is practically as low as that for R 0 /L T = 13.1.

Summary
In the present work, linear and nonlinear simulation investigations of the standard and SWITG mode for experimental profiles and parameters of shot# 33 536 [but setting T i = T e ] of ADITYA-U tokamak have been carried out using the global nonlinear gyrokinetic PIC code ORB5 and global linear eigenvalue gyrokinetic code GLOGYSTO. The growth rate and real frequency are calculated using GLOGYSTO and ORB5 for different toroidal mode numbers. We see that there is a good agreement in growth rate and real frequency values between GLOGYSTO and ORB5. As a consequence of the linear analysis, we found two peaks for R 0 /L T = 26.8 case as opposed to the single peak that is often seen in the linear analysis of the conventional ITG modes. The SWITG mode is represented by the second peak, which appears around k θ ρ s ≃ 1.2. The SWITGs are found to be suppressed for low values of R 0 /L T , meaning that only the standard ITG mode remains unstable, using linear stability analysis.
Nonlinear global simulations using ORB5 for ADITYA-U tokamaks are also performed for R 0 /L T = 26.8 and R 0 /L T = 13.1 (to check what transport looks like if SWITG is suppressed) keeping R 0 /L n fixed. Higher heat fluxes are observed in nonlinear simulations for higher values of R 0 /L T , which is in accordance with the linear growth rate's trend with respect to R 0 /L T . From the time evolution of the heat flux, we observed that the contribution of the SWITG mode to the overall heat flux is very small when compared to the contribution of the standard ITG mode (the region 0.0 ⩽ k θ ρ s ⩽ 0.8, which mostly contributes to the transport). It is found that the E × B shearing rate is larger than the linear mode frequency and growth rate of the SWITG mode, indicating that zonal flows are the primary mechanisms causing the SWITG mode to be suppressed. This is consistent with the findings of [72], which showed that E × B velocity shear can stabilise the slab SWITG mode. In our situation, it is also found that the E × B shearing rate is larger for the SWITG modes with R 0 /L T = 26.8, which also displays higher growth rate linearly, as contrasted to the lower shearing rate for those with R 0 /L T = 13.1 and lower growth rate linearly. In spite of the temperature scale length, which linearly determines the strength of the SWITG mode in comparison to the conventional ITG, the higher k θ ρ s part of the spectrum relevant to the SWITG mode contributes very little to the thermal ion heat flux. This is because the higher shearing rate appears to compensate for the higher growth rate of the mode.
When combined with other modes, such as TEMs, which in the presence of kinetic electrons at higher k θ ρ s , SWITG mode can be significant for electron thermal transport even though it has been shown to have very little contribution to net ion heat flux [10]. Therefore, when considering kinetic electron dynamics and interaction with TEMs, the smallness of the contribution of SWITG modes to heat transport determined here with an adiabatic electron model may differ. This issue will be addressed in a future communication. Apart from the TEM simulation, it will be interesting to take into account electromagnetic effects (finite β effects) on the SWITG mode in the linear and nonlinear regimes, which will also be addressed in future work. This effort is a first step towards understanding the microturbulence and transport in ADITYA-U. It might be possible to plan future ADITYA-U experiments using the knowledge gathered from this investigation on electrostatic microturbulence.
the gradient profiles can be obtained easily: where L T , L n are the temperature and density characteristic lengths, and ∆ T , ∆ n are constants. T e /T i = 1, q(s) = 0.86 − 0.16s + 2.52s 2 , ∆ T = 0.3, ∆ n = 0.3 and η = L n /L T = 3.12.
These input profiles are shown in figure A.1.
In the radial domain from 0.01 to 1.0 (in terms of normalized minor radius), the linear benchmark simulations for ORB5 are performed. In the toroidal direction, simulations are run as a full torus. Additional parameters are: a/ρ s = 185, l x = 2a/ρ s = 370.0, total number of grid points is 4 × 10 8 and 12 particles/cell are used in the linear simulation. The ITG instability is calculated at s 0 = 0.5, where the density and temperature gradients have maximum value. Here, we specify nine toroidal mode numbers: n = 5, 10,15,20,21,22,25,30,35 which cover the poloidal wave numbers k θ ρ s from 0.07 to 0.55. Each linear simulation corresponds to a single toroidal mode number n.
As illustrated in figure A.2, reasonable agreement is obtained for the real frequencies ω r and growth rates γ with GENE code [73] from where the real frequencies ω r and growth rates γ data are taken. Ad-hoc MHD equilibrium is used for this run. The poloidal mode structure of most unstable mode n = 21 is also shown in figure A.3. Note that for the plots in figure A.2, we run the code ORB5 only; the results corresponding to the other code (GENE) in the same figure are taken from [73]. When comparing growth rates obtained using ORB5 and GENE, a good agreement is observed for the lower k θ ρ s values, but minor inconsistencies are seen for k θ ρ s ⩾ 0.3. The growth rates for the two curves agree within 20%. It has been reported in [22] about the discrepancies due to the differences in the model equations considered by the two codes. The growth rates for the largest k θ ρ s are lower in ORB5, as compared to GENE. There is also a good agreement in real frequency values between ORB5 and GENE, that is within 10% with some small deviations, which are probably due to other minor differences between the two codes [22].

Appendix C. The Rosenbluth-Hinton test
The zonal flow damping tests for ADITYA-U are carried out in the linear regime. Low frequency electrostatic modes known as zonal flows, which are produced by turbulence on their own and are crucial for controlling the turbulence [58]. The Rosenbluth and Hinton theory on the collisionless damping of zonal flow [59] is verified by the collisionless damping of the zonal electrostatic field to a non-zero steady-state value.
As zonal flows are considered as a crucial saturation mechanism in turbulent regimes, particularly for ITG turbulence, a correct prediction of residual level (a non-zero steady-state value) is an important test for gyrokinetic codes [13,61,74].
( v th,i R 0 q ) √ 1 + 2 (23 + 16τ e + 4τ 2 e ) q 2 (7 + 4τ e ) 2 (C.3) with s as the minor radius of the considered magnetic surface i.e. 0.6 in the present simulation, R 0 is the major radius, q(s) is the safety factor on the surface of interest, τ e = T e /T i , and ) .
Input parameters are hydrogen ions, a = 274ρ s , R 0 = 0.75 m, R 0 /a = 3 and B 0 = 1.0 T. The temperature and density profiles are taken as flat. The safety factor profile is a monotonic increasing and given as q(s) = 1.25 + 0.67 s 2 + 2.38 s 3 − 0.06 s 4 , at s 0 = 0.6, q(s 0 ) = 2.0. The grid numbers are N s = 256, N θ⋆ = 64, N φ = 32, the number of markers is N = 1.0 × 10 7 and the time step size is ∆t = 10Ω −1 ci . The initial condition is prepared to obtain an axisymmetric ion density perturbation. The density perturbation δn i (s) ∼ δn i,0 sin(k r s), where k r is the radial wave number in units of (1/a). The results of the ORB5 simulations are shown in figure C.1. In figure C.1, the potential Φ 00 (s 0 , t) normalized to the initial value Φ 00 (s 0 , t 0 ) is displayed as a function of time. The exponential decay predicted by equation (C.1) and the residual calculated from equation (C.2) are also plotted for reference. The result indicates that the real frequency and residual amplitude are in agreement with the Rosenbluth-Hinton theory [59].