* [11] J. Fröhlich and A. Pizzo, _Lie-Schwinger block-diagonalization and gapped quantum chains_ , Commun. Math. Phys. 375 (2020), 2039–2069, arXiv:1812.02457.
  * [12] David Gosset and Evgeny Mozgunov, _Local gap threshold for frustration-free spin systems_ , J. Math. Phys. 57 (2016), 091901.
  * [13] M. B. Hastings, _Lieb-Schultz-Mattis in higher dimensions_ , Phys. Rev. B 69 (2004), 104431.
  * [14] M. B. Hastings, _The stability of free Fermi Hamiltonians_ , J. Math. Phys. 60 (2019), 042201, arXiv:1706.02270.
  * [15] M. B. Hastings and X. G. Wen, _Quasi-adiabatic continuation of quantum states: The stability of topological ground-state degeneracy and emergent gauge invariance_ , Phys. Rev. B 72 (2005), 045141.
  * [16] J. Henheik and S. Teufel, _Adiabatic theorem in the thermodynamic limit. Part II: Systems with a gap in the bulk_ , arXiv:2012.15239, 2020.
  * [17] M. J. Kastoryano and A. Lucia, _Divide and conquer method for proving gaps of frustration free Hamiltonians_ , J. Stat. Mech.: Theor. Exp. (2018), 033105\.
  * [18] T. Kato, _Perturbation theory for linear operators_ , Springer-Verlag, 1995\.
  * [19] M. Lemm and E. Mozgunov, _Spectral gaps of frustration-free spin systems with boundary_ , J. Math. Phys. 60 (2019), 051901, arXiv:1801.08915.
  * [20] E.H. Lieb and D.W. Robinson, _The finite group velocity of quantum spin systems_ , Commun. Math. Phys. 28 (1972), 251–257.
  * [21] S. Michalakis and J.P. Zwolak, _Stability of frustration-free Hamiltonians_ , Commun. Math. Phys. 322 (2013), 277–302.
  * [22] A. Moon and Y. Ogata, _Automorphic equivalence within gapped phases in the bulk_ , J. Funct. Analysis 278 (2020), 108422, arXiv:1906.05479.
  * [23] B. Nachtergaele, Y. Ogata, and R. Sims, _Propagation of correlations in quantum lattice systems_ , J. Stat. Phys. 124 (2006), 1–13.
  * [24] B. Nachtergaele, V. Scholz, and R. F. Werner, _Local approximation of observables and commutator bounds_ , Operator Methods in Mathematical Physics (J. Janas, P. Kurasov, A. Laptev, and S. Naboko, eds.), Operator Theory: Advances and Applications, vol. 227, Birkäuser Verlag, 2013, pp. 143–149.
  * [25] B. Nachtergaele, R. Sims, and A. Young, _Quasi-locality bounds for quantum lattice systems and perturbations of gapped ground states. Part II. perturbations of frustration-free spin models with gapped ground states_ , arxiv:2010.15337.
  * [26] by same author, _Quasi-locality bounds for quantum lattice systems. I. Lieb-Robinson bounds, quasi-local maps, and spectral flow automorphisms_ , J. Math. Phys. 60 (2019), 061101.
  * [27] Y. Ogata, _A $\mathbb{Z}_{2}$-index of symmetry protected topological phases with reflection symmetry for quantum spin chains_, arXiv:1904.01669, 2019\.
  * [28] by same author, _A $\mathbb{Z}_{2}$-index of symmetry protected topological phases with time reversal symmetry for quantum spin chains_, Commun. Math. Phys. 374 (2020), 705–734.
  * [29] by same author, _A $H^{3}(G,T)$-valued index of symmetry protected topological phases with on-site finite group symmetry for two-dimensional quantum spin systems_, arXiv:2101.00426, 2021.
  * [30] Y. Ogata and H. Tasaki, _Lieb-Schultz-Mattis type theorems for quantum spin chains without continuous symmetry_ , Commun. Math. Phys. (2019), online, arXiv:1808.08740, published online, 5 Feb 2019.
  * [31] D.W. Robinson, _Properties of propagation of quantum spin systems_ , J. Australian Math. Soc. Series B. Applied Math. (1976), 387–399.
  * [32] D. Ruelle, _Statistical mechanics: Rigorous results_ , Benjamin, New York, 1969\.
  * [33] B. Simon, _The statistical mechanics of lattice gases_ , vol. 1, Princeton University Press, 1993.
  * [34] K. Yosida, _Functional analysis_ , 6th ed., Grundlehren der mathematishen Wissenschaften, Springer-Verlag, 1980.


# Effect of injection conditions on the non-linear behavior of the ECDI and
related turbulent transport

E. Bello-Benítez, A. Marín-Cebrián and E. Ahedo  Equipo de Propulsión Espacial
y Plasmas, Universidad Carlos III de Madrid, Leganés 28911, Spain

(August 28, 2024)

###### Abstract

The electron-cyclotron drift instability (ECDI) has been proposed as one of
the main actors behind the anomalous transport of electrons in Hall thruster
devices. In this work, we revisit the theory and perform two-dimensional PIC
simulations under several conditions to analyze the non-linear behavior and
the induced transport under several boundary conditions. Simulation results
with fully-periodic boundaries and conditions faithful to the linear theory
show the growth of ECDI modes, ion-wave trapping vortexes and agree with the
existing literature in early times. In the long term, however, we observe very
mild oscillations and null anomalous current. The evolution towards this new
equilibrium is coherent to what can be expected from energy conservation. The
quenching of the oscillations seem to be highly related with the distortion of
ion-trapping vortexes in phase space after a long-term interaction of ion
particles with the electrostatic wave. This result suggests that sustained
oscillations and turbulent current could benefit from the renewal of ions by,
e.g., removing and injecting particles through axial boundaries instead of
applying periodicity. This second type of simulations shows that injection
conditions highly impact the late simulation behavior of ECDI oscillations,
where we identify several regimes depending on the value of the ion residence
time compared to the characteristic saturation time in the fully periodic
case. The intermediate regime, where these two times are close, is the only
one providing sustained oscillations and electron transport and seems to be
the relevant one in Hall devices.

##  I Introduction

The problem of anomalous electron cross-field transport remains as one of big
open challenges for the community of $\bm{E}\times\bm{B}$ plasmas. In the
field of plasma propulsion, this problem has been mainly studied in the
context of Hall-thruster discharges and represent one big obstacle on the way
towards to predictive efficient numerical models. The large drift of electrons
in the azimuthal direction of the Hall thruster is a source of several
families of azimuthal oscillations that are potential candidates to explain
the anomalous transport and have been observed experimentally [1, 2, 3, 4, 5,
6]. The classical explanation [7] for the impact of oscillations on transport
relies on the correlation of oscillations in density and electric field in the
$\bm{E}\times\bm{B}$ direction under the presence of a magnetic field.

The amount of articles devoted to the analysis of instabilities and turbulence
in Hall thrusters is extensive. With a macroscopic description for ions and
electrons, some of the authors of this article have conducted global [8] and
local [9] linear stability analyses. In these references two-stream, drift-
gradient and drift-dissipative instabilities are discussed. Similar recent
studies by other authors are [10, 11, 12, 13, 14].

When using a kinetic formulation for the electrons, the analytical studies of
instabilities are usually limited to a homogeneous and collisionless plasma.
For the conditions of a Hall thruster, where electrons are magnetized but ions
are not, the dispersion relation of the electron cyclotron drift instability
(ECDI) is obtained[15, 16]. This classical instability have been revisited,
during the last two decades, by several authors [17, 18, 19, 20, 21], in the
context of Hall thrusters. Also from the point of view of kinetic simulations
[17, 22, 23, 24, 25, 26] and experiments [27, 4, 5, 6].

Kinetic models aimed to analyze electron turbulence and turbulent transport in
the plane perpendicular to the applied magnetic field $\bm{B}_{0}$ can be
classified in 1D azimuthal[25, 22, 26] and 2D axial-azimuthal [17, 23, 24].
The latter case has been the subject of a recent benchmark by several groups
[28]. Many 2D simulations include a number of phenomena that makes challenging
to compare the results with the ECDI linear theory; such as, inhomogeneous
magnetic field, collisions, ionization or electrical connection between anode
and cathode. The 1D-azimuthal simulations are closer to the linear theory of
the ECDI but some of them still add effects that are not considered in the
dispersion relation, such as refreshing of particle velocities or collisions.
Moreover, the 1D simulation results reported in the literature have been seen
to be sensitive to the axial treatment of particles; observing a sustained
saturated state only when using a virtual axial dimension and particle
velocities are refreshed after leaving the domain axially [22, 26]. If
particle are not refreshed an unbounded growth of oscillations and heating are
seen [25, 29].

In section II, we revisit the main results from the classical dispersion
relation of the ECDI theory. A recently developed in-house two-dimensional PIC
code is introduced in section III and used to analyze the nonlinear evolution
of the ECDI under several conditions. In every case, a 2D domain perpendicular
to $\bm{B}_{0}$ is considered; which we will refer to as the axial-azimuthal
plane, making the analogy with a Hall discharge geometry. Initially we
consider a fully periodic domain and simulate under the assumptions of the
classical ECDI theory. These results can be compared with the existing
literature [25, 22, 26] on 1D azimuthal simulations of the classical ECDI, but
open the possibility to the growth of modes with an axial component. Similar
results to those shown elsewhere are found during early simulation times; but
our long-term behavior shows no oscillation capable of driving an turbulent
electron transport, which disagrees with the literature. We try to give an
explanation from an energy conservation perspective.

Periodic results establish a clear relation between ion behavior and the
existence of electron transport in the long term, what suggests that renewal
of ion particles could help to develop sustained ECDI oscillations and axial
transport. This motivate us to replace axial boundary conditions from periodic
to removal/injection of particles, which mimics the generation and loss of
particles in a finite plasma such as the Hall discharge. Results are shown in
section IV, where we are able to establish a clear relation between the ion
velocity (related with the residence time) and the long-term behavior of
short-wavelength oscillations and transport in the plasma.

##  II The classical ECDI: linear analysis

We attempt to study the stability of an homogeneous, collisionless plasma at
equilibrium subjected to mutually perpendicular magnetic
$\bm{B}_{0}=B_{0}\bm{1}_{x}$ and electric $\bm{E}_{0}=E_{0}\bm{1}_{z}$ fields.
Throughout the article subindex ‘0’ and ‘1’ stand for equilibrium and
perturbed conditions. Consistent with Hall thruster discharges: the strength
of the stationary magnetic field keeps electrons well-magnetized and ions
unmagnetized, and the plasma currents are low enough to neglect the
contribution of the self-induced magnetic field $\bm{B}_{1}$. At equilibrium,
electrons have a drifted-Maxwellian velocity distribution function (VDF) with
density $n_{0}$, temperature $T_{e0}$, and an $\bm{E}_{0}\times\bm{B}_{0}$
drift velocity $\bm{u}_{e0}=u_{ye0}\bm{1}_{y}$, $u_{ye0}=E_{0}/B_{0}$.
Equilibrium ions are assumed cold, have a density $n_{0}$, and move with a
constant velocity $\bm{u}_{i0}=u_{zi0}\bm{1}_{z}$, parallel to $\bm{E}_{0}$.
This homogeneous equilibrium disregards the electrostatic acceleration of ions
due to $\bm{E}_{0}$. From analogy with a Hall-thruster geometry, let us refer
throughout the article to the directions of $\bm{B}_{0}$, $\bm{E}_{0}$ and
$\bm{E}_{0}\times\bm{B}_{0}$ as radial ($x$), axial ($z$) and azimuthal ($y$).

Next, we summarize the main aspects of the well-known linear stability
analysis [15, 16, 18, 19, 20] From the linearly perturbed Vlasov equation, the
relation between the electron density and the electric potential perturbations
of wavevector $\bm{k}$ and complex frequency $\omega$, for propagation
perpendicular to $\bm{B}_{0}$ (i.e., the component of $\bm{k}$ parallel to
$\bm{B}_{0}$ is zero), is

$\frac{n_{e1}}{n_{0}}=\Big{[}1-g(\omega_{e},b_{e};\omega_{ce})\Big{]}\frac{e\phi_{1}}{T_{e0}}$
(1)

with

$g(\omega_{e},b_{e};\omega_{ce})=\exp(-b_{e})\Big{[}I_{0}(b_{e})+2\sum_{m=1}^{\infty}\frac{\omega_{e}^{2}I_{m}(b_{e})}{\omega_{e}^{2}-m^{2}\omega_{ce}^{2}}\Big{]}.$
(2)

In these equations $\bm{E}=-\nabla\phi$, $c_{e0}=\sqrt{T_{e0}/m_{e}}$ is the
electron thermal velocity, $\omega_{ce}=eB_{0}/m_{e}$ is the cyclotron
frequency, $b_{e}=k^{2}\rho_{e0}^{2}$, $\rho_{e0}=c_{e0}/\omega_{ce}$ is the
electron Larmor radius, $\omega_{e}=\omega-k_{y}u_{ye0}$ is the electron
Doppler-shifted frequency, $\bm{k}=k_{y}\bm{1}_{y}+k_{z}\bm{1}_{z}$, and
$I_{m}$ are the modified Bessel functions of the first kind. According to
equation (2), the perturbed electron response has resonances at the electron
gyrofrequency and its harmonics, i.e. $\omega_{e}=m\omega_{ce}$.

Next, the density perturbations of cold ions follow

$\frac{n_{i1}}{n_{0}}=\frac{k^{2}c_{s0}^{2}}{\omega_{i}^{2}}\frac{e\phi_{1}}{T_{e0}},$
(3)

with $\omega_{i}=\omega-k_{z}u_{zi0}$ the ion Doppler-shifted frequency, and
$c_{s0}=\sqrt{T_{e0}/m_{i}}$ the ion sound speed. Using Eqs. (1) and (3),
together with the linearized Poisson equation, yields the two-dimensional
dispersion relation

$1+k^{2}\lambda_{D0}^{2}=\frac{k^{2}c_{s0}^{2}}{\omega_{i}^{2}}+g(\omega_{e},b_{e};\omega_{ce}).$
(4)

This equation is solved for the complex frequencies
$\omega=\omega_{r}+\text{i}\gamma$ with all other parameters fixed, including
$\bm{k}$.

