# Time-dependent electric transport in nodal loop semimetals

Zoltán Okvátovity<EMAIL_ADDRESS>Department of Theoretical
Physics Budapest University of Technology and Economics, 1521 Budapest,
Hungary  MTA-BME Lendület Topology and Correlation Research Group, Budapest
University of Technology and Economics, 1521 Budapest, Hungary  László
Oroszlány  Department of Physics of Complex Systems, ELTE Eötvös Loránd
University, 1117 Budapest, Hungary  MTA-BME Lendület Topology and Correlation
Research Group, Budapest University of Technology and Economics, 1521
Budapest, Hungary  Balázs Dóra  Department of Theoretical Physics Budapest
University of Technology and Economics, 1521 Budapest, Hungary  MTA-BME
Lendület Topology and Correlation Research Group, Budapest University of
Technology and Economics, 1521 Budapest, Hungary

###### Abstract

Close to the Fermi energy, nodal loop semimetals have a torus-shaped, strongly
anisotropic Fermi surface which affects their transport properties. Here we
investigate the non-equilibrium dynamics of nodal loop semimetals by going
beyond linear response and determine the time evolution of the current after
switching on a homogeneous electric field. The current grows monotonically
with time for electric fields perpendicular to the nodal loop plane however it
exhibits non-monotonical behavior for field orientations aligned within the
plane. After an initial non-universal growth $\sim Et$, the current first
reaches a plateau $\sim E$. Then, for perpendicular directions, it increases
while for in-plane directions it decreases with time to another plateau, still
$\sim E$. These features arise from interband processes. For long times or
strong electric fields, the current grows as $\sim E^{3/2}t$ or $\sim
E^{3}t^{2}$ for perpendicular or parallel electric fields, respectively. This
non-linear response represents an intraband effect where the large number of
excited quasiparticles respond to the electric field. Our analytical results
are benchmarked by the numerical evaluation of the current from continuum and
tight-binding models of nodal loop semimetals.

##  I Introduction

Recently, the investigation of properties of topologically non-trivial states
in solids has become one of the main focuses of condensed matter physics. The
theoretical prediction[1] and experimental realization[2] of topological
insulators have inspired the prediction of exotic phenomena such as the
topological magnetoelectric effect[3] and opened the door for exciting new
applications in tools for measuring fundamental constants [4], in
thermoelectric devices[5], and in architectural elements of spintronics
devices[6]. The bulk topological insulators, just as ordinary insulators, are
characterized by a gap separating the valence and conduction bands. However,
in these materials, the topological properties of bulk states, characterized
by the $\mathbb{Z}_{2}$ invariant[7], guarantee the presence of robust, spin
polarized states on the perimeter of samples.

Topology can still impact the properties of systems in the absence of a band
gap. The interplay of topology and symmetry can also stabilize robust features
in these so-called topological semimetals [8, 9, 10]. In Dirac and Weyl
semimetals[11, 8] band degeneracies near the Fermi-level occur at a discrete
set of points in the Brillouin zone. Recently these systems have been
intensively studied both theoretically [12] and experimentally culminating in
several interesting observations such as the chiral anomaly, anomalous Hall
conductivity, and Fermi arc surface states[13, 14, 15, 16].

Figure 1:  Schematic picture of the current in nodal loop semimetals when the
electric field perpendicular to the nodal loop (blue line) and parallel to it
(red line). The black dashed lines represents the timescales which split up
the time domain into four different region.

The story of topological semimetals does not end with the Dirac and Weyl
points. In certain materials, called nodal line semimetals, band crossings can
appear not only in selected points but through continuous lines which may form
closed loops or traverse the whole Brillouin zone[8, 9, 10]. Stabilizing these
nodal lines against gap opening requires the presence of some additional
symmetries[17, 10, 18, 19].

When the protecting symmetry is broken, either a finite band gap opens or the
nodal line breaks into several nodal points in the Brillouin zone. Nodal line
semimetals, in general, do not host protected edge states [9], however,
localized surface states can appear between the surface projection of the
nodal lines [18, 20] called drumhead states. These surface states owing to
their dispersionless nature may provide a fertile ground for correlation-
induced effects[21, 22].

Model systems and material realizations of nodal line semimetals have been
recently proposed in hyperhoneycomb structures[23, 24], superlattices made of
topological insulators[10, 25], alkaline-earth metal crystals[20, 26, 27] and
cold atomic systems [28, 29].

There has been also intense experimental progress to investigate the surface
properties using angle-resolved photoemission spectroscopy (ARPES) [30, 8, 31]
or magnetotransport experiments to reveal the bulk characteristics [32, 33,
34]. The characteristic features of nodal line semimetals have been identified
in various materials e.g., PbTaSe2 [35, 36], ZrSiTe and ZrSiSe [32, 37] or
Ca3P2 [20]. Although the lack of protected edge states makes nodal line
semimetals quite challenging to identify experimentally the peculiar nature of
their Fermi surface endows them with characteristic electronic and magnetic
properties[38, 39, 40].

In this paper, we investigate the non-equilibrium dynamics beyond the linear
response of nodal loop semimetals after a sudden switch of a homogeneous
electric field. We consider a simple continuum model with a single nodal loop
located in the $(p_{x},p_{y})$ plane in the momentum space at the Fermi level.
Following the works in Refs. [41, 42], we provide a detailed derivation of the
temporal behavior of the electric current for two distinct cases: when the
electric field is perpendicular to the plane of the nodal loop ($z$ direction)
or parallel to it ($x$ direction). The short time/weak electric field limit of
the current is obtained using first-order perturbation theory in the electric
field while for the long time/strong electric field limit, we use the
Dykhne–Davis–Pechukas (DDP) method. The main results are displayed in Fig. 1.
We find that in both cases, the time domain splits into four different
regions. In the ultrashort regime, the current is linear in time and electric
field, and the slope is determined by the high-energy cutoff. In the second
region, the current in the $z$ direction reaches a plateau then starts to
increase linearly while in the $x$ direction $1/t$ decay is observed. For the
third region, both currents are constant yet again. This tendency of the
current agrees with the suggested behavior from the frequency-dependent
optical conductivity from Ref. [43].

In the last temporal regime, the electric field dependence of the current
becomes non-linear. For the current, we obtained $\sim E^{3/2}t$ dependence in
$z$ direction and $\sim E^{3}t^{2}$ in $x$ direction. To test the
applicability of our result, we calculated the current by solving the
Schrödinger equation numerically and also compared it with the current
obtained from tight-binding calculations. The numerical results agree well
with the time and electric field dependence of the current obtained from
analytical calculations. These non-linear features of the electric response
are expected to be observable in transport measurements in nodal loop
semimetals.

The paper is organized as follows: in Sec. II, we briefly introduce the model
Hamiltonian and the current operators. In Secs. III and IV, the temporal and
electric field dependence of the current is investigated when the electric
field is perpendicular or parallel with the nodal loop, respectively. Then the
results of the tight-binding calculation are detailed in Sec. V. In Sec. VI,
the experimental possibilities are briefly discussed and in Sec. VII our main
results are summarized.

##  II  Model Hamiltonian and observables

We consider the effective low-energy Hamiltonian of a nodal loop semimetal
[39] as

$H=\left[\Delta-\frac{p^{2}_{x}+p^{2}_{y}}{2m}\right]\sigma_{x}+v_{\text{F}}p_{z}\sigma_{z}=P_{x}\sigma_{x}+P_{z}\sigma_{z},$
(1)

where $\sigma_{i}$’s ($i=x,z$) are the Pauli matrices, $m\approx
0.1-1~{}m_{e}$ ($m_{e}$ is the mass of an electron) is the effective mass [44,
45, 34], $\Delta\approx 0.1-1$ eV is the energy scale that defines the nodal
loop radius [46, 47] and $v_{\text{F}}\approx 10^{5}-10^{6}$ m/s is the Fermi-
velocity in the $z$ direction [44, 45]. Diagonalizing the Hamiltonian yields
the energy spectrum as $E_{\pm}({\bf p})=\pm\varepsilon_{\bf p}$ with $\pm$
the band index and $\varepsilon_{\bf
p}=\sqrt{(v_{\text{F}}p_{z})^{2}+(\Delta-(p^{2}_{x}+p^{2}_{y})/2m)^{2}}$. The
homogeneous electric field switched on at $t=0$ is introduced as a time-
dependent vector potential $\textbf{A}(t)$ through the Peierls substitution:
$\textbf{p}\rightarrow\textbf{p}-e\textbf{A}(t)$ at $t=0$. We are interested
in two different cases, when the electric field points to $x$ and $z$
direction which leads us to two different vector potentials
$\textbf{A}_{x}(t)=[Et\Theta(t),0,0]$ and
$\textbf{A}_{z}(t)=[0,0,Et\Theta(t)]$, respectively.

For each momentum $\bf{p}$ the Hamiltonian Eq. (1) with the time-dependent
vector potential represents a two-level system. Depending on the orientation
of the electric field, the instantaneous spectrum exhibits one or two
(avoided) level crossings, thus a distinct temporal behavior of the electric
current in the $x$ and $z$ directions is expected. We investigate the current
using the framework of the Landau–Zener dynamics [48], with the general time-
dependent Schrödinger equation given by

$\displaystyle H(t)$ $\displaystyle=$ $\displaystyle
P_{x}(t)\sigma_{x}+P_{z}(t)\sigma_{z},$ (2) $\displaystyle
i\hbar\partial_{t}\Psi_{p}(t)$ $\displaystyle=$ $\displaystyle
H(t)\Psi_{p}(t).$ (3)

It is convenient to perform a time-dependent unitary transformation first
[49], which diagonalizes $H(t)$, and brings us to the adiabatic basis [48]. In
the resulting equation, the positive and negative energy eigenstates are
readily distinguished, simplifying further analytic and numerical analysis.
The unitary transformation is given by

$U=\begin{bmatrix}\cos\left(\frac{\theta_{t}}{2}\right)&\sin\left(\frac{\theta_{t}}{2}\right)\\\
\sin\left(\frac{\theta_{t}}{2}\right)&-\cos\left(\frac{\theta_{t}}{2}\right)\end{bmatrix},$
(4)

where $\tan(\theta_{t})=P_{x}(t)/P_{z}(t)$. In the adiabatic basis, the
Schrödinger equation takes the form

$i\hbar\partial_{t}\Phi_{p}(t)=\left[\varepsilon_{\bf
p}(t)\sigma_{z}+F(t)\sigma_{y}\right]\Phi_{p}(t),$ (5)

where $\Psi_{p}(t)=U\Phi_{p}(t)$ and $F(t)\sigma_{y}=-i\hbar
U^{+}\partial_{t}U$ arise due to the explicit time-dependence of the unitary
transformation. $F(t)$ is referred to as diabatic coupling and is written as
[50]

$F(t)=\hbar\frac{P_{z}(t)\partial_{t}P_{x}(t)-P_{x}(t)\partial_{t}P_{z}(t)}{2\varepsilon^{2}_{\bf
p}(t)}.$ (6)

The initial condition of Eq. (5) corresponds to half filling at zero
temperature: $\Phi^{T}_{p}(t=0)=\left[0,1\right]$. The current operator for a
given momentum in the original basis is defined as $j_{p}(t)=\partial
H(t)/\partial{\bf A}(t)$ giving $j_{x}(t)=\frac{e}{m}(p_{x}-eEt)\sigma_{x}$
and $j_{z}(t)=-ev_{\text{F}}\sigma_{z}$. In the adiabatic basis, the
expression for the current operator is

$\displaystyle j_{x}(t)$ $\displaystyle=$
$\displaystyle\dfrac{e}{m}(p_{x}-eEt)\left(\sin(\theta_{t})\sigma_{z}-\cos(\theta_{t})\sigma_{x}\right),$
(7) $\displaystyle j_{z}(t)$ $\displaystyle=$ $\displaystyle-
ev_{\text{F}}\left(\cos(\theta_{t})\sigma_{z}+\sin(\theta_{t})\sigma_{x}\right).$
(8)

The contribution to the current from a given momentum mode is the expectation
value of these operators. By denoting
$\Phi^{T}_{p}(t)=\left[\alpha(t),\beta(t)\right]$, we introduce the transition
probability as $n_{\bf p}(t)=\left|\alpha(t)\right|^{2}$ which gives the
number of electrons excited from the lower to the upper band (and also number
of the holes remaining in the lower band). The contribution of a given
momentum state to the current expressed with the transition probability reads
as

