Abstract

This work examines the effect of electron pressure anisotropy on electromagnetic density fluctuations and the MHD stability parameter at the JET tokamak edge in the kinetic limit . Using a bi‑Maxwellian equilibrium and solving the linearized gyrokinetic Vlasov equation, we derive via the Fried‑Conte function. Accounting for density and magnetic field fluctuations yields. A critical anisotropy threshold is found. For (cases), the response is Landau‑damped () and ranges from 0.048% (weak turbulence,) to 4.370% (average turbulence,), staying below the 5% MHD ballooning threshold. The case gives the lowest fluctuations. For , becomes negative, signalling a kinetic mirror instability. In average turbulence, , exceeding the MHD threshold and indicating a pre‑ELM state. Comparison with Li‑BES measurements (Réfy et al. 2018) confirms moderate‑anisotropy cases, while extreme firehose () and mirror () match pre‑ELM triggering under NBI and ICRH heating, respectively.

1. Introduction

Despite considerable advances in the physics of confined plasmas, which have enabled modern tokamaks to come closer to fusion conditions ¹, a major limitation remains: the energy and particle transport measured experimentally often exceeds the predictions of classical or neoclassical models by several orders of magnitude, severely limiting confinement ². This anomalous transport is now attributed to small-scale electromagnetic turbulence, which is particularly active at the plasma edge, where strong pressure and current gradients sustain various microwave instabilities. In this peripheral region, electromagnetic fluctuations play a critical role in limiting the maximum value of the parameter, thereby compromising the achievement of high-performance fusion regimes.

The plasma edge region, comprising the pedestal and the scrape-off layer (SOL), is a particularly critical zone for the stability of tokamaks in H-mode, as it contains both the steepest pressure and temperature gradients and the highest levels of electromagnetic turbulence. The introduction of gyrokinetic theory in the 1980s led to a better understanding of the role of magnetic perturbations in electron dynamics at the edge, notably via the parallel streaming mechanism, which directly couples magnetic fluctuations to the electron distribution ³. At the same time, advances in sophisticated diagnostics such as lithium beam emission spectroscopy (Li-BES) on JET have enabled direct and quantitative measurement of electron density fluctuations at the plasma edge in H-mode, providing essential experimental data for the validation of kinetic models ⁴. This synergy between gyrokinetic theory, numerical simulations and experimental measurements has led to a better understanding of the role of electromagnetic instabilities, such as kinetic Alfvén waves (KAW), tearing modes and mirror and firehose instabilities, in the degradation of confinement at the tokamak edge ⁵.

In this article, we study the impact of electron density fluctuations induced by a fluctuating magnetic field on the stability parameter at the edge of the tokamak plasma, where the electromagnetic response becomes significant. Defined as the ratio of kinetic pressure to magnetic pressure, this parameter is crucial at the plasma edge: its fluctuations δβEM can indeed approach the thresholds of mirror () and firehose () instabilities during the pre-ELM phases ^{6, 7}. Within a rigorous kinetic formalism, the fluctuating electromagnetic electron density is calculated by integrating the perturbed bi-Maxwellian distribution function, obtained by linearizing the gyro-averaged Vlasov equation in the presence of the fluctuating vector potential, in the kinetic limit. The pressure anisotropy, induced by ICRH and NBI heating, plays a central role in modulating and, consequently, the fluctuations .

The paper is organised as follows. Section 2 discusses the effects of temperature anisotropy on the tokamak edge plasma. Section 3 presents the theoretical framework of the bi-Maxwellian and Maxwellian distribution functions, including the analytical derivation of the perturbed electromagnetic distribution function . Section 4 derives the analytical expression for the electromagnetic fluctuating density via the Fried-Conte function in the kinetic limit. Section 5 presents numerical applications at the edge of the tokamak for three cases of anisotropy and quantifies the impact on the parameter . Section 6 provides a comparative analysis with the experimental Li-BES measurements by Réfy et al. on JET. Section 7 concludes with a summary of the main results and future prospects.

2. The Effects of Thermal Anisotropy at the Edge of the Tokamak Plasma