Type | Description and symbol | Value and units
---|---|---
Fundamental plasma parameters | Ion mass, $m_{i}$ | 1 u
Electric field, $E_{0}$ |  $10^{4}$ V/m
Magnetic field, $B_{0}$ | 200 G
Plasma density, $n_{0}$ |  $10^{17}$ m-3
Ion axial velocity, $u_{zi0}$ | 2.5 km/s
Electron temperature, $T_{e0}$ | 6 eV
Derived plasma parameters | Electron azimuthal drift, $u_{ye0}$ | 500 km/s
Electron thermal speed, $c_{e0}$ | 1027 km/s
Sound speed, $c_{s0}$ | 23.97 km/s
Debye length, $\lambda_{D0}$ | 57.58 $\mu$m
Electron azimuthal-drift gyroradius, $\ell_{e0}$ | 142.1 $\mu$m
Electron Larmor radius, $\rho_{e0}$ | 292.0 $\mu$m
Electron plasma frequency, $\omega_{pe0}$ | 2.839 GHz
Electron gyrofrequency, $\omega_{ce}$ | 0.5600 GHz
Ion plasma frequency, $\omega_{pi0}$ | 66.26 MHz
| Lower-hybrid frequency, $\omega_{lh}$ | 13.07 MHz
Fundamental numerical parameters | Azimuthal domain length, $L_{y}$ | 5.359 mm
Axial domain length, $L_{z}$ | 2.679 mm
Number of azimuthal celss, $N_{y}$ | 100
Number of axial cells, $N_{z}$ | 50
Number of particles per cell, $N_{\textmd{ppc}}$ | 200
Time steps, $\Delta t$ |  $5\times 10^{-12}$ s
| Number of time steps, $N_{t}$ | $6\times 10^{5}$
| Number of time steps between print-outs, $N_{\textmd{print}}$ | 1000
Derived numerical parameters | Azimuthal cell size, $\Delta y$ | 53.59 $\mu$m
Axial cell size, $\Delta z$ | 53.59 $\mu$m
Table 1: Physical and numerical parameters of the reference simulation case.
The subscript ‘0’ stands for initial equilibrium conditions. Derived parameter
values are included for completeness but can be computed from fundamental
ones.

The solutions of the dispersion relation are pairs of modes of three types
[16]. First, there is a pair of ion-acoustic modes with $g\ll
k^{2}c_{s0}^{2}/\omega_{i}$ and real-valued frequencies
$\omega_{i}=\pm\omega_{\textmd{IA}}$; being

$\omega_{\textmd{IA}}=\frac{kc_{s0}}{\sqrt{1+k^{2}\lambda_{D0}^{2}}},$ (5)

which include a non-neutral term in the denominator that corrects the
quasineutral linear relation $\omega_{i}=\pm kc_{s0}$. Second, there are
electron Bernstein waves [30], with $g\gg k^{2}c_{s0}^{2}/\omega_{i}$ and
real-valued frequencies too, close to the resonances
$\omega_{e}=m\omega_{ce}$. Third, when $g\sim k^{2}c_{s0}^{2}/\omega_{i}$,
Bernstein waves become coupled with the ion acoustic modes, yielding one of
them the modified ion acoustic (MIA) pair with frequencies

$\omega_{i}=\pm\frac{kc_{s0}}{\sqrt{1+k^{2}\lambda_{D0}^{2}-g(\omega_{e})}}.$
(6)

For $g>1+k^{2}\lambda_{D0}^{2}$, one of MIA modes is unstable leading to the
so called ECDI. Since $g$ changes from $-\infty$ to $+\infty$ when crossing
the resonance $\omega_{e}=m\omega_{ce}$ with $\omega_{e}$ increasing, the MIA
mode starts becoming unstable at $\omega_{e}=(m\omega_{ce})^{+}$ and becomes
stable before reaching the $(m+1)$ resonance [18, 19, 20].

Figure 1 depicts, for the equilibrium solution of Table 1, the complex
frequency of the purely-azimuthal (i.e. $\bm{k}=k_{y}\bm{1}_{y}$) MIA mode,
and shows the growth rate of the ECDI linked to each $m\omega_{ce}$; as a
function of $k_{y}\ell_{e0}$, being $\ell_{e0}=u_{ye0}/\omega_{ce}$ the
azimuthal-drift gyroradius. The ion-acoustic frequency $\omega_{\textmd{IA}}$
is also plotted for comparison with the MIA mode. The same equilibrium
solution will be used in the kinetic simulations of the coming sections. For
$k_{y}>0$, the unstable mode has $\omega_{r}>0$, meaning propagation in the
$\bm{E}_{0}\times\bm{B}_{0}$ direction. Due to symmetry, for $k_{y}<0$, the
unstable mode has $\omega_{r}<0$ (i.e., moving again along
$\bm{E}_{0}\times\bm{B}_{0}$) and same $\gamma$. In the purely azimuthal ECDI,
one has $\omega\sim k_{y}c_{s}\ll\omega_{e}\simeq-
k_{y}u_{ye0}\approx-m\omega_{ce}$. Therefore, the electron (Bernstein-type)
response is quasi-steady and defines mostly the wavelength of the ECDI.
Regarding nonneutral effects the MIA is quasineutral as long as
$k_{y}^{2}\lambda_{D0}^{2}\ll 1$, i.e.

$\dfrac{\lambda_{D0}}{\rho_{e}}\ll\dfrac{u_{y_{e0}}}{mc_{e0}}$ (7)

and non-neutral effects tend to reduce the complex frequency. This effect is
very clearly seen in the curve $\omega_{\textmd{IA}}$ in figure 1(a), that
shows frequencies lowered with respect to a quasineutral ion-acoustic linear
relation. The comparison of the MIA and ion-acoustic frequencies demonstrates
that they follow the same trend but with significant deviations coming from
the coupling with the Bernstein terms.

(a)

(b)

Figure 1: ECDI dispersion relation for a hydrogen plasma with
$k_{\parallel}=0$, $k_{z}=0$ and the equilibrium conditions summarized in
table 1. The black solid line is the solution for an infinite plasma. Red
crosses stand for the discrete solutions in a finite plasma with $L_{y}=5.359$
mm. The blue solid line is the ion-acoustic frequency $\omega_{\textmd{IA}}$.
The dispersion relation has been solved with the same numerical approach than
reference [19].

If the plasma has a finite size $L_{y}$ along $y$, there is only a discrete
wave spectrum with $k_{y}=n2\pi/L_{y}$ and $n$ the number of wavelengths
fitting in the domain. Red crosses in Fig. 1 show that spectrum for
$L_{y}=12\pi\ell_{e0}=5.359$ mm, when resonances correspond to $n$ a multiple
of 6. That length has been chosen so that modes $n=1+6m$ capture approximately
the peaks in $\gamma$ associated to each resonance $m$. For the chosen
parameters, the fastest growing mode is $n=13$, close to the $m=2$ resonance,
with $\gamma/2\pi=34.7\times 10^{6}$ s-1, $k_{y}\ell_{e0}=2.167$ and
$\omega_{r}=46.6$ MHz. In terms of growth rate, the mode $n=7$, in the band
$m=1$, follows closely with $\gamma/2\pi=32.6\times 10^{6}$ s-1,
$k_{y}\ell_{e0}=1.167$ and $\omega_{r}=27.1$ MHz.

The effects of a nonzero (moderate) $k_{z}$ on the MIA modes are: introducing
an ion-Doppler shift $k_{z}u_{zi0}$ in the frequencies, shifting the unstable
bands in $k_{y}$, and changing mildly the growth rates (find a more exhaustive
analysis in [18]).

To reduce the computational cost, simulations here correspond to a hydrogen
plasma. The frequencies for hydrogen are, approximately, one order of
magnitude higher than those expected in xenon. This is a reasonable result,
since equation (4) shows that frequency and growth of the ECDI modes are
proportional to $kc_{s0}$ and, thus, scale with $1/\sqrt{m_{i}}$. This trend
is, indeed, retrieved in PIC simulations shown later in this paper. The use of
hydrogen instead of xenon in PIC simulations of the ECDI is a computational
advantage since it allows us to observe the same physical phenomena but in a
shorter time. With the same time step (still limited by the electron
dynamics), this means reducing the number of time steps by one order of
magnitude.

##  III The classical ECDI: nonlinear evolution

###  III.1 The numerical PIC model

Figure 2: Diagram summarizing the simulation axes, boundary conditions and
initial equilibrium state for simulations of the classical ECDI. A second type
of simulations is considered that substitutes left and right boundary
conditions to $\phi=0$ and injection/absorbing conditions for particles.

The nonlinear evolution and saturation of the classical ECDI is studied with a
2D axial azimuthal, 2D($z,\theta$), full PIC code developed in-house. The PIC
formulation follows ions and electrons. The electrostatic potential is
obtained from a Poisson solver. The numerical codes are implemented in Fortran
and use OpenMP shared-memory parallel computing.

The PIC code applies a standard Boris method to move electron and ion
macroparticles in the periodic domain, and it employs a particle-decomposition
strategy for parallel calculations. Macroparticles have equal and constant
weights (i.e., number of real particles per macroparticle). As already pointed
out, collisions between particles are totally disregarded.

The Poisson solver is able to use different schemes depending on boundary
conditions. When all boundaries are periodic (the case in the present section
III), spectral methods are well suited to solve the Poisson equation in the
Fourier complex space; here, the FFTW3 library [31] for Fourier and inverse
transform operations is used, with a zero average potential. If Neumann or
Dirichlet conditions are used in at least one boundary (the case in section
IV), the Poisson solver uses a second order finite difference scheme for the
Laplace operator and electric field, and the discrete linear system is solved
with the PARDISO direct-solver routines in the Math Kernel Library of INTEL.

Figure 3: Reference case with fully periodic conditions: time-evolution of
$E_{y}(y,z,t)$ in the $yz$-plane.

(a) $L_{y}=2.679\leavevmode\nobreak\ \textmd{mm}$

(b) $L_{y}=5.359\leavevmode\nobreak\ \textmd{mm}$

(c) $L_{y}=10.72\leavevmode\nobreak\ \textmd{mm}$

(d) $L_{y}=16.08\leavevmode\nobreak\ \textmd{mm}$

Figure 4: Reference case with fully periodic conditions: time-evolution of the
logarithm of normalized coefficients from the fast Fourier transform of
$E_{y}(t,y,z^{*})$ in coordinate $y$, at $z^{*}=5L_{z}/6$, for different
$L_{y}$.

(a) $0\leavevmode\nobreak\ \mu\textmd{s}<t<0.20\leavevmode\nobreak\
\mu\textmd{s}$

(b) $0.20\leavevmode\nobreak\ \mu\textmd{s}<t<0.75\leavevmode\nobreak\
\mu\textmd{s}$

(c) $2.5\leavevmode\nobreak\ \mu\textmd{s}<t<3.0\leavevmode\nobreak\
\mu\textmd{s}$

Figure 5: Reference case with fully periodic conditions: logarithm of
normalized coefficients from the 2D fast Fourier transform of
$E_{y}(t,y,z^{*})$ in coordinates $y$ and $t$, at $z^{*}=5L_{z}/6$, and for 3
different time windows. The absolute maxima are marked with a black cross.
Dashed magenta lines stand for ion-acoustic modes $\pm\omega_{\textmd{IA}}$
with $T_{e0}$ taken as the average value on the considered window: (a) 11.4
eV, (b) 40.1 eV and (c) 58.3 eV.

The electric and magentic fields felt by each species and their initial
macroscopic properties comply with the hypotheses and equilibrium state of
section II. That is to say, the initial electron population is randomly
sampled from a Maxwellian VDF with density $n_{0}$, velocity
$u_{ye0}\bm{1}_{y}$ and temperature $T_{e0}$. The initial ion population is
generated with $n_{0}$ and $u_{zi0}\bm{1}_{z}$ without dispersion. The values
of these properties and every other physical and numerical parameter are
gathered in Table 1. The macroparticle weights are kept constant throughout
the simulation and can be computed from initial conditions as $n_{0}\Delta
y\Delta z/N_{\textmd{ppc}}$ Electron and ion particles are moved in a periodic
domain with electric fields $\bm{E}_{e}=\bm{E}_{0}+\bm{E}_{1}$ and
$\bm{E}_{i}=\bm{E}_{1}$, respetively, being $\bm{E}_{1}=-\nabla\phi$ the local
fluctuation relative to $\bm{E}_{0}$ that comes as solution to the Poisson
equation with periodic boundary conditions. Figure 2(a) sketches the
simulation setup. At the initial equilibrium state, the axial current of
electrons, $j_{ze0}$ is zero. Any subsequent electron axial current is due to
the cross-field transport generated by plasma oscillations.

Regarding the macroscopic and kinetic results to be shown, they include a
moving-average during runtime on a window corresponding to
$N_{\textmd{print}}$ time steps. For an arbitrary variable $\varphi$, this is
defined as

$\tilde{\varphi}_{k}=\frac{\tilde{\varphi}_{k-1}(N_{\textmd{print}}-1)+\varphi_{k}}{N_{\textmd{print}}},$
(8)

being $k$ the time step index, $\varphi_{k}$ the instantaneous value of and
$\tilde{\varphi}_{k}$ its time-average value.

###  III.2 Onset, saturation, and vanishing of the ECDI

The time evolution in the $yz$-plane of $E_{y}(t,y,z)$ is represented in
figure 3. The initial equilibrium state is unstable and the ECDI start to grow
from any perturbation. The time-dependent solution is almost 1D in the
azimuthal direction, although an axial component is present in early times,
mainly. The oscillation amplitude gets a maximum around 0.2 $\mu$s and
decreases in later times until eventually a new equilibrium state, different
from the initial one, is reached. For $t<0.2\leavevmode\nobreak\ \mu$s, wave
modes are weakly mixed and the dominant monochromatic waves are easier to
observe. The top panels of figure 3 indicate that the dominant mode is $n=7$,
which is the closest one to the resonance $m=1$. For
$t>0.2\leavevmode\nobreak\ \mu$s, there is more mixing of modes; the bottom
panels of figure 3 show transitions to $n=5$ and 6 as dominant modes.

