1 Introduction
The prospect of creating electron–positron pair plasmas magnetically confined in dipole or stellarator geometries has been discussed since the early 2000s (Pedersen et al. Reference Pedersen, Boozer, Dorland, Kremer and Schmitt2003). In the near future, the first experiment aiming at this goal will be constructed (Pedersen et al. Reference Pedersen, Danielson, Hugenschmidt, Marx, Sarasola, Schauer, Schweikhard, Surko and Winkler2012). Recently, efficient injection and trapping of a cold positron beam in a dipole magnetic field configuration has been demonstrated by Saitoh et al. (Reference Saitoh, Stanja, Stenson, Hergenhahn, Niemann, Pedersen, Stoneking, Piochacz and Hugenschmidt2015) using a supported permanent magnet. This result is a key step towards further studies using the levitated magnetic coil with the ultimate aim of creating and studying the first manmade magnetically confined pair plasma in the laboratory.
It has been shown by Helander (Reference Helander2014) that pair plasmas possess unique gyrokinetic stability properties due to the mass symmetry between the particle species. For example, drift instabilities are completely absent in a straight geometry, e.g. in a slab, provided that the density and temperature profiles of the two species are identical (‘symmetric’ pair plasmas). The symmetry between the two species is broken if the temperature profiles of the electrons and positrons differ or there is an ion contamination. In these regimes, drift instabilities can be excited even in an unsheared slab geometry (Mishchenko et al. Reference Mishchenko, Zocco, Helander and Könies2018b ). In a sheared slab, pure pair plasmas are prone to current-driven reconnecting instabilities (Zocco Reference Zocco2017), but there are no drift waves. Note that asymmetry between the species is needed also in this case since the ambient electron flow velocity must differ from the positron one for the ambient current to be finite. In contrast to the slab geometry, a dipole magnetic field has finite curvature. In this case, the symmetry between the species is broken by curvature drifts and the plasma can be driven unstable by temperature and density gradients (Helander Reference Helander2014), even without ion contamination and for identical temperature profiles of the two species. This result also persists in the electromagnetic regime (Helander & Connor Reference Helander and Connor2016). The nonlinear stability of dipole pair plasmas has also been addressed by Helander (Reference Helander2017). More recently, Mishchenko, Plunk & Helander (Reference Mishchenko, Plunk and Helander2018a ) performed a detailed study of the gyrokinetic stability of a pure pair plasma in the dipole geometry, making use of both the Z-pinch and point-dipole limits. Again, it was found that such pair plasmas can be driven unstable by magnetic curvature, density and temperature gradients.
In this paper, we use the gyrokinetic code GENE (Jenko et al.
Reference Jenko, Dorland, Kotschenreuther and Rogers2000) to study the linear stability of plasmas with artificially varied mass ratios ranging from a conventional deuterium plasma to electron–positron plasmas in a variety of different magnetic geometries. In § 2 we introduce the GENE code and the assumptions and modifications required to run simulations for electron–positron plasmas. In § 3 we show the set-up and results for simulations in two different magnetic geometries. We employ the
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
model for an axisymmetric tokamak configuration and compare these with Pedersen et al. (Reference Pedersen, Boozer, Dorland, Kremer and Schmitt2003). We also run simulations for the W7-X stellarator geometry. We present the analytic theory of trapped-particle modes in electron–positron plasmas to explain the results of the numerical simulations. In § 4 we present our conclusions and outlook.
2 Physical assumptions
In this work GENE is operated in a linear mode, employing a local (‘flux-tube’) geometry. Throughout this paper we have considered only collisionless, electrostatic simulations, we have also assumed that the plasma density of species
$a$
denoted by
$n_{a}$
is sufficiently large when compared to the Brillouin density
$$\begin{eqnarray}n_{a}\gg n_{Ba}:=\frac{\unicode[STIX]{x1D716}_{0}m_{a}\unicode[STIX]{x1D6FA}_{ca}^{2}}{2\text{e}_{a}^{2}}.\end{eqnarray}$$
Here
$\unicode[STIX]{x1D716}_{0}$
is the vacuum permittivity,
$m_{a}$
the mass of the plasma species,
$\unicode[STIX]{x1D6FA}_{ca}=e_{a}B/m_{a}$
the gyrofrequency,
$B$
the magnitude of the magnetic field and
$e_{a}$
the charge of each plasma species. This assumption ensures that that Debye shielding can be neglected provided
$k_{\bot }\unicode[STIX]{x1D70C}\lesssim 1.$
Henceforth, we are working in a limit where
$k_{\bot }\unicode[STIX]{x1D706}_{D}\ll 1.$
Here, the following notation is employed
$$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D706}_{Da}^{2}}=\frac{\text{e}_{a}^{2}n_{a}}{\unicode[STIX]{x1D716}_{0}T_{a}},\quad \frac{1}{\unicode[STIX]{x1D706}_{D}^{2}}=\mathop{\sum }_{a}\frac{1}{\unicode[STIX]{x1D706}_{Da}^{2}},\end{eqnarray}$$
where
$T_{a}$
is the temperature of each species. It is important to note that this limit may not be satisfied in practice and this area warrants further study. Throughout this work,
$\unicode[STIX]{x1D6FD}_{a}=\unicode[STIX]{x1D707}_{0}n_{a}T_{a}/B^{2}$
. The electrostatic limit corresponds, as usual, to
$\unicode[STIX]{x1D6FD}_{a}=0$
.
The GENE code can handle plasmas of different mass ratios by changing the mass of the lighter species (electrons). An electron–positron plasma is obtained when the mass of the electrons is equal to the mass of the singly charged ion species – thus obtaining charge asymmetry and a mass ratio of unity.
The numerical resolution in the field-line-following direction is 128 grid points, while 16 and 8 grid points are employed for the parallel velocity and magnetic moment directions respectively. GENE employs numerical dissipation (hyperdiffusion). In this work, the strength of the hyperdiffusion is set to 0.2 for the parallel velocity and to 0.25 in the parallel direction.
3 Results
In toroidal fusion plasmas, transport is generally dominated by turbulence which is driven by microinstabilities. We have begun to investigate the microstability of quasineutral electron–positron plasmas in different magnetic configurations: a simple tokamak geometry and a stellarator geometry. Although drift waves and acoustic waves are not easily found in a pair plasma, it has been shown that the basic curvature-driven interchange instability can threaten the confinement of such plasmas in experimental devices. In the absence of magnetic shear, magnetic field curvature induces an interchange instability with a growth rate which scales like
$\unicode[STIX]{x1D6FE}\sim v_{t}/\sqrt{RL}$
, where
$v_{t}$
is the thermal velocity,
$R$
is the radius of curvature and
$L$
is the characteristic pressure gradient length scale of the plasma. In the absence of other microinstabilities, the interchange is thus a strong candidate for driving turbulence in a laboratory pair plasma.
3.1 Tokamak geometry
A simple test case for a magnetic geometry which is unstable to the collisionless interchange, is a circular tokamak with aspect ratio
$a/R$
. To this end we employ the well-known
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
geometry (Connor, Hastie & Taylor Reference Connor, Hastie and Taylor1971) model for which many studies have already been performed e.g. Dimits et al. (Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). This model geometry is parametrized by the magnetic shear
${\hat{s}}$
and the variable
$\unicode[STIX]{x1D6FC}$
introduced in the next section. Numerical simulations with varying plasma ratios, from deuterium plasmas to electron–positron plasmas, have also been addressed by Pedersen et al. (Reference Pedersen, Boozer, Dorland, Kremer and Schmitt2003) using the GS2 gyrokinetic code, this provides an excellent starting case for the work presented here and gives a basis for the comparison of different gyrokinetic codes.
The magnetic flux tube used for the simulations originates at radial position
$r/R=0.36.$
The normalised density and temperature gradients read
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{n,a}=-\frac{R}{n_{a}}\frac{\text{d}n_{a}}{\text{d}r},\quad \unicode[STIX]{x1D714}_{T,a}=-\frac{R}{T_{a}}\frac{\text{d}T_{a}}{\text{d}r},\end{eqnarray}$$
where the subscript
$a$
is the species index.
In this case, as in the aforementioned work of Pedersen et al. (Reference Pedersen, Boozer, Dorland, Kremer and Schmitt2003), we take the density gradient
$R/L_{n}=3$
, the temperature gradient
$R/L_{T}=3$
, the safety factor
$q=1.4$
, the magnetic shear
${\hat{s}}=0.8$
and equal temperatures
$T_{i}=T_{e}$
.
Growth rate
$\unicode[STIX]{x1D6FE}$
and corresponding real frequency
$\unicode[STIX]{x1D714}_{r}$
plotted as a function of the mass ratio. Simulations shown for an
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
geometry.
In figure 1 the growth rate
$\unicode[STIX]{x1D6FE}$
(normalised by
$c_{i}/a$
, where
$c_{i}=\sqrt{T_{i}/m_{i}}$
is the ion sound speed) and the corresponding real frequency
$\unicode[STIX]{x1D714}_{r}$
are plotted against the ratio of the electron mass to the ion mass. The sign of
$\unicode[STIX]{x1D714}_{r}$
is such that
$\unicode[STIX]{x1D714}_{r}>0$
corresponds to propagation in the ion/positron diamagnetic direction. The growth rate of the fastest growing mode for different values
$k_{y}\unicode[STIX]{x1D70C}_{i}$
is shown, where
$k_{y}$
is the binormal wavenumber and
$\unicode[STIX]{x1D70C}_{i}$
is the reference gyroradius,
$\unicode[STIX]{x1D70C}_{i}=c_{i}/\unicode[STIX]{x1D6FA}_{i}$
. Standard plasma mass ratios are shown on the left side of the figure, with the standard electron–hydrogen plasma occupying the left-most point. The pair plasma is at unity on the horizontal axis, where
$m_{e}=m_{i}$
.
When the mass ratio is varied to approach unity, we see an increase in the growth rate of the most unstable mode – the associated change in the sign of
$\unicode[STIX]{x1D714}_{r}$
corresponds to a regime switch from modes propagating in the electron diamagnetic direction, to modes propagating in the ion/positron diamagnetic direction.
In figures 2–5 we see how the modulus of the electrostatic potential
$\unicode[STIX]{x1D6F7}=|\unicode[STIX]{x1D719}|$
, the curvature
$\unicode[STIX]{x1D705}$
and the magnetic field strength
$B$
vary along the flux tube.
The
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
geometry.
$\unicode[STIX]{x1D6F7}$
,
$\unicode[STIX]{x1D705}$
versus
$\unicode[STIX]{x1D703}$
,
$m_{e}/m_{i}=1/1836$
.
The
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
geometry.
$\unicode[STIX]{x1D6F7}$
,
$B$
versus
$\unicode[STIX]{x1D703}$
,
$m_{e}/m_{i}=1/1836$
.
It is not possible to easily classify these modes as so-called temperature-gradient-driven modes and trapped-particle modes from the mode structure alone. Indeed, in this simple model tokamak both the magnitude of the magnetic field
$B$
and the curvature
$\unicode[STIX]{x1D705}$
have their minimum at the outboard midplane. Negative values of
$\unicode[STIX]{x1D705}$
correspond to ‘bad’ curvature. It is also important to note that both species are treated kinetically in the code and therefore particle trapping effects are present at all mass ratio values.
The
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
geometry.
$\unicode[STIX]{x1D6F7}$
,
$\unicode[STIX]{x1D705}$
versus
$\unicode[STIX]{x1D703}$
,
$m_{e}/m_{i}=1.0$
.
The growth rate is shown for different mass ratios in figure 6 as a function of
$k_{y}\unicode[STIX]{x1D70C}_{i}$
.
It is important to recall that we are considering plasmas with a sufficiently high density to neglect Debye shielding effects. For low enough density, the positron–electron interchange instability can be stabilised by the finite Debye length. The effects of Debye length have recently been investigated by Horn-Stanja et al. (Reference Horn-Stanja, Biancalani, Bottino, Mishchenko and Sunn Pedersen2018).
The
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
geometry.
$\unicode[STIX]{x1D6F7}$
,
$B$
versus
$\unicode[STIX]{x1D703}$
,
$m_{e}/m_{i}=1.0$
.
Growth rate
$\unicode[STIX]{x1D6FE}$
plotted as a function of
$k_{y}\unicode[STIX]{x1D70C}_{i}$
for different mass ratios. Simulations shown for an
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
geometry.
We note the excellent agreement of our results using the GENE gyrokinetic code operating in a linear flux-tube mode shown in figure 1 with those obtained by Pedersen et al. (Reference Pedersen, Boozer, Dorland, Kremer and Schmitt2003) using the GS2 gyrokinetic code (their figure 4). It is important to note that our growth rate is normalised using the ion sound speed,
$c_{i}$
, whereas the results of Pedersen et al. (Reference Pedersen, Boozer, Dorland, Kremer and Schmitt2003) are normalised to the ion thermal speed
$v_{ti}=\sqrt{2T_{i}/m_{i}}=\sqrt{2}c_{i},$
this accounts for the apparent discrepancy in the
$y$
-axis.
3.2 Wendelstein 7-X geometry
We further extend the results of Pedersen et al. (Reference Pedersen, Danielson, Hugenschmidt, Marx, Sarasola, Schauer, Schweikhard, Surko and Winkler2012) by also using the GENE code to study stability in a non-axisymmetric configuration. The stellarator has been proposed as a suitable trap for such a plasma. In this work we use the Wendelstein 7-X geometry to investigate plasma stability in a stellarator when the mass ratio is varied. Whilst this is unlikely to be experimentally relevant to the ongoing W7-X campaigns, the study of electron–positron plasmas in stellarators is nevertheless relevant since this type of device is one of the main candidates for a pair plasma trap.
For the stellarator geometry considered we will find it useful to employ the Boozer coordinate system
$(\sqrt{s},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D703})$
where
$\unicode[STIX]{x1D6FC}=\sqrt{s_{0}}/q_{0}(q\unicode[STIX]{x1D703}-\unicode[STIX]{x1D701})$
. Here,
$\unicode[STIX]{x1D703},\unicode[STIX]{x1D701}$
are the Boozer poloidal and toroidal angles, respectively,
$s_{0}$
is the radial location of the magnetic surface in terms of the normalised toroidal flux and
$q$
is the safety factor. The magnetic shear is defined as
$$\begin{eqnarray}{\hat{s}}=\frac{\text{d}\ln q}{\text{d}\ln r},\end{eqnarray}$$
where
$r$
is the minor radius of the flux surface.
We take the density gradient
$a/L_{n}=3$
and the temperature gradient
$a/L_{T}=3$
with
$$\begin{eqnarray}L_{n}=-\frac{1}{n}\frac{\text{d}n}{\text{d}x},\quad L_{T}=-\frac{1}{T}\frac{\text{d}T}{\text{d}x},\end{eqnarray}$$
where
$x=a\sqrt{s}$
. Our simulations are performed in a magnetic flux tube at the normalised radial position of
$r/a=0.49$
. Geometrically, this is the flux tube which originates at the centre of the bean-shaped cross-section in the device. This surface has a safety factor value of
$q=1.4$
, and a magnetic shear of
${\hat{s}}=-0.05$
.
Growth rate
$\unicode[STIX]{x1D6FE}$
and corresponding real frequency
$\unicode[STIX]{x1D714}_{r}$
plotted as a function of the mass ratio. Simulations shown for the Wendelstein 7-X geometry in the flux tube originating in the bean-shaped cross-section.
In figure 7 the normalised growth rate
$\unicode[STIX]{x1D6FE}$
(maximised over
$k_{y}\unicode[STIX]{x1D70C}_{i}$
) and the corresponding real frequency
$\unicode[STIX]{x1D714}_{r}$
are plotted against the ratio of the electron mass to the ion mass. (The peak value tends to be in the interval
$0.2<k_{y}\unicode[STIX]{x1D70C}_{i}<0.6$
.) Conventional plasma mass ratios correspond to the left side of the figure, with the standard electron–hydrogen plasma occupying the left-most point. The pair plasma is at unity on the
$x$
-axis, where
$m_{e}=m_{i}$
.
Growth rate
$\unicode[STIX]{x1D6FE}$
plotted as a function of
$k_{y}\unicode[STIX]{x1D70C}_{i}$
for different mass ratios. Simulations shown for the W7-X geometry in the flux tube originating in the bean-shaped cross-section.
As before, we see an increase in the growth rate of the instability as the mass ratio approaches unity. We can also once again identify two distinct regimes which are dominated by modes propagating in the electron diamagnetic direction
$(\unicode[STIX]{x1D714}_{r}<0)$
and modes propagating in the ion/positron diamagnetic direction
$(\unicode[STIX]{x1D714}_{r}>0)$
.
We note the interesting result that the real frequency is exactly zero for pair plasmas. In § 4 we show that this result can be derived analytically.
In figure 8 we see how the growth rate varies as a function of
$k_{y}\unicode[STIX]{x1D70C}_{i}$
. It is interesting to compare these results with those shown in figure 6 for the
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
geometry. We note in passing that in the stellarator, plasmas with a mass ratio close to unity tend to have their growth rates maximised for larger values of
$k_{y}\unicode[STIX]{x1D70C}_{i}$
than in the
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
geometry, this is due to the smaller connection length.
In the stellarator geometry considered, the regime change between the mode propagation direction appears to occur for a larger mass ratio than in the tokamak case. Indeed, the threshold mass ratio is an order of magnitude higher in the W7-X geometry than in the
${\hat{s}}{-}\unicode[STIX]{x1D6FC}$
tokamak case.
The W7-X geometry, bean flux tube
$\unicode[STIX]{x1D6F7}$
,
$\unicode[STIX]{x1D705}$
versus
$\unicode[STIX]{x1D703}$
,
$m_{e}/m_{i}=1/1836$
.
The W7-X geometry, bean flux tube.
$\unicode[STIX]{x1D6F7}$
,
$B$
versus
$\unicode[STIX]{x1D703}$
,
$m_{e}/m_{i}=1/1836$
.
The W7-X geometry, bean flux tube.
$\unicode[STIX]{x1D6F7}$
,
$\unicode[STIX]{x1D705}$
versus
$\unicode[STIX]{x1D703}$
,
$m_{e}/m_{i}=1.0$
.
The W7-X geometry, bean flux tube.
$\unicode[STIX]{x1D6F7}$
,
$B$
versus
$\unicode[STIX]{x1D703}$
,
$m_{e}/m_{i}=1.0$
.
In figures 9 and 10 we see how the modulus of the electrostatic potential
$\unicode[STIX]{x1D6F7}=|\unicode[STIX]{x1D719}|$
, the curvature
$\unicode[STIX]{x1D705}$
and the magnetic field strength
$B$
vary along the flux tube. We can also examine the mode structure for electron–positron plasmas (figures 11 and 12). In contrast to the tokamak case (figures 2–5), the mode structure is sensitive to the mass ratio
$m_{e}/m_{i}$
. Presumably, this is because the regions of bad curvature
$(\unicode[STIX]{x1D705}<0)$
and particle trapping (low
$B$
) do not coincide in the stellarator in contrast to the tokamak. In § 4 we show analytically that particle trapping plays a large role in the stability of pair plasmas.
We have also performed simulations using a variety of different temperature gradients and density gradients and using a different flux tube (not shown). In each case we observed similar trends in both the behaviour of the growth rate of the fastest growing mode and the structure of the modes along the flux tube for different mass ratios.
4 Analytic theory of trapped-particle modes for electron–positron plasmas
We note the interesting result that
$\unicode[STIX]{x1D714}_{r}\rightarrow 0$
as
$m_{e}\rightarrow m_{i}$
with
$\unicode[STIX]{x1D714}_{r}=0$
at a mass ratio of unity. This can be accounted for from the analytic theory of trapped-particle modes. The theory of trapped-electron modes in conventional plasmas has been studied extensively in Helander, Proll & Plunk (Reference Helander, Proll and Plunk2013). Here we follow a similar calculation but for plasmas of equal mass, which we will see plays an important role in the mode structure.
We consider an arbitrary configuration with nested magnetic flux surfaces, so that the magnetic field can be written as
$\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}$
, with
$\unicode[STIX]{x1D713}$
the toroidal magnetic flux and
$\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D703}-\unicode[STIX]{x1D704}\unicode[STIX]{x1D719}$
the Clebsch angle, constructed as usual from the poloidal and toroidal magnetic coordinates. In a gyroradius expansion, it follows from the zeroth-order drift kinetic equation that the equilibrium distribution function,
$f_{a0}$
, of each species
$a$
is a Maxwellian at rest whose density and temperature are constant on each flux surface. The linear stability of this equilibrium against drift-wave-ordered instabilities is governed by the gyrokinetic equation
$$\begin{eqnarray}\text{i}v_{\Vert }\unicode[STIX]{x1D735}_{\Vert }g_{a}+(\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{da})g_{a}=\frac{e_{a}\unicode[STIX]{x1D719}}{T_{a}}\text{J}_{0}\left(\frac{k_{\bot }v_{\bot }}{\unicode[STIX]{x1D6FA}_{a}}\right)(\unicode[STIX]{x1D714}-\unicode[STIX]{x1D714}_{\star a}^{\text{T}})f_{a0},\end{eqnarray}$$
in the collisionless and electrostatic approximation. Here
$\unicode[STIX]{x1D719}$
is the electrostatic potential perturbation,
$\text{J}_{0}$
is a zeroth-order Bessel function,
$\unicode[STIX]{x1D6FA}_{a}=e_{a}B/m_{a}$
the gyrofrequency and
$$\begin{eqnarray}g_{a}(\boldsymbol{R},v,\unicode[STIX]{x1D706},t)=f_{a1}(\boldsymbol{r},\boldsymbol{v},t)+\frac{e_{a}\unicode[STIX]{x1D719}(\boldsymbol{r},t)}{T_{a}}f_{a0}(v)\end{eqnarray}$$
denotes the non-adiabatic part of the perturbed distribution function, which in lowest order becomes independent of the gyroangle when written as a function of the guiding-centre position
$\boldsymbol{R}=\boldsymbol{r}-\boldsymbol{b}\times \boldsymbol{v}/\unicode[STIX]{x1D6FA}_{a}$
rather than the particle position
$\boldsymbol{r}$
. The parallel derivative is taken at constant magnetic moment
$\unicode[STIX]{x1D707}=m_{a}v_{\bot }^{2}/2B$
, and we shall use
$v$
and
$\unicode[STIX]{x1D706}=v_{\bot }^{2}/(v^{2}B)$
as our independent velocity-space variables. In addition to the mode frequency
$\unicode[STIX]{x1D714}=\unicode[STIX]{x1D714}_{r}+\text{i}\unicode[STIX]{x1D6FE}$
, two characteristic frequencies appear in equation (1), the drift frequency
$\unicode[STIX]{x1D714}_{da}=\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{v}_{da}$
, and the diamagnetic frequency
$\unicode[STIX]{x1D714}_{\star a}=(T_{a}k_{\unicode[STIX]{x1D6FC}}/e_{a})\,\text{d}\ln n_{a}/\text{d}\unicode[STIX]{x1D713}$
appearing in
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{\star a}^{\text{T}}=\unicode[STIX]{x1D714}_{\star a}\left[1+\unicode[STIX]{x1D702}_{a}\left(x^{2}-\frac{3}{2}\right)\right],\end{eqnarray}$$
with
$x^{2}=m_{a}v^{2}/2T_{a}$
. Here,
$n_{a}$
denotes the density,
$T_{a}$
the temperature,
$\unicode[STIX]{x1D702}_{a}=\text{d}\ln T_{a}/\text{d}\ln n_{a}$
and the wave vector has been written as
$\boldsymbol{k}_{\bot }=k_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}+k_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}$
. The system of equations is closed by the quasineutrality condition,
$$\begin{eqnarray}\mathop{\sum }_{a}\frac{n_{a}\text{e}_{a}^{2}}{T_{a}}\unicode[STIX]{x1D719}=\mathop{\sum }_{a}e_{a}\int g_{a}\text{J}_{0}\,\text{d}^{3}v.\end{eqnarray}$$
Apart from the usual assumptions in gyrokinetics, two approximations have been made in (4.1): electromagnetic effects and collisions have been neglected. The former are unimportant in the limit
$\unicode[STIX]{x1D6FD}\rightarrow 0$
, but it is in practice difficult to know a priori just how small
$\unicode[STIX]{x1D6FD}$
needs to be (typically below one or a few per cent, depending on the magnetic geometry). Collisions are negligible as long as the collision frequency is smaller than
$\unicode[STIX]{x1D714}f_{t}^{2}$
where
$f_{t}$
denotes the fraction of trapped particles.
From here we can obtain the dispersion relation obtained in Helander (Reference Helander2014) in the limit of zero Debye length and using the assumption that the mode frequency is much smaller than the bounce frequency, namely we obtain
$$\begin{eqnarray}\unicode[STIX]{x1D719}=\frac{1}{n}\int \bar{\unicode[STIX]{x1D719}}\left(1+\frac{\bar{\unicode[STIX]{x1D714}}_{d}(\bar{\unicode[STIX]{x1D714}}_{d}-\unicode[STIX]{x1D714}_{\star }^{\text{T}})}{\unicode[STIX]{x1D714}^{2}-\bar{\unicode[STIX]{x1D714}}_{d}^{2}}\right)f_{0}\,\text{d}^{3}v,\end{eqnarray}$$
where the overbar notation denotes the bounce average, defined for
$\unicode[STIX]{x1D719}$
as
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D719}}=\left.\int \frac{\unicode[STIX]{x1D719}(l)\,\text{d}l}{\sqrt{1-\unicode[STIX]{x1D706}B(l)}}\right/\!\int \frac{\text{d}l}{\sqrt{1-\unicode[STIX]{x1D706}B(l)}}\end{eqnarray}$$
the integrals which appear in this definition are taken along the magnetic field between two consecutive bounce points
$\unicode[STIX]{x1D706}B=1$
.
We can now declare an interest in modes satisfying
$\unicode[STIX]{x1D714}_{\star }^{\text{T}}\gg \bar{\unicode[STIX]{x1D714}}_{d}$
. This assumption is not satisfied in a general magnetic geometry but is well satisfied in the scenarios considered here. In this limit the dispersion relation reduces to
$$\begin{eqnarray}\unicode[STIX]{x1D719}=\frac{1}{n}\int \bar{\unicode[STIX]{x1D719}}\left(1-\frac{\bar{\unicode[STIX]{x1D714}}_{d}\unicode[STIX]{x1D714}_{\star }^{\text{T}}}{\unicode[STIX]{x1D714}^{2}}\right)f_{0}\,\text{d}^{3}v.\end{eqnarray}$$
From here we can make our first deduction. Namely, this dispersion relation is a linear function of the square of the complex frequency in which the coefficients are real variables. As such, it follows that
$\unicode[STIX]{x1D714}$
is either purely real (oscillatory modes) or purely imaginary (exponential growth). This explains why it is not unreasonable that we see the real frequency reaching exactly zero for electron–positron plasmas.
We can further explore this dispersion relation by transforming to pitch angle coordinates
$\text{d}^{3}v=\sum _{\unicode[STIX]{x1D70E}}(\unicode[STIX]{x03C0}v^{2}/\sqrt{1-\unicode[STIX]{x1D706}B})B\,\text{d}v\,\text{d}\unicode[STIX]{x1D706}$
where
$\unicode[STIX]{x1D70E}=\text{sgn}(v_{\Vert })$
and writing
$\bar{\unicode[STIX]{x1D714}}_{d}=x^{2}\tilde{\unicode[STIX]{x1D714}}_{d}(\unicode[STIX]{x1D706})$
. We can then reformulate (4.7) as a variational principle. This is accomplished by multiplying the equation by
$\unicode[STIX]{x1D719}^{\star }/B$
, where
$\unicode[STIX]{x1D719}^{\star }$
is the complex conjugate of
$\unicode[STIX]{x1D719}$
, and then integrating along the entire field line, using
$$\begin{eqnarray}\int _{-\infty }^{\infty }\unicode[STIX]{x1D719}^{\star }(l)\,\text{d}l\int _{1/B_{\text{max}}}^{1/B}\frac{\bar{\unicode[STIX]{x1D719}}_{j}\,\text{d}\unicode[STIX]{x1D706}}{\sqrt{1-\unicode[STIX]{x1D706}B}}=\int _{1/B_{\text{max}}}^{1/B_{\text{min}}}\mathop{\sum }_{j}\unicode[STIX]{x1D70F}_{j}|\bar{\unicode[STIX]{x1D719}}_{j}|^{2}\,\text{d}\unicode[STIX]{x1D706},\end{eqnarray}$$
where the sum is taken over all relevant magnetic wells, i.e. over all regions where the magnetic field strength
$B<1/\unicode[STIX]{x1D706}$
, and
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D719}}_{j}(l)=\frac{1}{\unicode[STIX]{x1D70F}_{j}(l)}\int \frac{\unicode[STIX]{x1D719}(l)\,\text{d}l}{\sqrt{1-\unicode[STIX]{x1D706}B}}\end{eqnarray}$$
denotes the bounce average of
$\unicode[STIX]{x1D719}$
over the
$j$
th such well, with
$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{j}(\unicode[STIX]{x1D706})=\int \frac{\text{d}l}{\sqrt{1-\unicode[STIX]{x1D706}B}}.\end{eqnarray}$$
Performing this calculation we obtain
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=S[\unicode[STIX]{x1D719}]=\frac{N[\unicode[STIX]{x1D719}]}{D[\unicode[STIX]{x1D719}]},\end{eqnarray}$$
where
$$\begin{eqnarray}N[\unicode[STIX]{x1D719}]=-\frac{3}{4}\int _{1/B_{\text{max}}}^{1/B_{\text{min}}}\mathop{\sum }_{j}\unicode[STIX]{x1D70F}_{j}(1+\unicode[STIX]{x1D702})\unicode[STIX]{x1D714}_{\star }\tilde{\unicode[STIX]{x1D714}}_{d}|\bar{\unicode[STIX]{x1D719}}_{j}|^{2}\,\text{d}\unicode[STIX]{x1D706}\end{eqnarray}$$
and
$$\begin{eqnarray}D[\unicode[STIX]{x1D719}]=\int _{-\infty }^{\infty }|\unicode[STIX]{x1D719}|^{2}\frac{\text{d}l}{B}-\frac{1}{2}\int _{1/B_{\text{max}}}^{1/B_{\text{min}}}\mathop{\sum }_{j}\unicode[STIX]{x1D70F}_{j}|\bar{\unicode[STIX]{x1D719}}_{j}|^{2}\,\text{d}\unicode[STIX]{x1D706}.\end{eqnarray}$$
The denominator of (4.11),
$D[\unicode[STIX]{x1D719}]$
, is positive since the Cauchy–Schwarz inequality,
$|\bar{\unicode[STIX]{x1D719}}|^{2}\leqslant \overline{|\unicode[STIX]{x1D719}|^{2}}$
, implies
$$\begin{eqnarray}\frac{1}{2}\int _{1/B_{\text{max}}}^{1/B_{\text{min}}}\mathop{\sum }_{j}\unicode[STIX]{x1D70F}_{j}|\bar{\unicode[STIX]{x1D719}}_{j}|^{2}\,\text{d}\unicode[STIX]{x1D706}\leqslant \frac{1}{2}\int |\unicode[STIX]{x1D719}^{2}|\,\text{d}l\int _{1/B_{\text{max}}}^{1/B_{\text{min}}}\frac{\text{d}\unicode[STIX]{x1D706}}{\sqrt{1-\unicode[STIX]{x1D706}B}}\leqslant \sqrt{1-\frac{B_{\text{min}}}{B_{\text{max}}}}\int |\unicode[STIX]{x1D719}|^{2}\frac{\text{d}l}{B}.\end{eqnarray}$$
Hence, we see that the stability of the modes is driven by the pressure gradient through the quantity
$(1+\unicode[STIX]{x1D702})\unicode[STIX]{x1D714}_{\star }\bar{\unicode[STIX]{x1D714}}_{d}.$
The sign of this quantity sets the stability.
Equation (4.11) is variational and assumes its minimum for the particular function
$\unicode[STIX]{x1D719}(l)$
which satisfies the integral equation (4.7), as follows from the fact that the variation vanishes,
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D714}^{2}}=\frac{\unicode[STIX]{x1D6FF}N}{N}-\frac{\unicode[STIX]{x1D6FF}D}{D}=\frac{1}{N}(\unicode[STIX]{x1D6FF}N-\unicode[STIX]{x1D714}^{2}\unicode[STIX]{x1D6FF}D)=0\end{eqnarray}$$
if, and only if, the integral equation (4.7) is satisfied. This follows directly from
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}N[\unicode[STIX]{x1D719}]=-\frac{3}{2}\int _{1/B_{\text{max}}}^{1/B_{\text{min}}}\mathop{\sum }_{j}\unicode[STIX]{x1D70F}_{j}(1+\unicode[STIX]{x1D702})\unicode[STIX]{x1D714}_{\star }\tilde{\unicode[STIX]{x1D714}}_{d}\bar{\unicode[STIX]{x1D719}}_{j}\overline{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}}_{j}\,\text{d}\unicode[STIX]{x1D706}\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}D[\unicode[STIX]{x1D719}]=2\int _{-\infty }^{\infty }\unicode[STIX]{x1D719}\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}\frac{\text{d}l}{B}-\int _{1/B_{\text{max}}}^{1/B_{\text{min}}}\mathop{\sum }_{j}\unicode[STIX]{x1D70F}_{j}\bar{\unicode[STIX]{x1D719}}_{j}\overline{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D719}}_{j}\,\text{d}\unicode[STIX]{x1D706}.\end{eqnarray}$$
We can use this variational principle to deduce that any unstable mode is purely growing i.e.
$\unicode[STIX]{x1D714}=\text{i}\unicode[STIX]{x1D6FE}$
provided there is at least one value of
$\unicode[STIX]{x1D706}$
for which the quantity
$(1+\unicode[STIX]{x1D702})\unicode[STIX]{x1D714}_{\star }\bar{\unicode[STIX]{x1D714}}_{d}$
is positive. We note that for the most unstable mode
$$\begin{eqnarray}\min _{\unicode[STIX]{x1D719}}S[\unicode[STIX]{x1D719}]=\unicode[STIX]{x1D714}^{2},\end{eqnarray}$$
where
$\unicode[STIX]{x1D714}^{2}$
is the frequency associated with the most unstable mode. Suppose now that
$(1+\unicode[STIX]{x1D702})\unicode[STIX]{x1D714}_{\star }\bar{\unicode[STIX]{x1D714}}_{d}$
is positive for at least one value of
$\unicode[STIX]{x1D706}$
, then we can construct a trial function
$\unicode[STIX]{x1D719}_{\text{trial}}$
, which is positive when
$(1+\unicode[STIX]{x1D702})\unicode[STIX]{x1D714}_{\star }\bar{\unicode[STIX]{x1D714}}_{d}$
is positive and zero otherwise. Note that
$\unicode[STIX]{x1D719}_{\text{trial}}$
need not (and indeed will not generally) satisfy the integral equation (4.7). It follows immediately that for such a trial function
$$\begin{eqnarray}S[\unicode[STIX]{x1D719}_{\text{trial}}]<0.\end{eqnarray}$$
As the functional
$S[\unicode[STIX]{x1D719}]$
attains its minimum for the value of
$\unicode[STIX]{x1D719}$
which does satisfy the integral equation, it follows that for the most unstable mode
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}=S[\unicode[STIX]{x1D719}]\leqslant S[\unicode[STIX]{x1D719}_{\text{trial}}]<0.\end{eqnarray}$$
Hence the most unstable mode is purely growing. This allows us to explain the trends shown in figures 1 and 7. Equation (4.11) tells us that
$\unicode[STIX]{x1D714}$
is purely real or purely imaginary. The curvature is bad for at least one value of
$\unicode[STIX]{x1D706}$
which implies that
$\unicode[STIX]{x1D714}^{2}<0$
for any unstable mode. Hence we deduce that
$\unicode[STIX]{x1D714}=\text{i}\unicode[STIX]{x1D6FE}$
for the any unstable mode in electron–positron plasmas.
One final result to which attention should be drawn is the following. If we take some estimate for the function
$\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{0}$
then provided
$$\begin{eqnarray}|\unicode[STIX]{x1D719}-\unicode[STIX]{x1D719}_{0}|<\unicode[STIX]{x1D716}\end{eqnarray}$$
the variational principle ensures
$$\begin{eqnarray}\unicode[STIX]{x1D714}^{2}-S[\unicode[STIX]{x1D719}_{0}]<\text{C}\unicode[STIX]{x1D716}^{2}\end{eqnarray}$$
for some constant
$\text{C}.$
That is, we obtain higher accuracy for the frequency than for the associated eigenfunction.
5 Conclusions
We have used the gyrokinetic code GENE to examine microinstabilities driven by temperature and density gradients in plasmas confined in different magnetic geometries for a range of different mass ratios between a conventional hydrogen plasma and electron–positron pair plasmas.
We considered an axisymmetric, circular flux surface, low
$\unicode[STIX]{x1D6FD}$
(tokamak) and were able to reproduce the results of Pedersen et al. (Reference Pedersen, Boozer, Dorland, Kremer and Schmitt2003) and also examine the behaviour of the real frequency of the modes as we varied the mass ratio. The change in regime from ion temperature-gradient-driven instabilities to trapped-particle-driven instabilities is clear as the plasma mass ratio approaches that of electron–positron plasmas. We considered a quasi-isodynamic stellarator configuration using the geometry of the W7-X device. We found similar behaviour of the growth rate and real frequency compared to the tokamak case. In particular we were able to identify two distinct regimes dominated by modes propagating in the electron diamagnetic direction (for mass ratios close to conventional values
$m_{e}/m_{i}=1/1836$
) and modes propagating in the ion/positron diamagnetic direction (for mass ratios close to a pair plasma
$m_{e}/m_{i}=1$
). In both the tokamak and the stellarator case we observed that the real frequency tends to zero as the mass ratio approaches unity.
Acknowledgements
We acknowledge T. S. Perdersen and the PAX/APEX experiment team for their interest in our work. We acknowledge T. Görler for helpful discussions and technical support with the GENE code.