$\displaystyle\left\langle j_{x}\right\rangle_{p}(t)$ $\displaystyle=$
$\displaystyle\dfrac{e}{m}\left(p_{x}-eEt\right)\sin(\theta_{t})\left(2n_{\bf
p}(t)-1\right)$ (9) $\displaystyle+$ $\displaystyle\dfrac{2\varepsilon_{\bf
p}(t)}{E}\partial_{t}n_{\bf p}(t)$ $\displaystyle\left\langle
j_{z}\right\rangle_{p}(t)$ $\displaystyle=$ $\displaystyle-
ev_{\text{F}}\cos(\theta_{t})\left(2n_{\bf p}(t)-1\right)$ (10)
$\displaystyle+$ $\displaystyle\dfrac{2\varepsilon_{\bf
p}(t)}{E}\partial_{t}n_{\bf p}(t)$

In both cases, the current consists of an intraband (first term) and interband
(second term) part which are also called conduction and polarization current
in QED terminology, respectively [51, 41, 42]. The total current is given by
the momentum integral of the momentum resolved contributions. We note that the
$n_{\bf p}(t)$ independent terms, corresponding to fully occupied or empty
states, give no contribution to the current and hence can be omitted. [41,
52]. The properties of $n_{\bf p}(t)$ and the electric current are discussed
in the following sections.

##  III  Current in the z direction

In this section, the constant electric field is aligned to the $z$ direction
i.e., it is perpendicular to the nodal loop. The time-dependent vector
potential is $\textbf{A}(t)=[0,0,Et\Theta(t)]$, and the variables in Eq. (2)
are $P_{x}=\Delta-(p^{2}_{x}+p^{2}_{y})/2m$ which remain time independent and
$P_{z}(t)=v_{\text{F}}(p_{z}-eEt)$. The evolution of the instantaneous
spectrum is plotted in Fig. 2.

Figure 2:  Visualization of the Landau–Zener dynamics when the $E$ points in
the $z$ direction. The gap between the two bands is determined by
$\hat{p}_{z}$ initially. During the time evolution, the gap starts to decrease
and when $\hat{p}_{z}-\hat{t}=0$ it closes at
$\hat{p}_{\perp}=\pm\sqrt{\hat{\Delta}}$, then start to increase. During the
time evolution, the system is driven through a quantum critical point [41].

The time-dependent Schrödinger equation and the current contribution in the
adiabatic basis for $t>0$ read as

$\displaystyle i\hbar\partial_{t}\Phi_{p}(t)$ $\displaystyle=$
$\displaystyle\left[\varepsilon_{\bf p}(t)\sigma_{z}+\dfrac{\hbar
v_{\text{F}}eE}{2\varepsilon^{2}_{\bf
p}(t)}\left(\Delta-\dfrac{p^{2}_{\perp}}{2m}\right)\sigma_{y}\right]\Phi_{p}(t),$
$\displaystyle\left\langle j_{z}\right\rangle_{p}(t)$ $\displaystyle=$
$\displaystyle-ev_{\text{F}}\dfrac{v_{\text{F}}(p_{z}-eEt)}{\varepsilon_{\bf
p}(t)}\left(2n_{\bf p}(t)-1\right)$ (12) $\displaystyle+$
$\displaystyle\dfrac{2\varepsilon_{\bf p}(t)}{E}\partial_{t}n_{\bf p}(t).$

Here, we introduced $p_{\perp}=\sqrt{p^{2}_{x}+p^{2}_{y}}$. The scaling
properties of the Schrödinger equation allow us to introduce dimensionless
variables as $\hat{t}=t/\tau_{z}$,
$\hat{\Delta}=\Delta\tau_{z}/\hbar,\hat{p}_{z}=v_{\text{F}}p_{z}\tau_{z}/\hbar$,
$\hat{p}_{\perp}=p_{\perp}\sqrt{\tau_{z}/2m\hbar}$ where the scaling factor
$\tau_{z}=\sqrt{\hbar/ev_{\text{F}}E}$ defines the natural time scale
connected to the electric field. The transition probability behaves
differently for $\hat{t}\ll 1$ and $\hat{t}\gg 1$, which also defines the
short- and long-time limits of the total current, respectively. As we show
below, for short times, the dominant contribution to the current is coming
from the polarization part while in the long-time limit, the current is
determined by the number of excited electrons in the conduction band [41, 53].

###  III.1  Short-time evolution of the current

The Schrödinger equation in Eq. (LABEL:eq:schadz) can be solved analytically
for arbitrary times and electric fields [51, 54], but it does not give an
immediately transparent solution for the transition probability. Therefore it
is more practical to obtain $n_{\bf p}(t)$ from approximate solutions in
different limits of $t$. In the short-time limit, employing first-order
perturbation theory in $E$ yields the transition probability as

$n_{\bf p}(t)=\frac{\left(\hbar v_{\text{F}}eE\right)^{2}}{4\varepsilon_{\bf
p}^{6}}\left(\Delta-\frac{p^{2}_{\perp}}{2m}\right)^{2}\sin^{2}\left(\frac{\varepsilon_{\bf
p}t}{\hbar}\right)$ (13)

which is valid except in the close vicinity of the nodal loop i.e.
$\varepsilon_{\bf p}\gg v_{\text{F}}eEt$ and resembles closely to the result
obtained for graphene [41]. To check the validity of our result, we calculated
the transition probability by solving the Schrödinger equation in Eq.
(LABEL:eq:schadz) numerically using the explicit Runge-Kutta method. We
obtained good agreement between the analytical and numerical results
visualized on the left side of Fig. 3.

Figure 3:  Transition probability for short (left) and long (right) times. The
$n_{\bf p}(t)$ resembles closely to the transition probability obtained for
graphene and Weyl semimetals [41, 42] i.e., dipolar for short times and
cylindrical for long times, but shifted to
$\hat{p}_{\perp}=\pm\sqrt{\hat{\Delta}}$.

Using this result in Eq. (12), the conduction part of the current contains
already higher (second) order terms in electric field and gives negligible
contribution to the total current. Only the polarization current contributes
to the linear order in the electric field as

$\displaystyle\left\langle j_{z}\right\rangle(t)$ $\displaystyle=$
$\displaystyle\frac{v_{\text{F}}^{2}e^{2}E}{(2\pi)^{2}\hbar^{2}}\int_{0}^{\Lambda_{z}}dp_{z}\int_{0}^{\Lambda_{\perp}}dp_{\perp}p_{\perp}\frac{(\Delta-p^{2}_{\perp}/2m)^{2}}{\varepsilon_{\bf
p}^{4}}$ (14) $\displaystyle\times$
$\displaystyle\sin\left(\frac{2\varepsilon_{\bf p}t}{\hbar}\right),$

$\Lambda_{z}$ and $\Lambda_{\perp}$ are the momentum cutoffs which arise from
the high-energy cutoff $W$ determined by the bandwidth as
$W=v_{\text{F}}\Lambda_{z}=\Lambda^{2}_{\perp}/2m$. In Eq. (14), two different
energy scales are present, $W$ and $\Delta$, which in turn determine three
different temporal regions: the ultrashort time transient response, when
$t\ll\hbar/W$, the second when $\hbar/W\ll t\ll\hbar/\Delta$ and the third
region when $\hbar/\Delta\ll t\ll\sqrt{\hbar/ev_{\text{F}}E}$.

For the ultrashort time transient response ($t\ll\hbar/W$), the current is
obtained as

$\left\langle
j_{z}\right\rangle_{p}(t)=\frac{v_{\text{F}}^{2}e^{2}Et}{(2\pi)^{2}\hbar^{3}}\frac{(\Delta-p^{2}_{\perp}/2m)^{2}}{\sqrt{(v_{\text{F}}p_{z})^{2}+(\Delta-p^{2}_{\perp}/2m)^{2}}^{3}}$
(15)

The momentum integral over $p_{z}$ and $p_{\perp}$ yields

$\left\langle
j_{z}\right\rangle(t)=\frac{me^{2}v_{\text{F}}E}{2\pi^{2}\hbar^{3}}Wt\ln\left[\sqrt{2}+1\right],$
(16)

in the $\Delta\ll W$ limit. This behavior has also been observed in Dirac and
Weyl fermions [41, 42] with the picture of classical particles accelerated by
an external electric field. These particles obey Newton’s equation with the
effective mass given by $m^{-1}_{zz}=\partial^{2}\varepsilon_{\bf
p}(t)/\partial p^{2}_{z}$.

In the second region, when $\hbar/W\ll t\ll\hbar/\Delta$, the current
saturates to a constant value similarly to graphene [41]. This can be
explained by symmetry considerations since for the electric field aligned to
the $z$ direction, the cylindrical symmetry of the system remains intact.
Therefore, the nodal loop can be thought of as two effective, graphene-like
systems with high-energy cutoff $W$ and $\Delta$, originating from states
outside or inside the nodal loop, respectively. Then, the current contribution
coming from the first graphene-like system saturates first when $\hbar/W\ll t$
[41, 42]. With increasing time, an additional linear term in time arises from
the second graphene-like system with cutoff $\Delta$. The total current can be
approximated as

$\left\langle
j_{z}\right\rangle(t)=\frac{me^{2}v_{\text{F}}E}{(2\pi)^{2}\hbar^{2}}\left[\frac{\pi^{2}}{8}+\frac{2t\Delta}{\hbar}\right].$
(17)

In the third temporal region, the additional $\Delta$ dependent part also
saturates, and the current reaches another constant value

$\left\langle j_{z}\right\rangle(t)=\frac{me^{2}v_{\text{F}}E}{16\hbar^{2}}.$
(18)

This result allows us to define a dc conductivity by taking the time-
independent current in Eq. (18) and divide it with the applied electric field
as $\sigma^{0}_{z}=j_{z}/E=me^{2}v_{\text{F}}/16\hbar^{2}$. This agrees with
the optical conductivity at $\omega\to 0$ in Ref. [43] up to a factor of two,
due to the spin degeneracy. As the frequency starts to increase, the optical
conductivity decreases with $1/\omega$ while for high frequencies, it tends to
a constant value $\sigma^{0}_{z}/2$. All these features in the optical
conductivity are in accord with our time-dependent current. Our analytical and
numerical results agree and are illustrated in Fig. 4.

Figure 4:  Temporal behavior of the current after switching on the electric
field. The blue and black dashed lines represent the polarization current in
Eq. (17) and Eq. (18), respectively. When $\hat{t}\gg 1$ the conductive
current becomes dominant, scaling as $\sim\hat{t}$.

###  III.2  Long-time, strong electric field limit

The long time, strong electric field limit, i.e.,
$t\gg\sqrt{\hbar/ev_{\text{F}}E}$ is out of the scope of perturbation theory.
Although the Schrödinger equation can be solved exactly for arbitrary times
and electric fields in terms of the parabolic cylinder functions [51], it does
not provide transparent expressions for the transition probability and the
current. Instead, we rely on the so-called WKB approach to obtain $n_{\bf
p}(t)$ which is often used to determine transition probability upon tunneling
through a barrier [55]. In practice, we use its temporal variant, the
Dykhne–Davis–Pechukas (DDP) method [56], which is also known as Landau–Dykhne
method for linear time dependence [53].

For long $t$ and large $E$, the dispersion relation has linear time dependence
with an (avoided) crossing visualized in Fig. 2. Then using the DDP method, we
obtain for the transition probability [41, 42, 53]

$n_{\bf p}(t)=\Theta(p_{z})\Theta(eEt-
p_{z})\exp\left[-\frac{\pi(\Delta-p^{2}_{\perp}/2m)^{2}}{\hbar
v_{\text{F}}eE}\right].$ (19)

The exponential term is also called the Schwinger pair production rate [57,
51]. The result in Eq. (19) agrees well with transition probability obtained
from numerical calculations, visualized in Fig. 3. Eq. (19) is applicable only
when $(p_{z},p_{z}-eEt)\gg|\Delta-p^{2}_{\perp}/2m|/v_{\text{F}}$ holds.

Using this, we calculate the total current, which is dominated by the
conductive part as

$\left\langle
j_{z}\right\rangle(t)=-\frac{2ev_{\text{F}}^{2}}{(2\pi)^{2}\hbar^{3}}\int_{-\infty}^{\infty}dp_{z}\int_{0}^{\infty}dp_{\perp}p_{\perp}\frac{p_{z}-eEt}{\varepsilon_{\bf
p}(t)}n_{\bf p}(t).$ (20)