The study of thermal anisotropy at the tokamak edge has highlighted its crucial role in transport and stability. Zhao et al. demonstrated using the UEDGE code that this anisotropy generates viscous forces that alter 2D transport and drift velocities in the boundary layer (SOL) ⁸. In terms of MHD stability, MISHKA-A simulations indicate that it strongly affects edge ballooning modes: the growth rate of instabilities increases for (typical of ICRH heating) and decreases for (typical of NBI neutral injection) ⁹. Theoretically, Ilgisonis et al. unified the description of the equilibrium and stability of anisotropic plasmas with macroscopic flows using the CGL approximation ¹⁰. The importance of rigorously modeling this anisotropy often induced by fast ions from auxiliary heating systems is confirmed by Z. Qu, who shows that a traditional isotropic model leads to major inconsistencies, particularly in the reconstruction of the safety profile q ¹¹. At the same time, Fielding and Haas laid the theoretical foundations for the dynamics of high-β tokamaks to define the stability limits under anisotropic pressure ¹². Building on this seminal work, the following sections will examine in detail the influence of these anisotropy effects on the fluctuating density by introducing a perturbed bi-Maxwellian distribution function.

3. Bi-MaxwellianDistributions as a Function of the Anisotropy Parameter A

In modern tokamaks, external heating methods (ICRH, NBI) generate distribution functions that deviate from isotropic equilibrium. Under the influence of a strong magnetic field, rapid cyclotron precession decouples particle dynamics, leading to the imposition of two distinct temperatures and ^{1, 2}. When collisions are insufficient to restore isotropy, which is common in the edge and pedestal regions, the plasma is described by a bi-Maxwellian distribution ¹³. This steady state, governed by the equilibrium between heating and wave-particle interactions, is characterised by the anisotropy parameter ³. In cylindrical coordinates (, this distribution is expressed as follows :

(1)

where is the thermal velocity, is the component parallel to, and is the perpendicular component. To obtain the bi-Maxwellian in modulus, we set and, the volume element becomes. After integration over the polar angle from , in the case where , we obtain the expression for the distribution of the modulus :

(2)

where is the imaginary error function ⁵. In the case where, the expression becomes:

(3)

where is the usual error function. The error function () physically represents the mixing of the two temperatures. It adjusts the shape of the curve to account for the fact that, for the same modulus , a particle can be either slow relative to or fast relative to . The increase in leads to an increase in the population of particles capable of giving energy to the wave (resonance), which explains why the instability thresholds often decrease when the perpendicular anisotropy is strong.

Figure 1(a) shows that ICRH heating injects energy into the gyration motion, which broadens the distribution at high velocities and lowers the peak amplitude to preserve the total number of particles. From a physical standpoint, this population of suprathermal particles perpendicular to the magnetic field promotes mirror instability and alters the local magnetic pressure, which directly impacts the stability of the β parameter. Figure 1(b) indicates that NBI heating, conversely, causes the distribution to narrow, characterized by low dispersion and a more intense peak shifted toward lower velocities. This concentration of energy in the parallel direction reduces stability against firehose instability while profoundly altering the plasma’s dielectric response and Landau damping.

4. Perturbation of the Bi-Maxwellian Distribution Function

At the edge of the H-mode tokamak plasma, unlike in the core where the electrostatic response dominates, electromagnetic fluctuations governed by the parallel vector potential play a predominant role in electron dynamics and confinement stability.The rotational approximation allows us to express the magnetic field as ), and since the turbulent perturbations at the edge of a tokamak vary much more abruptly in the direction perpendicular to the field lines than in the direction parallel to them ( or ), We find that ¹⁴. In this context, the anisotropic bi-Maxwellian distribution, constitutes the natural starting point for deriving the electromagnetic perturbation distribution function .We decompose the total distribution function into a bi-Maxwellian equilibrium part and a small perturbation :

(4)

In the presence of the fluctuating vector potential, the linearized gyrokineticVlasov equation for electrons is written:

(5)

At the electromagnetic limit ( dominant,), the parallel driving force is given by, the solution obtained by integration along the undisturbed trajectories gives:

(6)

Either

(7)

By introducing the anisotropic diamagnetic frequencies and , the complete perturbation including gradient corrections is written :

(8)

where is the density diamagnetic frequency, , the wave components, the thermal drift frequency related to the temperature gradient, and represents an attenuation factor handling singularities in plasma response calculations. The frequency of thermal drift with spatial gradient is defined by:

(9)

where represents the intensity of the temperature gradient relative to the density gradient ¹⁶. For a bimaxwellian distribution, the normal pressure modifies this frequency ¹⁵.Since the diamagnetic drift arises from the gyration motion of the particles, it is the perpendicular pressure that drives the diamagnetic flux; the term then becomes a function of:

(10)

where is the length scale of the density gradient. It is related to the speed at which the electron guidance centers move due to the macroscopic density gradient. This perturbation then becomes crucial for the study of electronic temperature gradient (ETG) modes; the expression for becomes:

(11)

This is an electromagnetic gyrokinetic response, within the limit . Within this electromagnetic formalism, we adopt the parallel component of the vector potential, denoted , as the fundamental variable describing the magnetic fluctuations; this accounts for the twisting of the field lines and induces, through coupling with the parallel velocity of the electrons.

Figure 2(a) shows that temperature anisotropy structurally alters the equilibrium distribution. For, the profile widens and shifts toward higher velocities, whereas for, it narrows around a more intense peak; these shifts directly shape the morphology of the electromagnetic perturbation around the resonant phase velocity. Analysis of Figure 2(b) reveals the impact of these changes on the perturbation. Its real part reflects an energy redistribution whose polarity reverses depending on the sign of the anisotropy, leading to a decrease in density for and the appearance of a positive peak for. At the same time, the imaginary part shows that parallel anisotropy () amplifies Landau damping through a more pronounced trough, whereas strong perpendicular anisotropy can reverse the sign of this component, marking the transition from a dissipative regime to unstable wave growth. Ultimately, deviating from isotropy not only changes the shape of the particle reservoir but also fundamentally redefines the energy exchange channels between the turbulence and the plasma at the tokamak’s periphery.

5. Analytical Expression of Fluctuating Electron Density

The electromagnetic perturbed density is obtained by integrating the perturbed distribution function over the entire velocity space, which naturally brings the Fried-Conte plasma dispersion function into being a central tool for dealing with the Landau resonance.

(12)

To describe a temperature anisotropy, it is natural to use cylindrical velocity coordinates . The integration is written as:

(13)

The equilibrium bi-Maxwellian distribution function is normalized with respect to such that the direct integration of its perpendicular moments gives: . By integrating the numerator with respect to , the term in brackets becomes:

(14)

The reduced integrated numerator is then written as:

(15)

If we define the dimensionless variable, then the numerator term n equation (15) becomes. By rearranging the terms in terms of the powers of, we obtain:

(16)

The numerator can therefore be written simply in the following form:. Introducing the normalized phase velocity, the normalized density then becomes ¹⁷:

(17)

Let us consider the integral part of equation (16), denoted by :

(18)

By factoring out the term from both terms, and:

(19)

To evaluate, we rewrite the numerator as follows . Then, by dividing it by, we get . The integral breaks down into a normalized Gaussian term and a term proportional to the plasma dispersion function.

(20)

The first term on the right-hand side is the normalized Gaussian integral, equal to unity. The second term defines the plasma dispersion function, , given by :

(21)

analytically extrapolated along the real axis 𝜁 in accordance with Landau's causality condition. If we incorporate these two results into the calculation of, we immediately arrive at the following expression:

(22)

The same procedure applies to ... Using the polynomial identity, we obtain, so that

(23)

The first integral is evaluated term by term using the standard Gaussian moments: and the latter term cancels out due to parity, as does the normalized Gaussian integral used above. The second integral is again equal to by definition. Combining these terms yields:

(24)

After combining the terms, the complete integral is written as:

(25)

In the kinetic limit, the Fried-Conte plasma dispersion function expands in first order as . From this, we obtain the following approximations and . By substituting these limits into the integral expression:

(26)

