# The physics of a small-scale tearing mode in collisionless slab plasmas

Chen Geng1, David Dickinson1 and Howard Wilson1,2 1York Plasma Institute,
Department of Physics, University of York, Heslington, York. YO10 5DD UK
2Culham Centre for Fusion Energy, Abingdon, Oxfordshire. OX14 3DB UK
<EMAIL_ADDRESS>
###### Abstract

Microtearing modes have been widely reported as a tearing parity electron
temperature gradient driven plasma instability, which leads to fine scale
tearing of the magnetic flux surfaces thereby resulting in reconnection of
magnetic field lines and formation of magnetic islands. In slab geometry it
has previously been shown that the drive mechanism requires a finite collision
frequency. However, we find in linear gyrokinetic simulations that a
collisionless fine-scale tearing parity instability exists even at low and
zero collision frequency. Detailed studies reveal that these slab modes are
also driven by electron temperature gradient but are sensitive to electron
finite Larmor radius effects, and have a radial wavenumber much smaller than
the binormal wavenumber, which is comparable to the ion Larmor radius.
Furthermore, they exist even in the electrostatic limit and electromagnetic
effects actually have a stabilising influence on this collisionless tearing
mode. An analytic model shows that this collisionless small scale tearing mode
is consistent with a tearing parity slab electron temperature gradient (ETG)
mode, which can be more unstable than the twisting parity ETG mode that is
often studied. This small-scale tearing parity mode can lead to magnetic
islands, which, in turn, can influence turbulent transport in magnetised
plasmas.

  * •

##  1 Introduction

In general, electromagnetic micro-instabilities in magnetised plasmas can be
categorised as tearing or twisting parity modes. Tearing parity modes, in
which the fluctuating parallel component of the magnetic potential is an even
function about the rational surface, perturb the magnetic field to form
magnetic islands. Twisting parity modes have a parallel component of the
magnetic potential which is odd about the rational surface, and cause a
rippling of the flux surface. Microtearing modes (MTMs) are a type of tearing
parity micro-instability.

MTMs are a candidate for anomalous electron heat transport in tokamak plasmas
[1, 2, 3]. They have been studied extensively since the 1960s. They are
characterised by large toroidal and poloidal wavenumbers, comparable to the
reciprocal of the ion Larmor radius. An early analytic linear model for MTMs
was developed by J. F. Drake et. al [4]. In that work, the main drive
mechanism is shown to come from the free energy in the electron temperature
gradient, contributing to electromagnetic fluctuations. The collision
frequency was shown to be a key factor in this drive mechanism. In the simple
two dimensional sheared slab geometry, they studied the impact of collision
frequency by dividing it into high-collisional, semi-collisional and
collisionless regimes, and predicted MTMs to be stable at both low and high
collision frequency. Numerical calculations by N. T. Gladd et. al [5]
confirmed these slab results, demonstrating that a velocity dependent
collision operator is essential for instability.

More recent studies [6, 7] have observed MTMs in simulations neglecting the
velocity dependence of the collision operator and even in the limit that the
collision frequency tends to zero. Furthermore, gyrokinetic simulations have
found microtearing modes can exist towards the edge of MAST tokamak plasmas
and that these modes can be unstable at low collision frequency in toroidal
geometry [7]. Gyrokinetic simulations [8, 9, 10] have also demonstrated
unstable MTMs in the complete collisionless limit in a range of scenarios.
These are at odds with the slab results presented in [5, 4] but the mechanism
is not yet fully understood. Understanding the collisionless drive mechanism
is vital for clarifying the impact for transport in tokamak plasmas,
especially those operating at higher temperature such as ITER.

In this article, we show that a fine scale collisionless tearing parity mode
can, in fact, be unstable even in slab geometry. To identify the key physics,
we develop an analytic model of this collisionless microtearing instability in
slab geometry and conclude that the drive mechanism persists even in the
electrostatic limit, with finite electron Larmor radius effects playing an
important role. Probing the model in more detail, we show that the instability
is the tearing parity branch of the electron temperature gradient (ETG) mode,
which can be more unstable than the usual twisting parity ETG mode.

The paper is laid out as follows. In the next section we describe the magnetic
geometry. In Section 3 we employ the GS2 code to demonstrate the existence of
a microtearing instability in a collisionless slab. In Section 4 we discuss
details of our analytic model, which we use to identify the main physics
mechanisms in Section 5. We close in Section 6 with conclusions.

##  2 Slab geometry

Figure 1: Illustration of the slab geometry.

It is convenient to define the slab geometry before our discussion of the
physical plasma instability. We consider a simple infinite slab of plasma with
magnetic field lines in the $y-z$ plane and with density and temperature
gradients in the $x$ direction. The scale lengths are
$L_{n}^{-1}=-\mathrm{d}\mathrm{ln}n/\mathrm{d}x$ and
$L_{T}^{-1}=-\mathrm{d}\mathrm{ln}T/\mathrm{d}x$, respectively. Using $L_{T}$
as the reference length, the normalised temperature gradient is defined as
$\eta=L_{n}/L_{T}$. An external magnetic field $\bm{B}$ and a current density
$\bm{J}$ are applied along the $z$ direction, resulting in
$\bm{B}=B_{0}(\bm{\hat{z}}+(x/L_{s})\bm{\hat{y}})$, where $L_{s}$ represents
the scale length of the magnetic field shearing. We assume that
$|B_{y}|\ll|B_{z}|$, so restrict consideration to $x\ll L_{s}$.

##  3 GS2 simulation

GS2 is an initial value simulation code solving the gyrokinetic Vlasov-Maxwell
equations using an implicit algorithm [11]. It employs local flux tubes and is
designed to operate in a range of magnetic geometries including general
tokamak, cylindrical and slab plasmas.

We first employ GS2 (version v8.0.1 [12]) to benchmark the numerical results
obtained by Gladd et al [5] in the slab geometry. Ion and electron
temperatures are equal at the centre of the slab; however, the ion temperature
gradient is zero while the electron temperature gradient is finite. There is
also a finite density gradient and sheared magnetic field applied as described
in Section 2. The scale lengths for the sheared magnetic field and for each
species’ temperature gradient and density gradient are $L_{s}$, $L_{T}$ and
$L_{n}$, respectively. The mode frequency and growth rate in this paper are
normalised to the electron diamagnetic frequency
$\omega_{*e}=k_{y}v_{e}\rho_{e}/2L_{n}$. Here, the wavenumber $k_{y}$ and
spatial coordinate $x$ will be normalised to the ion gyro radius,
$\rho_{i}=v_{i}/\omega_{ci}$, where $v_{j}=\sqrt{2T_{j}/m_{j}}$ is the thermal
speed of species $j$.

Figure 2: The growth rate (top) and mode frequency (bottom) for microtearing
modes as a function of collision frequency. The triangle and circular symbols
are GS2 simulation results, while the solid and dashed lines are numerical
solutions of eigenmode equations provided in reference [5]. Two electron
temperature gradients, $\eta_{e}=5.0$ and $\eta_{e}=7.0$ are shown. Other
physical parameters include $k_{y}\rho_{i}=0.3$, $\beta=8\pi
n_{0}T/B^{2}=0.005$, $m_{i}/m_{e}=1836$ and $T_{i}=T_{e}$. The numerical
parameters in GS2 simulations are set as $nperiod=128$ and $ntheta=8$ for
collision frequencies smaller than $10\omega_{*e}$, while $nperiod=32$ and
$ntheta=8$ are used for the other cases. Figure 3: The mode structures for
collisional (top) and collisionless (bottom) microtearing instabilities in GS2
simulations. Note the difference in the abscissa scale. The left two panels
are normalised parallel magnetic potential and the right two are normalised
electrostatic potential. The real and imaginary parts are shown with solid and
dashed lines respectively. The eigenmodes are normalised such that
$\hat{A}_{\|}(x=0)=1$. Here $\eta_{e}=5.0$; other parameters are kept the same
as for figure (2).

As shown in figure (2) the linear GS2 results match well with Gladd’s model
[5] in the collisional regime, and both demonstrate the drive from the
electron temperature gradient. However, in the very low collision frequency
regime, GS2 reveals an unexpected tearing parity instability in this slab
geometry. The frequency is not continuous between the collisional and
collisionless regions, indicating that they are different instability
branches. The real and imaginary parts of the normalised electrostatic
potential $\bar{\phi}$ and normalised parallel magnetic potential
$\bar{A}_{\|}$ are shown in figure (3). $\bar{A}_{\|}$ has an even symmetry
while $\bar{\phi}$ is odd, which is a defining feature of tearing parity
modes. The collisionless one has a narrower structure, thus the characteristic
radial wavenumber $k_{x}\rho_{i}$ is much larger for the collisionless branch
than for the collisional one (see figure (4)). In fact, in these GS2
simulations we note that capturing the collisionless instability requires
challenging numerical settings. The parallel grid extent and resolution needs
to be sufficiently high to capture the unstable mode accurately. In GS2, the
parallel flux tube extent is controlled by $nperiod$, while $ntheta$ defines
the grid resolution within each $2\pi$ period. In our simulations, the
collisionless branch requires $nperiod=128$ and $ntheta=8$, while the
collisional branch is well converged for $nperiod=32$ and $ntheta=8$.

Note that the condition $k_{x}\rho_{i}\ll 1$ is assumed in the derivation of
references [4, 5]. This enables a Gamma function expansion in the quasi-
neutrality equation, ignoring the finite Larmor radius effects from electrons.
Figure (4) tests the validity of this assumption for the range of collision
frequencies, with $\eta_{e}=5.0$ and $\eta_{e}=7.0$. It shows that electron
temperature gradient $\eta_{e}$ has very little influence on the value of
$k_{x}\rho_{i}$ but collision frequency $\nu/\omega_{*e}$ has a big impact. In
the collisional regime, $k_{x}\rho_{i}$ remains small and the approximation is
valid; however, this is not the case for the collisionless regime.
Specifically, in the collisionless regime, $k_{x}\rho_{i}$ is about 40 times
larger, which leads to $k_{x}\rho_{e}$ approaching 1. This gives the first
clue that the finite Larmor radius effects from electrons might be an
important factor for the collisionless mode seen here. This finding informs a
new reduced gyrokinetic model, to be derived in Section 4.

Figure 4: The values of $k_{x}\rho_{i}$ (left axis) and $k_{x}\rho_{e}$ (right
axis) as a function of collision frequency in the GS2 simulations. Parameters
are kept the same as for figure (2). Figure 5: The effect of the Bessel factor
$\alpha$ on the collisionless instability (left panel, $\nu/\omega_{*e}=0.1$)
and the collisional MTM (right panel, $\nu/\omega_{*e}=20$). Here
$\eta_{e}=5.0$; other parameters are kept the same as for figure (2).

To further test this point, we examined the influence of electron finite
Larmor radius effects directly by probing a Bessel function parameter in GS2.
This parameter, $\alpha$, enables a suppression of finite Larmor radius
effects in the gyro-averaging Bessel function $\mathrm{J_{0}}(\alpha k)$.
$\alpha=1$ is the default case capturing full gyrokinetic physics. Turning
$\alpha$ down towards zero is equivalent to turning off the finite Larmor
radius effects for the given species in GS2 simulations. Figure (5) shows the
effects of the Bessel factor on the collisionless and collisional branches in
GS2 simulation. Here the electron temperature gradient is $\eta_{e}=5.0$ and
the collision frequency is set to $\nu/\omega_{*e}=0.1$ for the collisionless
simulation. Note that the gyro-averaging provides velocity dependent
dissipation, which is in some sense similar to the collision operator. Whilst
the collisionality is an important factor in the collisional slab model, this
indicates that for the collisionless instability, the electron finite Larmor
radius effects are required to confine the mode; these are neglected in the
collisional model. On the contrary, the collisional instability is insensitive
to the electron finite Larmor radius effects.