We can estimate the overall time and field dependence by rescaling the
integral. The transition between bands occurs only if the system is driven
through the touching points as plotted in Fig. 2. This holds for $0\ll
p_{z}\ll eEt$, so the number of excited electrons can be characterized by a
longitudinally growing cylinder of length $\sim Et$ [42]. Using the scaling
parameter $\tau_{z}$, we rescale $p_{\perp}$ as
$\hat{p}_{\perp}=p_{\perp}\sqrt{\tau_{z}/(2m\hbar)}$ in the second integral of
Eq. (20) which brings out an additional $\sim E^{1/2}$ factor. Consequently,
the total current should scale with $\sim E^{3/2}t$.

The integrals over $p_{z}$ and $p_{\perp}$ yield

$\left\langle
j_{z}\right\rangle(t)=\frac{2me^{2}v_{\text{F}}E}{(2\pi)^{2}\hbar^{2}}\sqrt{\frac{v_{\text{F}}eE}{\hbar}}tf_{z}\left(\frac{\sqrt{\pi}\Delta}{\sqrt{\hbar
v_{\text{F}}eE}}\right)$ (21)

where $f_{z}(x)=(1+\text{erf}(x))/2$ with $\text{erf}(x)$, the error function
[58]. The time and electric field dependence agree with our estimation from
scaling and resembles closely to electric current in graphene [41]. The number
of excited particles is given by

$\displaystyle N(t)$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{3}\hbar^{3}}\int d^{3}p~{}n_{\bf p}(t)$ (22)
$\displaystyle=$
$\displaystyle\frac{2meE}{(2\pi)^{2}\hbar^{2}}\sqrt{\frac{v_{\text{F}}eE}{\hbar}}tf_{z}\left(\frac{\sqrt{\pi}\Delta}{\sqrt{\hbar
v_{\text{F}}eE}}\right),$

which leads to $\left\langle j_{z}\right\rangle(t)=ev_{\text{F}}N(t)$. The
current increases linearly with time due to the increasing number of electron-
hole pairs which propagate with a constant $v_{\text{F}}$ velocity in the
conduction and valence band. Due to the nodal loop, an additional $\Delta$
dependent part also arises which is responsible to a factor of $2$ enhancement
in the electric current. In the $\Delta\to\infty$ or $E\to 0$ limit, due to
the structure of the dispersion relation, the electrons can tunnel to twice as
many states as in the $\Delta\to 0$ or $E\to\infty$ limit which explains the
factor of $2$ difference in the electric current, visualized in Fig. 5.

Figure 5:  Schematic plot of the dispersion relation in the $\hat{\Delta}\to
0$ and $\hat{\Delta}\to\infty$ limit. For a given energy (green dashed line)
slightly above the gap edge, there are twice as many empty states in the
conduction band for $\hat{\Delta}$ large, therefore the current is $2$ times
larger.

##  IV  Current in the $x$ direction

To calculate the current in the $x$ direction, we use the vector potential
$\textbf{A}(t)=[Et\Theta(t),0,0]$. The variables defined in Eq. (2) are
$P_{z}=v_{\text{F}}p_{z}$ and $P_{x}=\Delta_{\text{eff}}-(p_{x}-eEt)^{2}/2m$
where $\Delta_{\text{eff}}=\Delta-p^{2}_{y}/2m$. In contrast with the previous
case, the energy-momentum dispersion relation has a $t^{2}$ temporal
dependence for $t\to\infty$, but the nodal loop is shifted in the $x$
direction during the time evolution. The temporal behavior of the
instantaneous spectrum is visualized in Fig. 6.

Figure 6:  Left panel: The time evolution of the dispersion relation is
illustrated. The time-dependent vector potential shifts the nodal loop in the
$x$ direction. Right panel: The time evolution of $\hat{p}_{x}$ at different
values of $\hat{\Delta}_{\text{eff}}$ and $\hat{p}_{z}$. When
$\hat{\Delta}_{\text{eff}}>0$, $\hat{p}_{z}$ determines the gap, and
$\hat{\Delta}_{\text{eff}}>0$ sets the location of the two minimum/maximum
points. On the other hand, when $\hat{\Delta}_{\text{eff}}<0$ the gap starts
to increase and only one minimum/maximum remains.

The time-dependent Schrödinger equation and the electric current contribution
for a given momentum mode in the adiabatic basis are given by

$\displaystyle i\hbar\partial_{t}\Phi_{p}(t)$ $\displaystyle=$
$\displaystyle\left[\varepsilon_{\bf p}(t)\sigma_{z}+\frac{\hbar
v_{\text{F}}eEp_{z}(p_{x}-eEt)}{2m\varepsilon^{2}_{\bf
p}(t)}\sigma_{y}\right]\Phi_{p}(t),$ $\displaystyle\left\langle
j_{x}\right\rangle_{p}(t)$ $\displaystyle=$
$\displaystyle\frac{e}{m}\frac{(p_{x}-eEt)}{\varepsilon_{\bf
p}(t)}\left(\Delta_{\text{eff}}-\frac{(p_{x}-eEt)^{2}}{2m}\right)\times$ (24)
$\displaystyle\times$ $\displaystyle\left(2n_{\bf
p}(t)-1\right)+\frac{2\varepsilon_{\bf p}(t)}{E}\partial_{t}n_{\bf p}(t).$

We again introduce dimensionless variables as $\hat{t}=t/\tau_{x}$, where
$\tau_{x}=\sqrt[3]{2m\hbar/(eE)^{2}}$ coming from the Schrödinger equation,
giving $\hat{\Delta}_{\text{eff}}=\Delta_{\text{eff}}\tau_{x}/\hbar$,
$\hat{p}_{x/y}=p_{x/y}\sqrt{\tau_{x}/2m\hbar}$ and
$\hat{p}_{z}=v_{\text{F}}p_{z}\tau_{x}/\hbar$. In order to analyze the
Schrödinger equation, we apply the same approximations as before, namely
first-order perturbation theory for $\hat{t}\ll 1$ and the DDP method for
$\hat{t}\gg 1$.

###  IV.1  Short-time and weak electric field limit

We again use lowest order perturbation theory to obtain the transition
probability in the $\hat{t}\ll 1$ limit.

After some straightforward algebra, we get

$n_{\bf p}(t)=\left(\frac{\hbar
v_{\text{F}}eE}{2m}\right)^{2}\frac{p^{2}_{x}p^{2}_{z}}{\varepsilon_{\bf
p}^{6}}\sin^{2}\left(\frac{\varepsilon_{\bf p}t}{\hbar}\right)$ (25)

which is valid for $\varepsilon_{\bf p}\gg(eEt)^{2}/2m$. The transition
probability is evaluated by solving the Schrödinger equation in Eq.
(LABEL:eq:schadx) numerically and is plotted in Fig. 7. For short times and
low electric fields, we obtained good agreement between the numerical and
analytical result.

Figure 7:  Transition probability for short (top) and long times when
$\Delta_{\text{eff}}>0$ (middle) and $\Delta_{\text{eff}}<0$ (bottom). For
short times, the numerical result agrees with the result from first-order
perturbation theory. For long times, when $\Delta_{\text{eff}}>0$, we have
three different regions (separated with white lines) depending on how far the
nodal loop is shifted from its original position. In this case, the modified
DDP formula works well for the transition probability. In the last case, when
$\Delta_{\text{eff}}<0$ the simple DDP formula is reliable only in the
adiabatic limit ($\hat{p}_{z}\to\infty$). However, the disagreement has only a
minor effect on the current since the number of excited particles decreases
exponentially as $\Delta_{\text{eff}}\to-\infty$.

For short times, the dominant contribution to the current comes from the
polarization term as

$\displaystyle\left\langle j_{x}\right\rangle(t)$ $\displaystyle=$
$\displaystyle\frac{e^{2}v_{\text{F}}^{2}E}{2m^{2}(2\pi)^{2}\hbar^{2}}\int_{0}^{\Lambda_{z}}dp_{z}\int_{0}^{\Lambda_{\perp}}dp_{\perp}\frac{p^{3}_{\perp}p^{2}_{z}}{\varepsilon_{\bf
p}^{4}}\times$ (26) $\displaystyle\times$
$\displaystyle\sin\left(\frac{2\varepsilon_{\bf p}t}{\hbar}\right),$

where $p_{\perp}=\sqrt{p^{2}_{x}+p^{2}_{y}}$ and the high-energy cutoff
defined similarly to $\left\langle j_{z}\right\rangle(t)$ as
$W=v_{\text{F}}\Lambda_{z}=\Lambda^{2}_{\perp}/2m$. Here again, we have two
competing energy scales which separate the time domain into three distinct
regions.

The ultrashort response i.e., $t\ll\hbar/W$, is obtained by expanding the
integral in time up to first order. For $W\gg\Delta$, we obtain

$\left\langle
j_{x}\right\rangle(t)=\frac{e^{2}EtW^{2}}{(2\pi)^{2}\hbar^{3}v_{\text{F}}}\left[1-\sqrt{2}+\ln(1+\sqrt{2})\right]$
(27)

which grows linearly with time and electric field. This result is understood
from fully classical consideration by applying Newton’s equation with
effective mass $m^{-1}_{xx}=\partial^{2}\varepsilon_{\bf p}(t)/\partial
p^{2}_{x}$. For $\hbar/W\ll t\ll\hbar/\Delta$, the current does not saturate
but starts to decay in time as $1/t$ and tends to a minimum value. In this
time interval the current reads as

$\left\langle j_{x}\right\rangle(t)=\frac{e^{2}E}{(2\pi)^{2}\hbar
v_{\text{F}}}\left[\frac{\Delta\pi^{2}}{8}+\frac{\hbar\sin^{2}(Wt/\hbar)}{3t}\right].$
(28)

The oscillating part in Eq. (28) is not universal and comes from the sharp
energy cutoff. By applying a smooth exponential cutoff, $\exp(-p/\Lambda)$
instead of the sharp one, oscillations are absent and the second term is
modified to $2\hbar W^{2}t/(3(4W^{2}t^{2}+1))$. This also decays as $\sim 1/t$
with increasing time, similarly to the sharp cutoff scheme which means that it
is a universal characteristic feature of nodal loop semimetals.

For $\hbar/\Delta\ll t\ll\sqrt[3]{2m\hbar/(eE)^{2}}$, the current tends to a
time independent constant as

$\left\langle j_{x}\right\rangle(t)=\frac{e^{2}E}{(2\pi)^{2}\hbar
v_{\text{F}}}\left[\frac{\Delta\pi^{2}}{4}+\frac{\hbar(\sin^{2}(Wt/\hbar)-1/2)}{3t}\right].$
(29)

Taking the $t\to\infty$ limit in Eq. (29), the dc response is
$\sigma^{0}_{x}=j_{x}/E=e^{2}\Delta/16\hbar^{2}v_{\text{F}}$ which agrees with
Ref. [43]. Moreover, the optical conductivity grows linearly with the
frequency with increasing frequency, which corresponds to the $1/t$ decay in
our calculation. The current is also calculated numerically and agrees with
our analytical findings in Fig. 8.

Figure 8:  Short and long time behavior of the current after switching on the
electric field in $x$ direction for large and small values of $\hat{\Delta}$.
The blue dashed line represents the $\sim 1/t$ decay of the polarization
current while the black dashed line shows the $t\to\infty$ limit of the
polarization current from Eq. (29) which is $\sim\sigma^{0}_{x}E$. The insets
show the long time behavior when the conductive current is dominant. For large
$\hat{\Delta}$ a linear term in time can arise, due to the single channel
transitions, but it is overwhelmed by the leading $\sim\hat{t}^{2}$ term with
increasing time.

###  IV.2  Long-time evolution of the current

Once more we employ the DDP method, using the results of Ref. [56], to
elucidate the long time temporal behavior of the transition probability.

For long $t$ and large $E$, the dispersion relation displays $\sim t^{2}$ time
dependence with two (avoided) crossings, visualized in Fig. 6. The diabatic
coupling reads as

$F({\bf p},t)=\frac{\hbar
v_{\text{F}}eEp_{z}(p_{x}-eEt)}{2m\varepsilon^{2}_{\bf p}(t)},$ (30)

which is an odd function in time for $eEt\gg p_{x}$ and $p_{x}$ is small. The
wave function of Eq. (LABEL:eq:schadx) is rewritten as