It should be noted that, although fluctuations in the electrostatic potential 𝜙 are often treated explicitly in turbulence models, our formalism incorporates them here through the plasma’s kinetic response function. Indeed, in the limit of electromagnetic drift waves at the JET edge, the density fluctuation is coupled to the potential via the quasi-neutrality condition. The approach adopted here, by focusing on the direct response to the perturbation of the bi-Maxwellian equilibrium, captures the self-coherent response of the plasma. The obtained density amplitude levels, calibrated using experimental data from JET (Ref. ¹⁸, take this coupling dynamics into account. Consequently, ignoring 𝜙 as an independent source term is not an oversight, but rather a closure of the kinetic system in which the potential is considered a component of the medium’s response on the same footing as the density.

The critical threshold is defined by the exact cancellation of the real part,; thus, we have:

(27)

Thus, we find the following:

(28)

This writing highlights a structure with two very distinct contributions, which the separation into real and imaginary parts allows us to clearly identify:

(29)

(30)

where is the amplitude constant of the turbulent regime. The graph below illustrates the behavior of these two components as well as the modulus of as a function of , in mean turbulence regime (.

Analysis of the real part (Figure 3a) highlights the in-phase response of the electrons, characterized by the classical transition between thermal screening and particle inertia at the phase velocity . Strong perpendicular anisotropy) disrupts this regularity by inducing strongly oscillatory behavior. For its part, the imaginary part (Figure 3b) quantifies irreversible energy exchanges: the isotropic, parallel, and moderately perpendicular configurations () reflect a dissipative regime dominated by Landau damping (). Conversely, for, the curve becomes positive () at low phase velocities (), indicating a net transfer of energy from the particles to the wave direct kinetic evidence of the onset of the mirror instability. The modulus (Fig.3c) confirms the accuracy of the model: at the resonance, the curves for realistic tokamak configurations (to) fall strictly within the experimental band of the JET Li-BES diagnostic ¹⁸, with only the limiting case standing out due to a non-monotonic profile.Finally, all calculated amplitudes remain at least one order of magnitude below the critical threshold for ideal MHD (). This demonstrates the significant value of this kinetic approach, which is capable of capturing local transitions and instabilities well before the plasma reaches its limits of macroscopic destruction ¹⁹.

6. Applications to the JET Tokamak

The tables below summarize all the numerical values of the characteristic electronic parameters estimated at the edge of the JET tokamak plasma. The parameters at the edge are representative of a plasma in H mode with an electronic temperature and a density ¹⁴.As reported by Horton, W. (1999), the normalized phase velocity is and a magnetic turbulence amplitude , which fixes the order of magnitude of ^{20, 21, 22}.

Table 2 presents the numerical values of the ratio calculated at the plasma boundary, a region where magnetic fluctuation, characterized by the vector potential, constitutes the dominant perturbation. In this region, strong pedestal pressure gradients and high levels of MHD activity generate significant magnetic fluctuations that modulate the electron density through parallel streaming and electromagnetic Landau resonance.

When A exceeds the critical threshold , the real part becomes negative, indicating a phase shift of π between and , which is characteristic of the kinetically unstable mirror regime. For all cases where , the real part is positive.

The parameter represents the ratio of kinetic pressure to magnetic pressure. In a plasma where fluctuations are primarily driven by electrons (neglecting the ionic contribution at the boundary), we define:

(31)

Considering as a constant, the fluctuation induced by the electromagnetic perturbation is given by:

(32)

Purely perpendicular magnetic fluctuations () affect the total magnetic field strength only to second order. In the first order, the change in the magnetic field strength is exactly zero.This gives us the following expression:

(33)

From the values of the perturbed electron density calculated previously, it is possible to estimate the fluctuation of the stability parameter induced by turbulent perturbations in the plasma of the JET tokamak see table below.

We observe that for , and , the real part is positive. Conversely, for , the relative correction is the largest (5.78%) because is small; the magnetic correction therefore reduces the dominant term to a greater extent.

Figure 4 shows that the dynamics of the pedestal are governed by a kinetic transition at . Below this value, resonance dominates the effective behavior, whereas above it, the response becomes more rigid. The imaginary part reveals that the anisotropy drives the energy transfer, whereas Landau damping is dominant for . Beyond the resonance, the energy contribution from anisotropy becomes crucial, potentially leading to plasma instability. The analysis of the amplitudes confirms that, in the pedestal region , perpendicular anisotropy significantly reduces the MHD stability margin. This drives the plasma toward the nonlinearity point, thereby promoting the activation of ballooning modes and confirming the crucial role of anisotropy in the deterioration of confinement.

7. Comparison with Experimental Results

To assess the physical validity of our kinetic model, the calculated fluctuations in the electromagnetic electron density and those of the beta parameter are compared with the Li-BES measurements by Réfy et al., ¹⁸ at the JET periphery and with the MHD stability criterion by Hender et al., ¹⁹.The reference conditions are, (average turbulence).In this framework, the wave frequency is treated parametrically, fixed to typical experimental and simulation values characteristic of dominant core modes . Consequently, the derived threshold should be interpreted as the local kinetic drive inversion threshold for a pre-existing turbulent carrier wave, rather than an autonomous global mode dispersion threshold. This parametric scan allows us to isolate how the kinetic dissipation channel alone responds to equilibrium anisotropy variations.

For, the model gives and , well below the 5% MHD threshold. The modulus exceeds the upper bound of the Li-BES inter-ELM band , as expected: at the pedestal, electromagnetic fluctuations driven by are inherently larger than their electrostatic counterparts. The imaginary-to-real ratio confirms that the response is predominantly adiabatic, with Landau damping contributing less than 3.2% of the total amplitude consistent with the kinetic limit.When , the anisotropic diamagnetic frequency rad/s approaches rad/s, reducing the driving term to rad/s. The fluctuation level drops to reduced by a factor of relative to, and falling within the Li-BES inter-ELM band ¹⁸. The complete delta-beta gives, the lowest value among all five cases. This demonstrates that moderate ICRH-driven perpendicular anisotropy acts as a stabilising mechanism for electromagnetic edge fluctuations, suppressing the diamagnetic driving term through the near-resonance condition.

For, gives, and stable, but 45% above the isotropic value. This moderate enhancement is consistent with experimental observations in NBI-heated discharges, where mild parallel anisotropy () amplifies edge electromagnetic turbulence without triggering ELM events.The case yieldsand, the largest among the linearly stable configurations and critically close to the 5% MHD threshold (87% of the critical value). This near-threshold behaviour aligns with the observations of Solano et al. ²³, who reported enhanced magnetic fluctuations in JET H-mode phases with intense NBIinjectionfavouring. According to our model, a turbulence amplitude increase fromto merely 0.01144 would bring to the 5% boundary, providing a quantitative estimate of the ELM trigger margin under strong parallel anisotropy.The most physically significant result arises for, representative of intense ICRH heating ²⁴. Here exceeds , inverting the driving term to. As a consequence, a sign reversal signalling the transition from a Landau-damped response to kinetically unstable mirror-mode growth. This inversion occurs at the critical threshold. With the complete formula, in the average turbulence regime more than twice the MHD critical threshold placing the plasma in a pre-ELM instability state consistent with ICRH-triggered ELM events reported by Hender et al. ¹⁹. In the weak turbulence limit (),remains stable, confirming the turbulence-amplitude dependence of the instability onset.

The anisotropy-driven amplitude ratios equal to 0.254 (), 1.449 (), 2.199 (), and 5.000 () are strictly independent of, as they depend solely on the ratio. This invariance confirms that the anisotropy effects identified here are an intrinsic property of the bi-Maxwellian diamagnetic structure, not an artefact of the chosen wave number.

8. Conclusion

This study quantifies the impact of electron pressure anisotropy on the electromagnetic stability of the JET tokamak edge using a linear gyrokinetic model governed by an intrinsic critical threshold (). In stabilizing regimes (), the isotropic case () and the moderate mirror case (.5) effectively reduce plasma fluctuations to a minimum of due to a resonance effect. Conversely, firehose regimes () widen the frequency spread and amplify these fluctuations (which rise from to ), bringing the system dangerously close to the MHD threshold.

The extreme mirror regime (), representative of intense localized ICRH heating, crosses this critical threshold and triggers a phase inversion. The plasma then transitions to a kinetically unstable mirror mode growth, driving the amplitude of the fluctuations to twice the MHD threshold which mechanically explains the ELM triggers observed during the experiment. These theoretical predictions show excellent quantitative agreement with the JET Li-BES experimental data: the moderate cases align perfectly with the inter-ELM transition phases, while the extreme cases reproduce the pre-ELM behavior, with the remaining discrepancy factor explained by the action of nonlinear effects.

To overcome the current limitations of this study (related to the linear approximation, the slab geometry, and frozen ion dynamics), future work will focus on developing a nonlinear gyrokinetic framework. This extension will allow for the inclusion of turbulent saturation as well as full electromagnetic coupling. Finally, this new model will be rigorously compared and validated through dedicated benchmarking against the reference simulation codes GENE and GKW, using the precise parameters of the JET pedestal.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The data that support the findings of this study are available from the corresponding author, Said Mouigni Abdallah, upon reasonable request.

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors. This work was conducted within the framework of doctoral research at the Laboratory of Plasma Physics and Interdisciplinary Research (LPPRI), Cheikh Anta Diop University of Dakar (UCAD), Senegal.

References