Figure 6: The effect of plasma beta $\beta$ on the collisionless instability.
Here $\nu/\omega_{*e}=0.0$, $\eta_{e}=5.0$; other parameters are kept the same
as for figure (2).

Focusing on the zero collision frequency limit and varying the ratio of plasma
pressure to magnetic pressure, $\beta$, we found that this mode persists in
the electrostatic limit, as shown in figure (6). Indeed, it is more unstable
at lower $\beta$, and still exists when $\beta=0$. This means that
electromagnetic effects are stabilising for this mode, which is fundamentally
electrostatic in nature. Meanwhile, GS2 simulations with kinetic and adiabatic
ions demonstrate that the ion treatment has little impact on the collisionless
mode. This provides conclusive evidence that the main drive comes from
electrostatic electron physics. To capture the physics of the collisionless
mode and provide an interpretation of the GS2 results, we develop a new model
in the following section.

##  4 Modelling in slab geometry

In order to provide a physics interpretation of the GS2 simulation at low
collision frequency, we derive eigenmode equations valid in this limit. Here
we present two models derived from gyrokinetic theory. Section 4.1 describes a
simple case, focusing on zero collision frequency and zero $\beta$, which
demonstrates just the essential physics. Section 4.2 considers more
complicated factors including finite but small collision frequency and
electromagnetic effects, which can be compared in more detail with GS2 results
and help give a good foundation for future work.

###  4.1 Electrostatic model at zero collision frequency

Informed by the earlier GS2 results, we adopt an adiabatic ion response but
treat electrons kinetically, retaining finite Larmor radius effects. We first
consider the electrostatic limit, with perturbations only in the electrostatic
potential $\phi$. In Fourier space the gyrokinetic equation for electrons
yields:

$\left(\omega-i\frac{k_{y}}{L_{s}}v_{\|}\frac{\partial}{\partial
k}\right)\hat{g}(k)=-\frac{e}{T}\frac{n_{0}}{\pi^{3/2}v_{e}^{3}}\mathrm{e}^{-v^{2}/v_{e}^{2}}\left(\omega-\omega_{*e}^{T}\right)\mathrm{J_{0}}(k_{\perp}\rho_{\perp})\hat{\phi}(k)$
(1)

in which $\hat{g}(k)$ is the electron distribution perturbation in Fourier
space, $v_{e}$ is the electron thermal velocity,
$\rho_{\perp}=v_{\perp}/\omega_{ce}$ is the perpendicular velocity-dependent
electron gyro radius, $T_{e}=T_{i}=T$,
$\omega_{*e}^{T}=\omega_{*e}(1+\eta_{e}(\frac{v^{2}}{v_{e}^{2}}-\frac{3}{2}))$,
$k_{\perp}^{2}=k_{x}^{2}+k_{y}^{2}$ and $\mathrm{J_{0}}$ is the Bessel
function. The quasi-neutrality equation is

$n_{0}\frac{e\hat{\phi}(k)}{T}+\int_{-\infty}^{\infty}\mathrm{d}^{3}\bm{v}\cdot\hat{g}(k)\mathrm{J_{0}}(k_{\perp}\rho_{\perp})=-n_{0}\frac{e\hat{\phi}(k)}{T}$
(2)

Assuming small $k_{\perp}\rho_{\perp}$, expanding Bessel functions to second
order and conducting an inverse Fourier transform to real space, equations (1)
and (2) become:

$\displaystyle\left(\omega-\frac{k_{y}x}{L_{s}}v_{\|}\right)g(x)=-\frac{e}{T}\frac{n_{0}}{\pi^{3/2}v_{e}^{3}}\mathrm{e}^{-v^{2}/v_{e}^{2}}\left(\omega-\omega_{*e}^{T}\right)\left(1-\frac{k_{y}^{2}v_{\perp}^{2}}{4\omega_{ce}^{2}}+\frac{v_{\perp}^{2}}{4\omega_{ce}^{2}}\frac{\partial^{2}}{\partial
x^{2}}\right)\phi(x)$ (3) $\displaystyle
2n_{0}\frac{e\phi(x)}{T}=-\int_{-\infty}^{\infty}\mathrm{d}^{3}\bm{v}\cdot\left(1-\frac{k_{y}^{2}v_{\perp}^{2}}{4\omega_{ce}^{2}}+\frac{v_{\perp}^{2}}{4\omega_{ce}^{2}}\frac{\partial^{2}}{\partial
x^{2}}\right)g(x)$ (4)

Substituting (3) into (4) and normalising the variables as
$\bar{\omega}=\omega/\omega_{*e}$, $\bar{k_{y}}=k_{y}\rho_{e}$,
$\bar{x}=x/\rho_{e}$, $\bar{\phi}=e\phi/T$, we have

$\begin{split}\sqrt{\pi}\bar{\phi}=&\int_{-\infty}^{\infty}\mathrm{d}s\int_{0}^{\infty}t\mathrm{d}t\cdot\mathrm{e}^{-(s^{2}+t^{2})}\left(\bar{\omega}-1-\eta\left(s^{2}+t^{2}-\frac{3}{2}\right)\right)\times\\\
&\left[\left(1-\frac{\bar{k_{y}}^{2}t^{2}}{4}\right)\left(\frac{\left(1-\frac{\bar{k_{y}}^{2}t^{2}}{4}\right)\bar{\phi}+\frac{t^{2}}{4}\frac{\partial^{2}}{\partial\bar{x}^{2}}\bar{\phi}}{\bar{\omega}-2\epsilon\bar{x}s}\right)+\frac{t^{2}}{4}\frac{\partial^{2}}{\partial\bar{x}^{2}}\left(\frac{\left(1-\frac{\bar{k_{y}}^{2}t^{2}}{4}\right)\bar{\phi}+\frac{t^{2}}{4}\frac{\partial^{2}}{\partial\bar{x}^{2}}\bar{\phi}}{\bar{\omega}-2\epsilon\bar{x}s}\right)\right]\end{split}$
(5)

in which $s=v_{\|}/v_{e}$, $t=v_{\perp}/v_{e}$, $\epsilon=L_{n}/L_{s}$ and
$\eta=L_{n}/L_{T}$. Please note that we have normalised lengths to the
electron Larmor radius rather than the ion Larmor radius in the previous
sections. We consider $\bar{k_{y}}\ll 1$ in which case it can be neglected.
Neglecting third and fourth orders of the expansion in $k_{x}\rho_{e}$,
simplification of this equation yields a second order differential equation
for the electrostatic potential of the form
$\mathrm{C_{0}}\bar{\phi}+\mathrm{C_{1}}\bar{\phi}^{\prime}+\mathrm{C_{2}}\bar{\phi}^{\prime\prime}=0$,
where primes denote the differential with respect to $\bar{x}$ and the
coefficients $\mathrm{C_{0}}$, $\mathrm{C_{1}}$ and $\mathrm{C_{2}}$ are

$\displaystyle\mathrm{C_{0}}$
$\displaystyle=-\sqrt{\pi}-\frac{1}{4\epsilon\bar{x}}\left[\left(\bar{\omega}-1+\frac{1}{2}\eta\right)\mathrm{Z_{0,0}}-\eta\mathrm{Z_{2,0}}\right]-\frac{1}{8\epsilon\bar{x}^{3}}\left[\left(\bar{\omega}-1-\frac{1}{2}\eta\right)\mathrm{Z_{2,2}}-\eta\mathrm{Z_{4,2}}\right]$
(6) $\displaystyle\mathrm{C_{1}}$
$\displaystyle=\frac{1}{8\epsilon\bar{x}^{2}}\left[\left(\bar{\omega}-1-\frac{1}{2}\eta\right)\mathrm{Z_{1,1}}-\eta\mathrm{Z_{3,1}}\right]$
(7) $\displaystyle\mathrm{C_{2}}$
$\displaystyle=-\frac{1}{8\epsilon\bar{x}}\left[\left(\bar{\omega}-1-\frac{1}{2}\eta\right)\mathrm{Z_{0,0}}-\eta\mathrm{Z_{2,0}}\right]-\frac{1}{16\epsilon\bar{x}^{3}}\left[\left(\bar{\omega}-1-\frac{3}{2}\eta\right)\mathrm{Z_{2,2}}-\eta\mathrm{Z_{4,2}}\right]$
(8)

Here, $\mathrm{Z_{m,n}}=\mathrm{Z_{m,n}}(\bar{\omega}/2\epsilon\bar{x})$ is a
generalised plasma dispersion function

$\mathrm{Z_{m,n}}(\alpha)=\int_{-\infty}^{\infty}\frac{\mathrm{e}^{-s^{2}}s^{m}}{(s-\alpha)^{n+1}}\mathrm{d}s\quad,\quad\alpha\in\mathcal{C}\,,\,s\in\mathcal{R}\,\,\mathrm{and}\,\,m,n\in\mathcal{N}$
(9)

It can be shown that $\mathrm{Z_{0,0}}(\alpha)=i\pi\mathrm{W}(\alpha)$ where
$\mathrm{W}(\alpha)$ is the Faddeeva function. When $mn\neq 0$, there is a
pair of recurrence relations which can be used to relate
$\mathrm{Z_{m,n}}(\alpha)$ to $\mathrm{Z_{0,0}}(\alpha)$:

$\displaystyle\mathrm{Z_{m,n}}(\alpha)=\frac{m}{n}\mathrm{Z_{m-1,n-1}}(\alpha)-\frac{2}{n}\mathrm{Z_{m+1,n-1}}(\alpha)\qquad,\qquad
n\geq 1$ (10)
$\displaystyle\mathrm{Z_{m+1,0}}(\alpha)=\alpha\mathrm{Z_{m,0}}(\alpha)+\frac{(-1)^{m}+1}{2}\mathrm{\Gamma}(\frac{m+1}{2})$
(11)

When near the centre of the slab, where $\bar{x}=0$, the above coefficients
are well-defined with the limit of

$\displaystyle\mathrm{C_{0}}(\bar{x}=0)$
$\displaystyle=\sqrt{\pi}\left[-1+\frac{\bar{\omega}-1}{2\omega}+\frac{\epsilon^{2}\left(\bar{\omega}-1-2\eta\right)}{2\bar{\omega}^{3}}\right]$
(12) $\displaystyle\mathrm{C_{1}}(\bar{x}=0)$ $\displaystyle=0$ (13)
$\displaystyle\mathrm{C_{2}}(\bar{x}=0)$
$\displaystyle=\sqrt{\pi}\left[\frac{\bar{\omega}-1-\eta}{4\omega}+\frac{\epsilon^{2}\left(\bar{\omega}-1-3\eta\right)}{4\bar{\omega}^{3}}\right]$
(14)

The forms of these coefficients show that, when normalising to the diamagnetic
frequency $\omega_{*e}$, the mode frequency and growth rate are mostly
sensitive to magnetic shear scale length and electron temperature gradient.
Numerical solutions of this second order differential equation for
$\bar{\phi}$ are presented in Section 5.

###  4.2 Electromagnetic model with finite Lorentz collision operator

The above electrostatic model is valid only when the collision frequency is
zero. To compare with the GS2 results along the low collision frequency range,
we consider a classic Lorentz collision operator consisting of pitch-angle
scattering
$\mathrm{C}(\nu)=-\frac{i\nu}{2}\frac{\partial}{\partial\xi}(1-\xi^{2})\frac{\partial}{\partial\xi}$,
where $\nu$ is the collision frequency and $\xi=v_{\|}/v$ is the pitch angle.
Furthermore, it is useful to explore the influence of electromagnetic effects
including $\beta$, which will show the tendency of this mode to form magnetic
islands. Before we address this, we first update the above model. When
including the parallel magnetic potential and the Lorentz collision operator,
the gyrokinetic equation can be rewritten as