The linear theory showed that mode $m=2$ has a (slightly) higher growth rate
than mode $m=1$. In fact the Fourier analysis of $E_{y}(t,y,z)$ in $y$
[depicted in figure 4(b) for our simulation with $L_{y}=5.359$mm] shows some
contribution of mode $m=2$ to the early time spectrum. However, the fast
growth of several modes makes nonlinear effects important soon in the
simulation, which, together with the noise intrinsic to the PIC formulation,
makes tough the exact comparison of early PIC results with the linear results
from Vlasov equation [29, 32]. In addition, the long-term dominant modes in
the nonlinear stage may not coincide with the most unstable modes in the
linear dispersion relation [21]. The use of quiet start techniques by other
authors [32, 22] to minimize the noise of the initial population has not been
seen to be completely satisfactory.

Figure 6: Reference case with fully periodic conditions: time evolution, in
the phase space $(y,v_{y})$, of ions contained in an axial slab of width
$0.1L_{z}$. The formation and blurring of the vortex-like structure
characteristic of ion-wave trapping is observed.

The determination of the dominant frequency is more difficult since it depends
on time itself. In Figure 5, results of the 2D fast Fourier transform in $y$
and $t$ are shown for three different time windows, together with the ion-
acoustic curves $\pm\omega_{\textmd{IA}}$ for the average $T_{e}$ within the
corresponding time interval. Each window represents different stages in the
evolution of the ECDI:

  * (i)

$0<t<0.20\leavevmode\nobreak\ \mu\textmd{s}$. The peaks seem to concentrate in
bands near $k_{y}\ell_{e0}=m$, similarly to the theoretical dispersion
relation in figure 1. The maximum Fourier coefficient is located at
$k_{y}\ell_{e0}=1.166$ and $\omega_{r}=55.6$ MHz, with phase speed
$\omega_{r}/k_{y}=42.6$ km/s. This is mode number $n=7$, near the resonance
$m=1$. There are secondary bands close to $m=2$ and $m=3$. The results in this
early stage are qualitatively aligned with figure 1, but the dominant
frequency $\omega_{r}$ is larger than predicted by the linear theory.

  * (ii)

$0.2\leavevmode\nobreak\ \mu\textmd{s}<t<0.75\leavevmode\nobreak\
\mu\textmd{s}$. The bands of the upper spectrum near the resonances have been
blurred and there is an approximately linear relation between $k_{y}$ and
$\omega_{r}$, resembling a linear ion-acoustic relation. However, some parts
of the upper spectrum seem to follow the non-neutral acoustic frequency
$\omega_{\textmd{IA}}$. This is aligned with previous PIC simulations[23, 28]
and experiments[4]. The peak in the spectrum is at $k_{y}\ell_{e0}=0.8334$;
with $n=5$, $\omega_{r}=49.1$ MHz and $\omega_{r}/k_{y}=52.6$ km/s and matches
the ion-acoustic curve. Let us note that the lower part of the spectrum
(propagation in the $-\bm{E}_{0}\times\bm{B}_{0}$ direction), shows some
remnants of the counter-propagating ion-acoustic wave with frequency
$-\omega_{\textmd{IA}}$.

  * (iii)

$2.5\leavevmode\nobreak\ \mu\textmd{s}<t<3\leavevmode\nobreak\ \mu\textmd{s}$.
Even if a peak is identifiable at $k_{y}\ell_{e0}=1.166$ ($n=7$), the mixing
of different temporal and azimuthal scales result in a messy spectrum without
a clear dominant mode. The ion-acoustic behavior is only seen in the lower
part of the spectrum.

At $t\sim 3\leavevmode\nobreak\ \mu\textmd{s}$, the plasma seems to tend to a
new equilibrium with distorted ion and electron VDFs. Figure 6 plots the
evolution of ion particles, contained within an axial slab of width
$0.1L_{z}$, in phase space $(y,v_{y})$. During the growth and saturation of
the ECDI (this is $t<0.2\leavevmode\nobreak\ \mu$s), vortexes are formed in
phase space showing the characteristic behavior of ion trapping in the
electrostatic wave. The vortexes are shifted towards positive velocities,
matching the $\bm{E}_{0}\times\bm{B}_{0}$ and dominant mode propagation
directions. Later times reveal the distortion of those vortex structures until
they are fully blurred into a strongly one-sided distribution with a long tail
into positive velocities. This process coincides with the quenching of $E_{y}$
oscillations. Electrons one-dimensional VDFs at the end of the simulation are
shown in figure 7 to be fairly isotropic and flatter close to the average
velocity compared to a Maxwellian, which coincides with [25, 18].

Figure 7: Reference case with fully periodic conditions: normalized electron
distribution functions $f_{je}$ in the azimuthal ($j=y$, red dashed) and axial
($j=z$, blue dash-dotted) directions, being the abscissa coordinate $w_{je}$
the random velocity of electrons . The Maxwellian curve (black solid) is
included for comparison, using local final values of $n_{e}=10^{17}$ and
$T_{e}=58.8$ eV from the PIC simulation.

###  III.3 Evolution of the plasma energy

(a) Energy

(b) Energy derivative

(c) Effect of $L_{y}$ on $\bar{\mathcal{E}}$

(d) Effect of $L_{y}$ on $\left\langle j_{ze}\right\rangle$