$\Phi_{p}(t)=\begin{bmatrix}a_{1}(t)\exp\left[i\int^{t}_{0}\varepsilon_{\bf
p}(t^{\prime})dt^{\prime}\right]\\\
a_{2}(t)\exp\left[-i\int^{t}_{0}\varepsilon_{\bf
p}(t^{\prime})dt^{\prime}\right]\end{bmatrix}$ (31)

with initial conditions $a^{i}_{1}=a_{1}(0)=0$ and $a^{i}_{2}=a_{2}(0)=1$ and
the transition probability is $n_{\bf
p}(t)=\left|a^{f}_{1}/a^{f}_{2}\right|^{2}$, where $a^{f}_{1,2}$ denotes the
final states in the $t\to\infty$ limit. In the adiabatic limit, for a single
crossing point the connection between the initial and final states is given by

$\begin{bmatrix}a^{f}_{1}\\\ a^{f}_{2}\\\ \end{bmatrix}=\begin{bmatrix}1&0\\\
e^{iD_{c}}&1\end{bmatrix}\begin{bmatrix}a^{i}_{1}\\\ a^{i}_{2}\\\
\end{bmatrix},$ (32)

where $D_{c}$ is the time integral over the classically forbidden region where
$\varepsilon_{\bf p}(t)$ is imaginary. The limits of the integral are given by
the complex crossing points where $\varepsilon(t_{c})=0$ as

$t_{c}=\pm\frac{\sqrt{2m}}{eE}\sqrt{\Delta_{\text{eff}}\pm
i|v_{\text{F}}p_{z}|},$ (33)

and $D_{c}$ is

$D_{c}=\frac{2}{\sqrt{\beta}}\int^{\sqrt{\mu+i}}_{0}dz\sqrt{1+(\mu-z^{2})^{2}},$
(34)

where $\beta=(\hbar eE)^{2}/(2m|v_{\text{F}}p_{z}|^{3})$ and
$\mu=\Delta_{\text{eff}}/|v_{\text{F}}p_{z}|$. We can identify $\sqrt{\beta}$
as an adiabaticity parameter since as $\sqrt{\beta}\to 0$ the transition
probability also tends to $0$ [56]. The integral yields [56]

$D_{c}=\frac{\pi}{2\sqrt{\beta}}\sqrt{(\mu+i)(\mu^{2}+1)}_{2}F_{1}\left[\frac{1}{2},-\frac{1}{2};2;\frac{\mu+i}{\mu-i}\right],$
(35)

where ${}_{2}F_{1}(a,b;c;x)$ is the hypergeometric function[58]. When only one
channel is available to tunnel through the barrier i.e., when
$\Delta_{\text{eff}}>0$ and $|p_{x}|\ll\sqrt{2m\Delta_{\text{eff}}}$ or
$|p_{x}-eEt|\ll\sqrt{2m\Delta_{\text{eff}}}$, the transition probability is

$n^{1}_{\bf
p}(t)=\Theta\left(\sqrt{2m\Delta_{\text{eff}}}-|p_{x}|\right)\exp\left(-2\text{Im}\left[D_{c}\right]\right).$
(36)

To take into account both tunneling channels when
$\sqrt{2m\Delta_{\text{eff}}}\ll p_{x}\ll eEt-\sqrt{2m\Delta_{\text{eff}}}$
for $\Delta_{\text{eff}}>0$, we apply the matrix from Eq. (32) twice for the
two crossing points, but for the second time, the off-diagonal elements pick
up an extra minus sign. This extra minus sign arises due to the diabatic
coupling in Eq. (30), which has a freedom in its sign [59]. By fixing $F({\bf
p},t)$ to be positive for the first crossing point, then we have to insert a
negative sign for the off-diagonal terms since $F({\bf p},t)$ is an odd
function in time for $eEt\gg p_{x}$ [59]. The final states are given by

$\begin{bmatrix}a^{f}_{1}\\\ a^{f}_{2}\\\ \end{bmatrix}=\begin{bmatrix}1&0\\\
e^{iD_{c}}-e^{-iD^{*}_{c}}&1\end{bmatrix}\begin{bmatrix}a^{i}_{1}\\\
a^{i}_{2}\\\ \end{bmatrix},$ (37)

For the two channel tunneling case, the transition probability yields

$\displaystyle n^{2+}_{\bf p}(t)$ $\displaystyle=$
$\displaystyle\Theta\left(p_{x}-\sqrt{2m\Delta_{\text{eff}}}\right)\Theta\left(eEt-\sqrt{2m\Delta_{\text{eff}}}-p_{x}\right)$
(38) $\displaystyle\times$ $\displaystyle
4\sin^{2}\left(\text{Re}\left[D_{c}\right]\right)e^{-2\text{Im}\left[D_{c}\right]}.$

This method is also applicable for the opposite, $\Delta_{\text{eff}}<0$ case,
though with a modified time-dependent part due to the lack of the single
channel transition regions. In this case, the transition probability reads as

$n^{2-}_{\bf p}(t)=\Theta\left(p_{x}\right)\Theta\left(eEt-
p_{x}\right)4\sin^{2}\left(\text{Re}\left[D_{c}\right]\right)e^{-2\text{Im}\left[D_{c}\right]}.$
(39)

These transition probabilities are expected to work in principle only in the
adiabatic limit i.e., $\sqrt{\beta}\to 0$ [56, 59, 60]. However, for
$\Delta_{\text{eff}}>0$ case, we can indeed go beyond the adiabatic limit and
obtain results superior to $n^{2+}_{\bf p}(t)$ in Eq. (38). For large
$\Delta_{\text{eff}}$, we can treat our system as two, independent
Landau–Zener models. Its $S$ matrix is known exactly between the initial and
final states as[61]

$\begin{bmatrix}a^{f}_{1}\\\ a^{f}_{2}\\\
\end{bmatrix}=\begin{bmatrix}\sqrt{1-|R|^{2}}e^{-i\chi}&-R\\\
R&\sqrt{1-|R|^{2}}e^{+i\chi}\end{bmatrix}\begin{bmatrix}a^{i}_{1}\\\
a^{i}_{2}\\\ \end{bmatrix},$ (40)

where $R=\exp[iD_{c}]$ with the phase factor

$\chi=\frac{\pi}{4}+\frac{\lambda}{2}\ln\left(\frac{\lambda}{2}\right)-\frac{\lambda}{2}+\arg\left[\Gamma\left(1-\frac{i\lambda}{2}\right)\right]$
(41)

and $\lambda=1/(2\sqrt{\beta\mu})$. Applying the $S$ matrix for both crossing
points in the same way as in Eq. (37), we end up with

$\displaystyle n^{2\text{mod}}_{\bf p}(t)$ $\displaystyle=$
$\displaystyle\Theta$
$\displaystyle\left(p_{x}-\sqrt{2m\Delta_{\text{eff}}}\right)\Theta\left(eEt-\sqrt{2m\Delta_{\text{eff}}}-p_{x}\right)$
$\displaystyle\times$ $\displaystyle 4$
$\displaystyle\sin^{2}\left(\text{Re}\left[D_{c}\right]+\chi\right)e^{-2\text{Im}\left[D_{c}\right]}\left(1-e^{-2\text{Im}\left[D_{c}\right]}\right).$

This expression, though similar to $n^{2+}_{\bf p}(t)$ in Eq. (38), works much
better and is used instead of $n^{2+}_{\bf p}(t)$ in the following. We
compared the analytical transition probabilities from DDP and the numerical
result in the second and third row of Fig. 7. For $\Delta_{\text{eff}}>0$,
$n^{1}_{\bf p}(t)$ and $n^{2\text{mod}}_{\bf p}(t)$ agrees with the numerical
calculations even beyond the adiabatic limit (i.e. $p_{z}\to 0$) as advertised
above. Outside the nodal loop ($\Delta_{\text{eff}}<0$), the numerical and
analytical results differ. However this region gives a tiny contribution to
the current since as $\Delta_{\text{eff}}\to-\infty$ the transition
probability is exponentially suppressed due to the increasing gap. Finally,
the results for transition probability are summarized as

$\displaystyle n_{\bf p}(t)$ $\displaystyle=$
$\displaystyle\Theta(\Delta_{\text{eff}})\left(n^{1}_{\bf p}(t)+n^{1}_{eEt-
p_{x}}(t)+n^{2\text{mod}}_{\bf p}(t)\right)$ (43) $\displaystyle+$
$\displaystyle\Theta(-\Delta_{\text{eff}})n^{2-}_{\bf p}(t)$

The temporal and electric field scaling of the current can be estimated
similarly to the previous case. Rescaling the $p_{y}$ and $p_{z}$ variables
with $\tau_{x}$ gives a factor of $E$. The excitation to the upper band occurs
only upon complete nonadiabatic passage through the touching points which
holds for $0\ll p_{x}\ll eEt$ so the number of excited electrons scales with
$\sim Et$. The velocity itself is explicitly time-dependent which also brings
a factor of $\sim Et$. Combining these, we predict $E^{3}t^{2}$ scaling for
the total current.

This prediction is reproduced by using (LABEL:eq:nptxl4) for the conductive
current contribution in Eq. (24) as

$\left\langle j_{x}\right\rangle(t)=-\frac{2}{(2\pi)^{3}\hbar^{3}}\int
d^{3}p(\partial_{p_{x}}\varepsilon_{\bf p}(t))n_{\bf p}(t),$ (44)

while the polarization part gives only a subleading contribution $\sim E$. The
one channel tunnelings contribute with $\sim t$ terms to the current for
$t\to\infty$ while the other part gives the dominant, $\sim t^{2}$
contribution as

$\left\langle
j_{x}\right\rangle(t)=\frac{e^{4}E^{3}}{(2\pi)^{3}m\hbar^{2}v_{\text{F}}}t^{2}f_{x}\left(\sqrt[3]{\frac{2m}{(\hbar
eE)^{2}}}\Delta\right),$ (45)

where $f_{x}(x)$ contains the result of $p_{y}$ and $p_{z}$ integrals and
satisfies $f_{x}(0)=\textmd{const.}$ and $f_{x}(x\to\infty)\sim x^{2/3}$. This
means that for small and large $\hat{\Delta}$, the current scales as
$E^{3}t^{2}$ and $\Delta^{2/3}E^{23/9}t^{2}$, respectively. Given the large
electric field exponent close to 3, these can be written to a good
approximation as $E^{3}t^{2}$.

The field and time dependence of the electric current agree with our previous
estimation. The number of the excited particles are estimated as

$N(t)=\int d^{3}pn_{\bf p}(t)\approx\frac{e^{2}E^{2}t}{(2\pi)^{3}\hbar
v_{\text{F}}}f_{x}\left(\sqrt[3]{\frac{2m}{(\hbar eE)^{2}}}\Delta\right)$ (46)

which leads to $\left\langle j_{x}\right\rangle(t)=ev(t)N(t)$ where
$v(t)=eEt/m$ is the time-dependent part of the velocity operator. Therefore,
the current comes from the increasing number of particles excited to the upper
band, but also these particles are accelerated by the electric field.

##  V  Tight-binding model

Figure 9:  The electric current from the tight-binding model for $\hat{E}$ in
the $z$ direction. Top panel: Short-time response with $\delta/\gamma=1.5$.
The brown and light blue dashed lines represent the two time scales $1/\delta$
and $1/\Delta=1/(2\gamma-\delta)$, respectively. The polarization current
dominates, represented by the black dashed line. Bottom panel: Long-time
response with with $\delta/\gamma=1.5$. The orange dashed line represents the
$\sim E^{3/2}t$ prediction. For large electric fields and long times, Bloch
oscillations start to kick in. Figure 10:  The electric current from the
tight-binding model for $\hat{E}$ pointing to the $x$ direction. Top panel:
Short-time response for $\delta/\gamma=1.5$. The brown dashed line represents
the border of the ultrashort time domain at $1/\delta$. The black dashed line
represents the polarization current. Bottom panel: Long-time response with
$\delta/\gamma=1.97$. The orange dashed line represents our prediction $\sim
E^{3}t^{2}$. The additional constant term, $j^{0}_{x}$ is the dimensionless
static current coming from the polarization current. For longer times, the
Bloch oscillations would kick in.

To validate the results obtained from the continuum model, we performed tight-
binding calculations [29, 62] based on the lattice Hamiltonian defined as

$H_{\text{TB}}=\left[\delta-\gamma\cos(k_{x}a)-\gamma\cos(k_{y}a)\right]\sigma_{x}-\gamma\sin(k_{z}a)\sigma_{z},$
(47)