$\left(\omega-i\frac{k_{y}}{L_{s}}v\xi\frac{\partial}{\partial
k}-\frac{i\nu}{2}\frac{\partial}{\partial\xi}\left(1-\xi^{2}\right)\frac{\partial}{\partial\xi}\right)\hat{g}(k)=-\frac{e}{T}\frac{n_{0}}{\pi^{3/2}v_{e}^{3}}\mathrm{e}^{-v^{2}/v_{e}^{2}}\left(\omega-\omega_{*e}^{T}\right)\mathrm{J_{0}}(k_{\perp}\rho_{\perp})\left(\hat{\phi}(k)-v\xi\hat{A}_{\|}(k)\right)$
(15)

Again, expanding the Bessel function and conducting an inverse Fourier
transform results in

$\begin{split}\left(\omega-\frac{k_{y}x}{L_{s}}v\xi-\frac{i\nu}{2}\frac{\partial}{\partial\xi}\left(1-\xi^{2}\right)\frac{\partial}{\partial\xi}\right)g(x)=&-\frac{e}{T}\frac{n_{0}}{\pi^{3/2}v_{e}^{3}}\mathrm{e}^{-v^{2}/v_{e}^{2}}\left(\omega-\omega_{*e}^{T}\right)\times\\\
&\left(1-\frac{k_{y}^{2}v^{2}(1-\xi^{2})}{4\omega_{ce}^{2}}+\frac{v^{2}(1-\xi^{2})}{4\omega_{ce}^{2}}\frac{\partial^{2}}{\partial
x^{2}}\right)\left(\phi(x)-v\xi A_{\|}(x)\right)\\\ \end{split}$ (16)

Note that the perturbation of the electron distribution function depends on
both space and velocity $g(x)=g(x,v,\xi)$. Expanding the distribution function
in an orthogonal polynomial series,
$g(x,v,\xi)=\sum_{n=0}^{\infty}h_{n}(x,v)\mathrm{P}_{n}(\xi)$ in which
$\mathrm{P}_{n}$ is the Legendre polynomial of order $n$, we have

$\begin{split}\sum_{n=0}^{\infty}h_{n}&\left[\left(\omega+\frac{i\nu}{2}n(n+1)\right)\mathrm{P}_{n}(\xi)-\frac{k_{y}x}{L_{s}}v\frac{(n+1)\mathrm{P}_{n}(\xi)+n\mathrm{P}_{n+1}(\xi)}{2n+1}\right]=\\\
&-\frac{e}{T}\frac{n_{0}}{\pi^{3/2}v_{e}^{3}}\mathrm{e}^{-v^{2}/v_{e}^{2}}\left(\omega-\omega_{*e}^{T}\right)\left(1-\frac{k_{y}^{2}v^{2}(1-\xi^{2})}{4\omega_{ce}^{2}}+\frac{v^{2}(1-\xi^{2})}{4\omega_{ce}^{2}}\frac{\partial^{2}}{\partial
x^{2}}\right)\left(\phi(x)-v\xi A_{\|}(x)\right)\\\ \end{split}$ (17)

Applying the orthogonality relations for Legendre polynomials and integrating
over pitch angle $\xi$ from $-1$ to $1$ on both sides yields a set of
equations

$\begin{split}&\frac{2}{n+1}\left[\left(\omega+\frac{i\nu}{2}n(n+1)\right)h_{n}-k_{\|}v\left(\frac{n}{2n-1}h_{n-1}+\frac{n+1}{2n+3}h_{n+1}\right)\right]=\\\
&-\frac{e}{T}\frac{n_{0}}{\pi^{3/2}v_{e}^{3}}\mathrm{e}^{-v^{2}/v_{e}^{2}}\left(\omega-\omega_{*e}^{T}\right)\times\begin{cases}\left[\left(2-\frac{k_{y}^{2}v^{2}}{3\omega_{ce}^{2}}\right)\phi+\frac{v^{2}}{3\omega_{ce}^{2}}\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}\phi\right]&\text{if
$n=0$}\\\
\left[-\left(\frac{2}{3}v-\frac{4k_{y}^{2}v^{3}}{15\omega_{ce}^{2}}\right)A_{\|}-\frac{4v^{3}}{15\omega_{ce}^{2}}\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}A_{\|}\right]&\text{if
$n=1$}\\\
\left[\frac{k_{y}^{2}v^{2}}{15\omega_{ce}^{2}}\phi-\frac{v^{2}}{15\omega_{ce}^{2}}\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}\phi\right]&\text{if
$n=2$}\\\
\left[-\frac{k_{y}^{2}v^{3}}{35\omega_{ce}^{2}}A_{\|}+\frac{v^{3}}{35\omega_{ce}^{2}}\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}A_{\|}\right]&\text{if
$n=3$}\\\ 0&\text{if $n\geq 4$}\end{cases}\end{split}$ (18)

To derive a tractable model from the above, we adopt a matrix approach. The
equation (18) can be written in the matrix form as $\bm{M}\cdot\bm{h}=\bm{D}$.
The two dimensional matrix $\bm{M}$ is an infinite tridiagonal matrix, of
which the $n$-th row ($n$ starts from 1) is

$\begin{Bmatrix}\cdots&-\frac{(n-1)k_{y}xv}{(2n-3)(2n-1)L_{s}}&\frac{1}{2n-1}\left(\omega+\frac{n(n-1)}{2}i\nu\right)&-\frac{nk_{y}xv}{(2n-1)(2n+1)L_{s}}&\cdots\end{Bmatrix}_{\text{$n$-th
row}}$

The column vector $\bm{h}$ starts from the $h_{0}$ term and the column vector
$\bm{D}$ represents the driving terms in the right hand side of equation (18).
Note that the main difference between this model and collisional model in
reference [5] lays in $\bm{D}$. Without finite Larmor radius effects, $\bm{D}$
becomes a scalar thus it is possible to present $\bm{h}$ terms in a continued
fraction as in their model. In our model, re-writing as

$\bm{h}=\bm{M}^{-1}\cdot\bm{D}$ (19)

and noting that $\mathrm{d}^{2}\bm{M}/\mathrm{d}x^{2}=0$, we have

$\frac{\mathrm{d}^{2}\bm{h}}{\mathrm{d}x^{2}}=\bm{M}^{-1}\cdot\frac{\mathrm{d}^{2}\bm{D}}{\mathrm{d}x^{2}}-2\bm{M}^{-1}\cdot\frac{\mathrm{d}\bm{M}}{\mathrm{d}x}\cdot\bm{M}^{-1}\cdot\frac{\mathrm{d}\bm{D}}{\mathrm{d}x}+2\bm{M}^{-1}\cdot\frac{\mathrm{d}\bm{M}}{\mathrm{d}x}\cdot\bm{M}^{-1}\cdot\frac{\mathrm{d}\bm{M}}{\mathrm{d}x}\cdot\bm{M}^{-1}\cdot\bm{D}$
(20)

Substituting the Legendre series of $g(x,v,\xi)$ into the quasi-neutrality
equation (2) and expanding the Bessel function as before, we have

$n_{0}\frac{e\phi(x)}{T}=-2\pi\int_{0}^{\infty}v^{2}\mathrm{d}v\cdot\left[\left(1-\frac{k_{y}^{2}v^{2}}{6\omega_{ce}^{2}}\right)h_{0}+\frac{k_{y}^{2}v^{2}}{30\omega_{ce}^{2}}h_{2}+\frac{v^{2}}{6\omega_{ce}^{2}}\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}h_{0}-\frac{v^{2}}{30\omega_{ce}^{2}}\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}h_{2}\right]$
(21)

We define two row vectors as

$\displaystyle\bm{R_{\phi}}=\begin{pmatrix}1-\frac{k_{y}^{2}v^{2}}{6\omega_{ce}^{2}}&0&\frac{k_{y}^{2}v^{2}}{30\omega_{ce}^{2}}&0&0&\cdots\end{pmatrix}$
(22)
$\displaystyle\bm{Q_{\phi}}=\begin{pmatrix}\frac{v^{2}}{6\omega_{ce}^{2}}&0&-\frac{v^{2}}{30\omega_{ce}^{2}}&0&0&\cdots\end{pmatrix}$
(23)

Therefore the quasi-neutrality equation (21) becomes

$n_{0}\frac{e\phi(x)}{T}=-2\pi\int_{0}^{\infty}v^{2}\mathrm{d}v\cdot\left(\bm{R_{\phi}}\cdot\bm{h}+\bm{Q_{\phi}}\cdot\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}\bm{h}\right)$
(24)

To complete the set of electromagnetic equations, Ampere’s Law is also needed.
The parallel perturbed current density in Fourier space is

$\hat{j}_{e\|}=-\int_{-\infty}^{\infty}\mathrm{d}^{3}\bm{v}\cdot\hat{g}(k)\mathrm{J_{0}}(k_{\perp}\rho_{\perp})v_{\|}e$
(25)

Similarly, we expand the Bessel function, define another two row vectors as

$\displaystyle\bm{R_{A}}=\begin{pmatrix}0&\frac{2}{3}-\frac{k_{y}^{2}v^{2}}{15\omega_{ce}^{2}}&0&\frac{k_{y}^{2}v^{2}}{35\omega_{ce}^{2}}&0&0&\cdots\end{pmatrix}$
(26)
$\displaystyle\bm{Q_{A}}=\begin{pmatrix}0&\frac{v^{2}}{6\omega_{ce}^{2}}&0&-\frac{v^{2}}{35\omega_{ce}^{2}}&0&0&\cdots\end{pmatrix}$
(27)

and follow the same procedure above. The parallel perturbed current density in
real space becomes

$j_{e\|}=-\int_{0}^{\infty}v^{3}\mathrm{d}v\cdot\left(\bm{R_{A}}\cdot\bm{h}+\bm{Q_{A}}\cdot\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}\bm{h}\right)$
(28)

Therefore,

$\left(k_{y}^{2}-\frac{\partial^{2}}{\partial
x^{2}}\right)A_{\|}(x)=-\mu_{0}\int_{0}^{\infty}v^{3}\mathrm{d}v\cdot\left(\bm{R_{A}}\cdot\bm{h}+\bm{Q_{A}}\cdot\frac{\mathrm{d}^{2}}{\mathrm{d}x^{2}}\bm{h}\right)$
(29)

Substituting the expressions for $\bm{h}$ in equation (19) and
$\mathrm{d}^{2}\bm{h}/\mathrm{d}x^{2}$ in equation (20) into equations (24)
and (29), yields the final expressions incorporating both electromagnetic
effects and a Lorentz collision operator.

Simplifying the parameters in the drive term $\bm{D}$ and normalising the
results using the same parameters as in the electrostatic model above, as well
as $u=v/v_{e}$ and $\bar{A}_{\|}=A_{\|}/\rho_{e}B$. The normalised matrix
equations become

$\displaystyle\bar{\phi}-\frac{2}{\sqrt{\pi}}\int_{0}^{\infty}\mathrm{d}u\cdot
u^{2}\mathrm{e}^{-u^{2}}\left(\bar{\omega}-1-\eta\left(u^{2}-\frac{3}{2}\right)\right)\cdot\left(\bm{\bar{R_{n}}}\cdot\bm{\bar{h}}+\bm{\bar{Q_{n}}}\cdot\frac{\mathrm{d}^{2}}{\mathrm{d}\bar{x}^{2}}\bm{\bar{h}}\right)=0$
(30)
$\displaystyle\left(\bar{k_{y}}^{2}-\frac{\partial^{2}}{\partial\bar{x}^{2}}\right)\bar{A}_{\|}-\frac{2}{\sqrt{\pi}}\beta\int_{0}^{\infty}\mathrm{d}u\cdot
u^{3}\mathrm{e}^{-u^{2}}\left(\bar{\omega}-1-\eta\left(u^{2}-\frac{3}{2}\right)\right)\cdot\left(\bm{\bar{R_{j}}}\cdot\bm{\bar{h}}+\bm{\bar{Q_{j}}}\cdot\frac{\mathrm{d}^{2}}{\mathrm{d}\bar{x}^{2}}\bm{\bar{h}}\right)=0$
(31)