Figure 8: Fully-periodic case. Time evolution of (a) the energy time-
derivative and (b) the electrostatic, species and total energies per particle
(i.e., $\bar{\mathcal{E}}=\mathcal{E}/\left\langle n\right\rangle V$. Panels
(c) and (d) plot energy per particle and volume-averaged electron axial
current for different $L_{y}$.

Let us get a further insight on the ECDI by analyzing the plasma energy stored
in the plasma. The total energy in the plasma domain has contributions from
the electrostatic field, electrons and ions, according to

$\mathcal{E}=\mathcal{E}_{E}+\mathcal{E}_{e}+\mathcal{E}_{i},$ (9)

with

$\mathcal{E}_{E}=\frac{\varepsilon_{0}}{2}\int_{\mathcal{V}}E_{1}^{2}\leavevmode\nobreak\
\text{d}\mathcal{V}.$ (10)

the energy of the electromagnetic oscillations,

$\mathcal{E}_{s}=\int_{\mathcal{V}}\left[\frac{1}{2}m_{s}n_{s}u_{s}^{2}+\frac{3}{2}p_{s}\right]\text{d}\mathcal{V},\quad
s=e,i$ (11)

the total energy of electrons and ions; and $\mathcal{V}$ is the volume of the
domain. For convenience, the equilibrium field $\bm{E}_{0}$ has been omitted
from the definition of $\mathcal{E}_{E}$, since it would only introduce a
constant offset in the energy. The energy of each species is approximated in
the PIC formulation as

$\mathcal{E}_{s}=\sum_{p}\frac{1}{2}m_{s}v_{p}^{2}W_{p},$ (12)

where the sum is on every particle in the domain with $v_{p}$ and $W_{p}$ the
particle speed and weight, respectively.

In a consistent situation, the work done by the electric field should act as a
mechanism that converts species energy on electric-field energy, and the other
way around. However, the assumptions behind the linear theory of the ECDI
forced us to let ions and electrons feel different electric fields. Because of
this non-conventional feature, it can be proved that the total energy changes
according to

$\frac{\partial\mathcal{E}}{\partial
t}=\int_{\mathcal{V}}(\bm{j}_{e}\cdot\bm{E}_{0})\leavevmode\nobreak\
\text{d}\mathcal{V}=E_{0}\int_{\mathcal{V}}j_{ze}\leavevmode\nobreak\
\text{d}\mathcal{V}=E_{0}\langle j_{ze}\rangle\mathcal{V},$ (13)

with $\bm{j}_{e}$ the electron current density and $\langle j_{ze}\rangle$ the
volume-averaged $j_{ze}$. Therefore, this kind of simulation will not show a
proper conservation of energy and the equilibrium electric field will pump
energy into the isolated system. However, this source of energy requires also
an axial electron current $j_{ze}$ to be developed. This means that the energy
is conserved initially until the instability is triggered, and any other
stationary energy state should satisfy a null $\langle j_{ze}\rangle$, i.e.,
no turbulent electron transport.

Figure 8(a) shows the evolution of electrostatic, species, and total energies
in the domain per real particle. Here, average energies per particle
$\bar{\mathcal{E}}=\mathcal{E}/\left\langle n\right\rangle\mathcal{V}$ are
used, being $\left\langle n\right\rangle$ the volume-averaged density and
$\left\langle n\right\rangle\mathcal{V}$ the real number particles in the
domain. Initially $\bar{\mathcal{E}}_{E}=0$,
$\bar{\mathcal{E}}_{i}=m_{i}u_{zi0}^{2}/2=0.0324$ eV and
$\bar{\mathcal{E}}_{e}=3T_{e0}/2+m_{e}u_{ye0}^{2}/2=9.71$ eV. Eventually, the
total energy saturates and becomes stationary. There is significant heating of
electrons and, to a lesser extent, of ions; which is coherent with other works
[29, 25, 22, 26]. The electrostatic energy is much lower than those of
electrons and ions, and approaches zero for late simulation times, when
oscillations are very weak.

The balance (13) is evaluated in figure 8(b). The rate
$\textmd{d}\bar{\mathcal{E}}/\textmd{d}t$ is approximated numerically and
compared with the source term $E_{0}\langle j_{ze}\rangle/n_{0}$. The two
curves show an excellent match apart from the noise inherent to the PIC
approach, specially problematic in the calculation of $j_{ze}$. The PIC
simulation is able to approximately replicate the theoretical energy balance
(13), which is a good sign of the validity of such simulations. The change in
energy shows a peak and then decreases tending to zero for late simulation
times, meaning that a new stationary equilibrium state of the plasma is
reached. As already said, the energy source is related with the development of
a net axial current, so the new equilibrium holds an average $j_{ze}$ equal to
zero. However, at mid simulation times, saturated ECDI modes are effective in
inducing an axial transport of electrons.

###  III.4 Effect of the domain’s azimuthal length

As pointed out in section II, the azimuthal domain length $L_{y}$ is an
important parameter that determines the unstable ECDI modes that can be
excited in a finite domain. On the contrary, the effect of $L_{z}$ has been
seen to be small, although the chosen axial length is large enough to allow
the formation of waves with axial component fitting $L_{z}$. Increasing
$L_{y}$ allows larger scales to develop and increases the spectral resolution
in the $k_{y}$-space, better capturing the continuous dispersion relation for
an infinite plasma. In figure 4, the time evolution of azimuthal Fourier
coefficients for $E_{y1}$ are represented for several $L_{y}$ multiples of
2.679 mm, the reference case corresponding to plot 4(b). The main trends
identified in the shortest case are also recovered in larger domains: modes
close to resonant bands $m=1$ and $m=2$ are excited in early times, modes
close to $m=1$ dominate at mid times, and oscillations quench after the
spectrum peaks are passed. These trends become clearer as longer domains are
used. Some differences worth mentioning among cases are: (i) a mode with
azimuthal wavelength equal to $L_{y}$ (i.e. $n=1$) is present only at long
$L_{y}$ and (ii) the enhanced spectral resolution of longer $L_{y}$ allows to
better capture shorter resonant bands, such as the modes close to $m=2$.
Similar conclusions on these two points are reached in reference [29].

Figure 8(c) shows that the energy per real particle increases with parameter
$L_{y}$, probably due to the greater number of unstable excited modes [29].
With the exception of the shortest case $L_{y}=2.679$ mm, which shows a slower
energy increase and lower $\langle j_{ze}\rangle$ [most probably due to a poor
spectral resolution, unable to properly capture the peaks of the growth rate
of Fig 1] the rates of energy increase during the growth period are similar
for $L_{y}\geq 5.359$ mm. This suggests that the maximum $\langle
j_{ze}\rangle$ and the time to reach it are numerically robust and thus
physical. Therefore, results with $L_{y}=5.359$ mm are representative of
simulations with longer $L_{y}$ multiples of 5.359 mm. This behavior disagrees
with reference [29], where changes in $L_{y}$ can drastically change the
transients of $\langle j_{ze}\rangle$ and, thus, the energy growth rates.

Since the long term behavior of $\langle j_{ze}\rangle$ tends to zero in every
case, we can conclude that the long-wavelength mode with $n=1$ that develops
with increasing $L_{y}$ is not effectively producing an axial electron
transport. This is consistent with the late evolution of $n_{e}$ and $E_{y}$
(not included), which shows completely out of phase oscillations.

###  III.5 Effect of the ion mass

(a) $m_{i}=1\textmd{u}$

(b) $m_{i}=4\textmd{u}$

(c) $m_{i}=9\textmd{u}$

Figure 9: Fully-periodic case, Time-evolution of the logarithm of normalized
azimuthal Fourier coefficients for different values of $m_{i}$.  Figure 10:
Scaling factor, with the ion mass $m_{i}$, of several variables: dominant
$\omega_{r}$ in the interval $t\omega_{Pi}<83.15$, peak averaged $j_{ze}$,
final $\mathcal{E}$ and $\Delta t_{\textmd{sat}}$. The curves $\sqrt{m_{i}}$
and $1/\sqrt{m_{i}}$ are also included for comparison.

In this subsection, the simulation is repeated for increasing values of the
ion mass. While hydrogen mass, $m_{H}$ has been used in order to speed up the
dynamics and minimize the computational cost, the practical cases of interest
deal with heavier ions, e.g., xenon or krypton. In addition, an increased
azimuthal length $L_{y}=16.08$ mm is used to have a proper spectral resolution
and a fair comparison, since an increased $m_{i}$ narrows the unstable bands
in the dispersion relation near the resonances what needs an increased
spectral resolution.

In the early evolution of the ECDI, the dispersion relation (4) shows that the
frequency should scale inversely proportional to $\sqrt{m_{i}}$. However, this
does not have to be the case of the nonlinear dynamics or other magnitudes.
The scaling factors for several quantities, defined as the ratio with respect
to the value in the reference case $m_{i}=1$ u are computed in figure 10:
dominant $\omega_{r}$ in the interval $t\omega_{Pi}<83.15$, peak of
$\left\langle j_{ze}\right\rangle$, final $\mathcal{E}$ and $\Delta
t_{\textmd{sat}}$. The observed trends are coherent with results reported in
reference [33]. If the dominant frequency is measured in early simulation
times, the observed scaling perfectly matches the linear theoretical result.
The saturation time $\Delta t_{\textmd{sat}}$, defined as the time of maximum
average $\left\langle j_{ze}\right\rangle$, seems to scale proportional to
$\sqrt{m_{i}}$. The 1D azimuthal Fourier spectra are compared in figure 9,
where the time axis is normalized with the inverse of the ion plasma frequency
$\omega_{Pi}$. While the saturation of the early oscillations produced by the
ECDI happens at similar $t\omega_{Pi}$ (as shown by the scaling of $\Delta
t_{\textmd{sat}}$), the quenching of short wavelength $m=1$ and $\left\langle
j_{ze}\right\rangle$ seems to be slightly slower with increasing $m_{i}$.

The level of turbulent current, measured by the peak $\left\langle
j_{ze}\right\rangle$, initially decreases with $m_{i}$, which could be
explained by the weaker oscillations observed in $n_{e}$ (not included here),
but seems to saturate. The current $\left\langle j_{ze}\right\rangle$
decreases, however, slower than $1/\sqrt{m_{i}}$ which leads to a final
$\mathcal{E}$ increasing with $m_{i}$, in accordance with equation (13).
Extrapolating these results, with xenon we can expect somehow smaller level of
transport but dilated in time, what would cause a significant heating such as
that observed in similar studies in the literature [22, 21, 29].

###  III.6 Saturation behavior in previous literature

Figure 11: Finite plasma with axial injection. Time-evolution of $E_{y}$ in
the $yz$-plane.

The results we have shown in this section pursue the simulation of scenarios
as close as possible to the classical ECDI. While we observe the growth and
saturation of oscillations due to an instability, modes that induce an axial
electron transport seem to vanish at long times. The late simulation behavior
differ from other 1D [25, 22, 26, 29] and 2D [23, 28, 17, 24] simulations on
similar configurations, which report sustained oscillations not vanishing with
time. Even if our model is 2D, we disregard many of the effects included in
other 2D works (e.g., inhomogeneous magnetic field, anode-cathode circuit or
collisions), which are closer to Hall-discharge-like configurations. A fairer
comparison is that with the existing 1D azimuthal models [25, 22, 26, 29],
whose sustained oscillation for late simulation times either saturate[22, 26]
or forever grow [25, 29], depending on the axial treatment of particles.

(a) $u_{zi0}=1\leavevmode\nobreak\ \textmd{km/s}$

(b) $u_{zi0}=10\leavevmode\nobreak\ \textmd{km/s}$

(c) $u_{zi0}=100\leavevmode\nobreak\ \textmd{km/s}$

Figure 12: Finite plasma with axial injection. Late behavior (at $t=5\mu$s) of
$E_{y}$ in the $yz$-plane, for different $u_{zi0}$.

As seen in section III.2, the stages observed in the evolution of oscillations
are greatly related to the interaction of ions with the electrostatic wave and
the distribution of ion particles in phase space. Our approach, similarly to
references [25, 29], disregard $\bm{E}_{0}$ on the ion motion to prevent the
consequent axial inhomogeneity and be able to apply periodic conditions, not
only azimuthally, but also axially. Because there is no change in velocity for
particles going through boundaries, the change in the VDF of ions and
electrons from the initial one is solely a result of the ECDI. On the other
hand, references [22, 26] account for a virtual axial dimension, so that
particles that leave the domain through axial boundaries after interacting
with the electrostatic wave are re-injected with a refreshed velocity sampled
from a prescribed VDF. This is, the refreshing axial conditions (and the
refreshing rate determined but the axial length) modify the VDF of particles
and can have a significant impact on the particle-wave interaction and overall
simulation behavior. This point is extensively analyzed in the next section.

When no virtual axial dimension or refreshing are considered, previous
literature [29, 25, 22] observe an unlimited growth of the oscillations and
heating. The reasons that explain why they do not see the quenching of short-
wavelength waves are still to be investigated, but let us speculate about a
couple of possibilities. An explanation could be that the evolution towards a
new equilibrium is not inherent to all simulations of this type and depends on
the selection of parameters (although this is the behavior we have seen in all
of our simulations). Another possibility, in line with the results in section
III.5, could be that more simulation time is needed with xenon than considered
in [29, 25, 22] to observe our late simulation stage.

##  IV THE ECDI IN AN AXIALLY-INJECTED FINITE-PLASMA

In the previous section, we pointed out the possible interaction of particle
injection with the ECDI. Fully periodic boundary conditions preserve the VDF
of the particles resultant of interacting with the instability. When there is
particle removal and injection through boundaries, particles that have already
interacted with the wave are removed from the simulation and new particles are
injected having a different VDF that could possibly modify the behavior of the
instability.

Following this reasoning, a new simulation setup is presented here that
replaces axial periodic conditions by injection ones, while periodic
conditions are kept in the azimuthal direction. Any particle leaving the
domain though axial boundaries is removed from the simulation. A constant flux
of ions $n_{0}u_{zi0}$ is injected through the left boundary, with zero
temperature and velocity $u_{zi0}\bm{1}_{z}$. Constant fluxes of electrons
$\pm n_{0}c_{e0}/\sqrt{2\pi}$, corresponding to a half-Maxwellian, are
injected through left and right boundaries, whose velocities are sampled from
a Maxwellian VDF with temperature $T_{e0}$ and velocity $u_{ye0}\bm{1}_{y}$.
The injection fluxes are chosen such that they match the amount of ions and
electrons leaving the domain in equilibrium conditions, making this
equilibrium equivalent to the fully periodic domain. In contrast with the
particle refreshing approach used in [22, 26], once the instability arises,
the injected fluxes do not need to coincide with fluxes leaving the domain and
the number of particles and electric charge are not conserved anymore. For
that reason, axial conditions on potential are also changed from periodic to
fixed potential $\phi=0$ and the finite-difference PARDISO version of the
Poisson solver is used.

(a) Energy balance

(b) Total and species energies

Figure 13: Finite plasma with axial injection. Time evolution of (a) energy
balance and (b) electrostatic (c), species and total energies; for an axial-
injection simulation with $u_{zi0}=10$ km/s. Values are given per real
particle.

###  IV.1 Dependence on the ion residence time

In section III, we showed the usual stages that we observe in the evolution of
the ECDI, being the final one the quenching of short-wavelength oscillations.
These stages are related with the ion velocity distribution and the formation
and blurring of vortex-like structures in phase space. In this sense, adding
axial boundaries that inject and absorb particles can have a major impact on
the late simulation behavior since old particles that have interacted with the
wave are removed and new, non-trapped particles are injected with the original
distribution.

In this scenario with axial injection, more similar to a Hall discharge, the
ion residence time

$\Delta t_{i0}=L_{z}/u_{zi0}$ (14)

(infinite in a fully periodic domain independently of $L_{z}$) is the key
parameter in maintaining a saturated ECDI, with sustained short-wavelength
oscillations and a nonzero axial electron transport. Three regimes will be
distinguished depending on $\Delta t_{i0}$ being much smaller, of the same
order, or much higher than the saturation time $\Delta t_{\textmd{sat}}$ of
the full-periodic configuration for the same equilibrium plasma, where
$\langle j_{ze}\rangle$ peaks.

(a) $u_{zi0}=1$ km/s

(b) $u_{zi0}=5$ km/s

(c) $u_{zi0}=10$ km/s

(d) $u_{zi0}=15$ km/s

(e) $u_{zi0}=25$ km/s

(f) $u_{zi0}=100$ km/s

Figure 14: Finite plasma with axial injection. Time-evolution of the logarithm
of normalized azimuthal Fourier coefficients for different values of
$u_{zi0}$.

We consider again a hydrogen plasma and the reference case of Table 1. The
saturation time was $\Delta t_{\textmd{sat}}\sim 0.16$-0.18 $\mu$s. We run
cases with $u_{zi0}=1$, 5, 10, 15, 20, 25 and 100 km/s, yielding, for
$L_{z}=2.679$ mm, residence times $\Delta t_{i0}=2.68$, 0.536, 0.268, 0.179,
0.134, 0.107 and 0.0268 $\mu$s, covering the three expected regimes.

The evolution of azimuthal Fourier coefficients in figure 14 and the average
current $j_{ze}$ in figure 15 give an idea of the time evolution, where the
transition from one limit regime to the other can be appreciated. The final
$E_{y}$ in the axial-azimuthal plane after $5\leavevmode\nobreak\ \mu$s for
cases $u_{zi0}=1$, 10 and 100 km/s are shown in figure 12.

The case having $u_{zi0}=1$ km/s corresponds to regime $\Delta t_{i0}\gg\Delta
t_{\textmd{sat}}$ and the final $E_{y}$ resembles the periodic case. Looking
at figure 14(a), we see that the time evolution of Fourier coefficients is
only similar at very early times when the onset of the ECDI on the initial
population; after which oscillations are mild. From the energy point of view,
the different axial boundary conditions here limit the electron heating and
total energy in the domain, which is probably limiting the energy and
contributing to the quenching of oscillations after the first ECDI onset. The
evolution of $\langle j_{ze}\rangle$ is similar to the periodic case, where
the current induced by the onset of the ECDI at early simulation times fades
away at large enough times.

On the opposite limit we have the case $u_{zi0}=100$ km/s with $\Delta
t_{i0}\ll\Delta t_{\textmd{sat}}$, where the final $E_{y}$ shows a sustained
but very weak short-wavelength oscillation. The azimuthal spectrum shows that
dominant modes concentrate around resonancea $m=1$ and, to a lesser extend,
$m=2$. In this regime ECDI modes cannot develop completely since most of the
ions leave the domain before full trapping occurs, so that we get an early
stage of the ECDI observed with periodic conditions. According to figure 15,
these mild oscillatory modes induce a weak $j_{ze}$ in the electrons, which is
very small compared with values produced by a fully developed ECDI.

In the intermediate regime $\Delta t_{i0}\sim\Delta t_{\textmd{sat}}$ we
consider cases with $u_{zi0}$ between 5 and 25 km/s (always for hydrogen). The
plasma behavior resembles the saturated behavior of periodic simulation,
because the removal and injection of particles happen at the proper rate that
keeps feeding ECDI modes and allows them to fully develop, preventing the
quenching of the oscillations observed in the regime $\Delta t_{i0}\gg\Delta
t_{\textmd{sat}}$. From these cases, those with $u_{zi0}=5$ and 25 km/s are
halfway between regimes and show features of both limit and intermediate
regimes. In every case, the electric field shows sustained oscillations with a
clear dominant mode close to the resonance $m=1$, that effectively induces a
significant $j_{ze}$ in the long term. The magnitude of the induced transport
depends on the case, being more significant, and comparable to the peak values
of the periodic case, for the cases with $u_{zi0}=10$ and 15 km/s. In these
two cases the dominance of modes close to $m=1$ does not happen at every time
and there is some kind of intermittent or pulsed behavior of long and short
wavelength modes, which results in an oscillatory evolution of the average
$j_{ze}$. This intermittency is also observed in the time evolution of figure
11 for the case $u_{zi0}=10$ km/s. For greater and smaller values of ion
velocity, $j_{ze}$ gets diminished, which is expected after our observations
of null $j_{ze}$ in upper and lower limit regimes.

From the point of view of ion particles in phase space, the final pictures for
cases $u_{zi0}=1$, 10 and 25 km/s are shown in an axial slab
$0.7<z/L_{z}<0.75$. Here, it is confirmed that only when $\Delta
t_{i0}\sim\Delta t_{\textmd{sat}}$ the vortex-like structure characteristic of
ion-wave trapping is preserved and long-term axial transport exists. The
existence of different regimes depending on the ion residence time could also
explain why, with a virtual axial length, some groups observe a transition to
an ion-acoustic mode [20] while others do not [26, 34].

The sensitivity of these results to $L_{y}$ has been tested to ensure that the
same conclusions apply to larger domains. The evolution of electron currents
are plotted in figure 17 varying $L_{y}$ while fixing $u_{zi0}=10$ km/s. As
with fully periodic simulations, the shortest case gives the most different
transient but in all cases a large $j_{ze}$ develops and they all show the
characteristics of the intermediate regime $\Delta t_{i0}\sim\Delta
t_{\textmd{sat}}$. In light of the fully periodic results, where larger
$L_{y}$ seemed to favor the formation of long domain-size modes, it is
surprising that, here, increasing $L_{y}$ mitigates the intermittency of short
and long scales and yields less oscillatory currents. A possible explanation
could be that the saturation time is slightly affected by $L_{y}$.

Previous results can be extrapolated to xenon in a Hall discharge using the
scaling laws suggested in section III.5 and a typical $\Delta t_{i0}$. Using
an average ion velocity $u_{zi0}=17$ km/s and anode-to-cathode length
$L_{z}=3.35$ cm (e.g., from results in [35]), $\Delta t_{i0}$ is estimated to
be 1.95 $\mu$s. The value of $\Delta t_{\textmd{sat}}$ for xenon mass is
estimated to be 1.83 to 2.06 $\mu s$. This similarity suggests that the
intermediate regime of the ECDI in a finite plasma could develop in a
conventional Hall discharge and supports the idea of the ECDI possibly being
an important actor in the anomalous electron transport.

###  IV.2 Energy balance

(a) Axial electron current

(b) Time-averaged $\left\langle j_{ze}\right\rangle$

Figure 15: Finite plasma with axial injection. (a) Time evolution of the
volume-averaged electron axial current $\left\langle j_{ze}\right\rangle$ for
several $u_{zi0}$ and (b) its time average for $t\geq 2$ $\mu$s.

(a) $u_{zi0}=1$ km/s

(b) $u_{zi0}=10$ km/s

(c) $u_{zi0}=25$ km/s

Figure 16: Finite plasma with axial injection. Late behavior (at $t=5\mu$s) in
the phase space plane $(y,v_{y})$ of ion particles contained in the axial slab
$0.7\leq z/L_{z}\leq 0.75$.

For the fully periodic case, equation (13) shows that the evolution of total
energy in the domain is fed by the equilibrium electric field $E_{0}$ and is
tied to the presence of an average $j_{ze}$. Therefore, an energetically
stationary state does not allow for an axial electron current. This
theoretical conclusion is retrieved in our periodic simulations.

Axial injection and removal of particles involve energy inputs and losses
through axial boundaries that have to be accounted for in the energy balance,
yielding

$\frac{\textmd{d}\mathcal{E}}{\textmd{d}t}=E_{0}\langle
j_{ze}\rangle\mathcal{V}-P_{\textmd{out}},$ (15)

where $P_{\textmd{out}}$ gathers the net energy outflow through axial
boundaries and can be computed from the energies of removed and injected
particles. The new term opens the possibility to have a balance between the
energy input by $E_{0}$ and boundary losses, allowing for an energetically
stationary behavior and an average $j_{ze}$ at the same time. Any energy loss
term, such as inelastic collisions, could play a similar role in the balance.

The different terms in the balance equation (15) are computed and represented
in figure 13 for the case with $u_{zi0}=10$ km/s, in the intermediate regime
with a net $j_{ze}$. The balance is approximately fulfilled and the energy is
close to stationary. There is, however, small changes in the total energy as
well as in the partial energies of ions and electrons. The losses introduced
by the axial boundary conditions limit enormously the heating that was
observed in the fully periodic cases without losses. Actually, the total and
species energies remain within levels close to the initial population. The
electrostatic energy is, again, a minor contribution to the total energy.

(a)

(b)

Figure 17: Finite plasma with axial injection. (a) Time evolution of the
volume-averaged electron axial current $\left\langle j_{ze}\right\rangle$ for
several $L_{y}$ and (b) its time average for $t\geq 2$ $\mu$s.

##  V Conclusions

The first part of this article is focused on the simulation of the classical
ECDI with the PIC formulation and settings as close as possible to the linear
kinetic theory. Unstable short-wavelength modes are seen to grow in the
initially homogeneous plasma that fit qualitatively well the features of the
theoretical ECDI dispersion relation in early simulation times. Close to
saturation, some parts of the Fourier spectra show similarities with ion-
acoustic modes. After saturation, the short scale modes vanish and the plasma
tends to a new equilibrium with mild or long-wavelength oscillations and much
more mode mixing. It is only during the growth and saturation of the ECDI
modes that a turbulence-based axial current is induced in the electrons. The
initial growth and saturation are related with the ion distribution that
yields the characteristic vortex-like structure in phase space. The vanishing
of oscillations seemed related with the blurring of those vortexes.

This behavior differs from the unlimited growth reported in the literature of
1D ECDI simulations with no virtual axial length [25, 22, 29], where the
quenching of short modes was not seen. It is possible that the evolution we
observe is not inherent to every ECDI simulation and depends on the choice of
parameters. Other possibility is that more simulation time is needed in these
works to reproduce the full behavior reported here, which is in line with our
parametric analysis on the ion mass.

When axial boundaries are replaced by removal/injection surfaces, the
refreshing of particles (mimicking the finite plasma in a Hall discharge) can
yield completely different behaviors. These boundary conditions imply that
particles that have interacted with the electrostatic wave are removed from
the simulation eventually and new particles with the original distribution are
injected. A key parameter is the ion injection velocity $u_{zi0}$, while its
effect is negligible in the ECDI dispersion relation and periodic simulations.
The intermediate regime where saturation and ion-residence times are similar
is the only one yielding long-term short-scale oscillations and a turbulent-
based axial current, and it is the most interesting one in the context of Hall
discharges. A rough estimation with typical magnitudes shows that the
condition for this regime is likely to happen in a conventional Hall thruster
discharge.

## Acknowledgments

This work has been supported by the R&D project PID2022-140035OB-I00 (HEEP)
funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”.
E. Bello-Benítez thanks the financial support by the European Research Council
under the European Union’s Horizon 2020 research and innovation programme
(project ZARATHUSTRA, grant agreement No 950466). A. Marín-Cebrián thanks the
financial support by MCIN/AEI/10.13039/501100011033 under grant FPU20/02901.
The authors acknowledge useful discussions on this work with J. J. Ramos. This
work is part of the PhD thesis of E. Bello-Benítez, under development, and it
will be included as a chapter within it.

## References

  * Choueiri and Ziemer [2001] E. Choueiri and J. Ziemer, Quasi-steady magnetoplasmadynamic thruster performance database, Journal of Propulsion and Power 17, 967 (2001).
  * Ellison  _et al._ [2012] C. Ellison, Y. Raitses, and N. Fisch, Cross-field electron transport induced by a rotating spoke in a cylindrical Hall thruster, Physics of Plasmas 19, 013503 (2012).
  * MacDonald  _et al._ [2011] N. A. MacDonald, M. A. Cappelli, S. R. Gildea, M. Martinez-Sanchez, and W. A. Hargus Jr, Laser-induced fluorescence velocity measurements of a diverging cusped-field thruster, Journal of Physics D: Applied Physics 44, 295203 (2011).
  * Tsikata  _et al._ [2009] S. Tsikata, N. Lemoine, V. Pisarev, and D. Gresillon, Dispersion relations of electron density fluctuations in a Hall thruster plasma, observed by collective light scattering, Physics of Plasmas 16, 033506 (2009).
  * Tsikata and Minea [2015] S. Tsikata and T. Minea, Modulated electron cyclotron drift instability in a high-power pulsed magnetron discharge, Physical Review Letters 114, 185001 (2015).
  * Tsikata  _et al._ [2017] S. Tsikata, A. Héron, and C. Honoré, Hall thruster microturbulence under conditions of modified electron wall emission, Physics of Plasmas 24, 053519 (2017).
  * Janes and Lowder [1966] G. Janes and R. Lowder, Anomalous electron diffusion and ion acceleration in a low-density plasma, Physics of Fluids 9, 1115 (1966).
  * Bello-Benítez and Ahedo [2021] E. Bello-Benítez and E. Ahedo, Axial-azimuthal, high-frequency modes from global linear-stability model of a Hall thruster, Plasma Sources Science and Technology 30, 035003 (2021).
  * Ramos  _et al._ [2021] J. J. Ramos, E. Bello-Benítez, and E. Ahedo, Local analysis of electrostatic modes in a two-fluid E x B plasma, Physics of Plasmas 28, 052115 (2021).
  * Koshkarov  _et al._ [2019] O. Koshkarov, A. Smolyakov, Y. Raitses, and I. Kaganovich, Self-organization, structures, and anomalous transport in turbulent partially magnetized plasmas with crossed electric and magnetic fields, Physical Review Letters 122, 185001 (2019).
  * Frias  _et al._ [2012] W. Frias, A. Smolyakov, I. Kaganovich, and Y. Raitses, Long wavelength gradient drift instability in Hall plasma devices. I. Fluid theory, Physics of Plasmas 19, 072112 (2012).
  * Nikitin  _et al._ [2017] V. Nikitin, D. Tomilin, A. Lovtsov, and A. Tarasov, Gradient-drift and resistive mechanisms of the anomalous electron transport in Hall effect thrusters, Europhysics Letters 117, 45001 (2017).
  * Sorokina  _et al._ [2019] E. A. Sorokina, N. A. Marusov, V. P. Lakhin, and V. I. Ilgisonis, Discharge oscillations in Morozov’s stationary plasma thruster as a manifestation of large-scale modes of gradient drift instability, Plasma Physics Reports 45, 1 (2019).
  * Romadanov  _et al._ [2016] I. Romadanov, A. Smolyakov, Y. Raitses, I. Kaganovich, T. Tian, and S. Ryzhkov, Structure of nonlocal gradient-drift instabilities in hall exb discharges, Physics of Plasmas 23, 122111 (2016).
  * Forslund  _et al._ [1970] D. Forslund, R. Morse, and C. Nielson, Electron cyclotron drift instability, Physical Review Letters 25, 1266 (1970).
  * Wong [1970] H. Wong, Electrostatic electron-ion streaming instability, Physics of Fluids 13, 757 (1970).
  * Adam  _et al._ [2004] J. Adam, A. Herón, and G. Laval, Study of stationary plasma thrusters using two-dimensional fully kinetic simulations, Physics of Plasmas 11, 295 (2004).
  * Ducrocq  _et al._ [2006] A. Ducrocq, J. Adam, A. Héron, and G. Laval, High-frequency electron drift instability in the cross-field configuration of Hall thrusters, Physics of Plasmas 13, 102111 (2006).
  * Cavalier  _et al._ [2013] J. Cavalier, N. Lemoine, G. Bonhomme, S. Tsikata, C. Honoré, and D. Grésillon, Hall thruster plasma fluctuations identified as the exb electron drift instability: Modeling and fitting on experimental data, PoP 20, 082107 (2013).
  * Lafleur  _et al._ [2016a] T. Lafleur, S. Baalrud, and P. Chabert, Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model, Physics of Plasmas 23, 053503 (2016a).
  * Janhunen  _et al._ [2018a] S. Janhunen, A. Smolyakov, D. Sydorenko, M. Jimenez, I. Kaganovich, and Y. Raitses, Evolution of the electron cyclotron drift instability in two-dimensions, Physics of Plasmas 25, 082308 (2018a).
  * Lafleur  _et al._ [2016b] T. Lafleur, S. Baalrud, and P. Chabert, Theory for the anomalous electron transport in hall effect thrusters. i. insights from particle-in-cell simulations, Physics of Plasmas 23, 053502 (2016b).
  * Lafleur and Chabert [2017] T. Lafleur and P. Chabert, The role of instability-enhanced friction on ‘anomalous’ electron and ion transport in Hall-effect thrusters, Plasma Sources Science and Technology 27, 015003 (2017).
  * Coche and Garrigues [2014] P. Coche and L. Garrigues, A two-dimensional (azimuthal-axial) particle-in-cell model of a hall thruster, Physics of Plasmas 21, 023503 (2014).
  * Janhunen  _et al._ [2018b] S. Janhunen, A. Smolyakov, O. Chapurin, D. Sydorenko, I. Kaganovich, and Y. Raitses, Nonlinear structures and anomalous transport in partially magnetized ExB plasmas, Physics of Plasmas 25, 11608 (2018b).
  * Asadi  _et al._ [2019] Z. Asadi, F. Taccogna, and M. Sharifian, Numerical study of electron cyclotron drift instability: Application to hall thruster, Frontiers in Physics 7, 10.3389/fphy.2019.00140 (2019).
  * Adam  _et al._ [2008] J. Adam, J. Boeuf, N. Dubuit, M. Dudeck, L. Garrigues, D. Gresillon, A. Heron, G. Hagelaar, V. Kulaev, N. Lemoine, _et al._ , Physics, simulation and diagnostics of Hall effect thrusters, Plasma Physics and Controlled Fusion 50, 124041 (2008).
  * Charoy  _et al._ [2019] T. Charoy, J. P. Boeuf, A. Bourdon, J. A. Carlsson, P. Chabert, B. Cuenot, D. Eremin, L. Garrigues, K. Hara, I. D. Kaganovich, A. T. Powis, A. Smolyakov, D. Sydorenko, A. Tavant, O. Vermorel, and W. Villafana, 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas, Plasma Sources Science and Technology 28, 105010 (2019).
  * Tavassoli  _et al._ [2023] A. Tavassoli, M. Papahn Zadeh, A. Smolyakov, M. Shoucri, and R. J. Spiteri, The electron cyclotron drift instability: A comparison of particle-in-cell and continuum Vlasov simulations, Physics of Plasmas 30, 10.1063/5.0134457 (2023), 033905.
  * Swanson [2003] D. Swanson,  _Plasma Waves, 2nd Edition_ (IOP Publishing, Bristol, UK, 2003).
  * Frigo and Johnson [2005] M. Frigo and S. G. Johnson, The design and implementation of FFTW3, Proceedings of the IEEE 93, 216 (2005).
  * Tavassoli  _et al._ [2021] A. Tavassoli, O. Chapurin, M. Jimenez, M. Papahn Zadeh, T. Zintel, M. Sengupta, L. Couëdel, R. J. Spiteri, M. Shoucri, and A. Smolyakov, The role of noise in pic and Vlasov simulations of the Buneman instability, Physics of Plasmas 28, 10.1063/5.0070482 (2021).
  * Chen  _et al._ [2024] L. Chen, Z.-C. Kan, W.-F. Gao, P. Duan, J.-Y. Chen, C.-Q. Tan, and Z.-J. Cui, Growth mechanism and characteristics of electron drift instability in Hall thruster with different propellant types, Chinese Physics B 33, 015203 (2024).
  * Smolyakov  _et al._ [2019] A. Smolyakov, T. Zintel, L. Couedel, D. Sydorenko, A. Umnov, E. Sorokina, and N. Marusov, Anomalous electron transport in one-dimensional electron cyclotron drift turbulence, Plasma Physics Reports 46 (2019).
  * Bello-Benítez and Ahedo [2023] E. Bello-Benítez and E. Ahedo, Stationary axial model of the Hall thruster plasma discharge: electron azimuthal inertia and far plume effects, Plasma Sources Science and Technology 32, 115011 (2023).


#  Massless Dirac Fermions in a Square Optical Lattice

Jing-Min Hou1<EMAIL_ADDRESS>Wen-Xing Yang1,3 Xiong-Jun Liu2 1Department of
Physics, Southeast University, Nanjing, 211189, China
2 Department of Physics, Texas A&M University, College Station, Texas
77843-4242, USA
3 Institute of Photonics Technologies, National Tsing-Hua University, Hsinchu
300, Taiwan

(January 5, 2009 )

###### Abstract

We propose a novel scheme to simulate and observe massless Dirac fermions with
cold atoms in a square optical lattice. A $U(1)$ adiabatic phase is created by
two laser beams for the tunneling of atoms between neighbor lattice sites.
Properly adjusting the tunneling phase, we find that the energy spectrum has
conical points in per Brillouin zone where band crossing occurs. Near these
crossing points the quasiparticles and quasiholes can be considered as
massless Dirac fermions. Furthermore, the anisotropic effects of massless
Dirac fermions are obtained in the present square lattice model. The Dirac
fermions as well as the anisotropic behaviors realizeded in our system can be
experimentally detected with the Bragg spectroscopy technique.

###### pacs:

37.10.Jk, 03.75.Ss, 05.30.Fk

##  I Introduction

Realization of two-dimensional (2D) systems of massless Dirac fermions is of
great fundamental importance, in the light of many exotic phenomena obtained
in such systems, such as zero modes, fractional statistics, unconventional
Landau levels, parity anomaly, chirality, and anomalous quantum Hall effects
Semenoff ; Jackiw ; Haldane . However, two-dimensional massless Dirac field
have not been observed untill the creation of graphene, a monolayer of
graphite Novoselov2 ; Novoselov3 . Electrons in graphene, obeying a linear
dispersion relation, behave like massless Dirac fermions Novoselov2 ;
Novoselov3 ; Zhang ; Li ; Zheng ; Gusynin ; Hou ; Jachiw2 ; Pachos .

Besides graphene, physicists also make efforts to search for other physical
systems, e.g. patterned 2D electron gases Park  and ultracold atoms in the
honeycomb optical lattice Zhu ; Zhao ; Shao ; Wu  , to simulate massless Dirac
fermions. Realization of honeycomb optical lattice opens new possibility of
studying Dirac fermions in cold atoms which provide an extremely clean
environment and controllable fashion unique access to the study of complex
physics Jaksch ; Greiner ; Lewenstein . Nevertheless, all of the above systems
require the hexagonal symmetry. Then, it is very attractive to find a system
without the hexagonal symmetry to observe massless Dirac fermions.

Ultracold atom systems provide an ideal platform to study many interesting
physics in condensed matters. To investigate the effects of gauge fields with
ultracold atoms, several schemes have been proposed to create an artificial
Abelian gauge field Jaksch2 ; Dum ; Juzeliunas1 ; Juzeliunas2 ; Juzeliunas3 ;
Gunter  or a non-Abelian gauge field Osterloh ; Ruseckas ; Lu  for neutral
atoms with laser fields. Many effects have been studied for cold atoms in an
effective gauge field, e.g., Stern-Gerlach effect for chiral moleculesLi2 ,
Double and negative reflectionJuzeliunas4 , Landau levelsJacob , spin Hall
effect Liu ; Zhu2 , induced spin-orbit couplingLiu2 , magnetic monopolePietila
, spin field effect transistorsVaishnav . Furthermore, some groups have
realized the light-induced gauge fields in experimentsDutta ; Lin .

In this paper, we propose a scheme to generate a staggered gauge field with
laser fields. A 2D square lattice model under this artificial gauge field has
a spectrum behaving like massless Dirac fermions. Furthermore, our lattice
model does not have the hexagonal symmetry. In our scheme, the energy bands of
the system exhibit degeneracy points where the conduction and valence bands
intersect. Near the these crossing points the dispersion relation is linearly
dependent on the momentum, say, is of the Dirac type. The present scheme
suggests a new direction to study Dirac fermions in the optical lattice
without the hexagonal symmetry.

Figure 1: (a) The atomic levels and the interactions between atoms and laser
fields. (b) Schematic representation of the experimental setup with the two
laser beams incident on the cloud of atoms. (c) Schematic of the square
optical lattice and the designed phase factor (denoted by arrows). (d) The
scheme of overlapping the two state-selective optical lattices.

##  II Model

We consider a system of ultracold fermionic atoms with four levels shown in
FIG.1 (a). This atomic level configuration can be experimentally realized with
alkali atom 6Li Fuchs . We choose the atomic states
$2S_{1/2}(F=1/2,m_{F}=1/2)$, $2S_{1/2}(F=3/2,m_{F}=3/2)$,
$2P_{1/2}(F=1/2,m_{F}=1/2)$ and $2P_{1/2}(F=1/2,m_{F}=-1/2)$ as
$|1\rangle,|2\rangle,|3\rangle$ and $|4\rangle$, respectively. The cold atoms
are trapped in two state-selective optical potentials as shown in FIG.1 (c)
and (d). We assume that the states $|1\rangle$ and $|2\rangle$ have the same
the state-selective optical potential, say sublattice $A$, and $|4\rangle$
only perceives the other state-selective optical potential, say sublattice
$B$. Here, for convenience, we assume that atoms in state $|3\rangle$ also
perceive sublattice $A$. However, this is unnecessary in our scheme, for the
population of the quantum state $|3\rangle$ is finally eliminated. The two
sublattices have the lattice spacings $2l_{x}$ and $l_{z}$ in the $x$ and $z$
directions, respectively. The two sublattices make up a 2D rectangular lattice
with the lattice spacings $l_{x}$ and $l_{z}$, especially a 2D square lattice
for $l_{x}=l_{z}$, when overlapping together as shown in FIG.1 (c) and (d).
Without loss of generality, we suppose that atoms with internal states
$|1\rangle$ and $|2\rangle$ are trapped in odd columns and ones with internal
states $|4\rangle$ in even columns in the whole overlapped lattice. For
convenience, we assume that the 2D square lattice considered here is in the
$x-z$ plane as shown in FIG.1 (c). Two additional laser beams along the $y$
direction are added. When the potential barrier of the optical lattice along
the $y$ direction is high enough, the tunneling along this direction between
different planes is suppressed seriously, then every layer is an independent
2D lattice in $x-z$ plane.

Using $\\{|1\rangle,|2\rangle,|3\rangle,|4\rangle\\}$ as the basis, the
Hamiltonian of free ultracold fermions in the optical lattice can be written
in the second quantized form as follows,

$\displaystyle\hat{H}_{0}$ $\displaystyle=$ $\displaystyle\int
d^{2}r\hat{\Psi}^{\dagger}\left(-\frac{\hbar^{2}}{2m}\nabla^{2}+V({\bf
r})\right)\hat{\Psi},$ (1)

where
$\hat{\Psi}^{\dagger}=(\hat{\Psi}_{1}^{\dagger},\hat{\Psi}_{2}^{\dagger},\hat{\Psi}_{3}^{\dagger},\hat{\Psi}_{4}^{\dagger})$
and
$\hat{\Psi}=(\hat{\Psi}_{1},\hat{\Psi}_{2},\hat{\Psi}_{3},\hat{\Psi}_{4})^{T}$
($T$ denotes the matrix transposition) with $\hat{\Psi}_{i}({\bf r})$ and
$\hat{\Psi}_{i}^{\dagger}({\bf r})$ being field operators corresponding to
annihilating and creating an atom with the internal quantum state $|i\rangle\
(i=1,2,3,4)$ at coordinate position ${\bf r}$ respectively. Here, $V({\bf r})$
is the trap potential matrix as

$\displaystyle V({\bf r})=\left(\matrix{V_{A}({\bf r})&0&0&0\cr 0&V_{A}({\bf
r})&0&0\cr 0&0&V_{A}({\bf r})&0\cr 0&0&0&V_{B}({\bf r})}\right),$ (2)

where $V_{X}({\bf r})(X=A,B)$ are the two state-selective periodic potentials.
The ground state $|1\rangle$ is coupled to the excited state $|3\rangle$ via a
laser field with the corresponding Rabi frequency $\Omega_{1}e^{-iq_{1}z}$ and
the state $|2\rangle$ is coupled to the excited state $|3\rangle$ via a laser
field with the corresponding Rabi frequency $\Omega_{2}e^{iq_{2}z}$ as shown
in FIG.1 (a) and (b). The corresponding light-atom interaction Hamiltonian is,

$\displaystyle\hat{H}_{1}$ $\displaystyle=$ $\displaystyle\int
d^{2}r\hat{\Psi}^{\dagger}M\hat{\Psi},$ (3)

with

$\displaystyle M=\hbar\left(\matrix{0&0&\Omega_{1}e^{iq_{1}z}&0\cr
0&0&\Omega_{2}e^{-iq_{2}z}&0\cr\Omega_{1}e^{-iq_{1}z}&\Omega_{2}e^{iq_{2}z}&0&0\cr
0&0&0&0}\right),$ (4)

where $\Omega_{j}(j=1,2)$ are the Rabi frequencies. Additionally, the quantum
state $|1\rangle$ is coupled to the quantum state $|4\rangle$ via a laser
field propagating in the $y$ direction with Rabi frequency
$\Omega_{3}e^{iq_{3}y}$. Because $e^{iq_{3}y}$ is a constant in the $x-z$
plane, we can omit this phase factor by supposing the two-dimensional lattice
on the $y=0$ plane. The corresponding interaction Hamiltonian is

$\displaystyle\hat{H}_{2}=\int d^{2}r\hat{\Psi}^{\dagger}N\hat{\Psi},$ (5)

with

$\displaystyle N=\hbar\left(\matrix{0&0&0&\Omega_{3}\cr 0&0&0&0\cr
0&0&0&0\cr\Omega_{3}&0&0&0}\right).$ (6)

The total Hamiltonian can be written as
$\hat{H}=\hat{H}_{0}+\hat{H}_{1}+\hat{H}_{2}$.

The Hamiltonian (3) can be diagonalized by the matrix,

$\displaystyle U=\left(\matrix{\cos\theta&-\sin\theta
e^{iqz}&0&0\cr\frac{\sqrt{2}}{2}\sin\theta
e^{-iqz}&\frac{\sqrt{2}}{2}\cos\theta&-\frac{\sqrt{2}}{2}e^{-iq_{2}z}&0\cr\frac{\sqrt{2}}{2}\sin\theta
e^{-iqz}&\frac{\sqrt{2}}{2}\cos\theta&\frac{\sqrt{2}}{2}e^{-iq_{2}z}&0\cr
0&0&0&1}\right),$ (7)

where $q=q_{1}+q_{2}$ and $\tan\theta=|\Omega_{1}|/|\Omega_{2}|$.
Correspondingly, we obtain the dressed states as

$\displaystyle|\chi_{1}\rangle$ $\displaystyle=$
$\displaystyle\cos\theta|1\rangle-\sin\theta e^{iqz}|2\rangle,$ (8)
$\displaystyle|\chi_{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}}{2}\sin\theta
e^{-iqz}|1\rangle+\frac{\sqrt{2}}{2}\cos\theta|2\rangle-\frac{\sqrt{2}}{2}e^{-iq_{2}z}|3\rangle,$
(9) $\displaystyle|\chi_{3}\rangle$ $\displaystyle=$
$\displaystyle\frac{\sqrt{2}}{2}\sin\theta
e^{-iqz}|1\rangle+\frac{\sqrt{2}}{2}\cos\theta|2\rangle+\frac{\sqrt{2}}{2}e^{-iq_{2}z}|3\rangle,$
(10) $\displaystyle|\chi_{4}\rangle$ $\displaystyle=$
$\displaystyle|4\rangle,$ (11)