where $k_{i}$’s ($i=x,y,z$) are the wave numbers in different directions,
$\gamma$ is the hopping integral, $a$ is the lattice constant and $\delta$
determines the radius of the nodal loop. The nodal loops are located in the
$k_{x}-k_{y}$ plane for $k_{z}a=0$ and $\pi$. In order to avoid their overlaps
and the concomitant Bloch oscillations, we use $aeEt/\hbar<\pi$ in the
numerics. We can identify the parameters of the continuum model by expanding
Eq. (47) in $k_{i}$ up to second order which gives $\Delta=2\gamma-\delta$,
$m=\hbar^{2}/a^{2}\gamma$ and $v_{\text{F}}=\gamma a/\hbar$. We solved the
time-dependent Schrödinger equation on a finite cubic lattice with $300$ unit
cells in each direction with an adaptive grid. We introduce the dimensionless
electric field as $\hat{E}=eEv_{\text{F}}\hbar/\gamma^{2}=eEa/\gamma$. The
obtained current for various electric fields are shown in Figs. 9 and 10,
displaying nice agreement with our analytical findings.

In the $z$ direction, the current has two plateaus for short times, while for
long times it scales with $\sim E^{3/2}t$. In the $x$ direction, the current
decays for short times and grows $\sim E^{3}t^{2}$ for later times. Due to
additional linear subleading terms scaling is adversely affected in the $x$
direction as compared to the $z$ orientation.

We calculated numerically the polarization current using first-order
perturbation theory in the electric field which is represented with the black
dashed lines in Fig. 9 and Fig. 10. Similarly to the continuum model, the
$t\to\infty$ limit defines the static current as $\sigma^{0}_{x}E$ and
$\sigma^{0}_{z}E$ for $x$ and $z$ direction, respectively. Overall, the long
and short-time behavior of the current from tight-binding calculations agrees
well with the analytical and numerical results from the continuum model.

##  VI  Experimental possibilities

So far, we discussed and evaluated the real evolution of the electric current,
after switching on an electric field. Here we briefly discuss possible
experiments related to measuring the current.

First, by following time dependence of the current, the observation of the
characteristic crossover timescales allows us to determine various parameters
of nodal loop semimetals. The electric field dependent timescales $\tau_{x/z}$
separate the short- and long-time behavior. By experimentally identifying
these timescales, the Fermi velocity and the effective mass can be obtained as
$v_{\text{F}}=\hbar/(eE\tau^{2}_{z})$ and $m=\tau^{3}_{x}(eE)^{2}/(2\hbar)$.
The short-time evolution of the current also contains useful information about
the nodal loop parameter $\Delta$ and the bandwidth. When the electric field
aligned to the $z$ direction, two plateaus are observable when $t>\hbar/W$ and
$t>\hbar/\Delta$. Identifying the temporal crossing points when the current
reaches these plateaus allows us to estimate the bandwidth and the nodal loop
parameter. Similarly, the peak in $j_{x}$ current before the $1/t$ decay also
determines the bandwidth.

One way to detect the nonlinear electric response is to experimentally realize
these nodal loop semimetals in ultracold atomic systems [63, 64]. Applying a
weak magnetic-field gradient [63] or tuning the frequency difference of laser
waves responsible for the optical lattice [65] create a constant force, which
is equivalent to switching on an electric field in solid-state systems. The
main advantages of these measurements are the absence of scattering and
dissipation and strong electric fields are not needed to obtain the nonlinear
response. According to our tight-binding calculations, the current from the
electron-hole pair creation is observable before the Bloch oscillations kick
in.

Another way to obtain the short and late time electric response of nodal loop
semimetals is to measure the current in a solid-state realization. However, in
these systems scattering processes appear due to the phonons, impurities, etc.
The Drude picture provides a simple way to interpret our results in the
presence of impurities: The charge carriers move ballistically until a
momentum transfer happens due to scattering processes. The average lifetime
can be characterized by the relaxation time $\tau_{sc}$, which introduces a
restricting time scale to our system. For $t>\tau_{sc}$, the current become
stationary, which is estimated by substituting $t\to\tau_{sc}$ in the
corresponding expressions for the current. Consequently, long-time features of
the electric current are only observable if $\tau_{x/z}<\tau_{sc}$, which is
equivalent to $E>E_{x/z,c}$, where $E_{x/z,c}$ is the critical electric field
which separates the linear from the non-linear regions. It is defined as
$E_{z,c}=\hbar/(ev_{\text{F}}\tau^{2}_{sc})$ for $z$ direction and
$E_{x,c}=\sqrt{2m\hbar/(e^{2}\tau^{3}_{sc})}$ for $x$ direction. The
scattering time is estimated as $\tau_{sc}\sim 10^{-2}-10^{-1}$ ps [66, 45],
which implies that the minimal electric field required to observe nonlinear
transport is $E_{x/z,c}\sim 10^{5}-10^{7}$ V/m. For $E_{E}>E_{x/z,c}$, the
current changes its slope as a function of the electric field, but an even
larger electric field window may be required to obtain the corresponding
exponents.

##  VII  Conclusion

Table 1:  Time evolution and electric field dependence of the current in nodal loop semimetals. | Ultrashort response | Kubo I | Kubo II | Long time response
---|---|---|---|---
| $t\ll\hbar/W$ | $\hbar/W\ll t\ll\hbar/\Delta$ | $\hbar/\Delta\ll t\ll\tau_{x/z}$ | $\tau_{x/z}\ll t$
$x$ direction | $j_{x}\sim Et$ | $j_{x}\sim E(const.+1/t)$ | $j_{x}\sim E$ | $j_{x}\sim E^{3}t^{2}$
$z$ direction | $j_{z}\sim Et$ | $j_{z}\sim E(const.+t)$ | $j_{z}\sim E$ | $j_{z}\sim E^{3/2}t$

In this paper, we have investigated the time evolution of the non-equilibrium
electric current of nodal loop semimetals after switching on a homogeneous
electric field. We considered the two characteristic cases, namely when the
electric field is within the plane of the loop or perpendicular to it. To
calculate the current, we determined the transition probabilities by using a
variety of techniques, including first-order perturbation theory for short
times and weak electric fields and the DDP method for long times and strong
electric fields. Based on this, the intra- and interband contributions to the
electric current are identified.

For short times and weak electric fields, the interband processes dominate the
current for both electric field orientations, and the ensuing time dependence
can also be formally understood from a Kubo formula calculation of the optical
conductivity. For long times and strong fields, on the other hand, the current
originates from intraband processes, namely by the increasing number of
excited quasiparticles in the initially empty upper band. In addition, we
benchmarked our analytical results by the numerical solution of the
Schrödinger equation both in the continuum limit and for tight-binding models.
Our results are summarized in Table 1.

###### Acknowledgements.

This research is supported by the National Research, Development and
Innovation Office - NKFIH within the Quantum Information National Laboratory
of Hungary and Quantum Technology National Excellence Program (Project No.
2017-1.2.1-NKP-2017-00001), K119442, K134437, K131938, FK124723, and by the
Romanian National Authority for Scientific Research and Innovation, UEFISCDI,
under Project No. PN-III-P4-ID-PCE-2020-0277. L. O. also acknowledges support
from of the NRDI Office of Hungary and the Hungarian Academy of Sciences
through the Bolyai and Bolyai+ scholarships.