in which

$\displaystyle\bm{\bar{R_{n}}}=$
$\displaystyle\begin{pmatrix}1-\frac{\bar{k_{y}}^{2}u^{2}}{6}&0&\frac{\bar{k_{y}}^{2}u^{2}}{30}&0&0&0&\cdots\end{pmatrix}$
(32) $\displaystyle\bm{\bar{Q_{n}}}=$
$\displaystyle\begin{pmatrix}\frac{u^{2}}{6}&0&-\frac{u^{2}}{30}&0&0&0&\cdots\end{pmatrix}$
(33) $\displaystyle\bm{\bar{R_{j}}}=$
$\displaystyle\begin{pmatrix}0&\frac{2}{3}-\frac{\bar{k_{y}}^{2}u^{2}}{15}&0&\frac{\bar{k_{y}}^{2}u^{2}}{35}&0&0&\cdots\end{pmatrix}$
(34) $\displaystyle\bm{\bar{Q_{j}}}=$
$\displaystyle\begin{pmatrix}0&\frac{v^{2}}{6\omega_{ce}^{2}}&0&-\frac{v^{2}}{35\omega_{ce}^{2}}&0&0&\cdots\end{pmatrix}$
(35) $\displaystyle\bm{\bar{h}}=$
$\displaystyle\bm{\bar{M}}^{-1}\cdot\bm{\bar{D}}$ (36)
$\displaystyle\begin{split}\frac{\mathrm{d}^{2}}{\mathrm{d}\bar{x}^{2}}\bm{\bar{h}}=&\bm{\bar{M}}^{-1}\cdot\frac{\mathrm{d}^{2}}{\mathrm{d}\bar{x}^{2}}\bm{\bar{D}}-2\bm{\bar{M}}^{-1}\cdot\frac{\mathrm{d}}{\mathrm{d}\bar{x}}\bm{\bar{M}}\cdot\bm{\bar{M}}^{-1}\cdot\frac{\mathrm{d}}{\mathrm{d}\bar{x}}\bm{\bar{D}}\\\
&+2\bm{\bar{M}}^{-1}\cdot\frac{\mathrm{d}}{\mathrm{d}\bar{x}}\bm{\bar{M}}\cdot\bm{\bar{M}}^{-1}\cdot\frac{\mathrm{d}}{\mathrm{d}\bar{x}}\bm{\bar{M}}\cdot\bm{\bar{M}}^{-1}\cdot\bm{\bar{D}}\end{split}$
(37)

with

$\bm{\bar{M}}=\begin{pmatrix}\bar{\omega}&-\frac{2}{3}\epsilon\bar{x}u&0&0&0&\cdots\\\
-\frac{2}{3}\epsilon\bar{x}u&\frac{1}{3}(\bar{\omega}+i\bar{\nu})&-\frac{4}{15}\epsilon\bar{x}u&0&0&\cdots\\\
0&-\frac{4}{15}\epsilon\bar{x}u&\frac{1}{5}(\bar{\omega}+3i\bar{\nu})&-\frac{6}{35}\epsilon\bar{x}u&0&\cdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&\cdots&-\frac{2(n-1)}{(2n-3)(2n-1)}\epsilon\bar{x}u&\frac{1}{2n-1}(\bar{\omega}+\frac{n(n-1)}{2}i\bar{\nu})&-\frac{2n}{(2n-1)(2n+1)}\epsilon\bar{x}u&\cdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots\end{pmatrix}$ (38)

and

$\bm{\bar{D}}=\begin{pmatrix}\bar{\phi}+\frac{1}{6}u^{2}\frac{\mathrm{d}^{2}}{\mathrm{d}\bar{x}^{2}}\bar{\phi}\\\
-\frac{2}{3}u\bar{A}_{\|}-\frac{4}{15}u^{3}\frac{\mathrm{d}^{2}}{\mathrm{d}\bar{x}^{2}}\bar{A}_{\|}\\\
-\frac{1}{30}u^{2}\frac{\mathrm{d}^{2}}{\mathrm{d}\bar{x}^{2}}\bar{\phi}\\\
\frac{1}{35}u^{3}\frac{\mathrm{d}^{2}}{\mathrm{d}\bar{x}^{2}}\bar{A}_{\|}\\\
0\\\ 0\\\ \vdots\end{pmatrix}$ (39)

These represent an infinite tridiagonal matrix and a column vector,
respectively. In practice, they will be truncated, albeit at a large size
during the calculation. Equations (30) and (31) will eventually lead to a
system of two simultaneous second order differential equations for both
electrostatic potential $\phi$ and parallel magnetic potential $A_{\|}$.
Generally, however, the coefficients for each of the terms are not easy to
simplify and reveal the insight of physics, except in some special cases. For
the collisional microtearing theory where the finite Larmor radius effects are
not considered, the drive term $\bm{D}$ has non-zero elements only in the
first two terms; thus only the left top four elements in the inverse matrix of
$\bm{M}$ will contribute to the results. In this case, the matrix products can
be simplified to produce the continued fraction and electron parallel
conductivity in reference [5]. In another simple case, when collision
frequency $\nu=0$, at the centre plane of the slab $x=0$, the tridiagonal
matrix $\bm{M}$ becomes diagonal, and the calculation will be significantly
simplified. Dropping electromagnetic terms, the equation becomes equivalent to
the electrostatic model discussed in the previous section.

##  5 Numerical results and discussion

We have established two reduced models for the collisionless micro-scale
tearing instability considering finite Larmor radius effects from electrons.
Based on the electrostatic eigenmode equations (6), (7) and (8), and
electromagnetic eigenmode equations (30) and (31), we have developed two codes
to calculate the complex eigenmode frequency $\omega$ for each of these
models. The algorithm of the codes is an eigensolver based on an iteration
method described in Chapter 5.4 of reference [13]. Here, for our
electromagnetic equations with our chosen parameters, the complex frequency
$\omega$ is typically converged to a relative tolerance of $10^{-6}$ when the
matrices are of size 30 by 30. The boundary conditions for the tearing
instability for both models are spatially localised tearing parity
requirements

$|\phi|,|A_{\|}|\to 0\text{ as
}|x|\to\infty\quad,\quad\phi(0)=0\quad,\quad\frac{\mathrm{d}A_{\|}}{\mathrm{d}x}\bigg{|}_{x=0}=0$
(40)

Please note that we did not assume any parity in the derivation; thus our
model is also capable of looking for twisting parity solutions, whose boundary
conditions are

$|\phi|,|A_{\|}|\to 0\text{ as }|x|\to\infty\quad,\quad
A_{\|}(0)=0\quad,\quad\frac{\mathrm{d}\phi}{\mathrm{d}x}\bigg{|}_{x=0}=0$ (41)

For the collisionless limit $\nu=0$, solving either electrostatic or
electromagnetic eigenmode equations gives results close to those of GS2 (see
figure (8) for the point at $\nu=0$). We know that the electrostatic model
describes the electron temperature gradient (ETG) mode. However, the ETG mode
is usually considered to be a twisting parity mode. Nevertheless, as with any
eigenmode problem, there is a family of solutions (harmonics) with alternating
parity, in which the twisting parity is the fundamental harmonic. In fact,
eigenmodes of different harmonics can co-exist and there is no physical reason
why the fundamental one should be the most unstable. Indeed, such phenomena,
where the higher harmonics are more unstable, were found previously in both
ion temperature gradient (ITG) mode [14, 15] and ETG mode calculations [16,
17]; studies [18, 19] reported the existence of unstable high order ballooning
(twisting parity) modes and that parity transition can happen under certain
parameters; studies [20, 21] demonstrated the excitation of higher harmonics
and the parity mixture under certain scenarios. For our electrostatic model,
figure (7) shows different eigenvalues of which the even and odd harmonics are
twisting and tearing modes, respectively, solved with corresponding boundary
conditions. Our results show that the most unstable mode in this case is the
third order harmonic, which corresponds to the collisionless tearing parity
instability we have found.

Figure 7: Comparison of harmonics of eigenmode solution for GS2 and for
solutions of equations (6), (7) and (8). Here $\nu/\omega_{*e}=0.0$ and
$\eta_{e}=5.0$; other parameters are kept the same as in figure (2). Figure 8:
Comparison of our electromagnetic model (solid lines) with GS2 results
(triangle symbols) as a function of collision frequency in the collisionless
regime. Parameters are kept the same as in figure (2).

When including the collision frequency, the results from our electromagnetic
model are similar to those obtained from GS2, as shown in figure (8). Both
models show that the this mode is driven by electron temperature gradient,
consistent with the identification as an ETG. The mode growth rate decreases
as the collision frequency rises, but the growth rate in GS2 has a stronger
variation and switches to a different harmonic at $\nu=1.0$ for $\eta_{e}=5.0$
and $\nu=0.7$ for $\eta_{e}=7.0$, as indicated by the jump in frequency. Note
that our electromagnetic model results are consistent with the third order
harmonic in figure (7). One reason for the difference between our model and
GS2 might be that the collision operator in GS2 differs from our model [22,
23]. Though quantitatively slightly different, both our model and GS2 show
that this collisionless tearing parity instability tends to be stabilised by
collisions. To summarise, we conclude that the collisionless tearing parity
instability found here in slab geometry is a form of ETG instability, with a
different drive mechanism to the standard collisional slab MTM.

Figure 9: The structure of magnetic islands at $\nu/\omega_{*e}=0.0$ growing
from small amplitude (left) to large amplitude (right) calculated from our
electromagnetic model. Here $\eta_{e}=5.0$; other related parameters are kept
the same as in figure (2).

Although the underlying mechanism behind the collisionless tearing parity
instability studied here is different from the collisional MTM, it still leads
to magnetic reconnection and the formation of magnetic islands. Both GS2 and
our electromagnetic model provide the mode structure for the collisionless
instability as shown in figure (3). We can calculate the flux surfaces of the
magnetic field from the magnetic potential $A_{\|}$. The structure of magnetic
field lines is given by contours of the flux

$\psi(x,y)=\frac{B_{0}x^{2}}{2L_{s}}+\mathrm{Re}(A_{\|}\cdot\mathrm{e}^{ik_{y}y})$
(42)

The contour plot of constant levels of $\psi(x,y)$ gives the flux surfaces and
hence the magnetic structure. Figure (9) shows the island structure for the
collisionless mode when $\nu=0$. However, note that the amplitude of $A_{\|}$
is arbitrary in our linear model, so the width of the island is not
determined. From left to right, the three panels of figure (9) are examples to
show that under the same parameters the island shape can become more contorted
as its size grows from the order of electron Larmor radius $\rho_{e}$ to ion
Larmor radius $\rho_{i}$. It can also be found that as the island width grows,
a secondary island arises in the vicinity of the X-point. We believe that the
inflection points of $A_{\|}$ will finally provide a limit for the maximum
island width. How this will affect the particle and heat transport is to be
answered in future work.

##  6 Conclusion

We have shown that there is a collisionless micro-scale tearing parity
instability that can drive reconnection even in the absence of collisions. We
have established two models considering electron finite Larmor radius effects
to interpret the physics of this mode, which is shown to be stabilised when
the collision frequency increases. We identify the collisionless mode as a
tearing parity harmonic of the conventional slab ETG mode, which is the most
unstable harmonic for our parameters. The electromagnetic component results in
magnetic islands.