with the energy eigenvalues $E_{i}=(0,-\hbar\Omega,\hbar\Omega,0)$ with
$\Omega=\sqrt{|\Omega_{1}|^{2}+|\Omega_{2}|^{2}}$. Here, the state
$|\chi_{1}\rangle$ is a so-called dark state, which does not contain the
component of the excited atomic state $|3\rangle$, and
$|\chi_{2}\rangle,|\chi_{3}\rangle$ are bright states. In the dressed state
basis
$\\{|\chi_{1}\rangle,|\chi_{2}\rangle,|\chi_{3}\rangle,|\chi_{4}\rangle\\}$,
the vector field operator can be written as
$\hat{\Phi}=({\hat{\Phi}_{1},\hat{\Phi}_{2}},\hat{\Phi}_{3},\hat{\Phi}_{4})^{T}=U({\hat{\Psi}_{1},\hat{\Psi}_{2}},\hat{\Psi}_{3},\hat{\Psi}_{4})^{T}$,
where $\hat{\Phi}_{j}(j=1,2,3,4)$ represent destructing an atom in the dressed
state $|\chi_{j}\rangle(j=1,2,3,4)$. Thus, the Hamiltonian can be rewritten as

$\displaystyle\hat{H}=\int
d^{2}r\hat{\Phi}^{\dagger}\left[\frac{1}{2m}(-i\hbar\nabla-\tilde{\bf
A})^{2}+\tilde{V}({\bf r})+\tilde{N}\right]\hat{\Phi},$ (12)

