# Effect of pseudospin polarization on wave packet dynamics in graphene
antidot lattices (GALs) in the presence of a normal magnetic field

R. A. W. Ayyubi  Department of Physics, Quaid-I-Azam University, Islamabad
45320, Pakistan  N. J. M. Horing  Department of Physics, Stevens Institute of
Technology, Hoboken, New Jersey 07030, USA  K. Sabeeh  Department of Physics,
Quaid-I-Azam University, Islamabad 45320, Pakistan

###### Abstract

We have investigated the role of pseudospin polarization in electron wave
packet dynamics in pristine graphene and in a graphene antidot lattice subject
to an external magnetic field. Employing a Green’s function formalism, we show
that the electron dynamics can be controlled by tuning pseudospin
polarization. We find that in Landau quantized pristine graphene both the
propagation of an electron wave packet and Zitterbewegung oscillations
strongly depend on pseudospin polarization. The electron wave packet is found
to propagate in the direction of initial pseudospin polarization. We also show
that, in this system, the propagation of an electron can be enhanced in any
desired direction by carving a one dimensional antidot lattice in that
direction. The study suggests that a graphene antidot lattice can serve as a
channel for electron transport with the possibility of tunability by means of
pseudospin polarization, antidot potential and applied normal magnetic field
strength.

##  I INTRODUCTION

It has been discovered that the charge carriers in two-dimensional materials
such as graphene have, in addition to spin, other binary degrees of freedom.
These are sublattice pseudospin and valley pseudospin.Xu  _et al._ (2014);
Mecklenburg and Regan (2011) The graphene honeycomb lattice is composed of two
triangular sublattices; this allows representing the wave function of
electrons as a two-component spinor, like a spin-$1/2$ particle.Mecklenburg
and Regan (2011) Similarly extrema in the band structure at $K$ and
$K^{\prime}$ points also define a binary degree of freedom which is known as
the valley pseudospin. The established field of spintronics exploits spin for
electronic applications.Pesin and Macdonald (2012); San-Jose  _et al._ (2009);
Song  _et al._ (2015); Kamalakar  _et al._ (2015); Loss and DiVincenzo (1998);
Shabbir  _et al._ (2018) Efforts to understand and control these additional
degrees of freedom with the aim of utilizing them in technological
applications has led to the emerging fields of pseudospintronics and
valleytronics.Wang  _et al._ (2012); Van Tuan and Roche (2016) Towards this
aim, in this work we focus on the sublattice pseudospin degree of freedom of
the prototypical 2D material, which is graphene. Specifically, we address the
role of pseudospin polarization in the propagation of electrons in graphene in
the presence of both an applied magnetic field and an antidot lattice.
We also investigate the role of pseudospin polarization on Zitterbewegung
(ZB), which is the oscillatory motion of the wave packet, another important
phenomenon observed in Dirac materials.Castro Neto  _et al._ (2009) These
oscillations have features which differ in the presence and in the absence of
a magnetic field. In monolayer graphene, in the absence of a magnetic field,
ZB oscillations are transient and disappear after some time, while for finite
B, ZB oscillations are permanent and do not die with time. Maksimova  _et al._
(2008); Rusin and Zawadzki (2008) A recent study suggests that ZB oscillations
in monolayer graphene in the absence of magnetic field strongly depends on
initial pseudospin polarization, which is named as pseudospinorial
Zitterbewegung (PZB). Serna  _et al._ (2019) These studies invite us to
explore the effect of initial pseudospin polarization on ZB in graphene in the
presence of a magnetic field, which we have investigated in this study.
Our main finding is that the interplay of pseudospin polarization, applied
magnetic field and antidot lattice leads to enhanced propagation in the
direction of the antidot lattice; this raises the possibility of creating one-
dimensional channels for electrons in Landau quantized graphene. The control
features for propagation in these channels are the applied magnetic field,
strength of antidot lattice and pseudospin polarization. Further, we also show
that there are persistent Zitterbewegung (ZB) Oscillations in the presence of
an applied magnetic field in pristine graphene, which are sensitive to
pseudospin polarization.

##  II Graphene in A Perpendicular Magnetic Field: GREEN’S FUNCTION

Graphene is a single layer of carbon atoms packed in a honeycomb lattice. The
effective Hamiltonian of graphene has the form Wallace (1947)

$H=\gamma\hskip 1.70717pt\bm{\sigma}_{\nu}\hskip 1.70717pt.\hskip
1.70717pt\bf{p}$ (1)

where $\sigma_{\nu}=[\sigma\textsubscript{x},1_{\nu}\hskip
2.84526pt\sigma\textsubscript{y}]$ and $\sigma\textsubscript{x}$,
$\sigma\textsubscript{y}$ are Pauli’s spin matrices, which act on the
sublattice/pseudospin space and represent the sublattice degree of freedom of
graphene’s honeycomb lattice structure, also $1_{\nu}=1$ or -1 for $K$ or
$K^{\prime}$ valleys, and $\gamma=\frac{3}{2}\alpha_{h}\textit{d}\approx
10^{6}ms^{-1}$ plays the role of a density-independent Fermi velocity
($\alpha_{h}$ is the hopping amplitude originating from the tight binding
approximation in which the lattice spacing is d).
We consider a graphene sheet placed on the $xy$-plane in a perpendicular and
uniform magnetic field B= $B\hat{z}$ with vector potential
$\textbf{A}=\frac{1}{2}(\textbf{B}\times\textbf{r})$. The magnetic field is
introduced by minimal substitution
$\textbf{p}\rightarrow\textbf{p}-e\textbf{A}$ in Eq. (1). The requirement of
gauge invariance leads to

$\mathcal{G}(\textbf{r},\textbf{r}^{\prime};t,t^{\prime})=C(\textbf{r},\textbf{r}^{\prime})\hskip
2.84526pt\mathcal{G}^{\prime}(\textbf{r}-\textbf{r}^{\prime};t-t^{\prime}),$
(2)

where the Green’s function
$\mathcal{G}^{\prime}(\textbf{r}-\textbf{r}^{\prime};t-t^{\prime})$ is
spatially translationally invariant and gauge invariant and it is related to
usual Green’s function
$\mathcal{G}(\textbf{r},\textbf{r}^{\prime};t,t^{\prime})$ with the help of
Peierls phase factor $C(\textbf{r},\textbf{r}^{\prime})$, which carries all
aspects of the lack of translational invariance in a magnetic field and all
gauge dependence as

$C(\textbf{r},\textbf{r}^{\prime})=exp\left(\frac{\textit{i}e}{2\hbar}\textbf{r}\hskip
1.42262pt\cdot\hskip
1.42262pt\textbf{B}\times\textbf{r}^{\prime}-\phi(\textbf{r})+\phi(\textbf{r}^{\prime})\right)$
(3)

($\phi(\textbf{r})$ is an arbitrary gauge function). The translationally
invariant Green’s function is given by the equation of motion

$\left[\textit{i}\frac{\partial}{\partial
T}-\gamma\sigma_{\nu}\cdot\left(\frac{1}{\textit{i}}\frac{\partial}{\partial\textbf{R}}-\frac{e}{2}\textbf{B}\times\textbf{R}\right)\right]\mathcal{G}^{\prime}(\textbf{R},T)=\\\
I\textsubscript{2}\delta\textsuperscript{2}(\textbf{R})\delta(T)$ (4)

where $\textbf{R}=\textbf{r}-\textbf{r}^{\prime}$ such that
$X=x-x^{\prime},Y=y-y^{\prime}$ and $T=t-t^{\prime}$.
The diagonal elements of the Green’s function matrix ($\omega$ representation)
for this system of monolayer graphene in a uniform, constant, perpendicular
magnetic field are ($\mathcal{G}^{\prime}_{11\atop 22}$ collectively
represents $\mathcal{G}^{\prime}_{11}$, $\mathcal{G}^{\prime}_{22}$)

$\mathcal{G}^{\prime}_{11\atop
22}(\textbf{R};\omega)=\frac{1}{4\pi\hbar\gamma^{2}}\hskip
2.27621pt\omega\hskip
2.27621pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}\left(X^{2}+Y^{2}\right)}\\\
\times\sum_{n=0}^{\infty}\frac{L_{n}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}\left(X^{2}+Y^{2}\right)\right]}{\frac{\omega^{2}}{\omega^{2}_{g}}-\left(n+\frac{1-1_{\nu}}{2}\right)}.$
(5)

We introduced the notation

$\omega_{g}=\gamma\sqrt{\frac{2eB}{\hbar}},$ (6)

which is the cyclotron frequency for Dirac fermions. The poles of
$\mathcal{G}^{\prime}_{11\atop 22}(\textbf{R};\omega)$ show that the energy
spectrum is Landau quantized, with Landau level index given by $n$.
The off-diagonal matrix elements can be calculated using the relation
($\gamma_{\nu}\equiv 1_{\nu}\gamma$) Horing and Liu (2009)

$\omega\mathcal{G}^{\prime}_{21\atop
12}=[\gamma\Pi\textsubscript{X}\textsubscript{Y}\pm\textit{i}\gamma_{\nu}\Pi\textsubscript{Y}\textsubscript{X}]\mathcal{G}^{\prime}_{11\atop
22}.$ (7)

In the above equation,
$\Pi\textsubscript{X}\textsubscript{Y}\equiv\frac{1}{\textit{i}}\frac{\partial}{\partial
X}+\frac{eB}{2}Y$ and
$\Pi\textsubscript{Y}\textsubscript{X}\equiv\frac{1}{\textit{i}}\frac{\partial}{\partial
Y}-\frac{eB}{2}X$ are momentum operators and $\mathcal{G}^{\prime}_{21\atop
12}$= $\mathcal{G}^{\prime}_{21}$ or $\mathcal{G}^{\prime}_{12}$ corresponds
to the upper or lower $\pm$, $\mp$ signs elsewhere in the equations. Making
use of Eq. (7), off-diagonal elements of Green’s function matrix for $K$ and
$K^{\prime}$ can be separately expressed as