## References

  * Bernevig et al. [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006), URL https://doi.org/10.1126/science.1133734.
  * Konig et al. [2007] M. Konig, S. Wiedmann, C. Brune, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007), URL https://doi.org/10.1126/science.1148047.
  * Qi et al. [2009] X.-L. Qi, R. Li, J. Zang, and S.-C. Zhang, Science 323, 1184 (2009), URL https://doi.org/10.1126/science.1167747.
  * Maciejko et al. [2010] J. Maciejko, X.-L. Qi, H. D. Drew, and S.-C. Zhang, Physical Review Letters 105, 166803 (2010).
  * Xu et al. [2017] N. Xu, Y. Xu, and J. Zhu, npj Quantum Materials 2, 51 (2017).
  * Brüne et al. [2012] C. Brüne, A. Roth, H. Buhmann, E. M. Hankiewicz, L. W. Molenkamp, J. Maciejko, X.-L. Qi, and S.-C. Zhang, Nature Physics 8, 485 (2012), ISSN 1745-2481, URL https://www.nature.com/articles/nphys2322.
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Reviews of Modern Physics 82, 3045 (2010), URL https://doi.org/10.1103/revmodphys.82.3045.
  * Yang et al. [2018] S.-Y. Yang, H. Yang, E. Derunova, S. S. P. Parkin, B. Yan, and M. N. Ali, Advances in Physics: X 3, 1414631 (2018), URL https://doi.org/10.1080/23746149.2017.1414631.
  * Fang et al. [2016] C. Fang, H. Weng, X. Dai, and Z. Fang, Chinese Physics B 25, 117106 (2016), URL https://doi.org/10.1088/1674-1056/25/11/117106.
  * Burkov et al. [2011] A. A. Burkov, M. D. Hook, and L. Balents, Physical Review B 84, 235126 (2011), URL https://doi.org/10.1103/physrevb.84.235126.
  * Armitage et al. [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Reviews of Modern Physics 90, 015001 (2018), URL https://doi.org/10.1103/revmodphys.90.015001.
  * Turner and Vishwanath [2013] A. Turner and A. Vishwanath, _Topological Insulators: Chapter 11. Beyond Band Insulators: Topology of Semimetals and Interacting Phases_ , Contemporary Concepts of Condensed Matter Science (Elsevier Science, 2013), ISBN 9780128086926.
  * Yasuoka et al. [2017] H. Yasuoka, T. Kubo, Y. Kishimoto, D. Kasinathan, M. Schmidt, B. Yan, Y. Zhang, H. Tou, C. Felser, A. P. Mackenzie, et al., Physical Review Letters 118, 236403 (2017), URL https://doi.org/10.1103/physrevlett.118.236403.
  * Burkov and Balents [2011] A. A. Burkov and L. Balents, Physical Review Letters 107, 127205 (2011), URL https://doi.org/10.1103/physrevlett.107.127205.
  * Xu et al. [2015] S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, et al., Science 349, 613 (2015), URL https://doi.org/10.1126/science.aaa9297.
  * Huang et al. [2015] S.-M. Huang, S.-Y. Xu, I. Belopolski, C.-C. Lee, G. Chang, B. Wang, N. Alidoust, G. Bian, M. Neupane, C. Zhang, et al., Nature Communications 6, 8373 (2015), URL https://doi.org/10.1038/ncomms8373.
  * Chiu et al. [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Reviews of Modern Physics 88, 035005 (2016), URL https://doi.org/10.1103/revmodphys.88.035005.
  * Fang et al. [2015] C. Fang, Y. Chen, H.-Y. Kee, and L. Fu, Physical Review B 92, 081201(R) (2015), URL https://doi.org/10.1103/physrevb.92.081201.
  * Chiu and Schnyder [2014] C.-K. Chiu and A. P. Schnyder, Physical Review B 90, 205136 (2014), URL https://doi.org/10.1103/physrevb.90.205136.
  * Chan et al. [2016] Y.-H. Chan, C.-K. Chiu, M. Y. Chou, and A. P. Schnyder, Physical Review B 93, 205132 (2016), URL https://doi.org/10.1103/physrevb.93.205132.
  * Liu and Balents [2017] J. Liu and L. Balents, Physical Review B 95, 075426 (2017), URL https://doi.org/10.1103/physrevb.95.075426.
  * Le et al. [2020] T. Le, Y. Sun, H.-K. Jin, L. Che, L. Yin, J. Li, G. Pang, C. Xu, L. Zhao, S. Kittaka, et al., Science Bulletin 65, 1349 (2020), URL https://doi.org/10.1016/j.scib.2020.04.039.
  * Mullen et al. [2015] K. Mullen, B. Uchoa, and D. T. Glatzhofer, Physical Review Letters 115, 026403 (2015), URL https://doi.org/10.1103/physrevlett.115.026403.
  * Ezawa [2016] M. Ezawa, Physical Review Letters 116, 127202 (2016), URL https://doi.org/10.1103/physrevlett.116.127202.
  * Phillips and Aji [2014] M. Phillips and V. Aji, Physical Review B 90, 115111 (2014), URL https://doi.org/10.1103/physrevb.90.115111.
  * Hirayama et al. [2017] M. Hirayama, R. Okugawa, T. Miyake, and S. Murakami, Nature Communications 8, 14022 (2017), URL https://doi.org/10.1038/ncomms14022.
  * Du et al. [2017] Y. Du, F. Tang, D. Wang, L. Sheng, E. jun Kan, C.-G. Duan, S. Y. Savrasov, and X. Wan, npj Quantum Materials 2, 41535 (2017), URL https://doi.org/10.1038/s41535-016-0005-4.
  * Yin et al. [2018] C. Yin, L. Li, and S. Chen, Physical Review A 97, 013604 (2018), URL https://doi.org/10.1103/physreva.97.013604.
  * Zhang et al. [2016] D.-W. Zhang, Y. X. Zhao, R.-B. Liu, Z.-Y. Xue, S.-L. Zhu, and Z. D. Wang, Physical Review A 93, 043617 (2016), URL https://doi.org/10.1103/physreva.93.043617.
  * Kim et al. [2015] Y. Kim, B. J. Wieder, C. L. Kane, and A. M. Rappe, Physical Review Letters 115, 036806 (2015), URL https://doi.org/10.1103/physrevlett.115.036806.
  * Bian et al. [2016a] G. Bian, T.-R. Chang, H. Zheng, S. Velury, S.-Y. Xu, T. Neupert, C.-K. Chiu, S.-M. Huang, D. S. Sanchez, I. Belopolski, et al., Physical Review B 93, 121113(R) (2016a), URL https://doi.org/10.1103/physrevb.93.121113.
  * Hu et al. [2016] J. Hu, Z. Tang, J. Liu, X. Liu, Y. Zhu, D. Graf, K. Myhro, S. Tran, C. N. Lau, J. Wei, et al., Physical Review Letters 117, 016602 (2016), URL https://doi.org/10.1103/physrevlett.117.016602.
  * Emmanouilidou et al. [2017] E. Emmanouilidou, B. Shen, X. Deng, T.-R. Chang, A. Shi, G. Kotliar, S.-Y. Xu, and N. Ni, Physical Review B 95, 245113 (2017), URL https://doi.org/10.1103/physrevb.95.245113.
  * Pezzini et al. [2017] S. Pezzini, M. R. van Delft, L. M. Schoop, B. V. Lotsch, A. Carrington, M. I. Katsnelson, N. E. Hussey, and S. Wiedmann, Nature Physics 14, 178 (2017), URL https://doi.org/10.1038/nphys4306.
  * Bian et al. [2016b] G. Bian, T.-R. Chang, R. Sankar, S.-Y. Xu, H. Zheng, T. Neupert, C.-K. Chiu, S.-M. Huang, G. Chang, I. Belopolski, et al., Nature Communications 7, 10556 (2016b), URL https://doi.org/10.1038/ncomms10556.
  * Chang et al. [2016] T.-R. Chang, P.-J. Chen, G. Bian, S.-M. Huang, H. Zheng, T. Neupert, R. Sankar, S.-Y. Xu, I. Belopolski, G. Chang, et al., Physical Review B 93, 245130 (2016), URL https://doi.org/10.1103/physrevb.93.245130.
  * Muechler et al. [2020] L. Muechler, A. Topp, R. Queiroz, M. Krivenkov, A. Varykhalov, J. Cano, C. R. Ast, and L. M. Schoop, Physical Review X 10, 011026 (2020), URL https://doi.org/10.1103/physrevx.10.011026.
  * Tai and Lee [2021] T. Tai and C. H. Lee, Physical Review B 103, 195125 (2021), URL https://doi.org/10.1103/physrevb.103.195125.
  * Oroszlány et al. [2018] L. Oroszlány, B. Dóra, J. Cserti, and A. Cortijo, Physical Review B 97, 205107 (2018), URL https://doi.org/10.1103/physrevb.97.205107.
  * Martín-Ruiz and Cortijo [2018] A. Martín-Ruiz and A. Cortijo, Physical Review B 98, 155125 (2018), URL https://doi.org/10.1103/physrevb.98.155125.
  * Dóra and Moessner [2010] B. Dóra and R. Moessner, Physical Review B 81, 165431 (2010), URL https://doi.org/10.1103/physrevb.81.165431.
  * Vajna et al. [2015] S. Vajna, B. Dóra, and R. Moessner, Physical Review B 92, 085122 (2015), URL https://doi.org/10.1103/physrevb.92.085122.
  * Barati and Abedinpour [2017] S. Barati and S. H. Abedinpour, Physical Review B 96, 155150 (2017), URL https://doi.org/10.1103/physrevb.96.155150.
  * Singha et al. [2017] R. Singha, A. K. Pariari, B. Satpati, and P. Mandal, Proceedings of the National Academy of Sciences 114, 2468 (2017), URL https://doi.org/10.1073/pnas.1618004114.
  * Rudenko and Yuan [2020] A. N. Rudenko and S. Yuan, Physical Review B 101, 115127 (2020), URL https://doi.org/10.1103/physrevb.101.115127.
  * Fu et al. [2019] B.-B. Fu, C.-J. Yi, T.-T. Zhang, M. Caputo, J.-Z. Ma, X. Gao, B. Q. Lv, L.-Y. Kong, Y.-B. Huang, P. Richard, et al., Science Advances 5, 1126 (2019), URL https://doi.org/10.1126/sciadv.aau6459.
  * Hosen et al. [2017] M. M. Hosen, K. Dimitri, I. Belopolski, P. Maldonado, R. Sankar, N. Dhakal, G. Dhakal, T. Cole, P. M. Oppeneer, D. Kaczorowski, et al., Physical Review B 95, 161101(R) (2017), URL https://doi.org/10.1103/physrevb.95.161101.
  * Vitanov and Garraway [1996] N. V. Vitanov and B. M. Garraway, Physical Review A 53, 4288 (1996), URL https://doi.org/10.1103/physreva.53.4288.
  * Cohen and McGady [2008] T. D. Cohen and D. A. McGady, Physical Review D 78, 036008 (2008), URL https://doi.org/10.1103/physrevd.78.036008.
  * Lehto and Suominen [2012] J. Lehto and K.-A. Suominen, Physical Review A 86, 033415 (2012), URL https://doi.org/10.1103/physreva.86.033415.
  * Tanji [2009] N. Tanji, Annals of Physics 324, 1691 (2009), URL https://doi.org/10.1016/j.aop.2009.03.012.
  * Ashcroft [1976] N. W. Ashcroft, _Solid State Physics_ (Cengage Learning, Inc, 1976), ISBN 0030839939.
  * Boross et al. [2011] P. Boross, B. Dóra, and R. Moessner, physica status solidi (b) 248, 2627 (2011), URL https://doi.org/10.1002/pssb.201100184.
  * Gavrilov and Gitman [1996] S. P. Gavrilov and D. M. Gitman, Physical Review D 53, 7162 (1996), URL https://doi.org/10.1103/physrevd.53.7162.
  * Casher et al. [1979] A. Casher, H. Neuberger, and S. Nussinov, Physical Review D 20, 179 (1979), URL https://doi.org/10.1103/physrevd.20.179.
  * Suominen [1992] K.-A. Suominen, Optics Communications 93, 126 (1992), URL https://doi.org/10.1016/0030-4018(92)90140-m.
  * Schwinger [1951] J. Schwinger, Physical Review 82, 664 (1951), URL https://doi.org/10.1103/physrev.82.664.
  * Zwillinger [2014] D. Zwillinger, _Table of Integrals, Series, and Products_ (Academic Press, 2014), eighth ed.
  * Davis [1976] J. P. Davis, The Journal of Chemical Physics 64, 3129 (1976), URL https://doi.org/10.1063/1.432648.
  * Dykhne [1962] A. Dykhne, JTEP 41, 1324 (1962).
  * Kazantsev et al. [1990] A. P. Kazantsev, G. I. Surdutovich, and V. P. Yakovlev, _Mechanical Action of Light on Atoms_ (WORLD SCIENTIFIC, 1990), URL https://doi.org/10.1142/0585.
  * Gu et al. [2019] F.-L. Gu, J. Liu, F. Mei, S. Jia, D.-W. Zhang, and Z.-Y. Xue, npj Quantum Information 5, 36 (2019), URL https://doi.org/10.1038/s41534-019-0148-9.
  * Tarruell et al. [2012] L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, Nature 483, 7389 (2012), URL https://doi.org/10.1038/nature10871.
  * Song et al. [2019] B. Song, C. He, S. Niu, L. Zhang, Z. Ren, X.-J. Liu, and G.-B. Jo, Nature Physics 15, 911 (2019), URL https://doi.org/10.1038/s41567-019-0564-y.
  * BenDahan et al. [1996] M. BenDahan, E. Peik, J. Reichel, Y. Castin, and C. Salomon, Physical Review Letters 76, 4508 (1996), URL https://doi.org/10.1103/physrevlett.76.4508.
  * Novak et al. [2019] M. Novak, S. N. Zhang, F. Orbanić, N. Biliškov, G. Eguchi, S. Paschen, A. Kimura, X. X. Wang, T. Osada, K. Uchida, et al., Physical Review B 100, 085137 (2019), URL https://doi.org/10.1103/physrevb.100.085137.


# First results from the FPGA/NIOS Adaptive FIR Filter Using Linear Prediction
Implemented in the AERA Radio Stations to Reduce Narrow Band RFI for Radio
Detection of Cosmic Rays

Zbigniew Szadkowski, , D. Głas, C. Timmermans, T. Wijnen,
for the Pierre Auger Collaboration Manuscript received June 03, 2014.Zbigniew
Szadkowski is with the University of Łódź, Department of Physics and Applied
Informatics, Faculty of High-Energy Astrophysics, 90-236 Łódź, Poland,
(e-mail: zszadkow @kfd2.phys.uni.lodz.pl, phone: +48 42 635 56 59).Dariusz
Głas is the PhD student with the University of Łódź, Department of Physics and
Applied Informatics, Faculty of High-Energy Astrophysics, 90-236 Łódź,
Poland,Charles Timmermans, Thei Wijnen are with Nikhef, Science Park,
Amsterdam, Netherlands, and IMAPP, Radboud University Nijmegen, Netherlands

###### Abstract

The FPGA/NIOS® FIR filter based on linear prediction (LP) to suppress radio
frequency interference (RFI) has been installed in several radio stations in
the Auger Engineering Radio Array (AERA) experiment. AERA observes coherent
radio emission from extensive air showers induced by ultra-high-energy cosmic
rays to make a detailed study of the development of the electromagnetic part
of air showers. Radio signals provide complementary information to that
obtained from Auger surface detectors, which are predominantly sensitive to
the particle content of an air shower at the surface. The radio signals from
air showers are caused by the coherent emission due to geomagnetic and charge-
excess processes. These emissions can be observed in the frequency band
between 10 - 100 MHz. However, this frequency range is significantly
contaminated by narrow-band RFI and other human-made distortions. A FIR filter
implemented in the FPGA logic segment of the front-end electronics of a radio
sensor significantly improves the signal-to-noise ratio. Theoretical
calculations show a high efficiency of this filter for mono-carrier as well as
for standard FM radio contaminations. The laboratory tests, performed on the
Altera® Cyclone® V DK-DEV-5CEA7N development kit confirmed the theoretical
expectations.

In this paper we present first results of the efficiency of the adaptive LP
FIR filter, deployed in real AERA station on pampas, with a comparison to the
currently used IIR notch filter with constant coefficients. The coefficients
for the linear predictor are dynamically refreshed and calculated in a Voipac
PXA270M ARM processor, which is implemented on a daughter-board placed in the
same digital unit as the FPGA. The laboratory tests confirms the stability of
the filter. Using constant LP coefficients the suppression efficiency remains
the same for hours, which corresponds to more than $\bf 10^{12}$ clock cycles.
We compared in real conditions several variants of the LP FIR filter with
various lengths and various coefficients widths (due to fixed-point
representations in the FPGA logic) with the aim to minimize the power
consumption for the radio station while keeping sufficient accuracy for noise
reduction.

The laboratory and real condition tests provide data to optimize the RFI
cleaning for the next generation of the AERA Front-End based on Cyclone® V
with Hardcore Processor System (HPS) and System on Chip (SoC).

##  I Introduction

Recently radio detection of cosmic-ray air showers relives a renaissance,
mainly thanks to a huge progress of the powerful digital signal processing
techniques in experiments such as LOPES [1], CODALEMA [2] or the Auger
Engineering Radio Array (AERA) [3], which is situated within the Pierre Auger
Observatory (PAO) [4].

Figure 1: A diagram showing a (FFT + Median filter + iFFT) chain cleaning the
signal from the RFI contamination. The 1st graph shows the ADC input as
unsigned data with an offset of ca. 2300 ADC-counts, the 2nd - the absolute
values of FFT coefficients in the frequency domain, the 3rd - FFT coefficients
“decontamined” by the median filter and 4th - signal converted back to the
time domain. Additionally, the 0th FFT coefficient has been zeroed. Thus, the
cleaned signal in the time domain is represented as signed data without the
offset. The amplitude of the signal remains roughly the same and the noise is
considerably reduced.

Results from the Pierre Auger Observatory, point to the need for very large
aperture detection systems for ultra-high energy cosmic rays. With its nearly
100% duty cycle, its high angular resolution, and its sensitivity to the
longitudinal air-shower evolution, the radio technique is particularly well-
suited for detection of ultra-high energy cosmic rays (UHECRs) in large-scale
arrays. AERA has been enlarged to 124 radio detector stations (RDSs), covering
an area of 6.5 km2 therefore allowing the detection of UHECRs.

Since the 1960s we know that the radio emission from air showers is strongly
correlated with the local geomagnetic field [5]. In addition to the
geomagnetic effect that can be described macroscopically [6], Askaryan [7] [8]
predicted that there should be an emission component related to the time-
variation of the negative net charge excess in air showers. Both [1] and [2]
confirmed these effects. Using the AERA setup, the Auger collaboration
quantified the relative strength of these effects [9] The observation of air
showers with radio-detection techniques can be done at almost all times.
Moreover, radio signals are sensitive to the development of the
electromagnetic component of particle showers in the atmosphere of the Earth.
In the last 10 years the radio-detection technique in the MHz region has been
revived and the present radio-detector arrays for cosmic-ray research are
equipped with low-noise and high-rate digital samplers. Simultaneously, the
number of stations within these arrays has grown from less than ten to more
than one thousand. The question to be addressed in the VHF band (MHz-range) is
not whether extensive air showers emit radiation. At this moment the main
question to answer is: _can we use radio signals to determine the primary
energy, the arrival direction, and the mass of cosmic rays with accuracies
which are equal to or better than those obtained by other techniques? And if
yes, can we build for an affordable price a huge surface-detector array based
on the radio-detection technique?_

##  II RFI suppression for real AERA data

Triggering directly on the radio signal of the air showers (instead of using
particle detectors as a trigger) poses some challenges for the data
acquisition, due to man-made radio-frequency interference (RFI). The
continuous background level is set by the radio emission from the Galactic
plane, but any man-made narrow band transmitters add to the level above which
one must detect air-shower pulses. Additionally, man-made pulsed RFI (from
sparking electrical equipment, airplanes, etc.) can mimic the short pulsed
signal from cosmic rays. Since the bandwidth and computational resources at
each triggering level are limited, one of the technical focuses for the first
stage of the array has been to develop various methods to reject RFI in order
to minimize efficiency losses from bandwidth saturation.

The energy threshold of radio detection of cosmic rays is limited by the
considerable radio background and noise. The very high level of RFI in the FM
and short wave band has to be eliminated by a band pass filter. Within the
remaining receiver bandwidth of 30 to 80 MHz the noise at the quiet-rural
environment of the Pierre Auger Observatory is dominated by the frequency
dependent galactic noise [10] with noise temperatures of 5000 K at 60 MHz

In addition to the galactic noise, there is a human made background. This
background consists of continuous signals, as from a few radio and TV
stations, and transients produced by machines. Without an effective trigger, a
stable and low level energy threshold is not guaranteed. Furthermore, the data
rate for communication of the triggered data to the central DAQ would exceed
the available bandwidth.

###  II-A FFT + Median filter + iFFT

For self-triggered measurements, the data will be digitized and processed in
real time by a powerful FPGA chip. The narrow peaks in the frequency domain
due to radio frequency interferences have to be strongly suppressed before
building a trigger. These peaks are removed in several stations, digitizing at
180 MHz, using a median filter. The filter works in the frequency domain using
the Fast Fourier Transform (FFT) routine provided by Altera . Furthermore, the
phase of the signal deformed by the steep band pass filter is reconstructed by
a deconvolution in the frequency domain. The median FPGA filter eliminates
mono-frequent carriers, but broadband radio pulses from cosmic showers are not
affected. After a second inverse FFT, signals are converted back to the time
domain. This chain of the digital signal processing strongly enhances the
signal to noise ratio, and thus improves the radio pulse detection sensitivity
(Fig. 1).

In order to suppress the strong man-made radio carrier signals below 30 MHz
and above 80 MHz and to fulfill the Nyquist theorem, the signals of the
antennas are filtered using an analogue band filter before being sampled by
the ADCs. The filters with the constant transfer function $H(f)$ are of high
order with a non-constant group delay, thus leading to a dispersion and
decreasing the amplitude of the input signals $S(f)$. In a frequency domain,
the resulting output signal $P(f)$ can simply be calculated using (1) :

$P(f)=H(f)\cdot S(f)$ (1)

Knowing $H(f)$ from measurements, it is now easy to get back $S(f)$ by
inverting (1). This operation can be implemented into the FPGA by placing a
complex multiplication unit directly after the FFT engine which multiplies the
output data of the engine with the precomputed coefficients of $1/H(f)$ stored
in a RAM. After transforming the signal back into the time domain, the
amplitude of the de-convoluted signal increases by about 20% compared to the
input signal. Since the galactic and electronic noise is completely
uncorrelated, the S/N ratio increases by the same amount.

Aliasing appears when converted pulses are located close to a border of
converted blocks. It manifests by a spurious contribution in the opposite
border of the block and in the neighboring block as well. This effect may
cause spurious triggers and has to be eliminated.

The problems can only be solved, without introducing dead time between the
blocks, by using an overlapping routine. Therefore the filter engine must run
in another clock domain with higher frequency (210 MHz instead of 180MHz).
Nevertheless, it requires additional resources e.g. FILO (First In Last Out)
memory based procedure to inverse in time a sequence of samples. The FFT
approach is generally very power consuming. This is a factor for a system
supplied from solar panels.

###  II-B IIR notch filter

As stated earlier, before triggering on a radio pulse it is advantageous to
increase the signal-to-noise ratio by filtering out any narrow band
transmitters from the digitized antenna signals. In stations digitizing at 200
MHz, this is accomplished in a computationally efficient manner by using a
series of infinite-impulse response (IIR) notch filters in a Cyclone®-IV FPGA.
The IIR filters operate on the time-domain signal, and the output of the
filter $y_{i}$ is a linear combination of input samples $x_{j}$ and delayed
feedback output samples $y_{j}$ from the filter:

$\displaystyle y_{i}=x_{i}$ $\displaystyle-$
$\displaystyle(2\cos\omega_{N}\cdot x_{i-1})+x_{i-2}+$ (2) $\displaystyle+$
$\displaystyle(2r\cos\omega_{N}\cdot y_{i-1})-(r^{2}\cdot y_{i-2})$

The normalized filter frequency $\omega_{N}=2\pi f_{N}/f_{S}$ is given by the
notch frequency $f_{N}$ and the sampling frequency $f_{S}$. The width
parameter, r, is provided with a value strictly between 0 and 1, with higher
values giving a narrower response function. For a narrow transmitter a typical
value is r=0.99. A complication in the implementation of the IIR filters in
high-frequency FPGAs arises from the fact that one cannot arbitrarily pipeline
the feedback computation. We have resolved this by using the scattered look-
ahead pipelining technique [11], which increases the filter complexity but
allows more time for computation in the FPGA.

Figure 2: The data flow of the FIR filter based on the LP method

The coefficients of the $x_{j}$ and $y_{j}$ in Eq. 2 can be precomputed for
the desired notch frequencies, converted to a fixed-point representation and
loaded into the FPGA at run-time. The current design allows for four
independent tunable notch filters for each polarization direction and is
applied on a Cyclone® III, IV and V FPGA. It has been operating in the field
for several years and successfully decreased and stabilized the threshold
settings.

$\uparrow i-p-D$$\uparrow i-D-1$$\uparrow i$…..

Figure 3:  An illustration of the method. The sine wave represents the signal
that is ‘fitted’ (although in actuality no sinusoidal fit is performed like
this) where sample number $i$ is predicted by using the samples $i-p-D$ to
$i-D$. The predicted values are then subtracted from the original values as
illustrated by the green dots below the horizontal bar with the minus sign,
reducing the variance of the signal.

###  II-C Linear predictor

Linear prediction is a mathematical operation where future values of a
discrete-time signal are estimated as a linear function of previous samples
[12]. This method is widely used in audio signal processing and speech
processing for representing the spectral envelope of a digital signal of
speech in compressed form, using the information of a linear predictive model
[13]. With the advent of faster signal processing techniques in FPGAs it is
now possible to apply similar techniques to the real-time processing of radio
signals in the 10 - 100 MHz region [14].

In the LP method the covariances for 1024 ADC samples can be calculated in the
FPGA fast logic block. Either the NIOS® processor, or the external ARM-
processor, solves the matrix of 32 or 64 linear equations and provides
coefficients needed for the FIR filter. The calculated coefficients are next
transferred to the fast logic block, updating appropriate registers. They are
used as the FIR coefficients in the ADC data filtering. Finally, the predicted
and delayed data (expected background) are subtracted from the ADC data to
clean the signal from periodic contaminations (Fig. 2).

Comparison of graphs in [16] indicates that the LP approach can eliminate RFI
narrow frequency contaminations as efficient as notch and FFT filters.

Figure 4: FFT for ADC and FIR filter with various D = 128 and 1, respectively,
for long-term (7.5 s) data. Time 7.5 s corresponds to $1.5\cdot 10^{9}$ cycles
(sampling = 200 MSps). Registered events contain 1024 samples. The radio
station LS009NS has been selected for a strong (100000 a.u. - arbitrary units)
contamination by 4 carriers. The FIR64 filter shows a very good efficiency of
the RFI suppression for all D factors (128 and 1).

Figure 5: Here, the radio station LS122EW has been selected for a weak (6000
a.u. - arbitrary units) contamination, focused on two narrow bands around 12
and 62 MHz. The suppression for 62 MHz band is pretty good for all D factors,
however for the 12 MHz band, the RFI suppression is most efficient for D=128.

Figure 6: Example of a potential signal distortion for D=1. It is recommended
to introduce a delay for a signal modification. For D=128 signals are almost
unaffected.  Figure 7: Histogram of distortion factors (DF) for several
hundreds AERA event filtered by the FIR filter with D=1, 32 and 128,
respectively. It is visible that a configuration with D=1 distorts the signal
too strong and cannot be used in a final design.

##  III Analysis of AERA events

Figures 4 \- 5 show spectra of original ADC traces and cleaned by the LP FIR
filter with various values for parameter D of 128, 32 and 1, respectively
(Fig. 3). For the mono-carrier contamination (Fig. 4) the RFI suppression is
very good for any D parameter. For more structured spectral contaminations the
suppression depends on the D parameter , nevertheless, we observe a reduction
of periodic noise.

The delay-line $D$ implies that there is a gap between the samples that are
used for the prediction and the sample that is to be predicted (Fig. 3). This
delay-line is necessary to allow transient signals to pass through the filter
unaltered. For D=1 the efficiency of the RFI suppression is maximal, however,
the signal is significantly affected and a distortion factor (DF) (Fig. 6),

$DF=\sum\limits_{k=-16}^{16}\left(1-\frac{(x_{FIR})_{k}}{(x_{ADC})_{k}}\right)^{2}$
(3)

introduced to estimate a quality of filtering, reaches a large, unacceptable
value. For laboratory tests we selected D=128, to keep a reasonable safety
margin in measurements (Fig. 7).

Figure 8: Plots showing a suppression of the sine signal with a sweep from
50.0 MHz to 50.2 MHz in 2 minutes. For 8 such a cycles we observe a perfect
repetition. Plots obtained in real measurements on the laboratory setup for
FIR64, D=128 and 14-bit LP coefficients.

##  IV Laboratory tests

The LP FIR filter has been tested with the Altera® DK-DEV-5CEA7N Development
Kit with a Cyclone® V FPGA and Texas Instr. ADS4249EVM Evaluation Module with
2-channel 14-bits 250MSps ADC (ADS4249). Both modules were connected through
the Altera® HSMC-ADC-BRIDGE providing the LVDS data transmission (Fig. 17).

Figure 9: Plots showing a suppression of the 50 MHz signal frequency modulated
(FM) with a deviation of 75 kHz and a modulation of 15 kHz. The signal
corresponds to Hi-Fi radio transmission. Even for wide-band FM the RFI
suppression is relevant. Plots obtained in real measurements on the laboratory
setup for FIR64, D=128 and 14-bit LP coefficients.

At first the filter was tested by a mono-carrier drifting the frequency from
50.0 MHz to 50.2 MHz in 120 s. The LP coefficients were not calculated by
NIOS® in a correlation with a generator running. Thus, these LP coefficients
were used for the data cleaning in several tens of minutes (which corresponds
to more than $10^{12}$ clock cycles). Fig. 8 shows a perfect long-term
stability of the filter. If the generator driving the filter uses the
frequency for which the LP coefficients were originally calculated, the
suppression is almost total. Therefore, LP coefficients do not need to be
refreshed very frequently.

Figure 10: Suppression factors for two sine pure carriers with frequencies
27.12 and 57.9 MHz and $V_{p-p}$ 50 and 200 mV, respectively, mixed with 10,
50, 100, 200 and 500 mV noise.

Figure 11: Suppression factors for two sine pure carriers with frequencies
27.12 and 57.9 MHz and $V_{p-p}$ 50 and 200 mV, respectively, mixed with 10,
50, 100, 200 and 500 mV noise. Bottom graph shows real time traces for 200 mV
noise. A suppression in the time domain is almost invisible

Secondly, we checked the Hi-Fi FM configuration (75 kHz deviation of the 50
MHz carrier with the maximal acoustic 15 kHz modulation). This is very
restrictive condition, which actually should not appear in real conditions in
Argentina. The band 30 - 80 MHz is used rather by narrow-band transmitter,
while the FM Hi-Fi transmission is selected for the band of 88 - 108 MHz, cut-
off by the band-pass analog filter. Nevertheless, even these critical
contaminations were successfully suppressed (Fig. 9).

We also tested a suppression efficiency (in the Fourier space) in the presence
of relatively strong white noise (Fig. 11). Two pure carriers with 27.12 and
57.9 MHz were wired mixed with noise. Fig. 10 shows that even when the noise
level reaches the signal level the suppression factor remains on a level of
5-10 (for signals 50 mV and 200 mV contaminated by 50 mV and 200 mV of a
noise). When the noise is small, the suppression factor reaches values up to
35. However, a strong asymmetry is observed. The strong signal (i.e. 200 mV of
27.12 MHz) is suppressed with a very high factor (35) while a 4 times smaller
signal (57.9 MHz) is almost not suppressed at all. For vice versa
configuration a suppression structure remain the same.

##  V Comparison to currently used IIR-notch filter

We compared suppression characteristics of the FIR filter based on the linear
prediction with the currently used IIR-notch filter with 4 band-reject bands.
Fig. 12 shows amplitudes of the $138^{th}$ frequency bin in the Fourier space
corresponding to a contribution of $1^{st}$ carrier contamination with a
various frequencies in a range of 27.00 - 27.24 MHz for a small (10 mV - upper
graph) and significant (100 mV - lower graph) noise. It is well visible that
the efficiency of the IIR filter is very high and, as expected, only around
the reject frequency (27.12 MHz) of the filter. Higher noise (comparison 100
mV vs. 10 mV) reduces the efficiency only slightly, thereby extending the
width of the rejection band (from $\sim\pm$20 kHz to $\sim\pm$30 kHz).

Figure 12: Amplitudes of Fourier modules ($\bar{X}_{freq_{1}=27.12MHz}$) for
FIR filter based on the linear predictor and IIR-notch filter. The width of
frequency bin is $\sim$0.1 MHz (1024-point FFT), all $\bar{X}_{freq_{1}}$
contributions fall into $138^{th}$ frequency bin. The signal with various
frequencies $freq_{1}$ = 27.00 - 27.24 with a grid of 10 kHz is provided by
the $1^{st}$ channel of Tektronic AFG3252C, the $2^{nd}$ signal from the
$2^{nd}$ channel of AFG3252C has a fixed frequency = 55.2 MHz, exactly equals
the reject frequency of he $3^{rd}$ stage of the IIR-notch filter. The reject
frequency of the $1^{st}$ IIR-not filter is 27.12 MHz. Additional 10 mV or 100
mV noise is provided by Agilent 33250A generator.  Figure 13: Amplitudes of
Fourier modules ($\bar{X}_{freq_{2}=55.2/57.9MHz}$) for the $2^{nd}$ mono-
carrier contamination with 55.2 or 55.9 MHz, respectively.  Figure 14: The
current consumption for several tested variants of the FIR filter with a
comparison to the currently used IIR one.  Figure 15: A graph showing the
current consumption of the core (@1.1 V). A reduction of the current
consumption corresponds to temporary disabling of the NIOS processor

Figure 16:  FFT for ADC and FIR filter with D = 128, for long-term (7.5 s)
data for the radio station LS009 and NS and EW polarizations. As for Fig. 4 we
see 4 mono-carriers, however with additional low frequency non-stationary RFI.
The presence of the variable in time RFI reduces an efficiency of the LP
filter.

Fig. 13 shows the $283^{th}$ (55.2 MHz) and $296^{th}$ (57.9 MHz) Fourier
modules vs. various frequencies in the $1^{st}$ channel contamination (27.00 -
27.24 MHz). As expected, the suppression of the IIR filter outside the reject
band is negligible. Note that for structured contaminations (two mono-carriers
with frequencies equal exactly the reject frequencies of the IIR-filter) the
efficiency of the FIR filter is comparable with the efficiency of the IIR one,
provided the background noise level is low in comparison to the contamination.

Figures 12 and 13 justify testing the FIR filter based on the linear predictor
as an adaptive filter, adjusting the suppression characteristics to changing
RFI conditions in the field.

##  VI Power consumption

A power consumption is the factor for systems supplied from solar panels. More
sophisticated filter provides a better accuracy of data processing. However,
the power efficiency may significantly decreased. We measured a power
consumption for all developed FIR variants and compared with the currently
used IIR filter. Fig.14 shows results for the current consumption with the
NIOS processor (left panel) and with temporary blocked NIOS (right panel). The
NIOS processor was used for the calculation of LP coefficients for the FIR
filter and for data transmission via UART to the PC for both filters Fig.15.

From Fig.14 it is visible that the FIR filter wit 32 stages and only 14-bit
coefficients is at least power efficient in comparable to the IIR one. Longer
FIR filters consume much more power and could be used exceptionally in
extremely contaminated environment, where the suppression efficiency becomes
the more important factor than the power efficiency.

##  VII Data from pampas

We have implemented the 32-stage LP filter (with D=128) in the AERA radio
station LS009. Fig. 16 show a suppression of the RFI in the higher frequency
range. In comparison to Fig. 4b the suppression efficiency is lower, however,
in current environment the huge contribution of non-stationary RFI was
observed. In a presence of variable RFI the efficiency of the LP filter is not
so high as for mono-carriers contamination. Fig. 16 shows that the noise
suppression for the EW polarization is almost negligible.

Figure 17: The laboratory measurement setup with the Altera® DK-DEV-5CEA7N
Development Kit with Cyclone® V FPGA and Texas Instr. ADS4249EVM Evaluation
Module with 2-channel 14-bits 250MSps ADC (ADS4249). Both modules were
connected by the Altera® HSMC-ADC-BRIDGE providing the LVDS data transmission.
Tektronix AFG3252C 2-channel generator provides 2 sine waveforms, Agilent
33250A provides noise signal.

##  VIII Conclusions

Analytical calculations (Fig. 4-5) show that a very high efficiency of the LP
filter can be obtained. The laboratory measurements confirmed that when the
data is contaminated by mono-carriers, the suppression factor is very high.
However, the environmental RFI as encountered in the Pierre Auger Observatory
has a much more sophisticated structure. Therefore, the efficiency of the LP-
filter significantly decreases.

Nevertheless, in each radio station the filter setup can be calculated,
thereby optimizing the suppression factor, depending on the location of each
station and the type of RFI contamination.

In AERA the Cyclone® IV FPGA EP4CE75F29I7 is currently in use. The laboratory
tests provide input for an optimization of the RFI cleaning for the next
generation of the AERA Front-End based on Cyclone® V with Hardcore Processor
System (HPS) and System on Chip (SoC).

## Acknowledgment

This work was supported by the Polish National Center for Research and
Development under NCBiR Grant No. ERA/NET/ASPERA/02/11, by the National
Science Centre (Poland) under NCN Grant No. 2013/08/M/ST9/00322, , and by the
Ministerie van Onderwijs, Cultuur en Wetenschap, Nederlandse Organisatie voor
Wetenschappelijk Onderzoek (NWO), Stichting voor Fundamenteel Onderzoek der
Materie (FOM), the Netherlands.

## References

  * [1] H. Falcke,W. D. Apel, A. F. Badea, et al., ”Detection and imaging of atmospheric radio flashes from cosmic ray air showers”, _Nature_ , vol. 435, pp. 313-316, May 2005.
  * [2] D. Ardouin, A. Bell´etoile, D. Charrier, et al., ”Radioelectric field features of extensive air showers observed with CODALEMA” _Astropart. Phys._ , vol. 26, pp. 341-350, Dec. 2006.
  * [3] S. Fliescher for the Pierre Auger Collaboration, ”Radio detection of cosmic ray induced air showers at the Pierre Auger Observatory”, _Nucl. Instr. Meth._ , ser. A, vol. 662, pp. 124-129, Jan. 2012.
  * [4] J. Abraham et al., [Pierre Auger Collaboration], “Properties and Performance of the Prototype Instrument for the Pierre Auger Observatory”, _Nucl. Instr. Meth._ , ser. A, vol. 523, pp. 50-95, May 2004.
  * [5] D. J. Fegan, ”Detection of elusive radio and optical emission from cosmic-ray showers in the 1960s”, _Nucl. Instr. Meth._ , ser. A, vol. 662, pp. 2-11, Jan. 2012.
  * [6] F. D. Kahn and I. Lerche, ”Radiation from Cosmic Ray Air Showers”, _Proc. Roy. Soc. A_ , vol. 289, pp. 206-213, Jan. 1966.
  * [7] G. A. Askaryan, _Journal of Exp. and Theoretical Phys._ , vol. 14, pp. 441, 1962.
  * [8] G. A. Askaryan, ”Coherent Radio Emission from Cosmic Showers in Air and in Dense Media”, _Journal of Exp. and Theoretical Phys._ , vol. 21, pp. 658-659, Jan. 1965.
  * [9] A. Aab et al., ”Probing the radio emission from air showers with polarization measurements”, Phys. Rev. D89, 052002, Mar. 2014.
  * [10] G. A. Dulk, W. C. Erickson, R. Manning, and J.-L. Bougeret, ”Calibration of low-frequency radio telescopes using the galactic background radiation”, _A &A_, vol. 365, pp. 294-300, Jan. 2001
  * [11] K. K. Parhi and D. G. Messerschmitt, _IEEE Trans. on Acoustics, Speech, and Signal Processing_ , vol. 37, pp. 1099-1117, 1989.
  * [12] J. Makhoul, “Linear prediction: A tutorial review” _Proc. of the IEEE_ , vol. 63, no. 4, pp. 561-580, Apr. 1975.
  * [13] D. Li, D. O’Shaughnessy, ”Speech processing: a dynamic and optimization-oriented approach”, pp. 41–48. ISBN 0-8247-4040-8 (2003)
  * [14] Z. Szadkowski, E.D. Fraenkel, A. M. van den Berg, ”FPGA/NIOS Implementation of an Adaptive FIR Filter Using Linear Prediction to Reduce Narrow-Band RFI for Radio Detection of Cosmic Rays”, _IEEE Trans. on Nucl. Science_ , vol. 60, pp. 3483-3490, Oct. 2013.
  * [15] Z. Szadkowski, E.D. Fraenkel, D. Głas, R. Legumina, ”An optimization of the FPGA/NIOS adaptive FIR filter using linear prediction to reduce narrow band RFI for the next generation ground-based ultra-high energy cosmic-ray experiment”, _Nucl. Instr. Meth._ , ser. A, vol. 732, pp. 535-539, June 2013.
  * [16] Z. Szadkowski, Ad M. van den Berg, E.D. Fraenkel, D. Głas, J. Kelley, C. Timmermans, T. Wijnen, ”Analysis of the efficiency of the filters suppressing the RFI being developed for the extension of AERA”, _33nd International Cosmic Ray Conference - July 2013 - Rio de Janeiro, Brazil_.