Our result stands as an example to show that tearing parity modes can arise
from a whole range of different drives, and there may be other possible ways
to get small scale tearing parity modes. These can have an impact on the
transport and can be very challenging to resolve numerically, posing problems
for attempts to simulate them. On the other hand, even if the tearing parity
eigenmodes are not the most unstable harmonic linearly, it still may be
possible that tearing harmonics can play a role non-linearly, leading to a
background degradation to the confining magnetic field everywhere that such
instabilities exist.

The remaining questions in our research include why and in what parameter
range does the tearing harmonic become the most unstable ETG mode in the
electrostatic model, and how does this collisionless mode behave in toroidal
geometry. The electromagnetic model (30) and (31) we obtain is a rather
complicated expression. Further simplification may reveal more physical
insight.

This work is funded by China Scholarship Council and the University of York.
It is also part of TDoTP project funded by EPSRC (EP/R034737/1).

## References

  * Wong et al. [2007] K. L. Wong, S. Kaye, D. R. Mikkelsen, J. A. Krommes, K. Hill, R. Bell, and B. LeBlanc.  Microtearing instabilities and electron transport in the nstx spherical tokamak.  _Phys. Rev. Lett._ , 99:135003, 2007.  doi:10.1103/PhysRevLett.99.135003.
  * Guttenfelder et al. [2011] W. Guttenfelder, J. Candy, S. M. Kaye, W. M. Nevins, E. Wang, R. E. Bell, G. W. Hammett, B. P. LeBlanc, D. R. Mikkelsen, and H. Yuh.  Electromagnetic transport from microtearing mode turbulence.  _Phys. Rev. Lett._ , 106:155004, 2011.  doi:10.1103/PhysRevLett.106.155004.
  * Doerk et al. [2011] H. Doerk, F. Jenko, M. J. Pueschel, and D. R. Hatch.  Gyrokinetic microtearing turbulence.  _Phys. Rev. Lett._ , 106:155003, 2011.  doi:10.1103/PhysRevLett.106.155003.
  * Drake and Lee [1977] J. F. Drake and Y. C. Lee.  Kinetic theory of tearing instabilities.  _Phys. Fluids_ , 20:1341, 1977.  doi:10.1063/1.862017.
  * Gladd et al. [1980] N. T. Gladd, J. F. Drake, C. L. Chang, and C. S. Liu.  Electron temperature gradient driven microtearing mode.  _Phys. Fluids_ , 23:1182, 1980.  doi:10.1063/1.863119.
  * Applegate et al. [2007] D. J. Applegate, C. M. Roach, J. W. Connor, S. C. Cowley, W. Dorland, R. J. Hastie, and N. Joiner.  Micro-tearing modes in the mega ampere spherical tokamak.  _Plasma Physics and Controlled Fusion_ , 49(8):1113–1128, 2007.  doi:10.1088/0741-3335/49/8/001.
  * Dickinson et al. [2013] D. Dickinson, C. M. Roach, S. Saarelma, R. Scannell, A. Kirk, and H. R. Wilson.  Microtearing modes at the top of the pedestal.  _Plasma Phys. Control. Fusion_ , 55:074006, 2013.  doi:10.1088/0741-3335/55/7/074006.
  * Swamy et al. [2014] Aditya K. Swamy, R. Ganesh, J. Chowdhury, S. Brunner, J. Vaclavik, and L. Villard.  Global gyrokinetic stability of collisionless microtearing modes in large aspect ratio tokamaks.  _Physics of Plasmas_ , 21(8):082513, 2014.  doi:10.1063/1.4893314.
  * Swamy et al. [2015] Aditya K. Swamy, R. Ganesh, S. Brunner, J. Vaclavik, and L. Villard.  Collisionless microtearing modes in hot tokamaks: Effect of trapped electrons.  _Physics of Plasmas_ , 22(7):072512, 2015.  doi:10.1063/1.4927579.
  * Predebon and Sattin [2013] I. Predebon and F. Sattin.  On the linear stability of collisionless microtearing modes.  _Physics of Plasmas_ , 20(4):040701, 2013.  doi:10.1063/1.4799980.
  * Kotschenreuther et al. [1995] M. Kotschenreuther, G. Rewoldt, and W. M. Tang.  Comparison of initial value and eigenvalue codes for kinetic toroidal plasma instabilities.  _Computer Physics Communications_ , 88:128, 1995.  doi:10.1016/0010-4655(95)00035-E.
  * Barnes et al. [2019] Michael Barnes, David Dickinson, William Dorland, Peter Alec Hill, Joseph Thomas Parker, Colin Malcolm Roach, Stephen Biggs-Fox, Ryusuke Numata, George Wilkie, Lucian Anton, Justin Ball, Jessica Baumgaertel, Greg Colyer, Michael Hardman, Joachim Hein, Edmund Highcock, Gregory Howes, Adrian Jackson, Michael T. Kotschenreuther, Jungpyo Lee, Huw Leggate, Noah Mandell, Adwiteey Mauriya, Tomo Tatsuno, and Ferdinand Van Wyk.  GS2 (v8.0.1).  2019\.  doi:10.5281/zenodo.2551067.
  * Dickinson [2012] D. Dickinson.  _Effects of profiles on microinstabilities in tokamaks_.  PhD thesis, University of York, 2012.
  * Gao et al. [2002] Zhe Gao, J. Q. Dong, G. J. Liu, and C. T. Ying.  Electromagnetic ion temperature gradient modes of tearing mode parity in high $\beta$ sheared slab plasmas.  _Physics of Plasmas_ , 9(5):1692–1697, 2002.  doi:10.1063/1.1471516.
  * Plunk et al. [2014] G. G. Plunk, P. Helander, P. Xanthopoulos, and J. W. Connor.  Collisionless microinstabilities in stellarators. iii. the ion-temperature-gradient mode.  _Physics of Plasmas_ , 21(3):032112, 2014.  doi:10.1063/1.4868412.
  * Lee et al. [1987] Y. C. Lee, J. Q. Dong, P. N. Guzdar, and C. S. Liu.  Collisionless electron temperature gradient instability.  _Phys. Fluids_ , 30:1331, 1987.  doi:10.1063/1.866248.
  * Zocco et al. [2015] A. Zocco, N. F. Loureiro, D. Dickinson, R. Numata, and C. M. Roach.  Kinetic microtearing modes and reconnecting modes in strongly magnetised slab plasmas.  _Plasma Physics and Controlled Fusion_ , 57(6):065008, 2015.  doi:10.1088/0741-3335/57/6/065008.
  * Xie and Xiao [2015] Hua-sheng Xie and Yong Xiao.  Unconventional ballooning structures for toroidal drift waves.  _Physics of Plasmas_ , 22(9):090703, 2015.  doi:10.1063/1.4931072.
  * Xie et al. [2018] Hua-Sheng Xie, Zhi-Xin Lu, and Bo Li.  Kinetic ballooning mode under steep gradient: High order eigenstates and mode structure parity transition.  _Physics of Plasmas_ , 25(7):072106, 2018.  doi:10.1063/1.5025949.
  * Pueschel et al. [2019] M. J. Pueschel, D. R. Hatch, D. R. Ernst, W. Guttenfelder, P. W. Terry, J. Citrin, and J. W. Connor.  On microinstabilities and turbulence in steep-gradient regions of fusion devices.  _Plasma Physics and Controlled Fusion_ , 61(3), 2019.  doi:10.1088/1361-6587/aaf8c1.
  * Ishizawa et al. [2019] A. Ishizawa, Y. Kishimoto, and Y. Nakamura.  Multi-scale interactions between turbulence and magnetic islands and parity mixture - a review.  _Plasma Physics and Controlled Fusion_ , 61(5), 2019.  doi:10.1088/1361-6587/ab06a8.
  * Abel et al. [2008] I. G. Abel, M. Barnes, S. C. Cowley, W. Dorland, and A. A. Schekochihin.  Linearized model fokker–planck collision operators for gyrokinetic simulations. i. theory.  _Phys. of Plasmas_ , 15:122509, 2008.  doi:10.1063/1.3046067.
  * Barnes et al. [2009] M. Barnes, I. G. Abel, W. Dorland, D. R. Ernst, G. W. Hammett, P. Ricci, B. N. Rogers, A. A. Schekochihin, and T. Tatsuno.  Linearized model fokker–planck collision operators for gyrokinetic simulations. ii. numerical implementation and tests.  _Phys. of Plasmas_ , 16:072107, 2009.  doi:10.1063/1.3155085.


# Two High-performance Schemes of Transmit Antenna Selection for Secure
Spatial Modulation

Feng Shu, Zhengwang Wang, Riqing Chen, Yongpeng Wu, and Jiangzhou Wang This
work was in part supported by the National Natural Science Foundation of China
(Grant Nos. 61771244, 61702258, 61472190, 61501238 and 61602245), the China
Postdoctoral Science Foundation (2016M591852), the Postdoctoral research
funding program of Jiangsu Province(1601257C), the Natural Science Foundation
of Jiangsu Province (Grants No. BK20150791), and the open research fund of
National Mobile Communications Research Laboratory, Southeast University,
China (No.2013D02).Feng Shu, Zhengwang Wang are with School of Electronic and
Optical Engineering, Nanjing University of Science and Technology, Nanjing,
210094, China.Riqing Chen, and Feng Shu are with the College of Computer and
Information Sciences, Fujian Agriculture and Forestry University, Fuzhou
350002, China.Yongpeng Wu is with the Shanghai Key Laboratory of Navigation
and Location-Based Services, Shanghai Jiao Tong University, Minhang
200240Jiangzhou Wang is with the School of Engineering and Digital Arts,
University of Kent, Canterbury CT2 7NT, U.K. E-mail<EMAIL_ADDRESS>
###### Abstract

In this paper, a secure spatial modulation (SM) system with artificial noise
(AN)-aided is investigated. To achieve higher secrecy rate (SR) in such a
system, two high-performance schemes of transmit antenna selection (TAS),
leakage-based and maximum secrecy rate (Max-SR), are proposed and a
generalized Euclidean distance-optimized antenna selection (EDAS) method is
designed. From simulation results and analysis, the four TAS schemes have an
decreasing order: Max-SR, leakage-based, generalized EDAS, and random
(conventional), in terms of SR performance. However, the proposed Max-SR
method requires the exhaustive search to achieve the optimal SR performance,
thus its complexity is extremely high as the number of antennas tends to
medium and large scale. The proposed leakage-based method approaches the Max-
SR method with much lower complexity. Thus, it achieves a good balance between
complexity and SR performance. In terms of bit error rate (BER), their
performances are in an increasing order: random, leakage-based, Max-SR, and
generalized EDAS.

###### Index Terms:

MIMO, secure spatial modulation, transmit antenna selection, signal-to-leakage
noise ratio, artificial noise.

##  I Introduction

Spatial modulation (sm) is emerging as a promising multiple-input-multiple-
output (MIMO) communication technology, which exploits both the index of
activated transmit antennas and amplitude phase modulation (APM) signals to
carry bits information[1, 2]. Due to low complexity of transmitter and
receiver, SM technology is applicable to energy-efficient scenarios. However,
in such a SM system, it is very possible that the confidential messages are
intercepted by unintend receivers due to the broadcast nature of wireless
channel. How to achieve a secure transmission in SM systems is becoming a hot
research topic in wireless networks.