$K(1_{\nu}=+1):\hskip 3.41432pt\mathcal{G}^{\prime}_{21\atop
12}(\textbf{R};\omega)=\frac{\omega_{g}^{2}}{8\pi\hbar\gamma^{3}}\hskip
2.27621pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}\left(X^{2}+Y^{2}\right)}\\\
\times\left(\textit{i}X\mp
Y\right)\sum_{n=1}^{\infty}\frac{L_{n-1}^{1}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}\left(X^{2}+Y^{2}\right)\right]}{\frac{\omega^{2}}{\omega^{2}_{g}}-n}$
(8) $K^{\prime}(1_{\nu}=-1):\hskip 3.41432pt\mathcal{G}^{\prime}_{21\atop
12}(\textbf{R};\omega)=\frac{\omega_{g}^{2}}{8\pi\hbar\gamma^{3}}\hskip
2.27621pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}\left(X^{2}+Y^{2}\right)}\\\
\times\left(\textit{i}X\pm
Y\right)\sum_{n=1}^{\infty}\frac{L_{n-1}^{1}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}\left(X^{2}+Y^{2}\right)\right]}{\frac{\omega^{2}}{\omega^{2}_{g}}-\left(n+1\right)}.$
(9)

Clearly, the off-diagonal elements interchange with the interchange of Dirac
points ($K$ and $K^{\prime}$). Horing (2017, 2018)
The Green’s functions matrix in time representation can be obtained by Fourier
transform of Eqs. (5), (8), (9) as

$\mathcal{G}^{\prime}_{\mu\nu}(\textbf{R};t)=\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}\hskip
2.27621pte^{-\textit{i}\omega t}\hskip
2.27621pt\mathcal{G}^{\prime}_{\mu\nu}(\textbf{R};\omega)$ (10)

where $\mu,\nu=1,2$ denote matrix indices.
Noting that $\mathcal{G}^{\prime}_{\mu\nu}(\textbf{R};\omega)$ has real energy
poles at $\epsilon_{K}=\pm\omega_{g}\hskip 2.84526pt\sqrt{n}$ ;
$\epsilon_{K^{\prime}}=\pm\omega_{g}\hskip 2.84526pt\sqrt{n+1}$ we employ
contour integration with $\omega\rightarrow\omega+\textit{i}\hskip
1.70717pt0^{+}$ for the retarded Green’s function with the contour closed in
the lower half plane running clockwise from $-\infty$ to $+\infty$. For the
$K$-point, we have (for the $K^{\prime}$-point
$\sqrt{n}\rightarrow\sqrt{n+1}$)

$\mathcal{G}^{\prime}_{11\atop
22}\left(x,x^{\prime};y,y^{\prime};t\right)=-\textit{i}\hskip
1.99168pt\eta_{+}(t)\hskip 1.99168pt\frac{\omega_{g}^{2}}{4\pi\hskip
1.42262pt\gamma^{2}}\hskip 3.41432pte^{-\frac{1}{2}\hskip
1.99168pt\zeta}\hskip
3.41432pt\sum_{n=0}^{\infty}L_{n}\left[\zeta\right]\hskip 2.84526pt\\\
\times\cos\left(\frac{\omega_{g}\hskip 2.27621ptl}{\gamma}\hskip
1.42262pt\frac{t}{\tau_{o}}\hskip 1.13809pt\sqrt{n}\hskip 2.27621pt\right)$
(11)

$\mathcal{G}^{\prime}_{21\atop
12}\left(x,x^{\prime};y,y^{\prime};t\right)=-\eta_{+}(t)\hskip
1.99168pt\frac{\omega_{g}^{3}\hskip 2.27621ptl}{8\pi\hskip
1.42262pt\gamma^{3}}\hskip 3.41432pte^{-\frac{1}{2}\hskip
1.99168pt\zeta}\hskip 3.41432pt\\\
\times\left[\textit{i}\left(\frac{x-x^{\prime}}{\textit{l}}\right)\mp\left(\frac{y-y^{\prime}}{\textit{l}}\right)\right]\sum_{n=1}^{\infty}\frac{L_{n-1}^{1}\left[\zeta\right]}{\sqrt{n}}\
\hskip 2.84526pt\\\ \times\sin\left(\frac{\omega_{g}\hskip
2.27621ptl}{\gamma}\hskip 1.42262pt\frac{t}{\tau_{o}}\hskip
1.13809pt\sqrt{n}\hskip 2.27621pt\right):$ (12)

Here, we have introduced an arbitrary constant length, $l$, chosen for
convenience to be the width of an impressed wave packet, and $\eta_{+}(t)=0,1$
for $t\leq 0,t\textgreater 1$, respectively, is the Heaviside unit step
function; also, we have defined

$\zeta=\frac{1}{4}\left(\frac{\omega_{g}\hskip
2.27621pt\textit{l}}{\gamma}\right)^{2}\left[\left(\frac{x-x^{\prime}}{\textit{l}}\right)^{2}+\left(\frac{y-y^{\prime}}{\textit{l}}\right)^{2}\right]$

and $\tau_{o}=\textit{l}/\gamma$.
To study the wave packet dynamics, we take the initial wave function to be a
Gaussian wave packet having nonvanishing average momentum $p_{0x}=\hbar
k_{0x}$ and width $l$,

$\psi(\textbf{r},0)=\frac{f(\textbf{r})}{\sqrt{\mid c_{1}\mid^{2}+\mid
c_{2}\mid^{2}}}\begin{pmatrix}c_{1}\\\ c_{2}\end{pmatrix}$ (13)
$f(\textbf{r})=\frac{1}{l\hskip 1.13809pt\sqrt{\pi}}\hskip
1.42262ptexp\left(-\frac{x^{2}+y^{2}}{2\hskip
1.13809ptl^{2}}+\textit{i}k_{0x}x\right),$

where $c_{1}$ and $c_{2}$ are the coefficients which set the initial
pseudospin polarization. Also $\psi(\textbf{r},0)$ can be taken as a smooth
enveloping function considering that lattice period is much smaller than the
width $l$ of the initial wave packet. The Green’s function matrix elements
$\mathcal{G}_{\mu\nu}(\textbf{r},\textbf{r}^{\prime},t)$ determine the time
evolution of an arbitrary initial state $\psi(\textbf{r},0)$; in
Schr$\ddot{o}$dinger representation it is given by

$\psi_{\mu}(\textbf{r},t)=\int d\textbf{r}^{\prime}\hskip
2.27621pt\mathcal{G}_{\mu\nu}(\textbf{r},\textbf{r}^{\prime},t)\hskip
2.27621pt\psi_{\nu}(\textbf{r}^{\prime},0)\hskip 3.41432pt,$ (14)

where $\mu,\nu=1,2$ denote the matrix indices, which corresponds to upper
component $\psi_{1}(\textbf{r},t)$ and lower component
$\psi_{2}(\textbf{r},t)$ of state $\psi_{\mu}(\textbf{r},t)$. The probability
density will be

$\rho(\textbf{r},t)=\mid\psi(\textbf{r},t)\mid^{2}=\mid\psi_{1}(\textbf{r},t)\mid^{2}+\mid\psi_{2}(\textbf{r},t)\mid^{2}.$
(15)

And to study ZB, the average value of coordinates can be represented as

$\bar{x}_{j}=\int\psi_{1}^{\ast}(\textbf{r},t)\hskip 1.13809ptx_{j}\hskip
1.13809pt\psi_{1}(\textbf{r},t)\hskip
1.13809ptd\textbf{r}+\int\psi_{2}^{\ast}(\textbf{r},t)\hskip
1.13809ptx_{j}\hskip 1.13809pt\psi_{2}(\textbf{r},t)\hskip
1.13809ptd\textbf{r}$ (16)

where $j=1,2$ with $x_{1}$=$x$ and $x_{2}$=$y$. The Peierls phase factor
defined in Eq. (3) for the choice
$\phi(\textbf{r})=\phi(\textbf{r}^{\prime})\equiv 0$ has the form

$C(\textbf{r},\textbf{r}^{\prime})=exp\left[\frac{\textit{i}}{4}\left(\frac{\omega_{g}l}{\gamma}\right)^{2}\left(\frac{yx^{\prime}}{l^{2}}-\frac{xy^{\prime}}{l^{2}}\right)\right].$

##  III Dynamics of A Gaussian wave packet with different pseudospin
polarizations in Landau quantized graphene: Zitterbewegung

To obtain results for the temporal evolution of the initial Gaussian wave
packet, ZB oscillations and the effect of initial-pseudospin polarization, we
have performed numerical calculations. To facilitate it, we have introduced
following dimensionless variables: Maksimova  _et al._ (2008); Demikhovskii
_et al._ (2008)

  * •

A dimensionless parameter which is suitable to replace wave vector $k_{0x}$ is
$a_{0}=k_{0x}\hskip 0.85358ptl$.

  * •

Distance for propagation of the wave packet can be measured in units of
initial width $l$ of wave packet.

  * •

Time can be measured in $\tau_{o}$=$l/\gamma$ ($\gamma$ is Fermi velocity,
$10^{6}$ m/s) units.

  * •

Some other variables can be combined to produce dimensionless variables, e.g
$b=\frac{\omega_{g}\hskip 1.70717ptl}{\gamma}$.