where $\tilde{\bf A}=i\hbar U\nabla U^{\dagger}$, $\tilde{V}({\bf r})=UV({\bf
r})U^{\dagger}+UMU^{\dagger}+\frac{\hbar^{2}}{2m}[(U\nabla
U^{\dagger})^{2}+\nabla U\cdot\nabla U^{\dagger}]$ and
$\tilde{N}=UNU^{\dagger}$. We straightforwardly calculate these matrices and
obtain,

$\displaystyle\tilde{\bf A}$ $\displaystyle=$ $\displaystyle-\hbar{\bf
e}_{z}\left(\matrix{-q\sin^{2}\theta&\frac{\sqrt{2}}{2}q\sin\theta\cos\theta
e^{iqz}&\frac{\sqrt{2}}{2}q\sin\theta\cos\theta
e^{iqz}&0\cr\frac{\sqrt{2}}{2}q\sin\theta\cos\theta
e^{-iqz}&\frac{1}{2}q\sin^{2}\theta+\frac{1}{2}q_{2}&\frac{1}{2}q\sin^{2}\theta-\frac{1}{2}q_{2}&0\cr\frac{\sqrt{2}}{2}q\sin\theta\cos\theta
e^{-iqz}&\frac{1}{2}q\sin^{2}\theta-\frac{1}{2}q_{2}&\frac{1}{2}q\sin^{2}\theta+\frac{1}{2}q_{2}&0\cr
0&0&0&0}\right),$ (13)

and

$\displaystyle\tilde{V}({\bf r})$ $\displaystyle=$
$\displaystyle\left(\matrix{V_{A}({\bf r})&0&0&0\cr 0&V_{A}({\bf
r})-\hbar\Omega&0&0\cr 0&0&V_{A}({\bf r})+\hbar\Omega&0\cr 0&0&0&V_{B}({\bf
r})}\right),$ (14)

and

$\displaystyle\tilde{N}$ $\displaystyle=$
$\displaystyle\hbar\left(\matrix{0&0&0&\Omega_{3}\cos\theta\cr
0&0&0&\frac{\sqrt{2}}{2}\Omega_{3}\sin\theta e^{-iqz}\cr
0&0&0&\frac{\sqrt{2}}{2}\Omega_{3}\sin\theta
e^{-iqz}\cr\Omega_{3}\cos\theta&\frac{\sqrt{2}}{2}\Omega_{3}\sin\theta
e^{iqz}&\frac{\sqrt{2}}{2}\Omega_{3}\sin\theta e^{iqz}&0}\right).$ (15)