Recently, physical layer security in MIMO systems have been widely
investigated, which exploit the uniqueness and time-varying characteristics of
channel to obtain secure transmission against eavesdropper[3, 4, 5, 6, 7].
Some literature focuses on the investigation of SM systems for secure
transmission [8, 9, 10, 11, 12, 13]. The authors in [8] derived the secrecy
mutual information of SM with finite alphabet input and proposed a precoding
scheme to improve security for SM or space shift keying (SSK) systems in[12].
Both [9] and [13] compressed artificial noise (AN) onto the null-space of the
desired channel to interfere unknown eavesdropper. Via transmitting AN at
transmitter[9] or full-duplex receiver[13] may achieve a high secrecy rate
(SR) approaching to the spectral efficiency (SE). In[10] and [11], the authors
generalized precoding-aided spatial modulation (PSM) to secrecy PSM with
constructing time-varying precoder[10] or optimizing the precoder by jointly
minimizing the receive power at eavesdropper while maximizing the receive
power at desired user[11]. PSM systems utilize the indices of receive antennas
to carry bits information, which reduces the complexity of receiver. However,
all transmit antennas are activated in PSM systems, which results in the
problems of inter-channel interference (ICI) and inter-antenna synchronization
(IAS).

On the other hand, transmit antenna selection (TAS) for improving the
performance of SM systems was firstly investigated in[14]. The authors of [14]
proposed two TAS methods: capacity-optimized antenna selection (COAS) and
Euclidean distance-optimized antenna selection (EDAS). The EDAS scheme showed
a better bit error rate (BER) performance and higher complexity than COAS
scheme. Therefore, some literature was intended to reduce the computational
and search complexity of the EDAS scheme[15, 16, 17, 18].

To the best of our knowledge, there is no research work concerning how to
design TAS methods for secure SM system. In this paper, we will focus on the
aspect research and make our effort to address this problem, our main
contributions are summarized as follows:

  1. 1.

To reach the bound of SR, a maximum secrecy rate (Max-SR) method is proposed
with exhaustive search. In accordance with simulation results, it can show the
best SR performance in a large range of SNR.

  2. 2.

To reduce the leakage of the confidential messages power for desired user to
eavesdropper, a leakage-based TAS method is proposed for secure SM system.
Also, its low complexity version, i.e., sorting-based solution, is presented
to achieve the same SR performance. Simulation results show that the proposed
leakage-based method provides a SR performance being close to that of the Max-
SR method with far lower-complexity.

  3. 3.

Finally, as a performance benchmark, we also generalize the conventional EDAS
method in non-secure SM system and make it become a secure EDAS for secure SM
system. Meanwhile, we analyze and compare the complexity for the three TAS
methods.

The remainder is organized as follows. In Section II, we describe secure SM
system model and give a definition for its average secrecy rate. Subsequently,
we propose two TAS methods: leakage-based and Max-SR, and generalize the
conventional EDAS scheme to secure SM system in Section III. In Section IV,
numerical simulation results are presented. Finally, we make our conclusions
in Section V.

Notation: throughout the paper, matrices, vectors, and scalars are denoted by
letters of bold upper case, bold lower case, and lower case, respectively.
Sign $(\cdot)^{-1}$, $(\cdot)^{H}$ denote inverse, conjugate transpose,
respectively. Notation $\mathbb{E}\\{\cdot\\}$ stands for the expectation
operation. Matrices $\textbf{I}_{N}$ denotes the $N\times N$ identity matrix,
and $\mathrm{tr}(\cdot)$ denotes matrix trace.

##  II System model

###  II-A Secure Spatial Modulation

Consider a typical secure SM system as shown in Fig. 1. In the system, there
are a transmitter (Alice) equipped with $N$ transmit antennas, a desired user
(Bob) with $N_{b}$ receive antennas, and an eavesdropper (Eve) with $N_{e}$
receive antennas, respectively. Without loss of generality, it is assumed that
$N$ is not a power of two, thus we have to select
$N_{t}={2^{\left\lfloor{\log_{2}^{{N}}}\right\rfloor}}$ out of $N$ transmit
antennas for mapping the bits to the antenna index. Notice that there are a
total $Q=\left(\begin{array}[]{l}N\\\ {N_{t}}\end{array}\right)$ patterns,
represented as
$\boldsymbol{\Omega}=\\{\boldsymbol{\Omega}_{1},\cdots,\boldsymbol{\Omega}_{Q}\\}$,
where $\boldsymbol{\Omega}_{k}$ denotes the antennas set of the $k$th pattern.
After Alice chooses one pattern from pattern set $\boldsymbol{\Omega}$ and
shares it with Bob through a low speed forward link, she activates one of
$N_{t}$ transmit antennas to emit $M$-ary APM symbol and uses the index of
activated antenna to convey partial bits information. As a result, the SE is
${\log_{2}}{MN_{t}}$ bits per channel use (bpcu).

Figure 1: System model of secure SM with TAS scheme

Referring to the secure SM system model in[9], the transmit baseband signal
with the aid of AN can be expressed as

$\displaystyle\mathbf{x}=\beta_{1}\sqrt{P_{S}}\mathbf{e}_{n}s_{m}+\beta_{2}\sqrt{P_{S}}\mathbf{P}_{AN}\mathbf{n}$
(1)

where $\mathbf{e}_{n}$ is the $n$th column of $\mathbf{I}_{N_{t}}$ for
$n\in\left\\{1,2,\cdots,N_{t}\right\\}$, and $s_{m}$ denotes the input symbol
equiprobably drawn from discrete $M$-ary APM constellation for
$n\in\left\\{1,2,\cdots,M\right\\}$. We normalize the input symbol power to
unit, i.e., $\mathbb{E}{|s_{m}|^{2}}=1$. In (1), matrix
$\mathbf{P}_{AN}\in\mathbb{C}^{N_{t}\times N_{t}}$ is the AN projection
matrix, $\mathbf{n}\sim\mathcal{CN}(0,\mathbf{I}_{N_{t}})$ is the random AN
vector, $P_{S}$ is the total transmit power, $\beta_{1}$, and $\beta_{2}$ are
the power allocation factors with $\beta_{1}^{2}+\beta_{2}^{2}=1$. Here, it is
particularly noted that all selected $N_{t}$ transmit antennas are used to
emit AN but only one of them emits APM symbol. Therefore, such secure SM
system alleviates the problems of ICI and IAS compared to conventional secure
MIMO systems.

Accordingly, the signals observed at Bob and Eve can be respectively
formulated as follows

$\mathbf{y}_{b}=\beta_{1}\sqrt{P_{S}}\mathbf{H}\mathbf{T}_{k}\mathbf{e}_{n}s_{m}+\beta_{2}\sqrt{P_{S}}\mathbf{H}\mathbf{T}_{k}\mathbf{P}_{AN}\mathbf{n}+\mathbf{n}_{b},$
(2)

and

$\mathbf{y}_{e}=\beta_{1}\sqrt{P_{S}}\mathbf{G}\mathbf{T}_{k}\mathbf{e}_{n}s_{m}+\beta_{2}\sqrt{P_{S}}\mathbf{G}\mathbf{T}_{k}\mathbf{P}_{AN}\mathbf{n}+\mathbf{n}_{e},$
(3)

where $\mathbf{H}\in\mathbb{C}^{N_{b}\times N}$ and
$\mathbf{G}\in\mathbb{C}^{N_{e}\times N}$ are the channel state information
(CSI) of desired and eavesdropping channels, respectively. Without loss of
generality, we assume $\mathbf{H}$ and $\mathbf{G}$ are frequency-flat
Rayleigh fading channels with each element obeying the distribution
$\mathcal{CN}(0,1)$ and the CSIs are available at Alice, which might be true
for active eavesdropper. $\mathbf{T}_{k}\in\mathbb{R}^{N\times
N_{t}},k=1,2,...,Q$, is the TAS matrix constituted by the specifically
selected $N_{t}$ columns determined by $\boldsymbol{\Omega}_{k}$ of
$\mathbf{I}_{N}$. Meanwhile
${\textbf{n}_{b}}\sim\mathcal{CN}(0,{\sigma_{b}^{2}}\mathbf{I}_{N_{b}})$ and
${\textbf{n}_{e}}\sim\mathcal{CN}(0,{\sigma_{e}^{2}}\mathbf{I}_{N_{e}})$
denote the complex additive white gaussian noise (AWGN) vectors at Bob and
Eve, respectively.

Here, we project the AN lying in the null-space of desired channel, therefore
the $\mathbf{P}_{AN}$ is given by

$\mathbf{P}_{AN}=\frac{1}{\mu}\left[{{{\mathbf{I}_{{N_{t}}}}-{\mathbf{H}_{k}}^{H}{{\left({{\mathbf{H}_{k}}{\mathbf{H}_{k}}^{H}}\right)}^{-1}}{\mathbf{H}_{k}}}}\right]$
(4)

where
$\mu={\parallel{{{\mathbf{I}_{{N_{t}}}}-{\mathbf{H}_{k}}^{H}{{\left({{\mathbf{H}_{k}}{\mathbf{H}_{k}}^{H}}\right)}^{-1}}{\mathbf{H}_{k}}}}\parallel_{\text{F}}}$
is the normalized factor of making the identity
$\mathrm{tr}(\mathbf{P}_{AN}^{H}\mathbf{P}_{AN})=1$ holds, and
$\mathbf{H}_{k}=\mathbf{H}\mathbf{T}_{k}$. It is obvious that
$\mathbf{H}_{k}\mathbf{P}_{AN}=\mathbf{0}$, this means the AN has no impact on
Bob.

Assuming Bob has the knowledge of $\mathbf{H}$ and $\mathbf{T}_{k}$, the
maximum-likelihood detector (MLD) is given by[2]

$[\hat{n},\hat{m}]=\mathop{\arg\min}\limits_{n\in[1,{N_{t}}],m\in[1,M]}\parallel{\mathbf{y}_{b}-\beta_{1}\sqrt{P_{S}}\mathbf{H}{\mathbf{T}_{k}}{\mathbf{e}_{n}}{s_{m}}}\parallel^{2}$
(5)

###  II-B Average Secrecy Rate and Problem Formulation

Here, we characterize the security by evaluating average secrecy rate, which
is formulated as

$\bar{R}_{s}={\mathbb{E}_{\mathbf{H},\mathbf{G}}}{(R_{s})}$ (6)

with

$R_{s}=\left[{I\left({\mathbf{x};{\mathbf{y}_{b}}\left|\mathbf{H},\mathbf{T}_{k}\right.}\right)-I\left({\mathbf{x};{\mathbf{y}_{e}}\left|\mathbf{G},\mathbf{T}_{k}\right.}\right)}\right]^{+}$
(7)

where $R_{s}$ denotes the definition of SR in[3].
$I\left({\mathbf{x};{\mathbf{y}_{b}}\left|\mathbf{H},\mathbf{T}_{k}\right.}\right)$
and
$I\left({\mathbf{x};{\mathbf{y}_{e}}\left|\mathbf{G},\mathbf{T}_{k}\right.}\right)$
are the mutual information over desired and eavesdropping channels,
respectively. According to [9], they are formulated as follows

$\displaystyle
I\left({\mathbf{x};{\mathbf{y}_{b}}\left|\mathbf{H},\mathbf{T}_{k}\right.}\right)={\log_{2}}^{M{N_{t}}}-\frac{1}{{M{N_{t}}}}\sum\limits_{n=1}^{N_{t}}\sum\limits_{m=1}^{M}\mathbb{E}_{\mathbf{n}_{b}}$
(8)
$\displaystyle{{{\log}_{2}}\sum\limits_{n^{\prime}=1}^{N_{t}}\sum\limits_{m^{\prime}=1}^{M}{\exp\left({\frac{{{{\parallel{{\mathbf{n}_{b}}}\parallel}^{2}}-{{\parallel{\boldsymbol{\delta}_{k,n,m}^{n^{\prime},m^{\prime}}+{\mathbf{n}_{b}}}\parallel}^{2}}}}{{\sigma_{b}^{2}}}}\right)}},$

and