The Landau level summation is performed up to the $10th$ Landau level in all
calculations since the results become convergent in this limit.
We consider four cases with different initial pseudospin polarizations
$\\{c_{1},c_{2}\\}$ for the Gaussian wave packet given in Eq. (13).
Case-1: $\\{c_{1},c_{2}\\}$=$\\{1,0\\}$ which corresponds to initial electron
probability of one at the sites of sublattice A.
Case-2:$\\{c_{1},c_{2}\\}$=$\\{1,1\\}$ corresponds to the situation where the
pseudospin is directed along the $x$-axis.
Case-3: $\\{c_{1},c_{2}\\}$=$\\{1,i\\}$ this corresponds to the pseudospin
directed along $y$-axis.
Case-4: $\\{c_{1},c_{2}\\}$=$\\{1,e^{\textit{i}\pi/4}\\}$ this implies that at
time $t$=$0$, the pseudospin lies in $x$-$y$ plane making an angle of $45^{o}$
with $x$-axis.
The numerical results obtained from Eqs. (14), (15) are plotted in Fig.1.
In Fig.1, electronic probability density is plotted for parameters
$t$=$1,5,10\hskip 2.27621pt\tau_{o}$, with momentum $k_{0x}=0.6\hskip
2.27621ptnm^{-1}$, width $l$=$2\hskip 2.27621ptnm$ so that
$a_{0}$=$k_{0x}\hskip 2.27621ptl$=$1.2$ and $\tau_{o}$=$2\hskip
2.27621ptfemtosecond$ and $B$=$3.3\hskip 2.27621pt$ in units of $Tesla$ i.e.
$\omega_{g}\simeq 1\times 10^{14}\hskip 2.27621ptHz$. $B$=$3.3\hskip
2.27621ptT$ is chosen to facilitate comparison with work in the literature and
also because it facilitate numerical computation as $\omega_{g}$ is a round
figure at this value. Maksimova  _et al._ (2008) Left panel shows the initial
wave packet at a very small time $t$=$1\hskip 2.27621pt\tau_{o}$, and as we
move from left to right the time evolution of initial Gaussian packet can be
seen for different initial pseudospin polarizations. The strength of the
electron probability density $\rho(x,y,t)$ is given by the color bar on the
right side.
As one can see in Fig.1 $(a),(b),(c)$, the wave packet spreads and propagates
in the plane of graphene sheet in the form of rings for the pseudospin
$c_{1}$=$1$ and $c_{2}$=$0$. Initially in Fig.1 $(a)$, maximum probability of
the electron is located at the origin, but as time increases, the wave packet
propagates and the electron density can be found at a radius of $r\simeq
24\hskip 2.27621ptnm$ in $20\hskip 2.27621ptfs$. Similarly for other three
cases, the wave packet propagates with its maximum probability density in the
direction of initial pseudospin polarizations, but in the shape of incomplete
rings, this is because the probability density gradually decreases in the
directions away from the direction of pseudospin polarization. It propagates
in the $x$ direction when initial pseudospin polarization is along $x$ axis
(see Figs.1 $(d),(e),(f)$), in the $y$ direction when initial pseudospin
polarization is along $y$ axis (see Figs.1 $(g),(h),(i)$) and it propagates in
$(x,y)$=$(1,1)$ direction when initial pseudospin is polarized along
$(x,y)$=$(1,1)$ direction (see Figs.1 $(j),(k),(l)$). This is because of the
conservation of chirality, in which momentum gets aligned with pseudospin and
$\vec{\sigma}_{\nu}\hskip 1.70717pt.\hskip 1.70717pt\vec{p}$ remains
conserved. Also the distance covered by the wave packet in these three cases
is $r\simeq 24\hskip 2.27621ptnm$ in $20\hskip 2.27621ptfs$. Hence the
direction of propagation of a wave packet in Landau quantized graphene can be
controlled using pseudospin polarization. Also, the wave packet propagates
without any splitting; splitting was observed in the case of monolayer
graphene in the absence of magnetic field.Maksimova  _et al._ (2008) We
propose that, experimentally, this type of controlled propagation of a wave
packet in any direction can easily be obtained using photonic graphene as test
beds.Song  _et al._ (2015)

Figure 1:  (Color online) The electron probability density $\rho(x,y,t)$ for
initial Gaussian wave packet with $B=3.3\hskip 2.84526ptTesla$,
$a_{0}$=$k_{0x}\hskip 2.27621ptl$=$1.2$ at times $1\tau_{o}$, $5\hskip
1.42262pt\tau_{o}$ and $10\hskip 1.42262pt\tau_{o}$. Time increases as we move
from left to right. Color bar given at extreme right side shows the strength
of probability density from its minimum to its maximum. Up to down: Four rows
with initial pseudospin polarizations $\\{1,0\\}$,$\\{1,1\\}$,$\\{1,i\\}$ and
$\\{1,e^{i\frac{\pi}{4}}\\}$ with pseudospin and propagation directed in the
radial, $x$, $y$ and $(x=1,y=1)$ directions respectively.

Also note that the electronic probability densities plotted in Fig.1 are not
symmetric with respect to both $x$ and $y$ axes:
$\rho(x,y,t)_{\\{c_{1},c_{2}\\}}\neq\rho(-x,y,t)_{\\{c_{1},c_{2}\\}}$ and
$\rho(x,y,t)_{\\{c_{1},c_{2}\\}}\neq\rho(x,-y,t)_{\\{c_{1},c_{2}\\}}$ for any
pseudospin polarization (subscript $\\{c_{1},c_{2}\\}$ defines the
corresponding pseudospin polarization): This means that center of the wave
packet is oscillating along both the $x$ as well as the $y$ directions; these
oscillations can be readily recognized as Zitterbewegung oscillations.

Figure 2: Average coordinates $\bar{x}(t)$ (left column) and $\bar{y}(t)$
(right column) versus time $(\tau_{o}=l/\gamma)$ corresponding to four
different values of pseudospin polarization (changing from up to down), at two
values of momentum $a_{0}$ with $B=3.3\hskip 2.84526ptTesla$. (a),(c),(e),(g):
$\bar{x}(t)$ versus time corresponding to pseudospin polarizations
$\\{1,0\\}$,$\\{1,1\\}$, $\\{1,i\\}$ and $\\{1,e^{i\frac{\pi}{4}}\\}$
respectively. (b),(d),(f),(h): $\bar{y}(t)$ versus time corresponding to
pseudospin polarizations $\\{1,0\\}$,$\\{1,1\\}$, $\\{1,i\\}$ and
$\\{1,e^{i\frac{\pi}{4}}\\}$ respectively. Here $l^{\prime}$=$l\times
10^{-5}$.

To examine this trembling motion we use Eq. (16) and solve it numerically for
expectation values of both $x$ and $y$ coordinates. Figs.2 $(a),(c),(e),(g)$
and Figs.2 $(b),(d),(f),(h)$ show the oscillations in the wave packet’s center
along $x$ and $y$ directions respectively. Results are plotted for two values
of initial momentum $k_{0x}$ = 0.6 and 0.8 i.e $a_{0}=k_{0x}\hskip
1.42262pt\textit{l}$=$1.2$ and $1.6$, given by solid and dotted lines for four
different values of pseudospin polarization. It is clear from Fig.2 that an
increase in momentum results in a decrease in amplitude of Zitterbewegung
oscillations without any other change in the behaviour of ZB. In further
discussion, ZB oscillations corresponding to different pseudospin
polarizations will be referred to as $\bar{x}_{\\{c_{1},c_{2}\\}}$ and
$\bar{y}_{\\{c_{1},c_{2}\\}}$.
The ZB oscillations shown in Fig. 2 are of the order of $nanometer$ (easily
detectable), except $x(t)_{\\{1,0\\}}$ and $x(t)_{\\{1,i\\}}$ which are of the
order of $0.1\hskip 2.27621ptpm$. On comparing the ZB oscillations
$x(t)_{\\{1,0\\}}$, $x(t)_{\\{1,1\\}}$, $x(t)_{\\{1,i\\}}$, we have concluded
that when the direction of initial momentum ($x$-axis) is perpendicular to
initial pseudospin polarizations ($z,y$-axis) then the ZB oscillations carry
very small amplitude in the direction of initial momentum. Note that, ZB
Oscillations in Fig. 2 $(g)$ and Fig. 2 $(h)$ are very similar to those of
Fig. 2 $(c)$ and Fig. 2 $(f)$ respectively. This is because pseudospin
polarization $\\{1,e^{i\frac{\pi}{4}}\\}$ has projection on both $x$ and $y$
axes so ZB oscillations in this case have detectable amplitude in both
directions.
In Figs.2 $(a)$-$(h)$ we have seen that initially the amplitude of the
Zitterbewegung oscillations increases, then these oscillations seems to die
out but reappear for all pseudospin polarizations. For example in Fig.2 $(c)$
these oscillations reappear at $t\simeq 18,55,84\hskip 1.42262pt\tau_{o}...$ .
Hence when there is a magnetic field applied to the system, Zitterbewegung
oscillations are not transient; rather they are recurrent. Also, in the
presence of the magnetic field, several ZB frequencies appear (see Figs.2
$(a)$-$(h)$). This is different from the ZB phenomenon observed in monolayer
graphene without magnetic field, in which ZB oscillations are transient having
a single frequency.Maksimova  _et al._ (2008)
Hence, the results in hand are in agreement with the previous studies in both
respects i.e., (1) in the presence of a magnetic field ZB oscillations are
permanent and (2) ZB oscillations strongly depend on initial pseudospin
polarization. Rusin and Zawadzki (2008) Serna  _et al._ (2019)

Figure 3: Average coordinates $\bar{x}(t)$ versus $\bar{y}(t)$ corresponding
to pseudospin polarization $\\{c_{1},c_{2}\\}$=$\\{1,1\\}$ for momentum
$a_{0}=1.2$ with $B=3.3\hskip 2.84526ptTesla$ (a) ZB Trajectory for $t=0$ to
$t=100\tau_{o}(0.2\hskip 2.27621ptps)$. (b) Zoomed view of the same trajectory
for $t=7\tau_{o}$ to $t=100\tau_{o}$ .

Finally, in Fig. 3 we have plotted the average coordinates
$\bar{x}(t)_{\\{c_{1},c_{2}\\}}$ and $\bar{y}(t)_{\\{}c_{1},c_{2}\\}$ against
each other to study the ZB trajectory of the center of the wave packet
corresponding to the Figs.2 $(c),(d)$. Initially due to large ZB, the center
of the packet sweeps a large area. For a better understanding, we have shown a
zoomed view of the ZB trajectory in Fig. 3 $(b)$. Due to the presence of
recurrent ZB, these trajectories do not disappear with time (infinite
trajectories).
It is well known that the electronic band structure of graphene can be
modified by introducing external ID potentials using nanopatterning. Esaki and
Tsu (1970) These 1D potentials can not only modify the energy spectrum of the
system but they also help us to control the charge transport properties of the
system. Heide  _et al._ (2018); Castro Neto  _et al._ (2009); Park  _et al._
(2008a); Choi  _et al._ (2014); Park  _et al._ (2008b). Potentials such as
antidot lattices can be carved on graphene by various techniques and lattice
parameters can be tuned.Pedersen  _et al._ (2008); Giesbers  _et al._ (2012)
With the aim of controlling electron propagation in graphene along a 1D
channel we introduce a 1D antidot lattice in the following section.

##  IV Graphene Antidot lattice in the presence of a magnetic field: Landau
Minibands

We consider a two-dimensional graphene sheet having a one-dimensional lattice
of quantum antidot potential barriers, with a quantizing magnetic field B,
which is perpendicular to the plane of the graphene sheet.
We model the antidot array as a row of Dirac delta functions. Following the
Kronig-Penny model for a quantum antidot lattice, we introduce an infinite
array of identical quantum antidot potential barriers periodically spaced
along the $x$-axis at $x_{n^{\prime}}=n^{\prime}d,y\equiv 0$ as

$U(\textbf{r})=U(x,y)=\alpha\sum_{n^{\prime}=-\infty}^{\infty}\delta(x-n^{\prime}d)\delta(y)\hskip
2.84526pt,$ (17)

where $\alpha>0$ , $n^{\prime}=-\infty$ to $+\infty$ and

$\alpha=U_{0}\hskip 2.84526pta^{2}$ (18)