In our scheme, we only consider the atoms in the dressed states
$|\chi_{1}\rangle$ and $|\chi_{4}\rangle$. Thus, we have to adiabatically
eliminate the populations of the dressed states $|\chi_{2}\rangle$ and
$|\chi_{3}\rangle$ and to avoid the atoms decaying into these two dressed
states. This can be realized in the steps. First, we start with the atoms in
the atomic state $|1\rangle$ and $\Omega_{1}=0$, $\Omega_{3}=0$ with
$\Omega_{2}$ finite, then slowly turn $\Omega_{1}$, we will end up with the
atoms in the dressed state $|\chi_{1}\rangle$ Scully . During this process,
the variation of $\Omega_{1}$ is slow enough to satisfy the adiabatic
condition $|\langle\chi_{j}|\partial/\partial
t|\chi_{1}\rangle|\ll|E_{j}-E_{1}|/\hbar=\Omega$ with $j=2,3$ Messiah ; jmhou
. In the second step, we adiabatically turn the Rabi frequency $\Omega_{3}$
on, we will end up with atoms in the dressed states $|\chi_{1}\rangle$ and
$|\chi_{4}\rangle$. To avoid the atoms decaying into the dressed states
$|\chi_{2}\rangle$ and $|\chi_{3}\rangle$, the adiabatic conditions
$\frac{1}{2m}|\tilde{\bf
A}_{j1}|^{2}=\frac{\hbar^{2}}{2m}q^{2}\sin^{2}\theta\cos^{2}\theta\ll|E_{j}-E_{1}|=\hbar\Omega$
and
$|\tilde{N}_{j4}|=\frac{\sqrt{2}}{2}\hbar\Omega_{3}\sin\theta\ll|E_{j}-E_{4}|=\hbar\Omega$
for $j=2,3$ are satisfied. This is to say, the off-diagonal elements of the
Hamiltonian are small enough to avoid the atoms decaying into the dressed
states $|\chi_{2}\rangle$ and $|\chi_{3}\rangle$.

Since the atoms are only in the dressed states $|\chi_{1}\rangle$ and
$|\chi_{4}\rangle$, we consider the reduced space with the dressed state basis
$\\{|\chi_{1}\rangle,|\chi_{4}\rangle\\}$. Therefore, the total Hamiltonian
can be reduced to

$\displaystyle\hat{H}$ $\displaystyle=$ $\displaystyle\int
d^{2}r\hat{\Phi}_{1}^{\dagger}\left[\frac{1}{2m}(-i\hbar\nabla-{\bf{A}})^{2}+{V}_{A}({\bf
r})\right]\hat{\Phi}_{1}$ (16) $\displaystyle+$ $\displaystyle\int
d^{2}r\hat{\Phi}_{4}^{\dagger}\left[-\frac{\hbar^{2}}{2m}\nabla^{2}+{V}_{B}({\bf
r})\right]\hat{\Phi}_{4}$ $\displaystyle+$ $\displaystyle\hbar\Omega_{e}\int
d^{2}r\left(\hat{\Phi}_{4}^{\dagger}\hat{\Phi}_{1}+\hat{\Phi}_{1}^{\dagger}\hat{\Phi}_{4}\right),$

where $\Omega_{e}=\tilde{N}_{14}/\hbar=\Omega_{3}\cos\theta$ and the $U(1)$
adiabatic gauge potential ${\bf A}=\tilde{\bf A}_{11}=\hbar
q\sin^{2}\theta{\bf e}_{z}$.

##  III Massless Dirac fermions

Taking the tight-binding limit, we can superpose the Bloch states to get
Wannier functions $w_{a}({\bf r}-{\bf r}_{i})$ and $w_{b}({\bf r}-{\bf
r}_{j})$ for sublattice $A$ and $B$, respectively. In the present case, we can
expand the field operator in the lowest band Wannier functions as,
$\hat{\Phi}_{1}({\bf
r})=\sum_{m(odd),n}\hat{a}_{m,n}e^{\frac{i}{\hbar}\int_{0}^{{\bf r}_{mn}}{\bf
A}\cdot d{\bf r}}w_{a}({\bf r}-{\bf r}_{mn})$ and $\hat{\Phi}_{4}({\bf
r})=\sum_{m(even),n}\hat{b}_{m,n}w_{b}({\bf r}-{\bf r}_{mn})$. Substituting
the above expression into Eq.(16), we can rewrite the Hamiltonian as follows,

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle-\sum_{(m(odd),n)}[t_{b}\hat{b}^{\dagger}_{m+1,n+1}\hat{b}_{m+1}+t_{a}e^{i\gamma}\hat{a}^{\dagger}_{m,n+1}\hat{a}_{m,n}$
(17) $\displaystyle+2t_{1}\hat{a}^{\dagger}_{m,n}\hat{b}_{m+1,n}+{\rm
H.c.}]+\hat{H}_{K},$

with
$H_{K}=\epsilon_{a}\sum_{(m(odd),n)}\hat{a}_{m,n}^{\dagger}\hat{a}_{m,n}+\epsilon_{b}\sum_{(m(even),n)}\hat{b}_{m,n}^{\dagger}\hat{b}_{m,n}$.
Here, the parameters have the following forms: $t_{a}=\int
d^{2}rw^{*}_{a}({\bf r}-{\bf
r}_{m,n+1})(-\hbar^{2}\nabla^{2}/2m+V_{A})w_{a}({\bf r}-{\bf r}_{mn})$,
$t_{b}=\int d^{2}rw^{*}_{b}({\bf r}-{\bf
r}_{m,n+1})(-\hbar^{2}\nabla^{2}/2m+V_{B})w_{b}({\bf r}-{\bf r}_{mn})$,
$t_{1}=\Omega_{e}\int d^{2}rw^{*}_{b}({\bf r}-{\bf r}_{m+1,n})w_{a}({\bf
r}-{\bf r}_{mn})$, $\epsilon_{a}=\int d^{2}rw^{*}_{a}({\bf r}-{\bf
r}_{m,n})(-\hbar^{2}\nabla^{2}/2m+V_{A})w_{a}({\bf r}-{\bf r}_{mn})$,
$\epsilon_{b}=\int d^{2}rw^{*}_{b}({\bf r}-{\bf
r}_{m,n})(-\hbar^{2}\nabla^{2}/2m+V_{B})w_{b}({\bf r}-{\bf r}_{mn})$ and
$\gamma=2\pi\sin^{2}\theta\hbar ql_{z}$ is the phase resulted from the
adiabatic gauge potential. In our scheme, we consider
$\epsilon_{a}=\epsilon_{b}$, so $\hat{H}_{K}$ in Eq. (17) can be dropped out
as a constant term, which does not affect the physics considered here.

First, we consider that the ideal conditions $\gamma=\pi$,
$t_{a}=t_{b}=t_{1}=t$ and $l_{x}=l_{z}=l$ are satisfied. In experiments, these
conditions can be achieved. Taking the Fourier transformation, $\hat{a}({\bf
k})=\sum_{(m(odd),n)}\hat{a}_{m,n}\exp(-i{\bf k}\cdot{\bf r}_{m,n})$ and
$\hat{b}({\bf k})=\sum_{(m(even),n)}\hat{b}_{m,n}\exp(-i{\bf k}\cdot{\bf
r}_{m,n})$, we obtain the total Hamiltonian as

$\displaystyle\hat{H}$ $\displaystyle=$
$\displaystyle-2t\sum_{k}[\cos(k_{z}l)\hat{b}^{\dagger}({\bf k})\hat{b}({\bf
k})-\cos(k_{z}l)\hat{a}^{\dagger}({\bf k})\hat{a}({\bf k})$ (18)
$\displaystyle+\cos(k_{x}l)\hat{a}^{\dagger}({\bf k})\hat{b}({\bf
k})+\cos(k_{x}l)\hat{b}^{\dagger}({\bf k})\hat{a}({\bf k})].$

Diagonalizing the above Hamiltonian (18), we obtain the quasiparticle energy
spectrum $E({\bf k})=2st\sqrt{\cos^{2}(k_{x}l)+\cos^{2}(k_{z}l)}$ with $s=\pm
1$ being the band index, which is similar to the spectrum of $\pi$ flux states
in quantum spin liquids Wen . This energy spectrum has two energy bands and
contains four zero-energy Dirac points, where the conduction and valence bands
intersect, in the first Brillouin zone at ${\bf
K}_{1}=\left({\pi}/{2l},{\pi}/{2l}\right),{\bf
K}_{2}=\left(-{\pi}/{2l},{\pi}/{2l}\right),{\bf
K}_{3}=\left({-\pi}/{2l},-{\pi}{2l}\right),{\bf
K}_{4}=\left({\pi}/{2l},-{\pi}/{2l}\right)$. Near the Dirac points, the energy
dispersion has standard cone-like shape as shown in Fig.2 (a) and (b) and the
spectrum is linear. The low-energy state dynamics are described by linearizing
their spectrum about the degeneracy points and are modeled by massless
relativistic fermions.

For simplicity, we only consider the part around the Dirac points ${\bf
K}_{1}$, and the physics around the other Dirac points are similar. Setting
${\bf k}={\bf K}_{1}+{\bf p}$, we linearize the Hamiltonian around the Dirac
point ${\bf K}_{1}$ as, $\hat{H}=\hbar
v_{0}\sum_{p}[p_{z}\hat{b}^{\dagger}({\bf p})\hat{b}({\bf
p})-p_{z}\hat{a}^{\dagger}({\bf p})\hat{a}({\bf
p})+p_{x}\hat{a}^{\dagger}({\bf p})\hat{b}({\bf
p})+p_{x}\hat{b}^{\dagger}({\bf p})\hat{a}({\bf p})]$ with $v_{0}=2tl/\hbar$,
which can be rewritten in coordinate space as, $\hat{H}=\int
d^{2}r\hat{\eta}^{\dagger}({\bf r})\hat{\cal H}\hat{\eta}({\bf r})$, where
$\hat{\eta}=(\hat{\eta}_{b},\hat{\eta}_{a})^{T}$ with $\hat{\eta}_{b}({\bf
r})=\int d^{2}pe^{-i{\bf p}\cdot{\bf r}}\hat{b}({\bf p})$ and
$\hat{\eta}_{a}({\bf r})=\int d^{2}pe^{-i{\bf p}\cdot{\bf r}}\hat{a}({\bf
p})$. Here, $\hat{\cal H}$ is the single-particle Hamiltonian as $\hat{\cal
H}=\hbar v_{0}(\hat{p}_{x}\sigma_{x}+\hat{p}_{z}\sigma_{z})$, where
$\sigma_{x}$ and $\sigma_{z}$ are Pauli matrixes. We obtain the eigenstates

$\displaystyle\phi_{\bf
p}^{s}=\frac{1}{\sqrt{2}}\left(\matrix{\cos\frac{\alpha}{2}+s\sin\frac{\alpha}{2}\cr
s\cos\frac{\alpha}{2}-\sin\frac{\alpha}{2}}\right)e^{{i}{\bf p}\cdot{\bf r}},$
(19)

where $s=\pm 1$ and $\tan\alpha=p_{z}/p_{x}$. The corresponding eigenenergies
are $E^{s}({\bf p})=s\hbar v_{0}p$ with $p=\sqrt{p_{x}^{2}+p_{z}^{2}}$. When
the wave vector is ${\bf p}=p_{x}\hat{x}+p_{z}\hat{z}$, the corresponding
group velocity and pseudospin vector are ${\bf
v}_{g}=sv_{0}(p_{x}\hat{x}+p_{z}\hat{z})/p$ and ${\bf
c}=(p_{x}\hat{x}+p_{z}\hat{z})/p$, respectively. It is easy to find that the
three vectors ${\bf v}_{g}$, ${\bf c}$ and ${\bf p}$ are collinear, i.e., they
are parallel to each other. There is an intimate relation between the
pseudospin and motion of the quasiparticle or quasihole: pseudospin can only
be directed along the propagation direction (say, for quasiparticles) or only
opposite to it (for quasiholes). As a result, quasiparticles or quasiholes
exhibit a linear dispersion relation $E=\hbar v_{0}k$, as if they were
massless relativistic particles but the role of the speed of light is played
here by the Fermi velocity $v_{0}$.

Figure 2: Energy dispersion for cold fermionic atoms in a square optical
lattice. (a) shows the energy dispersion and (b) represents the profiles of
the energy dispersion with $k_{z}=\pi/2l$ (blue line) and $k_{x}=\pi/2l$(red
star), for the ideal case $t_{a}=t_{b}=t_{1}=t$, $l_{x}=l_{z}=l$,
$\gamma=\pi$. (c) shows the energy dispersion and (d) represents the profiles
of the energy dispersion with $k_{z}=\pi/2l$ (blue line) and $k_{x}=\pi/2l$
(red line), for the anisotropic case $t_{a}=t_{b}=2t_{1}/3=2t/3$,
$l_{x}=l_{z}=l$, $\gamma=\pi$. (e) shows the energy dispersion and (f)
represents the profiles of the energy dispersion with $k_{z}=2\pi/5l$ (blue
line) and $k_{x}=\pi/2l$(red line), for the anisotropic case
$t_{a}=t_{b}=t_{1}=t$, $l_{z}=5l_{x}/4=5l/4$, $\gamma=\pi$. (g) shows the
energy dispersion and (h) represents the profiles of the energy dispersion
with $k_{z}=2\pi/5l$ (blue line) and $k_{x}=\pi/2l$ (red line), for the
anisotropic case with $t_{a}=t_{b}=t_{1}=t$, $l_{x}=l_{z}=l$, $\gamma=6\pi/5$.

In practice, the parameters may have fluctuations around the ideal conditions
considered above. Fortunately, even the parameters deviate from the ideal
ones, the massless Dirac fermion spectrum persists and remarkably exhibit
anisotropic behaviors, which are just pursued in References Park1 by adding
external periodic potentials on graphene. Here, we provide alternative methods
to exhibit anisotropic behaviors of massless Dirac fermions in a square
optical lattice by setting the parameters deviated from the ideal situation.