$\displaystyle
I\left({\mathbf{x};{\mathbf{y}_{e}}\left|\mathbf{G},\mathbf{T}_{k}\right.}\right)={\log_{2}}^{M{N_{t}}}-\frac{1}{{M{N_{t}}}}\sum\limits_{n=1}^{N_{t}}\sum\limits_{m=1}^{M}\mathbb{E}_{\mathbf{n}_{e}^{{}^{\prime}}}$
(9)
$\displaystyle{{{\log}_{2}}\sum\limits_{n^{\prime}=1}^{N_{t}}\sum\limits_{m^{\prime}=1}^{M}{\exp\left({{{{\parallel{{\mathbf{n}_{e}^{{}^{\prime}}}}\parallel}^{2}}-{{\parallel{\mathbf{W}^{-1/2}\boldsymbol{\alpha}_{k,n,m}^{n^{\prime},m^{\prime}}+{\mathbf{n}_{e}^{{}^{\prime}}}}\parallel}^{2}}}}\right)}}$

where $\boldsymbol{\delta}_{n,m}^{n^{\prime},m^{\prime}}$,
$\boldsymbol{\alpha}_{n,m}^{n^{\prime},m^{\prime}}$ and
$\mathbf{n}_{e}^{{}^{\prime}}$ are given by

$\boldsymbol{\delta}_{k,n,m}^{n^{\prime},m^{\prime}}=\beta_{1}\sqrt{P_{S}}\mathbf{H}\mathbf{T}_{k}(\mathbf{e}_{n}s_{m}-\mathbf{e}_{n^{\prime}}s_{m^{\prime}}),$
(10)
$\boldsymbol{\alpha}_{k,n,m}^{n^{\prime},m^{\prime}}=\beta_{1}\sqrt{P_{S}}\mathbf{G}\mathbf{T}_{k}(\mathbf{e}_{n}s_{m}-\mathbf{e}_{n^{\prime}}s_{m^{\prime}}),$
(11)

and

$\mathbf{n}_{e}^{{}^{\prime}}=\mathbf{W}^{-1/2}(\beta_{2}\sqrt{P_{S}}\mathbf{G}\mathbf{T}_{k}\mathbf{P}_{AN}\mathbf{n}+\mathbf{n}_{e})$
(12)

where $\mathbf{n}_{e}^{{}^{\prime}}\sim\mathcal{CN}(0,\mathbf{I}_{N_{e}})$ is
the AWGN vector, and $\mathbf{W}$ is the covariance matrix of interference
plus noise at Eve

$\mathbf{W}=\beta_{2}^{2}P_{S}\mathbf{G}\mathbf{T}_{k}\mathbf{P}_{AN}\mathbf{P}_{AN}^{H}\mathbf{T}_{k}^{H}\mathbf{G}^{H}+\sigma_{e}^{2}\mathbf{I}_{N_{e}}$
(13)

Finally, our objective is to maximize the SR by TAS scheme, which is written
as the following optimization problem

$\displaystyle\mathop{\max}~{}~{}~{}~{}~{}~{}~{}{R_{s}}$ (14)
$\displaystyle\text{subject~{}to}~{}~{}\mathbf{T}_{k}\in\\{\mathbf{T}_{1},\mathbf{T}_{2},...,\mathbf{T}_{Q}\\}$

which is an integer or binary optimization problem and is NP-hard. In the next
section, we will present three TAS schemes for secure SM system.

##  III Proposed Transmit Antenna Selection Methods

In this section, we propose two new TAS schemes: leakage-based and Max-SR, and
generalize the conventional EDAS method to the secure SM scenario. Then, we
analyze the complexity for the three TAS methods.

###  III-A Proposed Leakage-Based Antenna Selection Method

####  III-A1 Problem Formulation

Similar to multi-user MIMO system in[19], we view the receive power of
confidential messages at Eve as the so-called leakage. Thus, the signal-to-
leakage-and-noise ratio (SLNR) for the $n$th channel of the $k$th pattern is
written as

$SLNR_{n}(\mathbf{T}_{k})=\frac{{\beta_{1}^{2}P_{S}\parallel{\mathbf{H}{\mathbf{T}_{k}}{\mathbf{e}_{n}}\parallel^{2}}}}{{\beta_{1}^{2}P_{S}\parallel{\mathbf{G}{\mathbf{T}_{k}}{\mathbf{e}_{n}}\parallel^{2}}+N_{b}\sigma_{b}^{2}}}$
(15)

It is assumed that all $N_{t}$ TAs of the selected pattern are activated with
equiprobability to transmit confidential messages, then the optimization
problem of maximizing SLNR (Max-SLNR) is formulated as

$\displaystyle\mathop{\max}~{}~{}~{}~{}~{}~{}\sum\limits_{n=1}^{{N_{t}}}{SLNR_{n}}(\mathbf{T}_{k})$
(16)
$\displaystyle\text{subject~{}to}~{}~{}\mathbf{T}_{k}\in\\{\mathbf{T}_{1},\mathbf{T}_{2},...,\mathbf{T}_{Q}\\}$

with $\mathbf{T}_{k}$ being the optimization variable. However, the above
optimization problem can be solved by exhaustive search with complexity
$\mathcal{O}\left(QN_{t}\right)$ floating-point operations (FLOPs). To lower
its complexity, the low-complexity implementation is necessary.

####  III-A2 Low-complexity sorting-based solution

Considering all transmit antennas are uncorrelated, that is, each transmit
antenna results in different SLNR, the optimization problem in (16) is
equivalent to choose the largest $N_{t}$ values from all $N$ SLNRs. Once we
calculate the SLNR values per transmit antenna, we simply arrange them in
descending order by sorting method as follow

$\underbrace{SLN{R_{{\pi_{\rm{1}}}}}\geq SLN{R_{{\pi_{\rm{2}}}}}\geq\cdots\geq
SLN{R_{{\pi_{N_{t}}}}}}_{{N_{t}}~{}{\text{selected~{}antennas}}}\geq\cdots\geq
SLN{R_{{\pi_{N}}}}$ (17)

where $\\{{\pi_{1}},{\pi_{2}},...,{\pi_{N}}\\}$ is an ordered permutation of
$\\{1,2,...,N\\}$ and $SLNR_{\pi_{n}}$ is written as

$SLNR_{\pi_{n}}=\frac{{\beta_{1}^{2}P_{S}\mathrm{tr}\left({\mathbf{h}_{\pi_{n}}}{\mathbf{h}_{\pi_{n}}}^{H}\right)}}{{\beta_{1}^{2}P_{S}\mathrm{tr}\left({\mathbf{g}_{\pi_{n}}}{\mathbf{g}_{\pi_{n}}}^{H}\right)+N_{b}\sigma_{b}^{2}}}$
(18)

where ${\mathbf{h}_{\pi_{n}}}$ and ${\mathbf{g}_{\pi_{n}}}$ are the
$\pi_{n}$th column of $\mathbf{H}$ and $\mathbf{G}$, respectively.

###  III-B Proposed Max-SR Antenna Selection Method

In the previous subsection, it is hard for the proposed leakage-based method
to maximize SR. In what follows, we will propose the Max-SR method, which aims
to maximize SR. For a specific channel realization and the $k$th selected
pattern, the transmission rates of desired and eavesdropping channels are
respectively bounded as[14]

$\alpha\leq{R_{b}}\leq\alpha+{\log_{2}}{N_{t}},~{}\beta\leq{R_{e}}\leq\beta+{\log_{2}}{N_{t}}$
(19)

where $\alpha$ and $\beta$ are given by

$\begin{array}[]{l}\alpha=\frac{1}{{{N_{t}}}}\sum\limits_{n=1}^{{N_{t}}}{\log_{2}}{\left(1+\frac{{{{\beta_{1}^{2}P_{S}\parallel{\mathbf{H}{\mathbf{T}_{k}}{\mathbf{e}_{n}}s_{m}}\parallel}^{\rm{2}}}}}{{{N_{b}}\sigma_{b}^{2}}}\right)}\\\
\beta=\frac{1}{{{N_{t}}}}\sum\limits_{n=1}^{{N_{t}}}{\log_{2}}\left({1+\frac{{{{\beta_{1}^{2}P_{S}\parallel{\mathbf{G}{\mathbf{T}_{k}}{\mathbf{e}_{n}}s_{m}}\parallel}^{\rm{2}}}}}{{\beta_{2}^{2}P_{S}{{\parallel{\mathbf{G}{\mathbf{T}_{k}}{\mathbf{P}_{AN}}\mathbf{n}}\parallel}^{2}}+{N_{e}}\sigma_{e}^{2}}}}\right)\end{array}$
(20)

Then, the corresponding SR is bounded as

$\alpha-\beta-{\log_{2}}{N_{t}}\leq{R_{s}}\leq\alpha-\beta+{\log_{2}}{N_{t}}$
(21)

Notice that $\alpha$ and $\beta$ are the functions of $\mathbf{T}_{k}$, thus
we can carefully select a pattern for maximizing $\alpha-\beta$, which can be
formulated as

$\displaystyle\max~{}~{}~{}~{}\sum\limits_{n=1}^{{N_{t}}}{\log_{2}}{\left({\frac{1+\frac{{{{\beta_{1}^{2}P_{S}\parallel{\mathbf{H}{\mathbf{T}_{k}}{\mathbf{e}_{n}}s_{m}}\parallel}^{\rm{2}}}}}{{{N_{b}}\sigma_{b}^{2}}}}{1+\frac{{\beta_{1}^{2}P_{S}\parallel{\mathbf{G}{\mathbf{T}_{k}}{\mathbf{e}_{n}}s_{m}}\parallel}^{\rm{2}}}{{\beta_{2}^{2}P_{S}{{\parallel{\mathbf{G}{\mathbf{T}_{k}}{\mathbf{P}_{AN}}\mathbf{n}}\parallel}^{2}}+{N_{e}}{\sigma_{e}^{2}}}}}}\right)}$
(22)
$\displaystyle\text{subject~{}to}~{}~{}\mathbf{T}_{k}\in\\{\mathbf{T}_{1},\mathbf{T}_{2},...,\mathbf{T}_{Q}\\}$

In general, the optimal TAS pattern for (22) is usually obtained by an
exhaustive search approach due to the fact that the $\mathbf{P}_{AN}$ in (22)
is associated with the selected pattern.

###  III-C Generalized Euclidean Distance Antenna Selection Method

Now, we develop the generalized EDAS method. It is well-known that the mutual
information lacks closed-form expression for discrete-input continuous-output
memoryless channels (DCMC). Here, we present a lower bound for (8)

$\displaystyle
I\left({\mathbf{x};{\mathbf{y}_{b}}\left|\mathbf{H},\mathbf{T}_{k}\right.}\right)_{LB}=\log_{2}^{M{N_{t}}}+{N_{b}}\left({1-\frac{1}{{{{\ln}2}}}}\right)-\frac{1}{{M{N_{t}}}}$
(23)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\sum\limits_{n=1}^{N_{t}}\sum\limits_{m=1}^{M}{{{\log}_{2}}\sum\limits_{n=1}^{N_{t}}\sum\limits_{m=1}^{M}{\exp\left({\frac{{-\parallel\boldsymbol{\delta}_{n,m}^{n^{\prime},m^{\prime}}\parallel^{2}}}{{2\sigma_{b}^{2}}}}\right)}},$

which is similar to that in [20]. Note that the lower bound is approximately
tight validated by Monte-Carlo simulation, so we may use it as an effective
design metric for TAS scheme.

Notice that the term $\boldsymbol{\delta}_{k,n,m}^{n^{\prime},m^{\prime}}$ in
(23) has an important impact on the lower bound and the minimum Euclidean
distance, i.e.,
$d_{\text{min}}=\mathop{\min}\limits_{(n,m)\neq(n^{\prime},m^{\prime})}\parallel\boldsymbol{\delta}_{k,n,m}^{n^{\prime},m^{\prime}}\parallel^{2}$,
is the dominant term in high signal-to-noise ratio (SNR) region[20]. In other
words, the bigger $d_{\text{min}}$ is expected to yield a higher mutual
information. Therefore, we aim to select a TAS pattern of maximizing the
minimum Euclidean distance over desired channel or minimizing the minimum
Euclidean distance over eavesdropping channel[12], which are formulated as
follows