where $U_{0}$ is the barrier height, $a^{2}$ is area and $d$ is the uniform
separation of the periodic quantum antidot potential barriers as shown in
Fig.4

Figure 4: Two dimensional graphene sheet having a one dimensional lattice of
quantum antidot barriers at $y=0$. Solid black circles represent periodically
placed quantum antidots along the $x$-axis. $a^{2}$ is the dot area and $d$ is
the uniform separation of the periodically placed quantum antidot potential
barriers.

The Green’s function $\mathcal{G}(x_{1};x_{2};0,0;\omega)$ for the two
dimensional Kronig-Penney-like model having a one dimensional antidot lattice,
for an electron propagating directly along the axis of the antidot lattice
($y\equiv 0\equiv y_{1}\equiv y_{2}$ and suppressing further reference to $y$)
is given by Morgenstern Horing (2020)

$\mathcal{G}(x_{1};x_{2};\omega)=\mathcal{G}^{0}(x_{1};x_{2};\omega)+\frac{\alpha
d}{2\pi}\sum_{n^{\prime}=-\infty}^{\infty}\mathcal{G}^{0}(x_{1};n^{\prime}d;\omega)\\\
\times\int_{-\frac{\pi}{d}}^{\frac{\pi}{d}}dp\hskip
5.69054pte^{-\textit{i}pn^{\prime}d}\left[1-\alpha\hskip
2.84526pt\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\omega)\right]^{-1}\tilde{\mathcal{G}}^{0}(p;x_{2};\omega).$
(19)

Here, $\mathcal{G}^{0}(\textbf{r}_{1},\textbf{r}_{1}^{\prime})$ is the Green’s
function for graphene in a perpendicular magnetic field in the complete
absence of quantum antidot potential barriers, which is given by Eqs. (5),(8)
and (9), while the overhead dot represents the spatially translationally
invariant Green’s function i.e.
$\mathcal{G}^{0}(m^{\prime}d,n^{\prime}d;\omega)=\dot{\mathcal{G}}^{0}([m^{\prime}-n^{\prime}]d;\omega)$
where $m^{\prime}$ and $n^{\prime}$ are integers. Also the Green’s function
$\tilde{\mathcal{G}}(p)$ can be expanded in a Fourier series due to the
periodicity of the lattice, given by Morgenstern Horing (2020)

$\tilde{\mathcal{G}}(p)=\sum_{r=-\infty}^{\infty}e^{\textit{i}prd}\mathcal{G}(rd)$
(20)

where $r=0,1,2,3$… are integers, with

$\mathcal{G}(m^{\prime}d)=\frac{d}{2\pi}\int_{-\frac{\pi}{d}}^{\frac{\pi}{d}}dp\hskip
2.84526pte^{-\textit{i}pdm^{\prime}}\tilde{\mathcal{G}}(p),$ (21)

and

$\dot{\mathcal{G}}^{0}([m^{\prime}-n^{\prime}]d)=\frac{d}{2\pi}\int_{-\frac{\pi}{d}}^{\frac{\pi}{d}}dp\hskip
2.84526pte^{-\textit{i}p\left(m^{\prime}-n^{\prime}\right)d}\hskip
5.69054pt\dot{\tilde{\mathcal{G}}}^{0}(p).$ (22)

It is important to note that the Peierls phase factor is
$C(\textbf{r},\textbf{r}^{\prime})$=$1$ for the case involving propagation
directly along the axis of the antidot lattice (our choice $y\equiv 0\equiv
y_{1}\equiv y_{2}$ results in $\textbf{r}\parallel\textbf{r}^{\prime}$).
Therefore $C(\textbf{r},\textbf{r}^{\prime})$ does not appear in Eq. (19) and
the eigen-energy spectrum given by poles of Eq. (19) is unaffected by
$C(\textbf{r},\textbf{r}^{\prime})$.
The energy spectrum of this system can be obtained from vanishing of the
frequency poles of the Green’s function of Eq. (19):

$Det(I_{2}-\alpha\hskip
2.84526pt\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\omega))=0.$ (23)

Eqs. (5) and (8) taken jointly with Eq. (20) yields ($Y=0$,
$X^{2}+Y^{2}=X^{2}=(rd)^{2}$)

$\dot{\tilde{\mathcal{G}}}^{0}_{11\atop
22}(p;0,0;\omega)_{K}=\frac{1}{4\pi\hbar\gamma^{2}}\hskip
2.27621pt\omega\hskip 2.27621pt\sum_{r=-\infty}^{\infty}e^{iprd}\hskip
2.27621pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}\left(rd\right)^{2}}\\\
\times\sum_{n=0}^{\infty}\frac{L_{n}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}\left(rd\right)^{2}\right]}{\frac{\omega^{2}}{\omega^{2}_{g}}-n}$
(24)

and

$\dot{\tilde{\mathcal{G}}}^{0}_{21\atop
12}(p;0,0;\omega)_{K}=\frac{\textit{i}\hskip
1.70717pt\omega_{g}^{2}}{8\pi\hbar\gamma^{3}}\hskip
2.27621pt\sum_{r=-\infty}^{\infty}\left(rd\right)\hskip
2.27621pte^{iprd}\hskip
2.84526pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}\left(rd\right)^{2}}\\\
\times\sum_{n=1}^{\infty}\frac{L_{n-1}^{1}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}\left(rd\right)^{2}\right]}{\frac{\omega^{2}}{\omega^{2}_{g}}-n}.$
(25)

for the $K$ point in graphene. It is clear from Eqs. (24) and (25) that
$\dot{\tilde{\mathcal{G}}}^{0}_{11}(p;0,0;\omega)_{K}=\dot{\tilde{\mathcal{G}}}^{0}_{22}(p;0,0;\omega)_{K}$
and
$\dot{\tilde{\mathcal{G}}}^{0}_{12}(p;0,0;\omega)_{K}=\dot{\tilde{\mathcal{G}}}^{0}_{21}(p;0,0;\omega)_{K}$.
Hence Eq. (23) can be written as

$1-2\alpha\dot{\tilde{\mathcal{G}}}^{0}_{11}+\alpha^{2}\left({\dot{\tilde{\mathcal{G}}}^{0}_{11}}^{2}-{\dot{\tilde{\mathcal{G}}}^{0}_{12}}^{2}\right)=0.$
(26)

Note that since the antidot radius is very small (i.e
$\frac{\alpha\omega_{g}}{4\pi\hbar\gamma^{2}}\ll 1$), a root of Eq. (26)
$\omega$ approaches the pole position i.e $\omega\rightarrow\omega_{n}$, so
that the $n$-$th$ pole has the primary influence in determining the eigen-
energy root $\omega_{n}$. Therefore, we can make a reasonable first
approximation by dropping all other terms of the $n$-sum. Also, for
$\frac{\omega_{g}d}{8\gamma}>1$ , it suffices to keep only $r=-1,0,1$ terms of
the $r$-sum in Eqs. (24) and (25). This imposes following condition on antidot
spacing

$d>\frac{145}{\sqrt{B\hskip 2.84526pt(Tesla)}}\hskip 2.84526pt(nm).$ (27)

Finally, Eq. (26) can be written as

$1-2\frac{\omega_{g}\hskip 1.70717pt\Omega\hskip
1.70717pt\omega_{n}}{\omega_{n}^{2}-n\omega_{g}^{2}}+\frac{\omega_{g}^{2}\hskip
1.70717pt\Omega^{2}\hskip
1.70717pt\omega_{n}^{2}}{\left(\omega_{n}^{2}-n\omega_{g}^{2}\right)^{2}}-\frac{\omega_{g}^{4}\hskip
1.70717pt\kappa^{2}}{\left(\omega_{n}^{2}-n\omega_{g}^{2}\right)^{2}}$ (28)

where we have defined

$\Omega=\frac{\alpha\hskip 2.27621pt\omega_{g}}{4\pi\hbar\gamma^{2}}\hskip
2.27621pt\left[1+2\cos{pd}\hskip
3.41432pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}d^{2}}L_{n}\left(\frac{\omega^{2}_{g}}{4\gamma^{2}}d^{2}\right)\right]$
(29)

and

$\kappa=\frac{\alpha\hskip 2.27621pt\omega_{g}}{4\pi\hbar\gamma^{2}}\hskip
2.27621pt\left(\frac{\omega_{g}\hskip 2.27621ptd}{\gamma}\right)\sin{pd}\hskip
3.41432pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}d^{2}}L_{n}\left(\frac{\omega^{2}_{g}}{4\gamma^{2}}d^{2}\right).$
(30)

The four roots of Eq. (28) describe the energy spectrum at the $K$ point of
graphene having a 1-D antidot lattice placed in a uniform normal magnetic
field:

$\omega_{n,K}=\frac{\Omega\pm\sqrt{\Omega^{2}+4(n\pm^{\prime}\kappa)}}{2}\hskip
2.27621pt\omega_{g}.$ (31)

Similarly, for $K^{\prime}$ point, $n$ will be replaced by $n+1$ on the right
hand side of Eq. (31).
In this, we have the energy spectrum composed of broadened Landau levels
(Landau minibands) for a graphene antidot lattice in a quantizing magnetic
field. Each Landau level has split into two branches and each branch has
broadened into a small continuous band (subband) of energy instead of a single
energy. Fig. 5 (a) and Fig. 5 (b) show the Landau minibands at the location of
$K$ and $K^{\prime}$ points. The broadening is so small that it can not be
observed with the naked eye, so we multiplied a broadening factor $\beta=200$
with the oscillatory terms $\cos(pd)$ and $\sin(pd)$. The parameter $\beta$
introduces an increase in the amplitude of the minibands to facilitate
observation of the broadening of the Landau minibands; i.e we used

$\Omega=\frac{\alpha\hskip 2.27621pt\omega^{2}_{g}}{4\pi\hbar\gamma^{2}}\hskip
2.27621pt\left(1+2\hskip 1.42262pt\beta\hskip 2.27621pt\cos{pd}\hskip
3.41432pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}d^{2}}L_{n}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}d^{2}\right]\right)$

and

$\kappa=\frac{\alpha\hskip 2.27621pt\omega_{g}}{4\pi\hbar\gamma^{2}}\hskip
2.27621pt\left(\frac{\omega_{g}\hskip 2.27621ptd}{\gamma}\right)\beta\hskip
1.70717pt\sin{pd}\hskip
3.41432pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}d^{2}}L_{n}\left(\frac{\omega^{2}_{g}}{4\gamma^{2}}d^{2}\right)$