For simplicity, we only consider three cases with the existence of parameter
deviation from the ideal situation as follow: (i) $t_{a}=t_{b}\neq t_{1}=t$,
$l_{z}=l_{x}=l$, $\gamma=\pi$; (ii) $t_{a}=t_{b}=t_{1}=t$, $l_{z}\neq
l_{x}=l$, $\gamma=\pi$; (iii) $t_{a}=t_{b}=t_{1}=t$, $l_{x}=l_{z}=l$,
$\gamma=\pi+\delta$ with $\delta\neq 0$. The corresponding dispersion
relations are $E_{\rm i}({\bf
k})=2st\sqrt{\cos^{2}(k_{x}l)+(t_{a}/t)^{2}\cos^{2}(k_{z}l)}$, $E_{\rm
ii}({\bf k})=2st\sqrt{\cos^{2}(k_{x}l_{x})+\cos^{2}(k_{z}l_{z})}$ and $E_{\rm
iii}({\bf
k})=t[\cos(k_{z}l+\delta)-\cos(k_{z}l)]+st\sqrt{4\cos^{2}(k_{x}l)+[\cos(k_{z}l+\delta)+\cos(k_{z}l)]^{2}}$
for cases (i), (ii) and (iii), respectively, which are shown in FIG.2 (c)-(h).
For case (ii), the four Dirac points are ${\bf
K}_{1}=\left({\pi}/{2l_{x}},{\pi}/{2l_{z}}\right),{\bf
K}_{2}=\left(-{\pi}/{2l_{x}},{\pi}/{2l_{z}}\right),{\bf
K}_{3}=\left({-\pi}/{2l_{x}},-{\pi}/{2l_{z}}\right),{\bf
K}_{4}=\left({\pi}/{2l_{x}},-{\pi}/{2l_{z}}\right)$, which are dependent on
the lattice spacing in the $x$ and $z$ direction, while the Dirac points for
cases (i) are the same as those of the ideal case. For case (iii), the four
Dirac points are ${\bf K}_{1}=\left({(\pi-\delta)}/{2l},{\pi}/{2l}\right),{\bf
K}_{2}=\left({(-\pi-\delta)}/{2l},{\pi}/{2l}\right),{\bf
K}_{3}=\left({(-\pi-\delta)}/{2l},-{\pi}/{2l}\right),{\bf
K}_{4}=\left({(\pi-\delta)}/{2l},-{\pi}/{2l}\right)$. Around the Dirac points,
these spectra can be linearized as $E_{\rm i}^{s}({\bf p})=s\hbar v_{0}p_{1}$
with $p_{1}=\sqrt{p_{x}^{2}+f_{1}^{2}p_{z}^{2}}$, $E_{\rm ii}^{s}({\bf
p})=s\hbar v_{0}p_{1}$ with $p_{2}=\sqrt{p_{x}^{2}+f_{2}^{2}p_{z}^{2}}$ and
$E_{\rm iii}^{s}({\bf p})=\pm 2t\sin(\delta/2)+s\hbar v_{0}p_{3}$ with
$p_{3}=\sqrt{p_{x}^{2}+f_{3}^{2}p_{z}^{2}}$, where $f_{1}=t_{a}/t=t_{b}/t$,
$f_{2}=l_{z}/l$ and $f_{3}=\cos(\delta/2)$. The corresponding single-particle
Hamiltonian can be written as $\hat{\cal H}_{\rm i}=\hbar
v_{0}(\hat{p}_{x}\sigma_{x}+f_{1}\hat{p}_{z}\sigma_{z})$, $\hat{\cal H}_{\rm
ii}=\hbar v_{0}(\hat{p}_{x}\sigma_{x}+f_{2}\hat{p}_{z}\sigma_{z})$ and
$\hat{\cal H}_{\rm iii}=\pm 2t\sin(\delta/2)+\hbar
v_{0}(\hat{p}_{x}\sigma_{x}+f_{3}\hat{p}_{z}\sigma_{z})$ for cases
(i),(ii),(iii), respectively. In all cases, the quasiparticles or quasiholes
are still massless Dirac fermions and show chiral behavior. For the wave
vector ${\bf p}=p_{x}\hat{x}+p_{z}\hat{z}$, the group velocity are pseudospin
vector are ${\bf v}_{g}=sv_{t}(p_{x}\hat{x}+f_{j}^{2}p_{z}\hat{z})/p_{j}$ and
${\bf c}=(p_{x}\hat{x}+f_{j}p_{z}\hat{z})/p_{j}$ for $j=1,2,3$. Here, the
three vectors ${\bf v}_{g},{\bf c}$ and ${\bf p}$ are not collinear and the
dispersion relations near the Dirac points show anisotropic behaviors.

Figure 3: The Bragg spectroscopies for the ideal case, (a)
$t_{a}=t_{b}=t_{1}=t$, $l_{x}=l_{z}=l$, $\gamma=\pi$, and the anisotropic
cases, (b) $t_{a}=t_{b}=2t_{1}/3=2t/3$, $l_{x}=l_{z}=l$, $\gamma=\pi$, (c)
$t_{a}=t_{b}=t_{1}=t$, $l_{z}=5l_{x}/4=5l/4$, $\gamma=\pi$, (d)
$t_{a}=t_{b}=t_{1}=t$, $l_{x}=l_{z}=l$, $\gamma=6\pi/5$ and $q=\pi/10l$. Here,
we represent the Bragg spectroscopies with blue lines for the mentum
difference ${\bf q}$ in the $x$ direction and with red stars or red lines for
${\bf q}$ in the $z$ direction. The frequency difference $\omega$ is expressed
in units of $qv_{0}$ and the dynamic structure factor $S({\bf q},\omega)$ is
expressed in units of $q/8\pi^{2}nv_{0}$ with $n$ being the number density of
atoms in the system.

##  IV Bragg spectroscopy

Here, we propose to identify massless Dirac fermionic quasiparticles with
Bragg spectroscopy Stamper-Kurn , which is extensively used to probe
excitation spectra in condensed matter physics. In Bragg scattering, the
atomic gas is exposed to two laser beams, with wavevectors ${\bf k}_{1}$ and
${\bf k}_{2}$ and a frequency difference $\omega$. The light-atom interaction
Hamiltonian for Bragg scattering can be written as, $\hat{H}_{B}=\sum_{{\bf
p}_{1},{\bf p}_{2}}\hbar\Omega_{B}e^{-i{\bf q}\cdot{\bf r}}|\phi_{{\bf
p}_{2}}^{f}\rangle\langle\phi_{{\bf p}_{1}}^{i}|+{\rm H.c.}$ with ${\bf
q}={\bf p}_{2}-{\bf p}_{1}$, where the initial state $|\phi_{{\bf
p}}^{i}\rangle$ is a filled state under Fermi surface and the final state
$|\phi_{{\bf p}}^{f}\rangle$ is an empty state above Fermi surface. From the
Fermi’s golden rule, we obtain the dynamic structure factor as follows,

$\displaystyle S({\bf q},\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{N\hbar^{2}\Omega_{B}^{2}}\sum_{{\bf
p}}|\langle\phi_{{\bf p}+{\bf q}}^{f}|\hat{H}_{B}|\phi_{{\bf
p}}^{i}\rangle|^{2}$ (20) $\displaystyle\times\delta(\hbar\omega-E_{{\bf
p}+{\bf q}}^{f}+E_{{\bf p}}^{i}),$

where $N$ is the total number of atoms in the system.

Here, we consider the case of half filling of cold fermions in the optical
lattice, i.e. the Fermi energy surface is at zero energy level, which is just
at Dirac points for the cases except anisotropic case (iii). The numerical
evaluation results of the dynamic structure factor is shown in FIG. 3. We note
that there are lower cutoff frequencies $\omega_{r}$ for the fixed momentum
difference $q$ and $S({\bf q},\omega)$ are approximately linear to the
frequency difference $\omega$ for large frequency difference $\omega$ for
FIG.3 (a),(b), (c). However, for FIG.3(d), the cutoff disappears for the Fermi
surface is not at Dirac points for this case. FIG.3 (a) show that the bragg
spectroscopy curves for the momentum difference ${\bf q}$ in the $x$ and $z$
directions are identical in the ideal case, which is just a consequence of the
isotropy of the energy spectrum. From FIG.3 (b), (c) and (d), we clearly see
that the Bragg spectroscopies are different for ${\bf q}$ in the $x$ and $z$
directions in the anisotropic cases, which features the anisotropic behaviors
of those spectra.

##  V Conclusion

In summary, we have proposed a novel scheme to realize massless Dirac fermions
in a 2D square optical lattice with assistance of laser fields. For massless
Dirac fermions, the gap is zero and the linear dispersion law holds. Our
scheme is very robust against perturbations. Even the experimental situation
deviates from the ideal conditions, massless Dirac fermions persist and,
furthermore, exhibit novel features, i.e., anisotropic behaviors. Due to the
absence of hexagonal symmetry, our scheme suggests a new direction to study
Dirac fermions in the optical lattice.

###### Acknowledgements.

This work was supported by the Teaching and Research Foundation for the
Outstanding Young Faculty of Southeast University. X. J. Liu acknowledges
support from US NSF Grant No. DMR-0547875 and ONR under Grant No.
ONR-N000140610122.

## References

  * (1) G. W. Semenoff, Phys. Rev. Lett. 53, 2449 (1984).
  * (2) R. Jackiw, Phys. Rev. D 29, 2375 (1984).
  * (3) F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
  * (4) K. S. Novoselov, et al., Science 306, 666 (2004).
  * (5) K. S. Novoselov et al., Nature 438, 197 (2005).
  * (6) Y. Zhang, et al., Nature 438, 201 (2005).
  * (7) G. Li and E. Y. Andrei, Nature Phys. 3, 623 (2007).
  * (8) Y. Zheng and T. Ando, Phys. Rev. B 65, 245420 (2002).
  * (9) V. P. Gusynin et al., Phys. Rev. Lett. 95, 146801 (2005).
  * (10) C. Y. Hou et al., Phys. Rev. Lett. 98, 186809 (2007).
  * (11) R. Jackiw et al., Phys. Rev. Lett. 98, 266402 (2007).
  * (12) J. K. Pachos et al., Int. J. Mod. Phys. B, 21, 5113 (2007).
  * (13) P. H. Park and S. G. Louie, arXiv:0808.2127.
  * (14) S. L. Zhu et al., Phys. Rev. Lett. 98, 260402 (2007).
  * (15) E. Zhao et al., Phys. Rev. Lett. 97, 230404 (2006).
  * (16) L. B. Shao et al., arXiv:0804.1850.
  * (17) C. Wu et al., Phys. Rev. Lett. 99, 070401 (2007); C. Wu and S. Das Sarma, Phys. Rev. B 77, 235107 (2008).
  * (18) D. Jaksch et al., Phys. Rev. Lett. 81, 3108(1998).
  * (19) M. Greiner et al., Nature 415, 39, (2002).
  * (20) M. Lewenstein et al., Adv. Phys. 56, 243 (2007), and references therein.
  * (21) D. Jaksch and P. Zoller, New J. Phys. 5, 56 (2003).
  * (22) G. Juzeliūnas and P. Öhberg, Phys. Rev. Lett. 93, 033602 (2004).
  * (23) G. Juzeliūnas, P. Öhberg, J. Ruseckas, and A. Klein, Phys. Rev. A 71, 053614 (2005).
  * (24) G. Juzeliūnas, J. Ruseckas, P. Öhberg, and M. Fleischhauer, Phys. Rev. A 73, 025602 (2006).
  * (25) K. J. Günter, M. Cheneau, T. Yefsah, S. P. Rath, and J. Dalibard, Phys. Rev. A 79, 011604(R) (2009).
  * (26) R. Dum and M. Olshanii, Phys. Rev. Lett. 76, 1788 (1996).
  * (27) K. Osterloh, M. Baig, L. Santos, P. Zoller, and M. Lewenstein, Phys. Rev. Lett. 95, 010403 (2005).
  * (28) J. Ruseckas, G. Juzeliūnas, P. Öhberg, and M. Fleischhauer, Phys. Rev. Lett. 95, 010404 (2005).
  * (29) L. H. Lu and Y. Q. Li, Phys. Rev. A 76, 023410 (2007).
  * (30) Y. Li, C. Bruder, and C. P. Sun, Phys. Rev. Lett. 99, 130403 (2007).
  * (31) G. Juzeliūnas, J. Ruseckas, A. Jacob, L. Santos, and P. Öhberg, Phys. Rev. Lett. 100, 200405 (2008).
  * (32) A. Jacob, P. Öhberg, G. Juzeliūnas, and L. Santos, New J. Phys. 10, 045022(2008)
  * (33) X. J. Liu, X. Liu, L. C. Kwek, and C. H. Oh, Phys. Rev. Lett. 98, 026602 (2007).
  * (34) S. L. Zhu, H. Fu, C. J. Wu, S. C. Zhang, and L. M. Duan, Phys. Rev. Lett. 97, 240401 (2006).
  * (35) X. J. Liu, M. F. Borunda, X. Liu, and J. Sinova, Phys. Rev. Lett. 102, 046402 (2009).
  * (36) V. Pietilä and M. Möttönen, Phys. Rev. Lett. 102, 080403 (2009)
  * (37) J. Y. Vaishnav, J. Rusechas, C. W. Clark, and G. Juzeliūnas, Phys. Rev. Lett. 101, 265302 (2008).
  * (38) S. K. Dutta, B. K. Teo, and G. Raithel, Phys. Rev. Lett. 83. 1934 (1999).
  * (39) Y. J. Lin, R. L. Compton, A. R. Perry, W. D. Phillips, J. V. Porto, and I. B. Spielman, arXiv:0809.2976 (2008).
  * (40) J. Fuchs, G. J. Duffy, W. J. Rowlands, A. Lezama, P. Hannaford, and A. M. Akulshin, J. Phys. B: At. Mol. Opt. Phys. 40, 1117(2007).
  * (41) M. O. Scully and M. S. Zubairy, Quantum Optics, (Cambridge University Press, Cambridge, 1997).
  * (42) A. Messiah, Quantum Mechanics (North-Holland/ Elsevier Science, New York, 1962).
  * (43) J. M. Hou, L. J. Tian, and S. Jin, Phys. Rev. B 73, 134425 (2006).
  * (44) X. G. Wen, Quantum Field Theory of Many-Body Systems, (Oxford University Press, Oxford, 2004).
  * (45) C. H. Park et al., Nature Phys. 4, 213 (2008); C. H. Park et al., Phys. Rev. Lett. 101, 126804 (2008).
  * (46) D. M. Stamper-Kurn et al., Phys. Rev. Lett. 83, 2876 (1999).