${\mathbf{T}_{{k^{*}}}}=\mathop{\arg\max}\limits_{\mathbf{T}_{k}\in\\{\mathbf{T}_{1},\mathbf{T}_{2},...,\mathbf{T}_{Q}\\}}~{}~{}\mathop{\min}\limits_{(n,m)\neq(n^{\prime},m^{\prime})}\parallel\boldsymbol{\delta}_{k,n,m}^{n^{\prime},m^{\prime}}\parallel^{2}$
(24)
${\mathbf{T}_{{k^{*}}}}=\mathop{\arg\min}\limits_{\mathbf{T}_{k}\in\\{\mathbf{T}_{1},\mathbf{T}_{2},...,\mathbf{T}_{Q}\\}}~{}~{}\mathop{\min}\limits_{(n,m)\neq(n^{\prime},m^{\prime})}\parallel\boldsymbol{\alpha}_{k,n,m}^{n^{\prime},m^{\prime}}\parallel^{2}$
(25)

Here, we obtain two TAS patterns by (24) and (25), then we select the one
which results in higher $R_{s}$. Note that some low complexity methods for
efficiently solving the problems in (24) and (25) can be found in[15, 16, 18].

###  III-D Complexity Analysis and Comparison

Here, we make a complexity comparison concerning the above three methods. It
is obvious that the proposed leakage-based method has the lowest complexity of
$N(2N_{b}+2N_{e})+\mathcal{O}(N^{2})\approx\mathcal{O}(N^{2})$ FLOPs, where
the term $\mathcal{O}(N^{2})$ is the search complexity of sorting operation.
For the proposed Max-SR method, the complexity is about
$2QN_{t}[N_{t}^{2}(N_{b}+N_{e})+N_{t}(2N_{b}^{2}+N_{e})+2N_{e}+N_{b}]\approx\mathcal{O}(QN_{t}^{3})$
FLOPs. The complexity of generalized EDAS method is made up of three
terms[18]: (a) the computational complexity of the upper triangle matrix of
minimum Euclidean distance, (b) the search complexity for optimal pattern, and
(c) the complexity of computing secrecy rate for two patterns. The
complexities of terms (a) and (b) are given by
$\mathcal{O}_{a}=16M({{N^{2}}-N})\left({\frac{{{N_{b}}+{N_{e}}-2}}{3}}\right)+N\left[{2\left({{N_{b}}+{N_{e}}}\right)-1}\right]$
and
$\mathcal{O}_{b}=\left({{N_{t}}^{2}+3{N_{t}}+3}\right){2^{N-{N_{t}}}}-\frac{{N+5N+8}}{2}$
in [18]. The complexity of term (c) is as follows
${\rm{2}}{M^{\rm{2}}}{N_{t}}^{2}{N_{noise}}[{2(N_{t}+1)({N_{b}}+{N_{e}})+{N_{e}}^{2}}]+2Q[N_{t}^{2}(N_{e}+N_{b})+2N_{t}(N_{b}^{2}+N_{e}^{2})]$,
where $N_{noise}$ is the number of noise samples for evaluating the secrecy
rate. In summary, it is evident that their complexities are in increasing
order: leakage-based, Max-SR, and generalized EDAS.

##  IV simulation and numerical results

In this section, numerical simulation results are presented to analyze and
compare the performance of three TAS methods from different aspects: average
secrecy rate, cumulative density function (CDF) of secrecy rate and bit error
rate (BER) of Bob. Simulation parameters are set as follows: $N_{b}=N_{e}=2$,
$N=15$, $N_{t}=8$, and $\sigma_{b}^{2}=\sigma_{e}^{2}$. The SNR is defined as
$\frac{\beta_{1}^{2}P_{S}}{\sigma_{b}^{2}}$, and Quadrature Phase Shift Keying
(QPSK) modulation is employed.

Figure 2: Average secrecy rate of four TAS methods for $\beta_{2}=0.25$,
$\beta_{2}=0$
Figure 3: CDF curves of secrecy rate for SNR=-5dB, 0dB, 5dB and
$\beta_{2}=0.25$
Figure 4: Average BER curves of Bob for four TAS methods

Fig. 2 demonstrates the average SR for the three TAS methods described in
Section-III with $\beta^{2}_{2}=0.25$ and $\beta^{2}_{2}=0$, where the random
TAS method is used for performance reference. For $\beta^{2}_{2}=0.25$, the SR
performance of the proposed leakage-based method is very close to the Max-SR
method and better than the generalized EDAS method and random method. When we
adjust the value of $\beta^{2}_{2}$ to zero, we find the SR performance of the
proposed leakage-based method is almost identical to that of Max-SR and much
better than that of generalized EDAS method. It is also seen that random
method fails to achieve positive SR.

Fig. 3 shows the CDF curves of SR for four TAS methods with three different
values of SNRs: -5dB, 0dB, and 5dB. It is seen from Fig. 3 that the SR
performance of the proposed leakage-based method is between the proposed Max-
SR method and the generalized EDAS method. With increasing in the value of
SNR, the performance of generalized EDAS method approaches the two proposed
methods due to the fact that the minimum Euclidean distance dominates the
mutual information term in the high SNR region in accordance with (23). As SNR
increases, the CDF curves of SR moves to the right-hand side. That is, the
probability that SR falls in the large value interval grows gradually.
However, for all three SNR values, the four TAS methods show the same SR
performance trend: Max-SR $\geq$ leakage-based $\geq$ generalized EDAS $\geq$
random.

Fig. 4 illustrates the curves of average BER versus SNR for the three TAS
methods described in Section-III with $\beta_{2}^{2}=0.25$, where the random
method is used for performance reference. From Fig. 4, it is seen that the
generalized EDAS method exceeds the remaining three methods in terms of BER
performance. This is due to the fact the minimum Euclidean distance makes an
important impact on pair-error probability (PEP) in high SNR region. The BER
performances of the proposed leakage-based and Max-SR are worse than that of
generalized EDAS and better than that of random method.

##  V conclusion

In this paper, we have made an extensive investigation of TAS methods in
secure SM systems. Then, two high-performance TAS schemes: leakage-based and
Max-SR, have been proposed to improve the SR performance, and the EDAS method
has been generalized to provide a secure transmission. From simulation results
and complexity analysis, the proposed Max-SR is the optimal TAS scheme among
the three TAS methods in terms of SR performance while the proposed leakage-
based method achieves a SR performance near to the proposed Max-SR method with
far low complexity. By means of BER performance, the generalized EDAS methods
substantially outperforms the remaining two methods due to its main goal to
maximize minimum Euclidean distance, which makes a direct improvement on BER
performance.

## References

  * [1] R. Y. Mesleh, H. Haas, S. Sinanovic, W. A. Chang, and S. Yun, “Spatial modulation,” _IEEE Trans. Veh. Technol._ , vol. 57, no. 4, pp. 2228–2241, 2008.
  * [2] J. Jeganathan, A. Ghrayeb, and L. Szczecinski, “Spatial modulation: optimal detection and performance analysis,” _IEEE Commun. Lett._ , vol. 12, no. 8, pp. 545–547, 2008.
  * [3] A. D. Wyner, “The wire-tap channel,” _Bell System Tech. J._ , vol. 54, no. 8, pp. 1355–1387, 1975.
  * [4] J. Hu, F. Shu, and J. Li, “Robust synthesis method for secure directional modulation with imperfect direction angle,” _IEEE Commun. Lett._ , vol. 20, no. 6, pp. 1084–1087, 2016.
  * [5] F. Shu, X. Wu, J. Li, R. Chen, and B. Vucetic, “Robust synthesis scheme for secure multi-beam directional modulation in broadcasting systems,” _IEEE Access_ , vol. 4, no. 99, pp. 6614–6623, 2016.
  * [6] N. Zhao, F. R. Yu, M. Li, Q. Yan, and V. C. M. Leung, “Physical layer security issues in interference-alignment-based wireless networks,” _IEEE Commun. Mag._ , vol. 54, no. 8, pp. 162–168, 2016.
  * [7] H. M. Wang, Q. Yin, and X. G. Xia, “Distributed beamforming for physical-layer security of two-way relay networks,” _IEEE Trans. Signal Process._ , vol. 60, no. 7, pp. 3532–3545, 2012.
  * [8] X. Guan, Y. Cai, and W. Yang, “On the secrecy mutual information of spatial modulation with finite alphabet,” _in Proc. International Conference on Wireless Communications and Signal Processing_ , pp. 1–4, 2013.
  * [9] L. Wang, S. Bashar, Y. Wei, and R. Li, “Secrecy enhancement analysis against unknown eavesdropping in spatial modulation,” _IEEE Commun. Lett._ , vol. 19, no. 8, pp. 1351–1354, 2015.
  * [10] F. Wu, L. L. Yang, W. Wang, and Z. Kong, “Secret precoding-aided spatial modulation,” _IEEE Commun. Lett._ , vol. 19, no. Science, pp. 1544–1547, 2015.
  * [11] F. Wu, R. Zhang, L. L. Yang, and W. Wang, “Transmitter precoding-aided spatial modulation for secrecy communications,” _IEEE Trans. Veh. Technol._ , vol. 65, no. 1, pp. 467–471, 2016.
  * [12] S. R. Aghdam and T. M. Duman, “Physical layer security for space shift keying transmission with precoding,” _IEEE Wireless Commun. Lett._ , vol. 5, no. 2, pp. 180–183, 2016.
  * [13] C. Liu, L. L. Yang, and W. Wang, “Secure spatial modulation with a full-duplex receiver,” _IEEE Wireless Commun. Lett._ , vol. PP, no. 99, pp. 1–1, 2017\.
  * [14] R. Rajashekar, K. V. S. Hari, and L. Hanzo, “Antenna selection in spatial modulation systems,” _IEEE Commun. Lett._ , vol. 17, no. 3, pp. 521–524, 2013.
  * [15] K. Ntontin, M. D. Renzo, A. I. Perez-Neira, and C. Verikoukis, “A low-complexity method for antenna selection in spatial modulation systems,” _IEEE Commun. Lett._ , vol. 17, no. 12, pp. 2312–2315, 2013.
  * [16] J. Zheng and J. Chen, “Further complexity reduction for antenna selection in spatial modulation systems,” _IEEE Commun. Lett._ , vol. 19, no. 6, pp. 937–940, 2015.
  * [17] P. Yang, Y. Xiao, Y. L. Guan, S. Li, and L. Hanzo, “Transmit antenna selection for multiple-input multiple-output spatial modulation systems,” _IEEE Trans. Commun._ , vol. 64, no. 5, pp. 2035–2048, 2016.
  * [18] Z. Sun, Y. Xiao, P. Yang, S. Li, and W. Xiang, “Transmit antenna selection schemes for spatial modulation systems: Search complexity reduction and large-scale mimo applications,” _IEEE Trans. Veh. Technol._ , vol. 66, no. 9, pp. 8010–8021, 2017.
  * [19] M. Sadek, A. Tarighat, and A. H. Sayed, “A leakage-based precoding scheme for downlink multi-user mimo channels,” _IEEE Trans. Wireless Commun._ , vol. 6, no. 5, pp. 1711–1721, 2007.
  * [20] X. Guan, Y. Cai, and W. Yang, “On the mutual information and precoding for spatial modulation with finite alphabet,” _IEEE Wireless Commun. Lett._ , vol. 2, no. 4, pp. 383–386, 2013.