to exhibit the Landau minibands plotted in Fig. 5.

Figure 5: Landau minibands ($n=0-3$) in Landau quantized graphene having a one
dimensional antidot lattice with antidot potential $U_{0}$=$100\hskip
2.27621ptmeV$, antidot radius $a$=$10\hskip 2.27621ptnm$ and spacing
$d$=$100\hskip 2.27621ptnm$ in the presence of magnetic field $B$=$3.3\hskip
2.84526ptTesla$. (a) $K$ point. (b) $K^{\prime}$ point.

The energy spectra of the $K$ and $K^{\prime}$ points differs by a unit shift
even in the presence of an antidot lattice. The presence of antidots induced a
gap between conduction and valence bands at both $K$ and $K^{\prime}$ points.
In the case of Fig. 5, if we consider $\beta=1$ with all other parameters
having same values, then the gap opened at the location of $K$ and
$K^{\prime}$ points is $7.96\hskip 2.27621ptmeV$ and $132\hskip 2.27621ptmeV$
respectively. This gap between conduction and valence bands increases with the
increase in strength of magnetic field and reaches $36.25\hskip 2.27621ptmeV$
and $283.4\hskip 2.27621ptmeV$ for $K$ and $K^{\prime}$ points respectively
when $B$=$15\hskip 2.84526ptTesla$. Similarly this gap shows an increase with
increase in strength of the antidot lattice.
Note that when the antidot strength $\alpha$ approaches zero (which means no
antidot lattice), $\Omega$ and $\kappa$ also approaches zero, i.e

$\lim_{\alpha\rightarrow 0}\Omega=0=\lim_{\alpha\rightarrow 0}\kappa.$

In this case Eq. (31) reduces to the case of a discrete eigen-energy spectrum
of graphene in a normal, uniform magnetic field. The width of Landau minibands
for the $K$ point of graphene is

$\Delta\omega_{n,K}=[\omega_{n,K}]_{p=\frac{2\pi}{d}}-[\omega_{n,K}]_{p=\frac{\pi}{d}}.$

Evaluation of this expression for the width of Landau minibands yields

$\frac{\Delta\omega_{n,K}}{\omega_{g}}=\alpha^{\prime}\hskip
2.27621pt\varrho\pm\frac{1}{2}\sqrt{\alpha^{\prime
2}(1+\varrho)^{2}+(2\sqrt{n})^{2}}\\\ \mp\frac{1}{2}\sqrt{\alpha^{\prime
2}(1-\varrho)^{2}+(2\sqrt{n})^{2}}\hskip 2.84526pt,$ (32)

where we have defined

$\alpha^{\prime}=\frac{\alpha\hskip
2.27621pt\omega_{g}}{4\pi\hbar\gamma^{2}}\qquad\varrho=2\hskip
1.70717pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}d^{2}}\hskip
3.69885ptL_{n}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}d^{2}\right].$

Upper and lower signs are used for conduction and valence bands respectively.
For the width of the Landau minibands for the $K^{\prime}$ point, we have
$(2\sqrt{n+1})^{2}$ in place of $(2\sqrt{n})^{2}$ in Eq. (32).
To investigate electron propagation in the system under study, first we have
to find the Green’s function matrix of the system in time representation. This
is evaluated in section V.

##  V  Green’s Functions: Frequency and Time Representation

To evaluate the full Green’s function matrix
$\mathcal{G}(x_{1};x_{2};\omega)_{K}$ given in Eq. (19), we have to find the
$\mathcal{G}^{0}(x_{1};x_{2};\omega)$, $\mathcal{G}^{0}(x_{1};nd;\omega)$,
$\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\omega)$ and
$\tilde{\mathcal{G}}^{0}(p;x_{2};\omega)$ matrices. These four matrix Green’s
functions can be easily determined using
$\mathcal{G}^{0}(x_{1};x_{2};\omega)$, which is the Green’s function in the
presence of perpendicular and uniform magnetic field in absence of an antidot
lattice given in Eqs. (5) and (8).
For propagation along the antidot lattice i.e the $x$ axis only
($y_{1}$=0=$y_{2}$, $X$=$x_{1}-x_{2}$ and
$C(\textbf{r}_{1};\textbf{r}_{2})$=$1$ ), Eqs. (5) and (8) reduce to (for
$K$-point )

$\mathcal{G}^{0}_{11\atop
22}(x_{1},x_{2};\omega)_{K}=\frac{1}{4\pi\hbar\gamma^{2}}\hskip
2.27621pt\omega\hskip
2.27621pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}\left(x_{1}-x_{2}\right)^{2}}\\\
\times\sum_{n=0}^{\infty}\frac{L_{n}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}\left(x_{1}-x_{2}\right)^{2}\right]}{\frac{\omega^{2}}{\omega^{2}_{g}}-n}$
(33)

and

$\mathcal{G}^{0}_{21\atop 12}(x_{1},x_{2};\omega)_{K}=\frac{\textit{i}\hskip
1.70717pt\omega_{g}^{2}}{8\pi\hbar\gamma^{3}}\hskip
2.27621pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}\left(x_{1}-x_{2}\right)^{2}}\left(x_{1}-x_{2}\right)\\\
\times\sum_{n=1}^{\infty}\frac{L_{n-1}^{1}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}\left(x_{1}-x_{2}\right)^{2}\right]}{\frac{\omega^{2}}{\omega^{2}_{g}}-n}$
(34)

respectively. Matrix elements of
$\mathcal{G}^{0}(x_{1};n^{\prime}d;\omega)_{K}$ can be obtained by taking
$x_{2}$=$n^{\prime}d$ in Eqs. (33) and (34), while matrix elements of
$\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\omega)_{K}$ are given in Eqs. (24) and
(25).
Similarly, the matrix $\tilde{\mathcal{G}}^{0}(p;x_{2};\omega)_{K}$ can be
obtained by using $x_{1}$=$rd$ in Eqs. (33) and (34) jointly with Eq. (20).
Note that for $\frac{\omega_{g}d}{8\gamma}>1$ , it suffices to keep only
$r=-1,0,1$ of the $r$-sum in the expressions for matrix elements of
$\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\omega)_{K}$ and
$\tilde{\mathcal{G}}^{0}(p;x_{2};\omega)_{K}$, as discussed earlier. The above
four matrix Green’s functions completely determine the full Green’s function
$\mathcal{G}(x_{1};x_{2};\omega)_{K}$ given in Eq. (19).
As discussed earlier, for temporal evolution of the wave packet we require the
time representation of the Green’s function. Hence, to find the time
representation of the full Green’s function, we have to take the Fourier
transform of $\mathcal{G}(x_{1};x_{2};\omega)_{K}$ matrix using Eqs. (10) and
(19):

$\mathcal{G}_{\mu\nu}(x_{1},x_{2};t)_{K}=\int_{-\infty}^{\infty}d\omega\hskip
1.99168pte^{-\textit{i}\omega t}\mathcal{G}^{0}(x_{1};x_{2};\omega)\\\
+\frac{\alpha d}{2\pi}\int_{-\infty}^{\infty}d\omega\hskip
1.99168pte^{-\textit{i}\omega
t}\sum_{n^{\prime}=-\infty}^{\infty}\mathcal{G}^{0}(x_{1};n^{\prime}d;\omega)\\\
\times\int_{-\frac{\pi}{d}}^{\frac{\pi}{d}}dp\hskip
5.69054pte^{-\textit{i}pn^{\prime}d}\left[1-\alpha\hskip
2.84526pt\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\omega)\right]^{-1}\tilde{\mathcal{G}}^{0}(p;x_{2};\omega).$
(35)

There are two integrals in the above equation. The first integral was
evaluated in Eqs. (11) and (12) as

$\mathcal{G}^{0}_{11\atop
22}\left(x_{1},x_{2};t\right)=-i\eta_{+}(t)\frac{\omega_{g}^{2}}{4\hskip
1.42262pt\pi\hskip 1.42262pt\gamma^{2}}\hskip
3.41432pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}\left(x_{1}-x_{2}\right)^{2}}\hskip
3.41432pt\\\
\times\sum_{n=0}^{\infty}L_{n}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}\left(x_{1}-x_{2}\right)^{2}\right]\hskip
2.84526pt\cos\left(\frac{\omega_{g}\hskip 2.27621ptl}{\gamma}\hskip
1.42262pt\frac{t}{\tau_{o}}\hskip 1.13809pt\sqrt{n}\hskip 2.27621pt\right)$
(36)

and

$\mathcal{G}^{0}_{21\atop
12}\left(x_{1},x_{2};t\right)=-i\eta_{+}(t)\frac{\omega_{g}^{3}\hskip
1.42262ptl}{8\hskip 0.56905pt\pi\hskip 0.56905pt\gamma^{3}}\hskip
0.85358pte^{-\frac{\omega^{2}_{g}}{8\gamma^{2}}\left(x_{1}-x_{2}\right)^{2}}\left(\frac{x_{1}-x_{2}}{\textit{l}}\right)\\\
\times\sum_{n=1}^{\infty}\frac{L_{n-1}^{1}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}}\left(x_{1}-x_{2}\right)^{2}\right]}{\sqrt{n}}\
\sin\left(\frac{\omega_{g}\hskip 2.27621ptl}{\gamma}\hskip
1.42262pt\frac{t}{\tau_{o}}\hskip 1.13809pt\sqrt{n}\hskip 2.27621pt\right).$
(37)

Evaluation of the second term of Eq. (35) is a lengthy process. The matrix
elements $V_{i,j}(t)$ ($i.j$=$1,2$) of the second term in time representation
are given by (see Appendix A)

$V_{11\atop
22}=-i\eta_{+}(t)\frac{\omega_{g}}{2\pi}\sum_{n^{\prime}=-\infty}^{\infty}\Bigg{\\{}\gamma_{1}\left[c_{11\atop
22}(\eta)\right]_{\eta=0}+\sum_{n=1}^{\infty}\\\
\times\cos\left(\frac{\omega_{g}\hskip 2.27621ptl}{\gamma}\hskip
1.42262pt\frac{t}{\tau_{o}}\hskip 1.13809pt\sqrt{n}\hskip
2.27621pt\right)\Bigg{[}\gamma_{1}L_{n}\left[\Upsilon\right]c_{11\atop
22}(\eta)\\\
+\frac{\gamma_{2}}{\sqrt{n}}L_{n-1}^{1}\left[\Upsilon\right]c_{21\atop
12}(\eta)\Bigg{]}_{\eta=\sqrt{n}}\Bigg{\\}}$ (38)

and

$V_{12\atop
21}=-\eta_{+}(t)\frac{\omega_{g}}{2\pi}\sum_{n^{\prime}=-\infty}^{\infty}\Bigg{\\{}i\gamma_{1}\left[c_{12\atop
21}(\eta)\right]_{\eta=0}+\sum_{n=1}^{\infty}\\\
\times\sin\left(\frac{\omega_{g}\hskip 2.27621ptl}{\gamma}\hskip
1.42262pt\frac{t}{\tau_{o}}\hskip 1.13809pt\sqrt{n}\hskip
2.27621pt\right)\Bigg{[}\gamma_{1}L_{n}\left[\Upsilon\right]c_{12\atop
21}(\eta)\\\
+\frac{\gamma_{2}}{\sqrt{n}}L_{n-1}^{1}\left[\Upsilon\right]c_{22\atop
11}(\eta)\Bigg{]}_{\eta=\sqrt{n}}\Bigg{\\}}$ (39)

where

$\Upsilon=\frac{\omega^{2}_{g}l^{2}}{4\gamma^{2}}\left(\frac{x_{1}-n^{\prime}d}{l}\right)^{2},\qquad\gamma_{1}=\frac{\alpha\omega_{g}}{4\pi\hbar\gamma^{2}}\hskip
2.27621pt\hskip 2.27621pt\hskip 2.27621pte^{-\frac{\Upsilon}{2}},$ (40)
$\gamma_{2}=\frac{\textit{i}\hskip
1.70717pt\alpha\omega_{g}^{2}l}{8\pi\hbar\gamma^{3}}\hskip
2.27621pt\left(\frac{x_{1}-n^{\prime}d}{l}\right)\hskip 2.27621pt\hskip
2.84526pte^{-\frac{\Upsilon}{2}},$ (41)

and $c_{ij}(\eta)$ (where $\eta$=$\frac{\omega}{\omega_{g}}$) are matrix
elements of (see Eq. (47): notation $q\equiv pd$)

$\int_{-\pi}^{\pi}dq\hskip
2.84526pte^{{-\textit{i}qn^{\prime}}}\left[I-\alpha\hskip
2.84526pt\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\eta)\right]^{-1}\tilde{\mathcal{G}}^{0}(p;x_{2};\eta)$

These considerations yield the time representation of the full Green’s
function, i.e. $\mathcal{G}(x_{1},x_{2};t)_{K}$ for a Landau-quantized
monolayer graphene having a one dimensional antidot lattice. In the section
VI, this Green’s function $\mathcal{G}(x_{1},x_{2};t)_{K}$ will be employed to
study the temporal dynamics of a wave packet in the lattice system.

##  VI  Wave packet dynamics along A one dimensional antidot lattice

With the solution of Eq. (35) in hand, in the form of Eqs. (36), (37), (38)
and (39), the temporal study of an electron wave packet propagating along the
axis of the antidot lattice given by

$\psi(\textbf{r},0)=\frac{f(\textbf{r})}{\sqrt{\mid c_{1}\mid^{2}+\mid
c_{2}\mid^{2}}}\begin{pmatrix}c_{1}\\\ c_{2}\end{pmatrix}$ (42)
$f(\textbf{r})=\frac{1}{l\hskip 1.13809pt\sqrt{\pi}}\hskip
1.42262ptexp\left(-\frac{x^{2}+y^{2}}{2\hskip
1.13809ptl^{2}}+\textit{i}k_{0x}x\right)\delta(y)$

can be made using Eqs. (14) and (15). In Fig. (6), results for the probability
density $\rho(x,0,t)$ of an electron along the axis of the antidot lattice
$(x,y=0)$ are plotted for four different cases to examine the effect of the
antidot lattice.
In Fig. 6, $\rho(x,0,t)$ is represented at times $t=1\tau_{o}$ and
$t=5\tau_{o}$ with solid and dashed lines, respectively. Moving from left to
right, the three columns correspond to three different pseudospin
polarizations $\\{c_{1},c_{2}\\}$, while from up to down the four rows
represent the increase in antidot strength $E_{A}=U_{0}$ from zero $meV$ to
$330\hskip 2.27621ptmeV$ with constant magnetic field strength $E_{B}=66\hskip
2.27621ptmeV$ ($E_{B}\equiv\hbar\omega_{g}$) at $B=3.3\hskip 2.27621ptTesla$.
All three columns represent the propagation of the wave packet along the $x$
axis ($y=0$). One can see that with the increase in antidot strength, the
probability density of finding the electron starts increasing along the axis
of antidot lattice. Row 1 corresponds to the situation when there is no
antidot lattice and the wave packet propagates in graphene under the effect of
the perpendicular magnetic field only. Similarly, rows 2, 3 and 4 correspond
to $E_{A}=U_{0}=13.2,66$ and $330\hskip 2.27621ptmeV$ respectively. Clearly,
the probability of finding an electron along the antidot direction increases
with the introduction of the antidot lattice on a Landau quantized graphene
sheet; this can be treated as the propagation of a wave packet through a
quantum antidot wire, which is clearly supporting the propagation through it.
Also, the probability density gets more confined with increase of antidot
strength; this can be seen in rows 2, 3 and 4 of Fig. 6. Moving from up to
down across the rows 2 to 4, the spread of the packet gets smaller with a
clear increase in magnitude of probability density, which means that the
probability of finding an electronic current along the antidot increases. This
may be referred to as collimation of the electronic beam along the axis of
antidot lattice.

Figure 6: Electron probability density $\rho(x,0,t)$ for $K$ point along the
$x$ axis in the presence of a uniform magnetic field $B=3.3\hskip
2.27621ptTesla$ having strength $E_{B}=66\hskip 2.27621ptmeV$ for
$a_{0}=k_{0x}l=1.2$ with $l=2\hskip 2.27621ptnm$; also antidot lattice
parameters, spacing $d=100\hskip 2.27621ptnm$ and radius $a=10\hskip
2.27621ptnm$. Left to right: Three columns with initial pseudospin
polarization $\\{1,0\\},\\{1,1\\}$ and $\\{1,i\\}$ respectively. Up to down:
Four rows with antidot strength $E_{A}=U_{0}=0,13,66$ and $330\hskip
2.27621ptmeV$ respectively.

Further, we explore the effect of initial pseudospin polarization on the wave
packet dynamics in the presence of an antidot lattice, as shown in Fig. 6. In
Fig. 6, it can be seen that propagation of the wave packet is strongly
affected by the change of pseudospin polarization across the rows. Columns 1
and 3 correspond to pseudospin $\\{1,0\\}$ and $\\{1,i\\}$ respectively and
the propagation corresponding to these two columns is similar because the
pseudospins are both perpendicular to the axis of the antidot lattice, and the
center of the wave packet does not propagate with an increase in time; it is
always located at $x=0$ without being affected by antidot strength, as shown
in Figs. 6 $(a),(d),(g),(h)$ and Figs. 6 $(c),(f),(i),(l)$. But, one can see
in Figs. 6 $(b),(e),(h),(k)$ (column 2) that the center of the wave packet
propagates along the axis of the antidot lattice when the initial pseudospin
polarization is along the axis of the antidot lattice i.e. $x$ axis, at
$t=5\tau_{o}$ the dashed lines clearly show the propagation of the wave packet
when compared with the solid lines at $t=1\tau_{o}$. Hence, the wave packet
only propagates along the one dimensional antidot lattice when the initial
pseudospin polarization is parallel to the axis of antidot lattice.
Experimental relevance: We now address the question of experimental relevance
of this work. In this regard, we note that it is possible to initialize the
pseudospin and study the dynamics of an electron wave packet in graphene by
means of pump-probe laser spectroscopy. Experiments have been performed where
pseudospin initialization and its subsequent relaxation have been probed on
the femtosecond scale in graphene. Danz  _et al._ (2017); Aeschlimann  _et
al._ (2017); Trushin  _et al._ (2015) Recently it has been shown that
pesudospin can be manipulated by a coupled waveguide-cavity configuration in
graphene. Han  _et al._ (2020)
Another possibility is the creation of artificial graphene (honeycomb
lattices) in cold atom systems. In these systems, it is possible to generate
effective fields and study transport in graphene-like structures Polini  _et
al._ (2013); Gomes  _et al._ (2012); Bahat-Treidel  _et al._ (2010); Soltan-
Panahi  _et al._ (2010)
Recently it has become possible to treat pseudospin as a real measurable
angular momentum at par with electron spin in photonic graphene (honeycomb
array of evanescently coupled waveguides). Song  _et al._ (2015); Gomes  _et
al._ (2012) These photonic systems allow a great degree of control of initial
conditions in the study of wave packet dynamics.

##  VII Conclusions

We have studied the evolution of a two dimensional Gaussian wave packet in a
graphene sheet placed in a uniform and perpendicular magnetic field. We have
observed that the temporal dynamics of the wave packet strongly depends on the
initial pseudospin polarization. We have shown that direction of propagation
of the wave packet can be controlled through pseudospin polarization.
Further, we have studied the effect of pseudospin on trembling motion (ZB) of
the Gaussian wave packet in a graphene sheet placed in a quantizing magnetic
field for non zero values of the initial momentum $k_{0x}$. Initially, these
ZB oscillations seemingly die out but reappear i.e. the amplitude of the ZB
oscillations start to grow again without any fundamental frequency. So, for
non-zero magnetic fields, ZB oscillations have a recurrent character and they
do not die with time. The quantized (discrete) energy spectrum, which is a
consequence of the magnetic field, is the main reason for the recurrent
character of ZB oscillations. This property is completely different from the
zero magnetic field case, in which the energy spectrum is not quantized
(discrete) and ZB of the wave packets has a transient character. Rusin and
Zawadzki (2007) Most importantly, we have found that in the presence of a
magnetic field, ZB has strong dependence on initial pseudospin polarization,
which can be regarded as pseudospinorial Zitterbewegung (PZB), a name given in
a recent study. Serna  _et al._ (2019)
We have also studied wave packet dynamics in a one dimensional antidot lattice
in graphene in the presence of the magnetic field. For this, we have
determined the Green’s function matrices of the system both in frequency and
time representations. From the poles of the Green’s function, we have found
that along the axis of antidot lattice, the energy spectrum is composed of
Landau minibands with a unit shift at $K$ and $K^{\prime}$ points of graphene.
In the time evolution of a wave packet, we find that wave packet dynamics are
highly dependent on the initial pseudospin polarization; and the center of the
Gaussian wave packet can be made to propagate along the antidot lattice by
tuning the pseudospin parallel to the axis of antidot lattice. Also, when the
strength of the antidot potential is greater than the strength of the magnetic
field, the wave packet becomes more confined in space and the probability of
finding it on the axis of antidot lattice significantly increases. Hence, we
propose that quantum antidot channels can be made such that propagation of the
wave packet can be controlled using pseudospin polarization and the strength
of the antidot potential.
Finally, we would like to point out that this work may lead to new insights
for controlling currents in both natural and artificial graphene systems by
tuning initial pseudospin polarization, both in Landau quantized graphene and
in a graphene antidot lattice in the presence of a magnetic field. This may
lead to the preparation of graphene based nano gates in which two different
pseudospin polarizations perpendicular to each other can be used to get “ON”
and “OFF” states.

##  Appendix A Contour Integral

To solve the second integral in Eq. (35), let us represent the integral by $I$

$I=\frac{\alpha
d\omega_{g}}{2\pi}\sum_{n^{\prime}=-\infty}^{\infty}\int_{-\infty}^{\infty}d\eta\hskip
1.99168pte^{-\textit{i}\omega_{g}\eta
t}\mathcal{G}^{0}(x_{1};n^{\prime}d;\eta)\\\
\int_{-\frac{\pi}{d}}^{\frac{\pi}{d}}dp\hskip
5.69054pte^{-\textit{i}pn^{\prime}d}\left[I_{2}-\alpha\hskip
2.84526pt\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\eta)\right]^{-1}\tilde{\mathcal{G}}^{0}(p;x_{2};\eta).$
(43)

($I_{2}$ is unit matrix of order 2) where

$\eta=\frac{\omega}{\omega_{g}}\qquad and\qquad d\omega=\omega_{g}d\eta$

Note that in the above integral, each Green’s function has a real pole at
$\eta$=$\pm\sqrt{n}$. The poles of the term

$T_{1}(\eta)=\left[I_{2}-\alpha\hskip
2.84526pt\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\eta)\right]^{-1}\tilde{\mathcal{G}}^{0}(p;x_{2};\eta)$
(44)

can be seen to cancel by reexpressing the Green’s functions
$\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\eta)$ and
$\tilde{\mathcal{G}}^{0}(p;x_{2};\eta)$. For this purpose, by using Eqs. (24)
and (25), one can easily write
$\alpha\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\eta)$ matrix as

$\alpha\dot{\tilde{\mathcal{G}}}^{0}(p;0,0;\eta)=\frac{1}{\prod_{n=0}^{\infty}(\eta^{2}-n)}\begin{pmatrix}a_{11}&a_{12}\\\
a_{12}&a_{11}\end{pmatrix}$ (45)

where we have defined

$a_{11}(\eta)=\alpha_{1}\eta\sum_{m=0}^{\infty}L_{m}\left[\frac{\omega^{2}_{g}l^{2}}{4\gamma^{2}}\left(j\frac{d}{l}\right)^{2}\right]\prod_{n=0\atop
n\neq m}^{\infty}(\eta^{2}-n)$

and

$a_{12}(\eta)=\alpha_{2}\sum_{m=1}^{\infty}L_{m-1}^{1}\left[\frac{\omega^{2}_{g}l^{2}}{4\gamma^{2}}\left(j\frac{d}{l}\right)^{2}\right]\prod_{n=0\atop
n\neq m+1}^{\infty}(\eta^{2}-n)$

with

$\alpha_{1}=\frac{\alpha\omega_{g}}{4\pi\hbar\gamma^{2}}\hskip 2.27621pt\hskip
2.27621pt\sum_{j=-\infty}^{\infty}e^{ipjd}\hskip
2.27621pte^{-\frac{\omega^{2}_{g}l^{2}}{8\gamma^{2}}\left(j\frac{d}{l}\right)^{2}}$

and

$\alpha_{2}=\frac{\textit{i}\hskip
1.70717pt\alpha\omega_{g}^{2}l}{8\pi\hbar\gamma^{3}}\hskip
2.27621pt\sum_{j=-\infty}^{\infty}\left(j\frac{d}{l}\right)\hskip
2.27621pte^{ipjd}\hskip
2.84526pte^{-\frac{\omega^{2}_{g}l^{2}}{8\gamma^{2}}\left(j\frac{d}{l}\right)^{2}}.$

Similarly, the $\tilde{\mathcal{G}}^{0}(p;x_{2};\eta)$ matrix can be written
as

$\tilde{\mathcal{G}}^{0}(p;x_{2};\eta)=\frac{1}{\prod_{n=0}^{\infty}(\eta^{2}-n)}\begin{pmatrix}b_{11}&b_{12}\\\
b_{12}&b_{11}\par\end{pmatrix}$ (46)

where we have defined

$b_{11}(\eta)=\beta_{1}\eta\sum_{m=0}^{\infty}L_{m}\left[\frac{\omega^{2}_{g}l^{2}}{4\gamma^{2}}\left(j\frac{d}{l}\right)^{2}\right]\prod_{n=0\atop
n\neq m}^{\infty}(\eta^{2}-n)$

and

$b_{12}(\eta)=\beta_{2}\sum_{m=1}^{\infty}L_{m-1}^{1}\left[\frac{\omega^{2}_{g}}{4\gamma^{2}l^{2}}\left(j\frac{d}{l}\right)^{2}\right]\prod_{n=0\atop
n\neq m+1}^{\infty}(\eta^{2}-n),$

with

$\beta_{1}=\frac{\omega_{g}}{4\pi\hbar\gamma^{2}}\hskip 2.27621pt\hskip
2.27621pt\sum_{j=-\infty}^{\infty}e^{ipjd}\hskip
2.27621pte^{-\frac{\omega^{2}_{g}l^{2}}{8\gamma^{2}}\left(\frac{jd-
x_{2}}{l}\right)^{2}}$

and

$\beta_{2}=\frac{\textit{i}\hskip
1.70717pt\omega_{g}^{2}l}{8\pi\hbar\gamma^{3}}\hskip
2.27621pt\sum_{j=-\infty}^{\infty}\left(\frac{jd-x_{2}}{l}\right)\hskip
2.27621pte^{ipjd}\hskip
2.84526pte^{-\frac{\omega^{2}_{g}l^{2}}{8\gamma^{2}}\left(\frac{jd-
x_{2}}{l}\right)^{2}}.$

Note that in the above equations, $n$ is Landau index and $m$ is a dummy index
for the Landau levels; the maximum value of $n$ and $m$ will be same.
Finally substituting Eqs. (45) and (46) in Eq. (44), the matrix $T_{1}$ will
become

$T_{1}(\eta)=\Bigg{[}\begin{pmatrix}\prod_{n=0}^{\infty}(\eta^{2}-n)&0\\\
0&\prod_{n=0}^{\infty}(\eta^{2}-n)\end{pmatrix}\\\
-\alpha\begin{pmatrix}a_{11}(\eta)&a_{12}(\eta)\\\
a_{12}(\eta)&a_{11}(\eta)\end{pmatrix}\Bigg{]}^{-1}\hskip
3.98337pt\begin{pmatrix}b_{11}(\eta)&b_{12}(\eta)\\\
b_{12}(\eta)&b_{21}(\eta)\end{pmatrix}.$

In the above expression, the real poles ($\eta$=$\pm\sqrt{n}$) cancel. The
above expression can be solved numerically. The resultant $2\times 2$ matrix
with $c_{ij}(\eta)$ (where $i,j=1,2$) as matrix elements can be written as
(substituting $q=pd$)

$Q(\eta)=\int_{-\pi}^{\pi}dqe^{{-\textit{i}qn^{\prime}}}T_{1}(\eta)=\begin{pmatrix}c_{11}(\eta)&c_{12}(\eta)\\\
c_{21}(\eta)&c_{22}(\eta)\end{pmatrix}.$ (47)

Note that, the above integral can be numerically calculated by applying
trapezoidal rule in the limits $-\pi$ to $\pi$ while keeping the trapezoidal
step equal to $\frac{\pi}{10}$; this gives an accuracy up to five decimal
points for each value of $n^{\prime}$.
Putting Eq. (47) in Eq. (43), the integral $I$ becomes

$I=\frac{\alpha\omega_{g}}{2\pi}\sum_{n^{\prime}=-\infty}^{\infty}\int_{-\infty}^{\infty}d\eta\hskip
1.99168pte^{-\textit{i}\omega_{g}\eta
t}\mathcal{G}^{0}(x_{1},n^{\prime}d;\eta)\hskip 1.42262pt.\hskip
1.42262ptQ(\eta).$ (48)

Using Eqs. (33) and (34), the matrix
$\alpha\mathcal{G}^{0}(x_{1},n^{\prime}d;\eta)$ can be written as

$\alpha\mathcal{G}^{0}(x_{1},n^{\prime}d;\eta)=\begin{pmatrix}\gamma_{1}\frac{1}{\eta}&0\\\
0&\gamma_{1}\frac{1}{\eta}\end{pmatrix}+\sum_{n=1}^{\infty}\\\
\times\begin{pmatrix}\gamma_{1}\eta\frac{L_{n}\left[\Upsilon\right]}{\eta^{2}-n}&\gamma_{2}\frac{L_{n-1}^{1}\left[\Upsilon\right]}{\eta^{2}-n}\\\
\gamma_{2}\frac{L_{n-1}^{1}\left[\Upsilon\right]}{\eta^{2}-n}&\gamma_{1}\eta\frac{L_{n}\left[\Upsilon\right]}{\eta^{2}-n}\end{pmatrix}$
(49)

where $\Upsilon$, $\gamma_{1}$ and $\gamma_{2}$ are defined in Eqs. (40) and
(41). Using the matrices $Q(\eta)$ and
$\alpha\mathcal{G}^{0}(x_{1},n^{\prime}d;\eta)$ in Eq. (48), and breaking the
matrix into two matrices, we get

$I=\frac{\omega_{g}}{2\pi}\sum_{n^{\prime}=-\infty}^{\infty}\int_{-\infty}^{\infty}d\eta\hskip
1.99168pte^{-\textit{i}\omega_{g}\eta
t}\begin{pmatrix}\gamma_{1}\frac{1}{\eta}c_{11}(\eta)&\gamma_{1}\frac{1}{\eta}c_{12}(\eta)\\\
\gamma_{1}\frac{1}{\eta}c_{21}(\eta)&\gamma_{1}\frac{1}{\eta}c_{22}(\eta)\end{pmatrix}\\\
+\frac{\omega_{g}}{2\pi}\sum_{n^{\prime}=-\infty}^{\infty}\int_{-\infty}^{\infty}d\eta\hskip
1.99168pte^{-\textit{i}\omega_{g}\eta
t}\begin{pmatrix}M_{11}(\eta)&M_{12}(\eta)\\\
M_{21}(\eta)&M_{22}(\eta)\end{pmatrix}$ (50)

(The first and second matrices correspond to $n$=$0$ and $n>0$ Landau
minibands, respectively) In above expression

$M_{11\atop
22}(\eta)=\sum_{n=1}^{\infty}\Bigg{(}\gamma_{1}\eta\frac{L_{n}\left[\Upsilon\right]}{\eta^{2}-n}c_{11\atop
22}(\eta)\\\
+\gamma_{2}\frac{L_{n-1}^{1}\left(\Upsilon\right)}{\eta^{2}-n}c_{21\atop
12}(\eta)\Bigg{)},$

and

$M_{12\atop
21}(\eta)=\sum_{n=1}^{\infty}\Bigg{(}\gamma_{1}\eta\frac{L_{n}\left[\Upsilon\right]}{\eta^{2}-n}c_{12\atop
21}(\eta)\\\
+\gamma_{2}\frac{L_{n-1}^{1}\left[\Upsilon\right]}{\eta^{2}-n}c_{22\atop
11}(\eta)\Bigg{)}.$

In Eq. (50), the first matrix has a pole at $\eta$=$0$, while the second
matrix has poles at $\eta$=$\pm\sqrt{n}$. We now use contour integration with
the Jordan lemma (closing the contour in the lower half plane for $t>0$) to
evaluate the integrals. Results for the two terms in Eq. (50) are

$\int_{-\infty}^{\infty}d\eta\hskip 1.99168pte^{-\textit{i}\omega_{g}\eta
t}\frac{c_{ij}(\eta)}{\eta}=-i\pi\eta_{+}(t)\left[c_{ij}(\eta)\right]_{\eta=0},$
(51) $\int_{-\infty}^{\infty}d\eta\hskip 0.56905pte^{-\textit{i}\omega_{g}\eta
t}M_{11\atop
22}(\eta)=-i\pi\eta_{+}(t)\sum_{n=1}^{\infty}\cos\left(\frac{\omega_{g}\hskip
0.56905ptl}{\gamma}\hskip 0.56905pt\frac{t}{\tau_{o}}\hskip
0.56905pt\sqrt{n}\right)\\\
\times\left[\gamma_{1}L_{n}\left[\Upsilon\right]c_{11\atop
22}(\eta)+\frac{\gamma_{2}}{\sqrt{n}}L_{n-1}^{1}\left[\Upsilon\right]c_{21\atop
12}(\eta)\right]_{\eta=\sqrt{n}}$ (52)

and

$\int_{-\infty}^{\infty}d\eta\hskip 0.56905pte^{-\textit{i}\omega_{g}\eta
t}M_{12\atop
21}(\eta)=-\pi\eta_{+}(t)\sum_{n=1}^{\infty}\sin\left(\frac{\omega_{g}\hskip
0.56905ptl}{\gamma}\hskip 0.56905pt\frac{t}{\tau_{o}}\hskip
0.56905pt\sqrt{n}\right)\\\
\times\left[\gamma_{1}L_{n}\left[\Upsilon\right]c_{12\atop
21}(\eta)+\frac{\gamma_{2}}{\sqrt{n}}L_{n-1}^{1}\left[\Upsilon\right]c_{22\atop
11}(\eta)\right]_{\eta=\sqrt{n}}.$ (53)

($\eta_{+}(t)$ is Heaviside unit step function and $i,j$=$1,2$). In the
calculation of the expressions given by Eqs. (51),(52) and (53), we have also
used

$\left[c_{11\atop 22}(\eta)\right]_{\eta=-\sqrt{n}}=\left[c_{11\atop
22}(\eta)\right]_{\eta=\sqrt{n}}$

and

$\left[c_{12\atop 21}(\eta)\right]_{\eta=-\sqrt{n}}=-\left[c_{12\atop
21}(\eta)\right]_{\eta=\sqrt{n}}$

which we found during the calculations.
Hence, Eq. (51) along with Eqs. (52) and (53) provide the complete solution
for the integral $I$, which is the time representation of second term of the
full Green’s function $\mathcal{G}(x_{1},x_{2};t)_{K}$ given in Eq. (35).

## Data Availability

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

## References

  * Xu  _et al._ (2014) X. Xu, W. Yao, X. di, and T. Heinz, Nature Physics 10, 343 (2014).
  * Mecklenburg and Regan (2011) M. Mecklenburg and B. C. Regan, Phys. Rev. Lett. 106, 116803 (2011).
  * Pesin and Macdonald (2012) D. Pesin and A. Macdonald, Nature Material 11, 409 (2012).
  * San-Jose  _et al._ (2009) P. San-Jose, E. Prada, E. McCann, and H. Schomerus, Phys. Rev. Lett. 102, 247204 (2009).
  * Song  _et al._ (2015) D. Song, V. Paltoglou, S. Liu, Y. Zhu, D. Gallardo, L. Tang, J. Xu, M. Ablowitz, N. K. Efremidis, and Z. Chen, Nature communications 6, 6272 (2015).
  * Kamalakar  _et al._ (2015) M. V. Kamalakar, C. Groenveld, A. Dankert, and S. P. Dash, in  _Nature communications_, Vol. 6 (2015) p. 6766.
  * Loss and DiVincenzo (1998) D. Loss and D. P. DiVincenzo, Phys. Rev. A 57, 120 (1998).
  * Shabbir  _et al._ (2018) B. Shabbir, M. Nadeem, Z. Dai, M. Fuhrer, Q.-K. Xue, X. Wang, and Q. Bao, Applied Physics Reviews 5, 041105 (2018).
  * Wang  _et al._ (2012) L. Wang, J. Zhu, L. Han, L. Jin, C. Zhu, E. Wang, and S. Dong, ACS Nano 6, 6659 (2012).
  * Van Tuan and Roche (2016) D. Van Tuan and S. Roche, Phys. Rev. Lett. 116, 106601 (2016).
  * Castro Neto  _et al._ (2009) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * Maksimova  _et al._ (2008) G. M. Maksimova, V. Y. Demikhovskii, and E. V. Frolova, Phys. Rev. B 78, 235321 (2008).
  * Rusin and Zawadzki (2008) T. M. Rusin and W. Zawadzki, Phys. Rev. B 78, 125419 (2008).
  * Serna  _et al._ (2019) E. Serna, I. Rodríguez Vargas, R. Pérez-Álvarez, and L. Diago-Cisneros, Journal of Applied Physics 125, 203902 (2019).
  * Wallace (1947) P. R. Wallace, Phys. Rev. 71, 622 (1947).
  * Horing and Liu (2009) N. J. M. Horing and S. Y. Liu, Journal of Physics A: Mathematical and Theoretical 42, 225301 (2009).
  * Horing (2017) N. J. M. Horing, AIP Advances 7, 065316 (2017).
  * Horing (2018) N. J. M. Horing, AIP Advances 8, 049901 (2018).
  * Demikhovskii  _et al._ (2008) V. Y. Demikhovskii, G. M. Maksimova, and E. V. Frolova, Phys. Rev. B 78, 115401 (2008).
  * Esaki and Tsu (1970) L. Esaki and R. Tsu, IBM Journal of Research and Development 14, 61 (1970).
  * Heide  _et al._ (2018) C. Heide, T. Higuchi, H. B. Weber, and P. Hommelhoff, Phys. Rev. Lett. 121, 207401 (2018).
  * Park  _et al._ (2008a) C.-H. Park, Y.-W. Son, L. Yang, M. L. Cohen, and S. G. Louie, Nano Letters 8, 2920 (2008a).
  * Choi  _et al._ (2014) S. Choi, C.-H. Park, and S. G. Louie, Phys. Rev. Lett. 113, 026802 (2014).
  * Park  _et al._ (2008b) C.-H. Park, L. Yang, Y.-W. Son, M. L. Cohen, and S. G. Louie (2008).
  * Pedersen  _et al._ (2008) T. G. Pedersen, C. Flindt, J. Pedersen, N. A. Mortensen, A.-P. Jauho, and K. Pedersen, Phys. Rev. Lett. 100, 136804 (2008).
  * Giesbers  _et al._ (2012) A. J. M. Giesbers, E. C. Peters, M. Burghard, and K. Kern, Phys. Rev. B 86, 045445 (2012).
  * Morgenstern Horing (2020) N. J. Morgenstern Horing, AIP Advances 10, 035203 (2020).
  * Danz  _et al._ (2017) T. Danz, A. Neff, J. H. Gaida, R. Bormann, C. Ropers, and S. Schäfer, Phys. Rev. B 95, 241412 (2017).
  * Aeschlimann  _et al._ (2017) S. Aeschlimann, R. Krause, M. Chávez-Cervantes, H. Bromberger, R. Jago, E. Malić, A. Al-Temimy, C. Coletti, A. Cavalleri, and I. Gierz, Phys. Rev. B 96, 020301 (2017).
  * Trushin  _et al._ (2015) M. Trushin, A. Grupp, G. Soavi, A. Budweg, D. De Fazio, U. Sassi, A. Lombardo, A. C. Ferrari, W. Belzig, A. Leitenstorfer, and D. Brida, Phys. Rev. B 92, 165429 (2015).
  * Han  _et al._ (2020) C.-D. Han, H.-Y. Xu, and Y.-C. Lai, Phys. Rev. Research 2, 033406 (2020).
  * Polini  _et al._ (2013) M. Polini, F. Guinea, M. Lewenstein, H. Manoharan, and V. Pellegrini, Nature nanotechnology 8, 625 (2013).
  * Gomes  _et al._ (2012) K. Gomes, W. Mar, W. Ko, F. Guinea, and H. Manoharan, Nature 483, 306 (2012).
  * Bahat-Treidel  _et al._ (2010) O. Bahat-Treidel, O. Peleg, M. Grobman, N. Shapira, M. Segev, and T. Pereg-Barnea, Phys. Rev. Lett. 104, 063901 (2010).
  * Soltan-Panahi  _et al._ (2010) P. Soltan-Panahi, J. Struck, P. Hauke, A. Bick, W. Plenkers, G. Meineke, C. Becker, P. Windpassinger, M. Lewenstein, and K. Sengstock, Nature Physics 7 (2010), 10.1038/nphys1916.
  * Rusin and Zawadzki (2007) T. M. Rusin and W. Zawadzki, Phys. Rev. B 76, 195439 (2007).