# Giant microwave absorption in the vortex lattice in $s$-wave superconductors

T. Liu  Department of Physics, University of Washington, Seattle, WA 98195,
USA  M. Smith  Materials Science Division, Argonne National Laboratory,
Lemont, Illinois 60439, USA  A. V. Andreev  Department of Physics, University
of Washington, Seattle, WA 98195, USA  B. Z. Spivak  Department of Physics,
University of Washington, Seattle, WA 98195, USA

(August 28, 2024)

###### Abstract

In this article we study microwave absorption in superconductors in the
presence of a vortex lattice. We show that in addition to the conventional
absorption mechanism associated with the vortex core motion, there is another
mechanism of microwave absorption, which is caused by the time-dependence of
the quasiparticle density of states outside the vortex cores. This mechanism
exists even in the absence of vortex motion and provides the dominant
contribution to microwave absorption in a broad interval of physical
parameters. At low frequencies, the dissipative part of the microwave
conductivity $\sigma(\omega)$ is proportional to the inelastic relaxation
time, $\tau_{\mathrm{in}}$, which is typically much larger than the elastic
relaxation time, $\tau_{\mathrm{el}}$. At high frequencies $\sigma(\omega)$ is
proportional to the quasiparticle diffusion time across the inter-vortex
distance, $\tau_{\mathrm{D}}$, which is still larger than
$\tau_{\mathrm{el}}$.

##  I Introduction

Type-II superconductors placed in a magnetic field normal to the plane, which
is weaker than the upper critical field $H_{c2}$, host Abrikosov vortices
whose density is set by the condition that the average flux of the magnetic
field per vortex is equal to the flux quantum $\Phi_{0}=\frac{\pi\hbar c}{e}$
[1]. In the absence of pinning, the Magnus force on the vortices induced by
the transport current through the sample causes vortex motion and dissipation.
The corresponding dissipative conductivity has been extensively studied since
the work of Bardeen and Stephen [2, 3, 4, 5, 6] and is described by the
formula,

$\displaystyle\sigma_{BS}\sim$ $\displaystyle\sigma_{n}\frac{H_{c2}}{H},$ (1)

where $\sigma_{n}$ is the normal state conductivity and $H_{c2}$ is the upper
critical field. In the Bardeen-Stephen theory, dissipation arises from the
friction force caused by vortex motion. Equation (1) may be obtained by
expressing the friction force per unit line length of the vortex as
$\bm{F}=-\eta_{BS}{\bm{v}}_{v}$, where
$\eta_{BS}=\Phi_{0}H_{c2}\sigma_{n}/c^{2}$ is the vortex viscosity, and
${\bm{v}}_{v}$ is the vortex velocity. In the flux flow regime the latter is
given by

$\displaystyle{\bm{v}}_{v}=$
$\displaystyle\,c\,\frac{{\bm{E}}\times{\bm{H}}}{H^{2}}.$ (2)

Equating the rate of viscous energy dissipation to Joule heat,
$\frac{H}{\Phi_{0}}\eta_{BS}{\bm{v}}_{v}^{2}=\sigma_{BS}E^{2}$, one obtains
Eq. (1).

In the presence of disorder the vortex lattice is pinned, and is capable of
supporting a dissipationless current density which is smaller than the
critical current, $J_{c}$. Therefore, Eq. (1) is relevant only in the non-
linear (flux flow) regime where the transport current significantly exceeds
$J_{c}$. On the other hand, _ac_ -electric fields
${\bm{E}}={\bm{E}}_{0}\cos(\omega t)$ induce dissipation in superconductors
even in the linear regime. In this case, the microwave absorption coefficient
is controlled by the dissipative _ac_ -conductivity $\sigma(\omega)$. The
latter can be evaluated using the relation

$T\overline{\dot{S}}=\frac{1}{2}\sigma(\omega){\bm{E}}_{0}^{2},$ (3)

where $T$ is the temperature, $\dot{S}$ is the entropy production rate per
unit volume, and the overline $\overline{(...)}$ indicates averaging over
time. In most articles on microwave absorption, the dissipative conductivity
in the linear regime is evaluated phenomenologically using the Bardeen-Stephen
expression for the vortex viscosity, $\eta_{BS}$, and the vortex velocity
${\bm{v}}^{\prime}_{v}(t)$, which is modified by the pinning forces (see, for
example, [7, 8, 9]). The corresponding result, denoted by
$\sigma^{\prime}_{BS}(\omega)$ below, is proportional to the elastic momentum
relaxation time $\tau_{el}$ in the normal state.

In this article we describe a new mechanism of microwave absorption in
superconductors in the mixed state. This mechanism is caused by the spectral
flow of the quasiparticle energy levels in the presence of an _ac_ -electric
filed, and exists even in the absence of vortex motion. We show that in a
broad range of physical parameters the dissipative part of the conductivity
caused by this mechanism can be parametrically larger than
$\sigma^{\prime}_{BS}(\omega)$.

Figure 1: The black line denotes the dependence of the modulus of the order
parameter $\Delta(r)$ on the distance from the vortex core $r$. The blue line
shows the $r$-dependence of the edge of the quasiparticle spectrum.

The origin of this mechanism can be traced to a fundamental difference between
the quasiparticle kinetics in superconductors and electron kinetics in normal
metals. In superconductors, the density of quasiparticle states
$\nu(\epsilon,t)$ may depend on time, which means that energies of individual
single particle states change in time. In the adiabatic approximation, the
quasiparticles occupying these levels also travel in energy space. As a
result, a non-equilibrium quasiparticle distribution is created, whose
relaxation generates entropy and provides a mechanism for microwave
absorption.

The time-dependence of $\nu(\epsilon,t)$ arises because the microwave electric
field induces oscillations of the superfluid momentum

$\displaystyle\bm{p}_{s}(\bm{r},t)=$
$\displaystyle\,\frac{\hbar}{2}\left[\bm{\nabla}\chi(\bm{r},t)-\frac{2e}{\hbar
c}\bm{A}(\bm{r},t)\right].$ (4)

Here $\chi(\bm{r},t)$ is the phase of the order parameter, $\bm{A}(\bm{r},t)$
\- the vector potential, and $e$, $\hbar$, and $c$ are, respectively, the
electron charge, Planck’s constant, and the speed of light. The microwave
field induces acceleration of the condensate,
$\dot{\bm{p}}_{s}(\bm{r},t)\propto e{\bm{E}}(t)$. Since the density of states
$\nu(\epsilon,t)$ is a scalar, its time derivative can have a linear
dependence on $\dot{\bm{p}}_{s}(\bm{r})$ only in the presence of a _dc_
-superfluid momentum $\bar{\bm{p}}_{s}({\bm{r}})$;
$\dot{\nu}(\bm{r})\propto\dot{\bm{p}}_{s}(\bm{r})\cdot\bar{\bm{p}}_{s}(\bm{r})$.
In the case of a flux lattice, the _dc_ -superfluid momentum is the
equilibrium superfluid momentum about the vortex cores.

The non-equilibrium quasiparticle distribution created by the spectral flow
can relax via two channels: inelastic scattering and quasiparticle diffusion.
The relative importance of these two channels depends on frequency of the
microwave field $\omega$. In typical superconductors the inelastic relaxation
time, $\tau_{\mathrm{in}}$, significantly exceeds the diffusion time across
the vortex plaquette, $\tau_{\mathrm{D}}$. Therefore inelastic relaxation is
caused by those quasiparticles for which diffusion cannot lead to full
equilibration. These quasiparticles have energies below the threshold of
percolation between different vortex plaquettes. The percolation threshold
$\epsilon^{*}$ arises because in the presence of a superfluid momentum
$\bar{\bm{p}}_{s}(\bm{r})$ the energy gap of the quasiparticle states is
shifted down from $\Delta$ by $\delta\epsilon(\bar{p}_{s}(\bm{r}))$, as
illustrated in Fig. 1. The non-equilibrium quasiparticles with energies below
$\epsilon^{*}$ are trapped inside the vortex plaquettes. Because of pinning
the vortex lattice is distorted, and the spatial distribution of the
superfluid momentum $\bar{\bm{p}}_{s}(\bm{r})$ is not symmetric about the
vortex core. Therefore, the non-equilibrium part of the of the distribution
function cannot relax completely by diffusion across the plaquettes, and
ultimately its relaxation is achieved by inelastic scattering processes. Thus,
the low frequency conductivity $\sigma\sim K\tau_{\mathrm{in}}$ is
proportional to the inelastic relaxation time and a parameter $K$
characterizing the degree of lattice distortion. Since $\tau_{in}$ may exceed
the elastic relaxation time $\tau_{\mathrm{el}}$ by many orders of magnitude,
this contribution typically is much larger than $\sigma^{\prime}_{BS}$. 111The
mentioned above mechanism of absorption has been discussed in the context of
microwave absorption in superconductors in the presence of spatially uniform
supercurrent [10, 11] and the resistance of superconductors in flux flow
regime [12].

At higher frequencies the value of the conductivity is controlled by the
quasiparticle diffusion time across a vortex lattice plaquette, $\tau_{D}$,
which is assumed to be much smaller that the inelastic relaxation time but
larger than the elastic relaxation time
$\tau_{\mathrm{el}}\ll\tau_{\mathrm{D}}\ll\tau_{\mathrm{in}}$. In this case
the value of the conductivity $\sigma$ is still larger than
$\sigma^{\prime}_{BS}$.

For simplicity, we focus on the case where the thickness $d$ of the
superconducting film is smaller than the skin length, and the microwave
electric field ${\bm{E}}(t)$ is spatially uniform in the film. We also assume
that the distance between vortices, which is of order of the magnetic length
$l_{H}=\hbar c/|e|H$, is smaller than the Pearl length [19], so that the
magnetic field is also uniform in the film.

The consideration below is organized as follows. In Sec. II we present a
general kinetic theory of quasiparticle dynamics in superconductors with time-
dependent density of states. In Sec. III we discuss the spatial distribution
of the condensate acceleration $\dot{\bm{p}}_{s}(\bm{r})$ and show that it has
a significant component outside the vortex cores, which exists even in the
absence of vortex motion. In Sec. IV we use our formalism to evaluate
microwave conductivity arising from this component, and show that it becomes
dominant in a broad interval of system parameters. We discuss the results in
Sec. V.

##  II Quasiparticle kinetics in the presence of spectral flow

In this section we present a general description of diffusive quasiparticle
kinetics in superconductors in the presence of spectral flow and inelastic
relaxation.

In the diffusive regime, where the characteristic spatial scales exceed the
elastic mean free path, $v_{F}\tau_{\mathrm{el}}$ , the quasiparticle
distribution function $n(\epsilon,{\bm{r}},t)$ depends only on energy
$\epsilon$, coordinate ${\bm{r}}$ and time $t$. In the presence of spectral
flow, its evolution equation has the form ,

$\displaystyle\partial_{t}n(\epsilon,\bm{r},t)$
$\displaystyle+v_{\nu}(\epsilon,\bm{r},t)\partial_{\epsilon}n(\epsilon,\bm{r},t)$
$\displaystyle+\nabla_{i}[D_{ij}(\epsilon,\bm{r},t)\nabla_{j}n(\epsilon,\bm{r},t)]=I_{\mathrm{in}}\\{n\\}.$
(5)

Here $D_{ij}(\epsilon,\bm{r},t)$ is the diffusion tensor222 Apart from the
crystalline anisotropy, it acquires additional anisotropy in the presence of a
superfluid momentum. In isotropic systems, it may be expressed in terms of the
longitudinal and transverse diffusion coefficients in the form
$D_{ij}(\epsilon,\bm{r},t)=n_{i}n_{j}D_{\parallel}(\epsilon,\bm{r},t)+(\delta_{ij}-n_{i}n_{j})D_{\perp}(\epsilon,\bm{r},t)$,
where $\bm{n}$ is a unit vector along $\bm{p}_{s}$. It can be shown that this
anisotropy can be substantial only in clean superconductors, whereas in the
dirty regime it is small., and $I_{\mathrm{in}}\\{n\\}$ is the collision
integral, for which we will use the relaxation time approximation,

$I_{\mathrm{in}}\\{n\\}=-\frac{n(\epsilon)-n_{F}(\epsilon)}{\tau_{\mathrm{in}}},$
(6)

with $n_{F}(\epsilon)$ being the Fermi distribution function.

The second term in Eq. (II) describes the motion of quasiparticles in energy
space caused by the spectral flow. For sufficiently slow spatial variations of
the superfluid momentum $\bm{p}_{s}$ and other system parameters, the “level
velocity” $v_{\nu}(\epsilon,\bm{r},t)$ can be expressed in terms of the time
derivative of the local density of states, $\nu(\epsilon,\bm{r},t)$. This
relation has the same form as that in uniform systems [10], and may be
obtained using conservation of the number of energy levels. From the
continuity equation for the spectral current in energy space,
$\partial_{t}\nu(\epsilon,\bm{r},t)+\partial_{\epsilon}[v_{\nu}(\epsilon,\bm{r},t)\nu(\epsilon,\bm{r},t)]=0,$
one gets

$\displaystyle v_{\nu}(\epsilon,\bm{r},t)$
$\displaystyle=\,-\frac{1}{\nu(\epsilon,\bm{r},t)}\int_{0}^{\epsilon}d\tilde{\epsilon}\partial_{t}\nu(\tilde{\epsilon},\bm{r},t).$
(7)

Similarly, the diffusion coefficient $D(\epsilon,\bm{r},t)$ in this
approximation is assumed to have the same dependence on $\bm{p}_{s}$ and
$\epsilon$ as in a uniform superconductor,
$D(\epsilon,\bar{\bm{p}}_{s}(\bm{r}))$.

Entropy production in the system is caused by diffusion of quasiparticles as
well as their inelastic relaxation. The dissipative part of the macroscopic
microwave conductivity $\sigma(\omega)$ can be obtained by equating the Joule
heating losses to the energy dissipation rate Eq. (3). Linearizing with
respect to small deviations from the equilibrium distribution function
$n=n_{F}(\epsilon)+\delta n$, the entropy production rate is given by

$\displaystyle T\dot{S}=$ $\displaystyle\,T\int d\epsilon
d\bm{r}\nu(\epsilon,\bm{r},t)\bigg{[}\nabla_{i}\delta
nD_{ij}(\epsilon,\bm{r},t)\nabla_{j}\delta n(\epsilon,\bm{r},t)$
$\displaystyle+\frac{\delta
n^{2}(\epsilon,\bm{r},t)}{\tau_{\mathrm{in}}n_{F}(\epsilon)(1-n_{F}(\epsilon))}\bigg{]}.$
(8)

Here $n_{F}(\epsilon)$ is the Fermi distribution function.

The description of quasiparticle kinetics in the presence of spectral flow
given by Eqs. (II)-(II) has a broad range of applicability. It describes both
clean and dirty superconductors, and does not assume a particular pairing
symmetry. On the other hand, the criteria of their applicability, and specific
values of $\nu(\epsilon,\bm{r},t)$, $v_{\nu}(\epsilon,\bm{r},t)$, and
$D_{ij}(\epsilon,\bm{r},t)$ depend on the pairing symmetry and other
parameters of the system. In the clean regime $\Delta\tau_{\mathrm{el}}\gg 1$,
Eqs. (II) and (7) have been derived in Ref. [10] from the conventional
Boltzmann kinetic equation for quasiparticle distribution function [13]. We
note that Eq. (II) has the same form as Eq. (A32) in Ref. [5], which describes
the time evolution of the quasiparticle distribution function $f$ in the
“dirty limit”, $\Delta\tau_{\mathrm{el}}\ll 1$, 333A more general formulation
of the kinetic equation for quasiparticles in dirty limit involves two
distribution functions, $f$ and $f_{1}$, where the distribution function
$f_{1}$ is responsible for the charge imbalance. We neglected the latter
because electromagnetic field absorption in uniform superconductors does not
create a charge imbalance. In the weakly inhomegeneous regime considered here
the generation of charge imbalance is small in the spatial gradients. At the
same time, in the presence of pair breaking, it relaxes relatively quickly by
local elastic scattering.  which was derived by Larkin and Ovchinnikov [14].
However, the local relation (7) between the level velocity
$v_{\nu}(\epsilon,\bm{r},t)$ and the time-derivative of the density of states
can be derived from the Larkin-Ovchinnikov equations only when the
characteristic length scale of spatial inhomogeneity exceeds the
superconducting coherence length, $\xi=\sqrt{D_{n}/\Delta}$, where
$D_{n}=v_{F}^{2}\tau_{\mathrm{el}}/3$, with $v_{F}$ being the Fermi velocity,
is the normal metal diffusion coefficient.

In situations where the spectral flow is caused by the pair-breaking effect of
the condensate momentum $\bm{p}_{s}$, the level velocity in Eq. (7) may be
expressed in the form

$v_{\nu}(\epsilon,\bm{p}_{\mathrm{s}})=\dot{\bm{p}}_{s}\cdot\bm{V}(\epsilon,\bm{p}_{\mathrm{s}})$
(9)

where $\bm{V}(\epsilon,\bm{p}_{s})$ denotes the level sensitivity to changes
in $\bm{p}_{s}$, and is given by

$\bm{V}(\epsilon,\bm{p}_{\mathrm{s}})=-\frac{1}{\nu(\epsilon,p_{\mathrm{s}})}\int_{0}^{\epsilon}d\epsilon\frac{\partial\nu(\epsilon,\bm{p}_{\mathrm{s}})}{\partial\bm{p}_{\mathrm{s}}}.$
(10)

For $s$-wave superconductors, the $\bm{p}_{s}$-dependence of the density of
states and level sensitivity was determined in Refs. [10, 11].

Below, we use Eqs. (II)-(10) to study microwave absorption in films of type-II
$s$-wave superconductors in the presence of a pinned vortex lattice. We show
that in the London regime, where the inter-vortex distance $\sim l_{H}$
exceeds the core radius $\sim\xi$, with the exception of small temperatures,
$T\ll\Delta$, the main contribution to the microwave absorption comes from
quasiparticles which reside at distances $\sim l_{H}$ from the vortex cores.

##  III Spatial distribution of the condensate acceleration

At low frequencies in the London regime, where the order parameter $\Delta$
outside the vortex cores is approximately uniform, the quasiparticle density
of states depends on the local instantaneous condensate momentum
$\bm{p}_{s}(\bm{r},t)$, and its rate of change is proportional to the local
condensate acceleration, $\dot{\bm{p}}_{s}(\bm{r},t)$. In thin films, the
Pearl length [19], which characterizes the screening of the magnetic field,
practically always exceeds the inter-vortex distance, $\sim l_{H}$. Therefore,
in this regime we can neglect the small inhomogeneity of the magnetic field
${\bm{H}}(\bm{r})$ caused by the diamagnetic currents.

In the presence of a microwave field ${\bm{E}}(t)$ the vortex positions
$\bm{r}_{a}(t)$ become time-dependent, and the condensate acceleration is
given by,

$\displaystyle\dot{\bm{p}}_{s}(\bm{r},t)=$
$\displaystyle\,\hbar\hat{z}\times\sum_{a}\left[\frac{2(\bm{r}-\bm{r}_{a}(t))(\bm{r}-\bm{r}_{a}(t))\cdot\dot{\bm{r}}_{a}(t)}{|\bm{r}-\bm{r}_{a}(t)|^{4}}\right.$
$\displaystyle\left.-\frac{\dot{\bm{r}}_{a}(t)}{|\bm{r}-\bm{r}_{a}(t)|^{2}}\right]+e{\bm{E}}(t).$
(11)

Thus, only in the absence of vortex displacement the condensate acceleration
is given by the second term, $e\bm{E}(t)$; the modification caused by the
motion of the vortices is described by the first term.

Let us consider the spatial distribution of $\dot{\bm{p}}_{s}(\bm{r},t)$
inside the plaquette of a given vortex $a$, that is at
$|\bm{r}-\bm{r}_{a}|\lesssim l_{H}$. The term $a$ in the sum in Eq. (III),
which corresponds to the motion of the native core, $\dot{\bm{r}}_{a}$,
produces a near-field contribution which decays rapidly with the distance from
the core,

$\displaystyle\dot{\bm{p}}^{(n)}_{s}(\bm{\rho}_{a})=$
$\displaystyle\,\hbar\hat{z}\left[\frac{2\bm{\rho}_{a}\,(\bm{\rho}_{a}\cdot\dot{\bm{r}}_{a})-|\bm{\rho}_{a}|^{2}\dot{\bm{r}}_{a}}{|\bm{\rho}_{a}|^{4}}\right].$
(12)

Here we introduce the notation $\bm{\rho}_{a}=\bm{r}-\bm{r}_{a}$. The
remaining contribution, which is produced by the motion of the other vortices
together with the second term, does not fall off with distance and may be
approximated by a constant for $|\bm{r}-\bm{r}_{a}(t)|\lesssim l_{H}$. We
refer to it as the far field contribution and denote it by
$\dot{\bm{p}}^{(f)}_{s}$. By order of magnitude, $\dot{\bm{p}}^{(f)}_{s}$
coincides with the spatial average of the condensate acceleration in the
system, $\langle\dot{\bm{p}}_{s}(t)\rangle$. In the linear approximation in
the microwave field and at sufficiently low frequencies, where the viscous
forces are negligible in comparison to pinning,
$\langle\dot{\bm{p}}_{s}(t)\rangle$ is related to the electric field by the
Campbell formula,

$\displaystyle\dot{\bm{p}}^{(f)}_{s}\sim$
$\displaystyle\,\langle\dot{\bm{p}}_{s}(t)\rangle=\frac{\lambda_{L}^{2}}{\lambda_{C}^{2}(H)}\,e{\bm{E}}(t).$
(13)

Here $\lambda_{L}$ is the London length, and $\lambda_{C}(H)$ is the Campbell
[17] length. The latter depends on the pinning strength and characterizes the
macroscopic superfluid density of the system444 Equation (13) reflects the
fact that the time-derivative of the superfluid transport current can be
expressed in two equivalent forms,
$\langle\dot{\bm{j}}\rangle=\frac{c^{2}}{4\pi
e\lambda_{L}^{2}}\langle\dot{\bm{p}}_{s}(t)\rangle=\frac{c^{2}}{4\pi
e\lambda_{C}^{2}(H)}e{\bm{E}}(t)$.. Similarly, the typical vortex velocity
$\dot{\bm{r}}_{a}$, which determines the magnitude of
$\dot{\bm{p}}^{(n)}_{s}(\bm{u}_{a})$ may be estimated as

$\displaystyle\dot{\bm{r}}_{a}$
$\displaystyle\sim\,{\bm{v}}^{\prime}_{v}(t)=c\,\frac{{\bm{E}}(t)\times{\bm{H}}}{H^{2}}\left(1-\frac{\lambda_{L}^{2}}{\lambda_{C}^{2}(H)}\right).$
(14)

The second term in the brackets on the right hand side describes the pinning-
induced reduction of the average vortex velocity from the value in Eq. (2).

Using the expression for the typical vortex velocity
${\bm{v}}^{\prime}_{v}(t)$ in Eq. (14) the Bardeen-Stephen contribution to the
conductivity for a pinned vortex lattice may be estimated as

$\displaystyle\frac{\sigma_{BS}^{\prime}}{\sigma_{BS}}$
$\displaystyle\sim\left(1-\frac{\lambda_{L}^{2}}{\lambda_{C}^{2}(H)}\right)^{2}.$
(15)

Assuming that the pinning strength is determined by the vortex cores and is
independent of the magnetic field, the ratio of the Campbell and London length
may be expressed in the form [17, 18]

$\displaystyle\frac{\lambda^{2}_{C}(H)}{\lambda^{2}_{L}}=$
$\displaystyle\,1+\frac{\Phi_{0}Hd}{8\pi\lambda_{L}^{2}k}.$ (16)

Here $k$ is the average “spring constant”, which relates the average pinning
force on the vortex $\bm{F}_{pin}=-k\delta\bm{r}$, to the average vortex
displacement $\delta\bm{r}$.

The difference $\lambda^{2}_{C}(H)/\lambda^{2}_{L}-1$ characterizes the
effectiveness of pinning. At perfect pinning, $k\to\infty$, we have
$\lambda^{2}_{C}(H)/\lambda^{2}_{L}-1\to 0$ and
$\langle\dot{\bm{p}}_{s}\rangle\to e{\bm{E}}$. At finite $k$, pinning becomes
more effective as the magnetic field is reduced;
$\lambda^{2}_{C}(H)/\lambda^{2}_{L}-1$ decreases. In particular at strong
pinning, where the condensation energy of the vortex core changes by a factor
of the order unity when the vortex is displaced by a distance $\xi$ from the
equilibrium position, we have $k\sim d\Delta^{2}\nu_{n}$. In this case the
pinning effectiveness from Eq. (16) may be estimated as

$\displaystyle\frac{\lambda^{2}_{C}(H)}{\lambda_{L}^{2}}-1\sim$
$\displaystyle\,\frac{H}{H_{c2}}.$ (17)

Here we used the fact that for dirty superconductors
$\lambda_{L}^{-2}\sim\Delta\sigma_{n}/\hbar c^{2}\sim
e^{2}\nu_{n}D_{n}\Delta/\hbar c^{2}$.

In the next section, we evaluate the dissipation arising from the two
contributions to the condensate acceleration. The dissipation caused by the
near field part of the condensate acceleration,
$\dot{\bm{p}}^{(n)}_{s}(\bm{\rho}_{a})$, is localized to the vortex cores and
corresponds to the Bardeen-Stephen contribution to the conductivity, Eq. (15).
The dissipation due to the far field part, $\dot{\bm{p}}^{(f)}_{s}$, arises
from quasiparticles outside the vortex cores. We will show that these
quasiparticles produce the dominant contribution to the conductivity in a wide
interval of physical parameters.

##  IV Estimates of microwave conductivity

We now apply Eqs. (II), (6), and (II) to evaluate the microwave conductivity
$\sigma(\omega)$ in the London regime, where the inter-vortex distance $\sim
l_{H}$ significantly exceeds the vortex core size $\sim\xi$ . We focus on the
contribution to $\sigma(\omega)$, which arises from quasiparticles residing
outside the vortex cores, where the level velocity is described by Eqs. (9)
and (10). For simplicity we assume $T\sim\Delta$ and focus on the dirty limit,
$\Delta\tau_{\mathrm{el}}\ll 1$. In this case the dependence of the
quasiparticle density of states on the local superfluid momentum
$\bar{p}_{s}(\bm{r})$ is characterized by the dimensionless parameter

$\displaystyle\eta(p_{s})$ $\displaystyle=D_{n}p_{s}^{2}/\Delta.$ (18)

The gap in the quasiparticle spectrum is lowered from $\Delta$ by the amount

$\displaystyle\delta\epsilon(p_{s})$
$\displaystyle\sim\Delta\eta^{2/3}({p_{s}}).$ (19)

For $\epsilon>\Delta+\delta\epsilon(p_{s})$, the sensitivity
$\bm{V}(\epsilon,\bm{p}_{s})$ and $\nu(\epsilon,p_{s})$ are rapidly decreasing
functions of energy $\epsilon$ [10]. Therefore, the dominant contribution to
the microwave conductivity arises from quasiparticles with energies
$|\epsilon-\Delta|\lesssim\delta\epsilon(p_{s})$. In this energy interval the
values of the level sensitivity $\bm{V}(\epsilon,\bm{p}_{s})$, and the density
of states $\nu(\epsilon,p_{s})$ may be estimated as [10],

$\displaystyle\bm{V}(\epsilon,\bm{p}_{s})$ $\displaystyle\sim
D_{n}\bm{p}_{s}\eta^{-1/3}(p_{s}),$ (20a) $\displaystyle\nu(\epsilon,p_{s})$
$\displaystyle\sim\nu_{n}\eta^{-1/3}(p_{s}),$ (20b)

where $\nu_{n}$ is the normal state density of states at the Fermi level.

Although in the presence of superfluid momentum the diffusion tensor is
anisotropic, for $\Delta\tau_{\mathrm{el}}\ll 1$ the anisotropy is negligible,
and we set $D_{ij}=\delta_{ij}D(\epsilon,\bar{p}_{s})$. The value of the
diffusion coefficient in the relevant energy interval may be estimated as

$\displaystyle D\left(\bar{p}_{s}\right)\sim D_{n}\eta^{1/3}.$ (21)

This estimate can be obtained555 In general the diffusion coefficient is
expressed in terms of the quasiclassical Green’s functions $g^{R}$ and $g^{A}$
(see the appendix in [5]), which must be calculated using Usadel’s equation
[15]. Since the Usadel’s equation has already been solved in the relevant
regime [11], we will omit the details of this calculation. using the Larkin-
Ovchinnikov equations [14].

###  IV.1 Low frequency regime

At the lowest frequencies, the microwave conductivity is dominated by
inelastic relaxation processes. The reason is that diffusion cannot lead to
full relaxation of the nonequilibrium quasiparticle density for quasiparticles
with energies below the percolation threshold $\epsilon^{*}$. Such
quasiparticles are trapped inside a plaquette of a particular vortex. The size
of the trapping region $r(\epsilon)$ for energy $\epsilon$ may be estimated
using Eqs. (18) and (19) and the the dependence of the superfluid momentum
$\bar{p}_{s}(r)$ on the distance $r$ to the vortex core. Since for $r\lesssim
l_{H}$

$\displaystyle\bar{p}_{s}(r)\sim$ $\displaystyle\,\frac{\hbar}{r},$ (22)

we get,

$\displaystyle r(\epsilon)$
$\displaystyle\sim\xi\left(\frac{\Delta}{\Delta-\epsilon}\right)^{3/4}.$ (23)

The size of the trapping region increases with the quasiparticle energy
$\epsilon$ and becomes of the order of the inter-vortex distance $\sim l_{H}$
as $\epsilon$ approaches the percolation threshold $\epsilon^{*}$ whose value
may be estimated as

$\displaystyle\Delta-\epsilon^{*}$
$\displaystyle\sim\Delta\left(\frac{\xi}{l_{H}}\right)^{4/3}.$ (24)

At frequencies smaller than the inverse diffusion time across the inter-vortex
distance, $\tau_{\mathrm{D}}^{-1}$ the nonequilibrium distribution function of
the trapped quasiparticles becomes spatially-uniform and depends only on the
energy. This part of the distribution (zero mode of diffusion) can relax only
via inelastic processes.

To describe this slow time relaxation of the zero mode we linearize Eq. (II)
with respect to $\delta n$ and average the result over the area of spatial
confinement at energy $\epsilon$ in the $a$-th vortex. This yields

$\left(\partial_{t}+\frac{1}{\tau_{\mathrm{in}}}\right)\langle\delta
n(\epsilon,t)\rangle_{a}=-\langle\dot{{\bm{p}}}_{s}(\bm{r})\cdot\bm{V}(\epsilon,\bar{\bm{p}}_{\mathrm{s}}({\bm{r}}))\,\rangle_{a}\partial_{\epsilon}n_{F}(\epsilon),$
(25)

where $\langle\ldots\rangle_{a}$ denotes averaging over the
$\epsilon$-dependent confinement region of the $a$-th vortex.

It is important to note that in the case of perfectly symmetric lattice the
right hand side of Eq. (25) vanishes. However in the presence of disorder, the
superfluid momentum around the vortices is asymmetric, and this term is
nonzero. Thus, $\delta\langle n(\epsilon)\rangle_{a}$ and
$\langle\dot{{\bm{p}}}_{s}(\bm{r})\cdot\bm{V}(\epsilon,\bar{\bm{p}}_{s}(\bm{r}))\,\rangle_{a}$
are random quantities which fluctuate from plaquette to plaquette. Equations
(20a) and (22) show that these quantities are dominated by distances from the
core, which are of the order of the radius of the trapping region in Eq. (23).
In this region $\dot{{\bm{p}}}_{s}(\bm{r})$ may be approximated by
$\dot{\bm{p}}^{(f)}_{s}(\bm{r})$ in Eq. (13). Making this approximation,
substituting $\delta\langle n(\epsilon)\rangle_{a}$ from Eq. (25) into Eq.
(II), averaging over plaquettes, and using Eq. (3), we obtain the following
result for the microwave conductivity $\sigma(\omega)$ at $T\gtrsim\Delta$,

$\displaystyle\frac{\sigma(\omega)}{\sigma^{\prime}_{BS}}\sim$
$\displaystyle\,\frac{\tau_{\mathrm{in}}}{\tau_{\mathrm{el}}}\frac{1}{\left[1+\left(\omega\tau_{\mathrm{in}}\right)^{2}\right]}\,\frac{H}{H_{c2}}\frac{\lambda_{L}^{4}(H)}{\left(\lambda_{C}^{2}(H)-\lambda_{L}^{2}\right)^{2}}$
$\displaystyle\times\int_{0}^{\epsilon^{*}}d\epsilon\frac{\left\langle\langle\nu(\epsilon,\bar{p}_{s}(\bm{r})\rangle_{a}\langle{\bm{n}}\cdot\bm{V}(\epsilon,\bm{r})\rangle_{a}^{2}\right\rangle}{T\nu_{n}v_{\mathrm{F}}^{2}},$
(26)

where the outer brackets indicate averaging over plaquettes, and we introduced
the unit vector ${\bm{n}}$ in the direction of ${\bm{E}}$. To arrive at this
expression we used Eqs. (1), (15), and the Einstein relation for the normal
state conductivity, $\sigma_{n}=e^{2}D_{n}\nu_{n}$. We note that for the zero
mode limit only the second term in the brackets in Eq. (II) contributes to the
entropy production rate.

Substituting the estimates (20), (22), and (23) into Eq. (IV.1) we find that
the dominant contribution to microwave conductivity arises from trapped
quasiparticles with energies near the percolation threshold $\epsilon^{*}$
given by Eq. (24). Using (16), and the relation $\xi^{2}/l_{H}^{2}\sim
H/H_{c2}$, we can express the result in the form

$\frac{\sigma(\omega)}{\sigma_{BS}^{\prime}}\sim
K\frac{\Delta}{T}\frac{\tau_{\mathrm{in}}\Delta}{\left[1+\left(\omega\tau_{\mathrm{in}}\right)^{2}\right]}\frac{\lambda_{L}^{4}(H)}{\left[\lambda_{C}^{2}(H)-\lambda_{L}^{2}\right]^{2}}\left(\frac{H}{H_{c2}}\right)^{5/3},$
(27)

where the parameter

$K=\frac{\left\langle\langle{\bm{n}}\cdot\bm{V}(\epsilon^{*},\bm{r})\rangle_{a}^{2}\right\rangle}{\langle\langle
V(\epsilon^{*},\bar{p}_{s}(\bm{r})\rangle_{a}^{2}\rangle}$ (28)

characterizes the degree of spatial asymmetry of the condensate momentum
$\bm{p}_{s}(\bm{r})$ inside the trapping areas in the pinned vortex lattice.
The value of $K$ is non-universal and depends on the details of the pinning
potential. If the relative amplitude of fluctuations of
$\bm{V}(\epsilon,\bm{r})$ is of order unity, and their correlation radius is
on the order of the plaquette size, then $K\approx 1$. Using the estimate (17)
we obtain for strong pinning,

$\frac{\sigma(\omega)}{\sigma_{BS}^{\prime}}\sim\frac{\Delta}{T}\frac{\tau_{\mathrm{in}}\Delta}{\left[1+\left(\omega\tau_{\mathrm{in}}\right)^{2}\right]}\left(\frac{H_{c2}}{H}\right)^{1/3}.$
(29)

Since in typical situations $\tau_{\mathrm{in}}\Delta\gg 1$, Eqs. (27) and
(29) show that the inelastic relaxation gives the main contribution to the
conductivity in the low frequency regime.

###  IV.2 High frequency regime

At $\omega\tau_{\mathrm{in}}>1$ the contribution of inelastic relaxation to
$\sigma(\omega)$ in Eq. (27) decreases as $1/\omega^{2}$, and at sufficiently
high frequencies the microwave conductivity is dominated by spatial diffusion
of nonequilibrium quasiparticles inside the trapping regions.

The diffusive contribution to dissipation is present even in a perfect vortex
lattice. Therefore, to estimate it we neglect the distortion of the vortex
lattice. We assume the axially symmetric distribution of the superfluid
momentum $\bar{\bm{p}}_{s}$ inside a given plaquette, Eq. (22), and write the
diffusion equation Eq. (II) in polar coordinates, $r$ and $\theta$. The
dependence of the level velocity $v_{\nu}(r,\theta)$ in Eq. (9) and
nonequilibrim quasiparticle density $\delta n$ on the azimuth angle
corresponds to the first angular harmonic,

$\delta n(\epsilon,r,\theta;t)=n_{1}(\epsilon,r;t)\cos\theta.$

Substituting this form into the linearized Eq. (II) we get

$\displaystyle\partial_{t}n_{1}(\epsilon,r;t)+\frac{1}{r}\partial_{r}\left(rD(\epsilon,\bar{p}_{s}(r))\partial_{r}n_{1}(\epsilon,r;t)\right)$
$\displaystyle-\frac{D(\epsilon,\bar{p}_{s}(r))n_{1}(\epsilon,r;t)}{r^{2}}=-\langle\dot{p}_{s}^{(f)}\rangle
V(\epsilon,\bar{p}_{s}(r))\partial_{\epsilon}n_{F}.$ (30)

Here we neglected the inelastic collision integral and focused on the long-
distance part of the condensate acceleration given by Eq. (13).

For the relevant quasiparticle energies, the diffusion coefficient
$D(\epsilon,\bar{p}_{s}(r))$, and the level sensitivity
$V(\epsilon,\bar{p}_{s}(r))$ have a power-law dependence of $r$, which is
given by Eqs. (18), (20a), (21), and (22). For a microwave field of frequency
$\omega$, the solution of this equation is characterized by the length scale
$L_{\omega}$ corresponding to the diffusion distance of the relevant
quasiparticles during the oscillation period,
$D\left(\bar{p}_{s}(L_{\omega})\right)/L_{\omega}^{2}=\omega$, and may be
estimated as

$\displaystyle L_{\omega}$
$\displaystyle\sim\xi\left(\frac{\Delta}{\omega}\right)^{3/8}.$ (31)

For $\xi\ll r\ll L_{\omega}$, and $L_{\omega}\ll r$ the solution to Eq. (IV.2)
has a power-law form,

$\displaystyle n_{1}\sim$
$\displaystyle\,eE\frac{\lambda_{L}^{2}}{\lambda^{2}_{C}(H)}\partial_{\epsilon}n_{F}\left\\{\begin{array}[]{cc}\frac{r^{7/3}}{\xi^{4/3}},&r\ll
L_{\omega},\\\ -i\frac{L_{\omega}^{8/3}}{r^{1/3}\xi^{4/3}},&r\gg
L_{\omega}.\end{array}\right.$ (34)

Substituting Eq. (34) into the first term of the RHS of Eq. (II) and using Eq.
(13) we get

$\displaystyle\frac{\sigma(\omega)}{\sigma_{BS}^{\prime}}\sim$
$\displaystyle\,\frac{\Delta}{T}\frac{\lambda_{L}^{4}}{(\lambda_{C}^{2}-\lambda_{L}^{2})^{2}}\\!\left(\frac{H}{H_{c2}}\right)^{1/3}\\!\\!\left\\{\begin{array}[]{cc}1,&\omega\ll\tau_{D}^{-1},\\\
(\omega\tau_{D})^{-5/4},&\omega\gg\tau_{D}^{-1}.\end{array}\right.$ (37)

Here $\tau_{D}$ the diffusion time of the relevant quasiparticles across the
plaquette, $\tau_{D}=l_{H}^{2}/D(\bar{p}_{s}(l_{H}))$. It may be estimated as

$\displaystyle\tau_{D}\sim$
$\displaystyle\,\frac{\hbar}{\Delta}\left(\frac{l_{H}}{\xi}\right)^{8/3}.$
(38)

Comparing Eq. (27) and Eq. (37) we conclude that, for
$\omega<\tau^{-1}_{\mathrm{D}}$, diffusive relaxation given the dominant
contribution to the conductivity for $\omega\gtrsim\omega^{*}$, where the
crossover frequency $\omega^{*}$ may be estimated as

$\displaystyle\omega^{*}\sim$
$\displaystyle\,\sqrt{\frac{K\Delta}{\tau_{\mathrm{in}}}}\left(\frac{H}{H_{c2}}\right)^{2/3}\sim(\tau_{\mathrm{in}}\tau_{\mathrm{D}})^{-1/2}.$
(39)

At strong pinning, see Eq. (17), the expression for the conductivity in Eq.
(37) simplifies to

$\displaystyle\frac{\sigma(\omega)}{\sigma_{BS}^{\prime}}\sim$
$\displaystyle\,\frac{\Delta}{T}\left(\frac{H_{c2}}{H}\right)^{5/3}\\!\\!\left\\{\begin{array}[]{cc}1,&\omega\ll\tau_{D}^{-1},\\\
(\omega\tau_{D})^{-5/4},&\omega\gg\tau_{D}^{-1}.\end{array}\right.$ (42)

In this case, in a broad frequency interval,
$\omega^{*}<\omega<\tau_{D}^{-1}$, the _ac_ conductivity is dominated by
diffusion of quasiparticles outside the vortex cores for sufficiently weak
magnetic fields $H\ll H_{c2}$.

Finally, we note that replacing $\langle\dot{p}_{s}^{(f)}(r)\rangle$ in Eq.
(IV.2) by $\dot{p}_{s}^{(n)}(r)$ from Eq. (12) (which has the same angular
dependence) gives a contribution to dissipation, which is dominated by short
distances from the core, and reproduces the Bardeen-Stephen result
$\sigma^{\prime}_{BS}$.

##  V discussion of the results

We developed a theory of microwave absorption of type-II superconductors in
the presence of a strongly pinned vortex lattice. In this case, in addition to
the Bardeen-Stephen contribution to the dissipative conductivity,
$\sigma^{\prime}_{BS}$, which is caused by the vortex motion and is described
by Eqs. (1), (15), and (16), there is an another contribution. This new
contribution is caused by the quasiparticle spectral flow induced by the
condensate acceleration in the presence of a microwave field. This
contribution exists even in the absence of vortex displacements and exceeds
$\sigma_{BS}^{\prime}$ in a wide interval of physical parameters. For
$T\gtrsim\Delta$ this contribution is dominated by the quasiparticles residing
outside the vortex cores.

At low frequencies the relaxation of the non-equilibrium distribution
quasiparticles induced by the spectral flow is mediated by inelastic
processes. As a result, the dissipative conductivity $\sigma(\omega)$ in this
regime is controlled by the inelastic relaxation time $\tau_{\mathrm{in}}$,
see Eq. (27). At $\omega\ll\tau_{\mathrm{in}}$ the lowest frequencies
$\sigma(\omega)$ becomes proportional to $\tau_{\mathrm{in}}$. Therefore it
can be parametrically larger than $\sigma^{\prime}_{BS}$, which is
proportional to the elastic mean free time $\tau_{\mathrm{el}}$. The mechanism
of absorption in this case is similar to the Debye mechanism of
electromagnetic wave absorption is molecules [20]. The difference however is
that the Debye absorption coefficient in molecules is quadratic in frequency,
while the dissipative part of conductivity in superconductors
$\sigma(\omega)$, is nonzero at $\omega\to 0$. The low frequency conductivity
is proportional to the degree of distortion of the spatial distribution of the
superfluid momentum in the lattice, which is caused by the pinning and
characterized by the parameter $K$ in Eq. (28).

At frequencies above $\omega^{*}$ given by Eq. (39) the conductivity is
dominated by diffusion of quasiparticles across the vortex. The characteristic
time scale $\tau_{D}$ for this process is given by Eq. (38), and is still much
larger than $\tau_{\mathrm{el}}$. In this regime the conductivity exhibits a
non-trivial frequency dependence given by Eq. (37). The distortion of the
vortex lattice induced by the pining becomes inessential in this regime.

At small temperatures, $T\ll\Delta$, the quasiparticles reside only in the
vortex cores. In this case the low frequency contribution to conductivity,
which is proportional to $\tau_{\mathrm{in}}$ still exists, but the
coefficient $K$ in Eq. (28) is determined by the deformation of the core.

We note that in the case of weak pinning there is another mechanism of low
frequency microwave absorption, which is proportional to $\tau_{\mathrm{in}}$,
[18]. It is related to vortex motion and exists also in the flux flow regime
[12]. The main contribution to $\sigma$ in this case comes from quasiparticles
in the vortex cores.

Although the estimates for the micirowave conductivity were obtained in the
dirty regime, $\Delta\tau_{\mathrm{el}}\ll 1$, qualitatively, our conclusions
remain valid for arbitrary value of $\Delta\tau_{\mathrm{el}}$, provided
$l_{H}\gg v_{F}\tau_{\mathrm{el}}$. In the clean case,
$\Delta\tau_{\mathrm{el}}\gg 1$ one would need to use different estimates for
the sensitivities and relevant intervals, which can be found in Ref. [10].

We would like to thank M. Feigel’man and V. Geshkenbein for useful
discussions. The work of T.L. and A.A. was supported by the National Science
Foundation through the MRSEC Grant No. DMR-1719797, the Thouless Institute for
Quantum Matter (TIQM), and the College of Arts and Sciences at the University
of Washington.

## References

  * [1] A.A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics (Courier Corporation, 1975) google-Books-ID: E 9NtwNY7UcC.
  * [2] J. Bardeen and M. J. Stephen, ”Theory of the Motion of Vortices in Superconductors”, Phys. Rev. 140, A1197 (1965).
  * [3] P. Nozi‘eres and W. F. Vinen, “The motion of flux lines in type II superconductors,” The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 14, 667–688 (1966), publisher: Taylor & Francis eprint: https://doi.org/10.1080/14786436608211964.
  * [4] L. P. Gorkov, and N.B. Kopnin, “Vortex motion and resistivity of type-II superconductors in magnetic field”, Sov. Phys. Uspehi, 18, 496, 1975.
  * [5] A.I. Larkin and Yu. N. Ovchinnikov, “Vortex motion in superconductors”, in Nonequilibrium superconductivity” ed by D.N. Langerberg and A.I. Larkin, Elsevier Science Publishers B.V., 1986.
  * [6] G. Blatter, M. V. Feigel’man, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, “Vortices in high-temperature superconductors,” Reviews of Modern Physics 66, 1125–1388 (1994), publisher: American Physical Society.
  * [7] J. I. Gittleman and B. Rosenblum, ”The Pinning Potential and High-Frequency Studies of Type-II Superconductors”, J. Appl. Phys. 39, 2617–2621 (1968)
  * [8] X.-Q. Zhou, C. J. S. Truncik, W. A. Huttema, N. C. Murphy, P. J. Turner, A. J. Koenig, R.-X. Liang, D. A. Bonn, W. N. Hardy, and D. M. Broun, ”Microwave spectroscopy of vortex dynamics in ortho-II YBa2Cu3O6.52”, Phys. Rev. B 87, 184512 – Published 28 May 2013
  * [9] Mark W. Coffey and John R. Clem, ”Unified theory of effects of vortex pinning and flux creep upon the rf surface impedance of type-II superconductors”, Phys. Rev. Lett. 67, 386 – Published 15 July 1991
  * [10] M. Smith, A.V. Andreev, B.Z. Spivak, ” Giant microwave absorption in s- and d- wave superconductors”, Annals of Physics, Volume 417, 2020, 168105, ISSN 0003-4916, https://doi.org/10.1016/j.aop.2020.168105.
  * [11] M. Smith, A. V. Andreev, B.Z. Spivak, Debye mechanism of giant microwave absorption in superconductors, Phys. Rev. B 101, 134508 (2020)
  * [12] M. Smith, A. V. Andreev, M. V. Feigel’man, and B. Z. Spivak,  ”Conductivity of superconductors in the flux flow regime”,  Phys. Rev. B 102, 180507(R) – Published 24 November 2020.
  * [13] A. G. Aronov, Yu M. Gal’perin, V. L. Gurevich, and V. I. Kozub.  The Boltzmann-equation description of transport in superconductors.  _Advances in Physics_ , 30(4):539–592, 1981.  doi: 10.1080/00018738100101407.  URL https://doi.org/10.1080/00018738100101407.
  * [14] A I Larkin and Yu N Ovchinnikov.  Nonlinear effects during the motion of vortices in superconductors.  _JETP_ , 46(1):155, 1977.
  * [15] K. D. Usadel,  Generalized Diffusion Equation for Superconducting Alloys.  Phys. Rev. Lett. 25, 507 (1970).
  * [16] Cambridge Mathematical Library: A course of modern analysis, ET Whittaker, GN Watson, Ch. 20 (2000).
  * [17] A. M. Campbell,  ”The response of pinned flux vortices to low-frequency fields”,  Journal of Physics C: Solid State Physics 2, 1492  doi: 10.1088/0022-3719/2/8/318
  * [18] B. V. Pashinsky, M. V. Feigel’man, A. V. Andreev,  ”Microwave response of type-II superconductors at weak pinning”,  SciPost Phys. 14, 096 (2023),  https://doi.org/10.21468/SciPostPhys.14.5.096
  * [19] J. Pearl,  ”Current distribution in superconducting films carrying quantized fluxoids”,  Applied Physics Letters. 5 (4): 65–66 (1964).  doi: 10.1063/1.1754056
  * [20] Peter Debye.  _Polar molecules_.  Dover Publ., 1970.  Google-Books-ID: f70ingEACAAJ.


# Floquet engineering of titled and gapped Dirac materials

Andrii Iurov1111E-mail contact<EMAIL_ADDRESS><EMAIL_ADDRESS>Liubov Zhemchuzhna1,2, Godfrey Gumbs2,3, Danhong
Huang4,5, Kathy Blaise1 and Chinedu Ejiogu1 1Department of Physics and
Computer Science, Medgar Evers College of City University of New York,
Brooklyn, NY 11225, USA
2Department of Physics and Astronomy, Hunter College of the City University of
New York, 695 Park Avenue, New York, New York 10065, USA
3Donostia International Physics Center (DIPC), P de Manuel Lardizabal, 4,
20018 San Sebastian, Basque Country, Spain
4US Air Force Research Laboratory, Space Vehicles Directorate, Kirtland Air
Force Base, New Mexico 87117, USA
5Center for High Technology Materials, University of New Mexico, 1313 Goddard
SE, Albuquerque, New Mexico, 87106, USA
6US Military Academy at West Point, 606 Thayer Road, West Point, New York
10996, USA

###### Abstract

We have established a rigorous theoretical formalism for Floquet engineering,
or investigating and eventually tailoring most crucial electronic properties
of tetragonal molybdenum disulfide (1T′-MoS2), by applying an external high-
frequency dressing field in the off-resonant regime. It was recently
demonstrated that monolayer semiconducting1T′-MoS2 may assume a distorted
tetragonal structure which exhibits tunable and gapped spin- and valley-
polarized tilted Dirac bandstructure. From the viewpoint of electronics,
1T′-MoS2 is one of the most technologically promising nanomaterials and a
novel representative of an already famous family of transition metal
dichalcogenides. The obtained dressed states strongly depend on the
polarization of the applied irradiation and reflect the full complexity of the
initial low-energy Hamiltonian of non-irradiated material. We have calculated
and analyzed the obtained electron dressed states for linear and circular
types of the polarization of the applied field focusing on their symmetrical
properties, anisotropy, tilting and bandgaps, as well as topological
signatures. Since a circularly polarized dressing field is also known to
induce a transition into a new state with broken time-reversal symmetry and a
non-zero Chern number, the combination of these topologically non-trivial
phases and transitions between them could reveal some truly unique and earlier
unknown phenomena.

##  I Introduction

Floquet theory, or even Floquet engineering, which describes electronic
behavior of a wide range of quantum-mechanical systems under a periodic field,
Kibis et al. (2020, 2018) has become extremely popular over the last several
years. It is important to note that recent experimental advances in optical
and microwave physics, laser technology and emergent technical applications
for condensed-matter quantum optics have enabled the experimental verification
of such theoretical predictions and their applications in actual
optoelectronic devices. Cheng et al. (2019); Weitenberg and Simonet (2021);
Wang et al. (2018); Dehghani et al. (2014); Nakagawa et al. (2020)

The effect of the off-resonance dressing field strongly depends on the
polarization of this field. Circularly polarized dressing field is known to
open a bandgap and break the time-reversal symmetry in graphene. An important
technical challenge addressed by Floquet engineering is the confining
ballistic electrons within a specific spatial region of an optoelectronic
device which is directly related to Klein paradox Iurov et al. (2020a) in our
material and formation of the localized states Pereira et al. (2006); Gumbs et
al. (2014a); Schnez et al. (2010); Castro et al. (2008); Gumbs et al. (2014b)
due to an induced energy bandgap. Kibis (2010); Iurov et al. (2011) In
contrast, a linearly polarized optical field induces an in-plane anisotropy
and could affect the sequential-tunneling current of doped electrons for a
non-zero polarization angle. Ibarra-Sierra et al. (2022); Iurov et al. (2022)
A gap could also be opened in high-intensity field regime of the circularly
polarized irradiation. citesandoval2020floquet

The electron-photon dressed states been studied in a variety of two-
dimensional materials: Castro et al. (2022); Oka and Kitamura (2018)
nanotubes, Kibis et al. (2021); Hsu and Reichl (2006); Sentef et al. (2015)
graphene, Kibis (2010); Calvo et al. (2011); Vogl et al. (2021); Dal Lago et
al. (2017) silicene, Iurov et al. (2017a) transitional metal dichalcogenides,
Kibis et al. (2017) dice lattice and $\alpha-\mathcal{T}_{3}$, Iurov et al.
(2019); Dey and Ghosh (2019, 2018); Weekes et al. (2021); Tamang et al. (2021)
various types of nanoribbons Wang et al. (2021); Tahir et al. (2016); Iurov et
al. (2021), anisotropic phosphorene Iurov et al. (2017b) and others. Calvo et
al. (2020, 2015); Iurov et al. (2013); Mojarro et al. (2020) The effect of
circularly polarized dressing field was examined in silicene, one of the
isotropic limiting cases for 1T′-MoS2. It was found that a circularly
polarized field breaks the equivalence of the two valleys and, most
surprisingly, could either increase or decrease the bandgap depending on its
initial value in a stark contrast to any known material. Kibis et al. (2017)
It was also demonstrated that the off-resonant field could substantially
modify transport properties Iurov et al. (2022); Kristinsson et al. (2016);
Iurov et al. (2020b), excitonic behavior Iorsh et al. (2022) and topological
signatures of a two-dimensional lattice. Dey et al. (2020); Dey and Ghosh
(2019); Dal Lago et al. (2015); Iorsh et al. (2017)

A molybdenum disulfide MoS2 is one of monolayer transition-metal
dichalcogenides with a distorted tetragonal structure. Tang et al. (2017);
Zheng et al. (2016); Fei et al. (2017) In general, the monolayer
semiconducting MoS2 structure may assume a trigonal prismatic coordination of
the metal atoms (2H with hexagonal symmetry), as well as octahedral
coordination (1T with tetragonal symmetry) or 1T′ with a distorted and the
most exotic Liu et al. (2018); Mehmood et al. (2021) which is exactly the
subject of our study.

This is the one of the most innovative two-dimensional Dirac materials which
is thermodynamically stable and hence easily synthesized in its semiconducting
phase. It was also theoretically predicted to demonstrate a strong quantum
spin Hall effect and a great potential for optoelectronic and other
applications.In the presence of an external transverse electric field, this
material exhibits valley-spin-polarized Dirac bands and a phase transition
between the topological insulator and a regular band insulator similarly to
silicene. The energy band structure also demonstrates a special type of
anisotropy which is referred to as tilted Dirac bands. Gomes and Ramos (2021);
Tan et al. (2021a)

1T′-MoS2 exhibits unique features, such as tunable anisotropy, bandgap, spin-
and valley-polarized states and coexistence of different topological phases
similarly to silicene. Tabert and Nicol (2013); Kara et al. (2012); Zhao et
al. (2016); Drummond et al. (2012); Ezawa (2012) Some of these properties are
determined by the value of external electric field and thus could be
controlled, while the others are fixed by an internal spin-orbit coupling gap
and other lattice parameters of 1T′-MoS2. The zero-gap limit of 1T′-MoS2 bands
is found in 8-Pmmn borophene which was intensively studied overt the recent
years because of a strong anisotropy in all of its crucial physical
properties. Mannix et al. (2015); Lopez-Bezanilla and Littlewood (2016); Islam
(2018); Paul et al. (2019); Tan et al. (2021b) Very recently, a plasmonic gain
was reported in current-biased tilted Dirac nodes Park et al. (2022) Their
optical properties and optical conductivity have been also thoroughly
investigated. Mojarro et al. (2021); Tan et al. (2021a) Other monolayer
materials with similar tilted band structure, such as TaIrTe4, TaCoTe2 or
$\alpha-SnS_{2}$ have been recently synthesized.

The remaining part of the present paper is organized as follows: in Sec. II we
review some crucial properties of the low-energy Hamiltonian, dispersions and
the corresponding electronic states of 1T′-MoS2: tilting and anisotropy,
direct and indirect bandgaps. Next section III is intended to derive and
analyze the electron-photon dressed states in the presence of electromagnetic
field with both linear and circular polarizations. We have derived the
electron-light interaction Hamiltonian, found and thoroughly analyzed the
energy bandstructure of the irradiated dressed states both analytical and
numerically. The final conclusions and outlook are discussed in Section IV.

##  II Model and electronic states in 1T′-MoS2

Let us first define our units and estimate the sought quantities related to a
real-life experiment. We define our units for the energy and momentum as the
Fermi energy in graphene for an experimentally accessible yet quite large
electron density $n_{e}^{(0)}=10^{11}\,cm^{-2}$. This corresponds to the Fermi
momentum $k_{F}^{(0)}=\sqrt{2\pi n^{(0)}}=7.92\cdot 10^{7}\,m^{-1}$ and the
Fermi energy $E_{F}^{(0)}=\hbar v_{F}k_{F}^{(0)}=8.32\cdot
10^{-21}\,J=52.02\,meV$. Thus, our unit of length is obtained as
$l^{(0)}=1/k_{F}^{(0)}=1.26\cdot 10^{-8}\,m\backsim 10\,nm$.

The logical starting point of building our model would be presenting the low-
energy Hamiltonian and the electronic states for 1T′-MoS2. In the vicinity of
two inequivalent $K$ and $K^{\prime}$ points, corresponding to $\xi=\pm 1$ the
main Hamiltonian

$\hat{\mathcal{H}}_{\xi=\pm 1}^{\,1T^{\prime}}({\bf
k})=V_{1}k_{x}\,\Gamma^{(2,0)}+\left\\{-\xi V_{-}\Gamma^{(0,0)}+\xi
V_{+}\,\Gamma^{(3,0)}+V_{2}\,\Gamma^{(1,1)}\right\\}\,k_{y}+\Delta_{0}\left[\xi\Gamma^{(2,0)}-ir_{E}\,\Gamma^{(2,0)}\cdot\Gamma^{(3,0)}\right]$
(1)

is linear which means the 1T′-MoS2 belong to Dirac materials. Here, the
components of anisotropic Fermi velocity $V_{=}0.286$, $V_{+}=0.721$,
$V_{1}=0.387$ and $V_{2}=0.046$ are given in terms of $v_{F}=10^{6}\,m/s$ for
graphene. The spin-orbit coupling gap is $\Delta_{0}=0.81\,E^{(0)}$.
$r_{E}=E_{z}/E_{c}$ is the relative value of the our-of-plane electric field
and $E_{c}$ is its critical value which closes the bandgap in 1T′-MoS2.

Following Ref. [Tan et al., 2021a], Hamiltonian (1) is written in terms of
$4\times 4$ gamma matrices $\Gamma^{(0,0)}=\tau_{0}\otimes\sigma_{0}$ is a
$4\times 4$ unit matrix, $\Gamma^{(1,1)}=\tau_{1}\otimes\sigma_{1}$,
$\Gamma^{(2,0)}=\tau_{2}\otimes\sigma_{0}$ and
$\Gamma^{(3,0)}=\tau_{3}\otimes\sigma_{0}$, where $\otimes$ means outer
product (or Kronecker product) and $\tau_{i}$ and $\sigma_{i}$ are regular
$2\times 2$ Pauli matrices acting in pseudospin and real-spin spaces,
correspondingly.

Figure 1: (Color online) Valley- and spin-polarized energy dispersions
$\epsilon_{\tau=\pm 1}({\bf k}\,|\,\xi=1,s)$ for 1T′-MoS2 with a finite spin-
orbit coupling gap $\Delta_{0}$ in the absence of dressing field. Panels
$(a)$, $(b)$ and $(c)$ correspond to different values of the external electric
field $r_{E}=0.0$, $0.5$ and $1.0$, as labeled. We have chosen valley index
$\xi=1$ for all the presented cases.

In the matrix form, we can rewrite Hamiltonian (1) as

$\hat{\mathcal{H}}_{\xi=\pm 1}^{\,1T^{\prime}}({\bf
k})=\left\\{\begin{array}[]{cccc}-\xi\,(V_{+}+V_{-})k_{y}&0&r_{E}\Delta_{0}-iV_{1}k_{x}&\xi\Delta+V_{2}k_{y}\\\
0&-\xi\,(V_{+}+V_{-})k_{y}&\xi\Delta+V_{2}k_{y}&r_{E}\Delta_{0}-iV_{1}k_{x}\\\
0&0&\xi\,(V_{+}-V_{-})k_{y}&0\\\
0&0&0&\xi\,(V_{+}-V_{-})k_{y}\end{array}\right\\}+h.c.\,,$ (2)

where $h.c.$ means a Hermitian conjugate matrix. The low-energy dispersions in
the vicinity of two Dirac points are immediately obtained as

$\displaystyle\epsilon_{\tau=\pm 1}({\bf
k}\,|\,\xi,s)=-\xi\,V_{-}k_{y}+\tau\,\mathbb{S}({\bf k}\,|\,\xi,s)\,,$ (3)
$\displaystyle\mathbb{S}({\bf
k}\,|\,\xi,s)=\sqrt{\left[(\xi-s\,r_{E})\Delta_{0}+V_{2}k_{y}\right]^{2}+(V_{+}k_{y})^{2}+(V_{1}k_{x})^{2}}\,,$
(4)

where $\tau=\pm 1$ described the electron/hole states related to the
conduction and valence bands and $s=\pm 1$ is the real spin index. We
immediately discern that the obtained energy spectrum is anisotropic and
tilted relatively to the $k_{y}$ axis ( $-\xi\,V_{-}k_{y}$ terms gives the
contributions equal and opposite sing around $k_{y}=0$). The energy bandgap is
indirect and closes if $\xi-r_{E}$. Similarly to silicene, $\xi>s\,r_{E}$ is a
topological insulator phase and $\xi<s\,r_{E}$ is a regular band insulator.

The corresponding wave functions are

$\Psi_{\xi=\pm 1}^{\,1T^{\prime}}({\bf
k})=\left\\{\begin{array}[]{c}s\,\mathcal{D}_{\,s,\,\xi}({\bf
k},\Delta_{0})\\\ -\mathcal{D}_{\,s,\,\xi}({\bf k},\Delta_{0})\\\ 1\\\
-s\end{array}\right\\}\,\texttt{e}^{i\xi k_{x}x}\,\texttt{e}^{ik_{y}y}\,,$ (5)

where

$\mathcal{D}_{\,s,\,\xi}({\bf k},\Delta_{0})=\left\\{(\xi-
sr_{E})\,\Delta_{0}+V_{2}k_{y}-isV_{1}k_{x}\right\\}^{-1}\,\left[\xi
V_{+}k_{y}-\tau\,\mathbb{S}({\bf k}\,|\,\xi,s)\right]\,$ (6)

and $\mathbb{S}({\bf k}\,|\,\xi,s)$ was defined in Eq. (3).

Figure 2: (Color online) $k_{y}$-dependence of the energy dispersions
$\epsilon_{\tau=\pm 1}(k_{x}=0,k_{y}\,|\,\xi,s)$ for non-irradiated 1T′-MoS2.
Each panel corresponds to a different values of the perpendicular electric
field $r_{E}=0.0$, $0.5$, $1.0$ (for which the smaller bandgap for spin $s=1$
is closed) and $r_{E}=1.2$, according to our labels. We take $\xi=1$ and
$k_{x}=0$ for all considered cases.  Figure 3: (Color online) Angular
dependence of the constant-energy $\varepsilon_{0}=2.0\,E^{(0)}$ cut of the
dispersions $\epsilon_{\tau=\pm 1}(k_{x}=0,k_{y}\,|\,\xi,s)=2.0\,E^{(0)}$ for
non-irradiated 1T′-MoS2. Each panel corresponds to a different values of the
perpendicular electric field $r_{E}=0.0$, $0.5$, $1.0$ and $1.2$. We have
chosen $\xi=1$ for all plots.

First, we should say that we encounter a highly unusual electronic states
because even though the Hamiltonian is linear to the components of the wave
vector $k_{x}$ and $k_{y}$, the Fermi velocities $V_{+}$ and $V_{-}$ are
present along the main diagonal of the Hamiltonian and the gap terms which
stay finite for ${\bf k}=0$ are located off-diagonal as well. This Hamiltonian
structure looks highly unusual and was not encountered in any of the
previously considered materials such as graphene, silicene, regular transition
metal dichalcogenides, 8-pmmn borophene or even phosphorenes.

The obtained energy dispersions, presented in Figs. 1 and 2 are tilted and
anisotropic. There is no mirror symmetry for $k_{y}\leftrightarrow-k_{y}$
while such symmetry is obviously present for $k_{x}$ component of the electron
momentum, which is very well seen from Fig. 4. In general, we have four non-
equivalent Fermi velocities $V_{+}$, $V_{-}$, $V_{1}$ and $V_{2}$ which affect
different properties of 1T′-MoS2 and contribute to its anisotropy.

The two subbands in each of the valence and conduction bands depend on the
spin index $s=\pm 1$. Thus, we observe two inequivalent subbands and two
different bandgaps unless there is no external perpendicular electrostatic
field and $r_{E}=0$. The smaller bandgap could be closed for $r_{E}=\xi=\pm 1$
and our lattice becomes semi metallic. For $|r_{E}|<1$, the band structure
represents a topological insulator, and for $|r_{E}|>1$ \- regular band
insulator. This situation is similar to silicene. All the dispersions and band
gaps depend on the valley Index $\xi=\pm 1$ which makes our electronic states
spin- and- valley- polarized.

##  III Electron dressed states

In this Section, we aim to calculating and analyze the new electronic states
originating as a result of the electron-photon interaction due to an applied
high-frequency external optical field. These modified electronic states are
often defined as electron dressed states and the field is referred to as
dressing field.

The vector potential of the dressing field enters effective Hamiltonian
through the canonical substitution of the wave vector ${\bf k}$ which enter
Hamiltonian (1) $k_{i}\longrightarrow k_{i}-eA_{i}^{(P)}(t)$, where the time-
dependent vector potential of applied field $A_{i}^{(P)}$ is mainly determined
by its polarization.

The driven-induced gauge fields substantially renormalize the band gaps and
the spin-orbit splitting of our investigated material. In order to solve the
obtained eigenvalue problem, we must rely on a perturbative Floquet-Magnus
expansion of the interaction Hamiltonian in powers of inverse frequency which
works particularly well for an off-resonant high-frequency irradiation.
Goldman and Dalibard (2014)

Most general type is elliptical polarization (clockwise):

$\mbox{\boldmath$A$}^{(E)}(t)=\left[\begin{array}[]{c}A^{(E)}_{x}(t)\\\
A^{(E)}_{y}(t)\end{array}\right]=\frac{E_{0}}{\omega}\left\\{\begin{array}[]{c}\cos\Theta^{(p)}\cos(\omega
t)-\beta\,\sin\Theta^{(p)}\sin(\omega t)\\\ \sin\Theta^{(p)}\cos(\omega
t)+\beta\,\cos\Theta^{(p)}\cos(\omega t)\end{array}\right\\}\,,$ (7)

where polarization angle is $\Theta^{(p)}$. The vector potential for a
linearly polarized light corresponds to $\beta=0$ in (7)

$\mbox{\boldmath$A$}^{(L)}(t)=\left[\begin{array}[]{c}A^{(L)}_{x}(t)\\\
A^{(L)}_{y}(t)\end{array}\right]=\frac{E_{0}}{\omega}\left\\{\begin{array}[]{c}\cos\Theta^{(p)}\\\
\sin\Theta^{(p)}\end{array}\right\\}\,\cos(\omega t)\ .$ (8)

Since the Hamiltonian in Eq.(1) is linear in $k_{x,y}$, in the presence of
$\mbox{\boldmath$A$}^{(L)}(t)$ it only acquires an additional interaction
term, yielding

$\hat{\mathcal{H}}_{1}(\mbox{\boldmath$k$}\,|\,\tau)\Longrightarrow\hat{\mathbb{H}}^{(L)}(\mbox{\boldmath$k$},t)=\hat{\mathcal{H}}_{1}(\mbox{\boldmath$k$}\,|\,\tau)+\hat{\mathcal{H}}_{A}^{(L)}(t)\
,$ (9)

where the $k$ independent interaction term is

$\displaystyle\hat{\mathcal{H}}_{A}^{(L)}(t)=c_{0}\,\cos(\omega t)\times$ (10)
$\displaystyle\times\left\\{\begin{array}[]{cccc}\xi\left[V_{+}+c_{-}\right]\sin\Theta^{(p)}&0&iV_{1}\cos\Theta^{(p)}&-V_{2}\sin\Theta^{(p)}\\\
0&\xi\left[V_{+}+V_{-}\right]\sin\Theta^{(p)}&-V_{2}\sin\Theta^{(p)}&iV_{1}\cos\Theta^{(p)}\\\
0&0&\xi\left[V_{-}-V_{+}\right]\sin\Theta^{(p)}&0\\\
0&0&0&\xi\left[V_{-}-V_{+}\right]\sin\Theta^{(p)}\end{array}\right\\}+h.c.\ $
(15)

Here we notice an important difference between tilted 1T′-MoS2 and all the
previously studied materials - the optical-coupling constants
$c_{0}^{(i)}=eV_{i}E_{0}/\omega$, corresponding to the different Fermi
velocities, are not the same for various elements of Hamiltonian (10). The
only time-dependent term found in matrix (10) is $\backsimeq\cos(\omega t)$
and is related to the linearly polarized type of the irradiation, not to a
specific material.

Figure 4: (Color online) Energy dispersions $\epsilon^{(L)}_{\tau=\pm
1}(\lambda_{0},{\bf k}\,|\,\xi,s)$ for the case of linearly polarized
irradiation applied to a 1T′-MoS2 lattice. The corresponding dispersions for
non-irradiated materials and the same values of the lattice parameters are
shown as dashed lines for comparison. Each of the panels corresponds to a
specific value of the electron-light coupling parameter $\lambda_{0}=0.3$ and
$0.5$ and the direction of the linear polarization of the dressing field
relative to the $x-$ axis $\Theta^{(P)}$ in accordance with our plot labels.
All the left panels $(a)$, $(c)$, $(e)$ and $(g)$ describe the energy
dispersions as a function of the $k_{y}$ component of the wavevector, while
the right ones $(b)$, $(d)$, $(f)$ and $(g)$ show their $k_{x}$-dependence.
The vertical electric field $r_{E}=0.5$ (when both energy bandgaps are open)
is selected for all plots.  Figure 5: (Color online) Energy dispersions
$\epsilon^{(L)}_{\tau=\pm 1}(\lambda_{0},{\bf k}\,|\,\xi,s)$ for the case of
linearly polarized irradiation applied to a 1T′-MoS2 lattice. The
corresponding dispersions for non-irradiated materials and the same values of
the lattice parameters are shown as dashed lines for comparison. Each of the
panels corresponds to a specific value of the electron-light coupling
parameter $\lambda_{0}=0.3$ and $0.5$ and the direction of the linear
polarization of the dressing field relative to the $x-$ axis $\Theta^{(P)}$ in
accordance with our plot labels. All the left panels $(a)$, $(c)$, $(e)$ and
$(g)$ describe the energy dispersions as a function of the $k_{y}$ component
of the wavevector, while the right ones $(b)$, $(d)$, $(f)$ and $(g)$ show
their $k_{x}$-dependence. The vertical electric field $r_{E}=1.0$ (when the
smaller energy bandgap is closed) is selected for all plots.

If the direction of the linear polarization is along the $x-$axis,
$\Theta^{(p)}=0$ and Hamiltonian (10) is reduced to

$\hat{\mathcal{H}}_{A}^{(L)}(t)=i\,c_{0}V_{1}\,\cos(\omega
t)\,\left[\begin{array}[]{cccc}0&0&1&0\\\ 0&0&0&1\\\ -1&0&0&0\\\
0&-1&0&0\end{array}\right]=-\frac{c_{0}\,V_{1}}{2}\cos(\omega
t)\,\Gamma^{(2,0)}\,,$ (16)

where we use the following notations: $c_{0}=v_{F}eE_{0}/\omega$ ($E_{0}$ is
the amplitude of the electric field of our dressing irradiation) and
$\lambda_{0}=c_{0}/(\hbar\omega)$.

The Floquet-Magnus perturbation approach which has been widely used to obtain
the electronic states in the presence of a off-resonance and high-frequency
dressing field Goldman and Dalibard (2014) requires that the only time-
dependent second term in Eq. (10) is rewritten as

$\hat{\mathcal{H}}_{A}^{(L)}(t)=\hat{\mathbb{O}}_{1}(c_{0},\tau)\,\texttt{e}^{i\omega
t}+\hat{\mathbb{O}}_{1}^{\dagger}(c_{0},\tau)\,\texttt{e}^{-i\omega t}\ ,$
(17)

in which matrix $\hat{\mathbb{O}}_{1}(c_{0},\tau)$ is free from time
dependence

$\displaystyle\hat{\mathbb{O}}_{1}(c_{0},\tau)=\frac{c_{0}}{2}\times$ (18)
$\displaystyle\times\left\\{\begin{array}[]{cccc}\xi\left[V_{+}+V_{-}\right]\sin\Theta^{(p)}&0&iV_{1}\cos\Theta^{(p)}&-V_{2}\sin\Theta^{(p)}\\\
0&\xi\left[V_{+}+V_{-}\right]\sin\Theta^{(p)}&-V_{2}\sin\Theta^{(p)}&iV_{1}\cos\Theta^{(p)}\\\
0&0&\xi\left[V_{-}-V_{+}\right]\sin\Theta^{(p)}&0\\\
0&0&0&\xi\left[V_{-}-V_{+}\right]\sin\Theta^{(p)}\end{array}\right\\}+h.c.=$
(23)
$\displaystyle=\xi\,c_{0}\,\sin\Theta^{(p)}\left[V_{+}\Gamma^{(3,0)}-V_{-}\Gamma^{(0,0)}\right]-c_{0}V_{1}\cos\Theta^{(p)}\,\Gamma^{(2,0)}-c_{0}V_{2}\sin\Theta^{(p)}\,\Gamma^{(1,1)}\,.$

Using the Floquet-Magnus series expansion over the powers of
$1/(\hbar\omega)$, the effective Hamiltonian for the dressed state is obtained
as

$\displaystyle\hat{\mathcal{H}}_{\text{eff}}^{\,(L)}(\mbox{\boldmath$k$}\,|\,\tau)$
$\displaystyle=$
$\displaystyle\hat{\mathcal{H}}_{1}(\mbox{\boldmath$k$}\,|\,\tau)+\frac{1}{\hbar\omega}\,\left[\,\hat{\mathbb{O}}_{1}(c_{0},\tau),\,\hat{\mathbb{O}}_{1}^{\dagger}(c_{0},\tau)\,\right]$
(24) $\displaystyle+$
$\displaystyle\frac{1}{2(\hbar\omega)^{2}}\left\\{\left[\left[\,\hat{\mathbb{O}}_{1}(c_{0},\tau),\,\hat{\mathcal{H}}_{\,1}(\mbox{\boldmath$k$}\,|\,\tau)\,\right],\,\hat{\mathbb{O}}_{1}^{\dagger}(c_{0},\tau)\,\right]\,\,+\,\,h.c.\right\\}\,\,+\cdots\
.$

As it should be the case with the perturbation expansions, the effective
Hamiltonian (24) is presented as the initial non-perturbed Hamiltonian (1) and
a small correction due to the dressing field.

The following term
$\left[\,\hat{\mathbb{O}}_{1}(c_{0},\tau),\,\hat{\mathbb{O}}_{1}^{\dagger}(c_{0},\tau)\right]$
is $k$-independent and, therefore, largely defines the irradiation-induced
bandgap of our dressed states for $k=0$. However, for 1T′ MoS2 the consequent
expansion terms can also affect the dispersions at ${\bf k}=0$. In the case of
linearly polarized light, it is obviously zero since
$\hat{\mathbb{O}}_{1}(c_{0},\tau)$ in Eq. (18) is Hermitian.

The remaining term in Eq. (24) is written as
$\hat{\mathbb{T}}_{2}(\lambda_{0}\,|\,k,\theta_{\bf k})$.For $\Theta^{(p)}=0$,
the 2nd order correction to interaction Hamiltonian is reduced to

$\displaystyle\hat{\mathbb{T}}_{2}(\lambda_{0}\,|\,k,\theta_{\bf
k})=\frac{\lambda_{1}^{2}}{4}\,\left\\{\begin{array}[]{cccc}\xi\,V_{+}k_{y}&0&-r_{E}\Delta_{0}&-(\xi\Delta_{0}+V_{-}k_{y})\\\
0&4\xi\,V_{+}k_{y}&-(\xi\Delta_{0}+V_{-}k_{y})&-r_{E}\Delta_{0}\\\
0&0&-\xi\,V_{+}k_{y}&0\\\ 0&0&0&-\xi\,V_{+}k_{y}\end{array}\right\\}+h.c.=$
(29)
$\displaystyle=\frac{\lambda_{1}^{2}}{4}\,\left\\{\xi\,V_{+}k_{y}\,\Gamma^{(3,0)}-(\xi\Delta_{0}+V_{-}k_{y})\,\Gamma^{(1,1)}+r_{E}\Delta_{0}\,\Gamma^{(2,0)}\right\\}\
,$

The energy dispersions under linearly polarized dressing field could be
obtained analytically but the expressions are way too long and complicated,
even for $\Theta^{(p)}=0$. Therefore, we will focus on the numerical
investigation of the obtained dressed states.

In the simplest case of $r_{E}=\xi$ (we also choose $\xi=1$ for the sake of
clarity), both non-irradiate and irradiated energy dispersions have no bandgap
and the latter is expressed as $(\lambda_{0}=v_{F}eE_{0}/(\hbar\omega^{2})\ll
1)$

$\epsilon^{(L)}_{\tau=\pm 1}(\lambda_{0},{\bf
k}\,|\,\xi=1,s)=-\xi\,V_{-}k_{y}\pm\sqrt{(V_{1}k_{x})^{2}+\left[1-2\lambda_{0}^{2}V_{1}^{2}\right]^{2}\,\left(V_{1}^{2}+V_{+}^{2}\right)\,k_{y}^{2}\,.}$
(30)

Please keep in mind that each of the Fermi velocities $V_{1}$, $V_{2}$,
$V_{+}$ and $V_{-}$ is dimensionless and is given as a ration to the Fermi
velocity in graphene $v_{F}=10^{6}\,m/s$.

In comparison, a linearly polarized dressing field with the same coupling
constant $\lambda_{0}$ applied to graphene and a dice lattice leads to
dispersions $\pm\hbar v_{F}\,\sqrt{k_{x}^{2}+a(\lambda_{0})^{2}k_{y}^{2}}$,
where $a(\lambda_{0})=1-\lambda_{0}^{2}/2$ for graphene and
$a(\lambda_{0})=1-\lambda_{0}^{2}/4$ for a dice lattice, correspondingly.

Figure 6: (Color online) Energy dispersions $\epsilon^{(C)}_{\tau=\pm
1}(\lambda_{0},{\bf k}\,|\,\xi,s)$ for the case of circularly polarized
dressing field applied to a 1T′-MoS2 lattice. The corresponding dispersions
for non-irradiated materials and the same values of the lattice parameters are
shown as dashed lines for comparison. Each of the panels corresponds to a
specific value of the electron-light coupling parameter $\lambda_{0}=0.5$ and
$0.7$, as labeled. The left panels $(a)$ and $(c)$, describe the energy
dispersions as a function of the $k_{y}$ component of the wavevector, while
the right ones $(b)$ and $(d)$ show their $k_{x}$-dependence. The vertical
electric field $r_{E}=0.5$ (both bandgaps are open) and $\xi=1$ are selected
for all plots. Figure 7: (Color online) Energy dispersions
$\epsilon^{(C)}_{\tau=\pm 1}(\lambda_{0},{\bf k}\,|\,\xi,s)$ for the case of
circularly polarized dressing field applied to a 1T′-MoS2 lattice. The
corresponding dispersions for non-irradiated materials and the same values of
the lattice parameters are shown as dashed lines for comparison. Each of the
panels corresponds to a specific value of the electron-light coupling
parameter $\lambda_{0}=0.5$ and $0.7$, as labeled. The left panels $(a)$ and
$(c)$, describe the energy dispersions as a function of the $k_{y}$ component
of the wavevector, while the right ones $(b)$ and $(d)$ show their
$k_{x}$-dependence. The vertical electric field $r_{E}=0.5$ (both bandgaps are
open) and $\xi=-1$ are selected for all plots. Figure 8: (Color online) Energy
dispersions $\epsilon^{(C)}_{\tau=\pm 1}(\lambda_{0},{\bf k}\,|\,\xi,s)$ for
the case of circularly polarized dressing field applied to a 1T′-MoS2 lattice.
The corresponding dispersions for non-irradiated materials and the same values
of the lattice parameters are shown as dashed lines for comparison. Each of
the panels corresponds to a specific value of the electron-light coupling
parameter $\lambda_{0}=0.5$ and $0.7$, as labeled. The left panels $(a)$ and
$(c)$, describe the energy dispersions as a function of the $k_{y}$ component
of the wavevector, while the right ones $(b)$ and $(d)$ show their
$k_{x}$-dependence. The vertical electric field $r_{E}=1.0$ (semi-metallic
dispersions with no bandgap) and $\xi=1$ are selected for all plots. Figure 9:
(Color online) Energy dispersions $\epsilon^{(C)}_{\tau=\pm
1}(\lambda_{0},{\bf k}\,|\,\xi,s)$ for the case of circularly polarized
dressing field applied to a 1T′-MoS2 lattice. The corresponding dispersions
for non-irradiated materials and the same values of the lattice parameters are
shown as dashed lines for comparison. Each of the panels corresponds to a
specific value of the electron-light coupling parameter $\lambda_{0}=0.5$ and
$0.7$, as labeled. The left panels $(a)$ and $(c)$, describe the energy
dispersions as a function of the $k_{y}$ component of the wavevector, while
the right ones $(b)$ and $(d)$ show their $k_{x}$-dependence. The vertical
electric field $r_{E}=1.0$ (semi-metallic dispersions with no bandgap) and
$\xi=1$ are selected for all plots.

Linearly polarized irradiation is generally known to induce or modify the
existing anisotropy in two-dimensional lattices. If the initial bandgap is
zero, it stays zero in graphene or an $\alpha-\mathcal{T}_{3}$. Since we are
dealing with initially anisotropic lattice, the direction of the light
polarization matters, and the obtained energy dispersions are essentially
different for the different directions $\Theta^{(p)}$ of the linear
polarization of the dressing field.

We consider the states with $r_{E}=0.5$ and a finite bandgap and present our
results in Fig. 4, as well as initially gapless bandstructure shown in Fig. 5.
Here, we compare how the energy dispersions are modified in the presence of a
dressing field (solid lines) compared to the non-irradiated states (dashed
lines).

We clearly see that all the important characteristics of the energy
dispersions are affected by the linearly polarized light. Unlike graphene,
both the direct and indirect energy bandgaps are decreased when the off-
resonant dressing field is applied.

Tilting, anisotropy and $k_{x}\leftrightarrow-k_{x}$ mirror symmetry of the
electronic states are affected are affected by the linearly polarized
irradiation. The locations of the critical points (minima and maxima) of each
subband are shifted as well. We also note that energy subbands corresponding
to the different values of spin index values are affected differently by
linearly polarized light.

For $\Theta^{(p)}$ the dispersions in both $x-$ and $y-$ in-plane directions
are modified and we see a lot richer physical picture which includes breaking
the $k_{x}\leftrightarrow-k_{x}$ mirror symmetry of the dressed states. We
also see that for the semi-metallic states with $r_{E}=\xi$ for which one of
the bandgaps is closed, it also remains zero under linearly polarized light.
The slopes of the dispersions, however, are affected differently for the $x-$
and $y-$ directions revealing field-induced anisotropy, just as we observed
for graphene. The $k_{x}\leftrightarrow-k_{x}$ mirror symmetry remains
unaffected by such field only for $\Theta^{(p)}=0$.

###  III.1 Circularly polarized light

The opposite limit $\beta=1$ is related to circular polarization

$\mbox{\boldmath$A$}^{(C)}(t)=\left[\begin{array}[]{c}A^{(C)}_{x}(t)\\\
A^{(C)}_{y}(t)\end{array}\right]=\frac{E_{0}}{\omega}\left\\{\begin{array}[]{c}\cos\left[\Theta^{(p)}+\omega
t\right]\\\ \sin\left[\Theta^{(p)}+\omega t\right]\end{array}\right\\}\ .$
(31)

Since we are looking for stationary and time-independent states, the initial
phase $\Theta^{(p)}$ for circularly polarized field could be removed with no
loss of generality.

Thus, the vector potential for circularly polarized irradiation is simplified
as

$\mbox{\boldmath$A$}^{(C)}(t)=\left[\begin{array}[]{c}A^{(C)}_{x}(t)\\\
A^{(C)}_{y}(t)\end{array}\right]=\frac{E_{0}}{\omega}\left\\{\begin{array}[]{c}\cos\left(\omega
t\right)\\\ \sin\left(\omega t\right)\end{array}\right\\}\ .$ (32)

The interaction Hamiltonian for circularly polarized light with vector
potential (32) becomes

$\displaystyle\hat{\mathcal{H}}_{A}^{(C)}(t)=$ (33)
$\displaystyle=c_{0}\,\left\\{\begin{array}[]{cccc}\xi\left[V_{+}+V_{-}\right]\sin(\omega
t)&0&iV_{1}\cos(\omega t)&-V_{2}\sin(\omega t)\\\
0&\xi\left[V_{+}+V_{-}\right]\sin(\omega t)&-V_{2}\sin(\omega
t)&iV_{1}\cos(\omega t)\\\ 0&0&\xi\left[V_{-}-V_{+}\right]\sin(\omega t)&0\\\
0&0&0&\xi\left[V_{-}-V_{+}\right]\sin(\omega t)\end{array}\right\\}+h.c.=$
(38) $\displaystyle=c_{0}\left\\{\xi\,\sin(\omega
t)\left[V_{-}\,\Gamma^{(0,0)}+V_{+}\,\Gamma^{(3,0)}\right]-V_{2}\,\sin(\omega
t)\,\Gamma^{(1,1)}-V_{1}\,\cos(\omega t)\,\Gamma^{(2,0)}\right\\}\,.$

Using Eq. (33), we can immediately obtain the time-independent interaction
matrix $\hat{\mathbb{O}}_{1}(c_{0},\tau)$ as

$\displaystyle\hat{\mathbb{O}}_{1}(c_{0},\xi)=c_{0}\left\\{-i\xi\,\left[V_{-}\,\Gamma^{(0,0)}+V_{+}\,\Gamma^{(3,0)}\right]+iV_{2}\,\Gamma^{(1,1)}-V_{1}\,\Gamma^{(2,0)}\right\\}=$
(39)
$\displaystyle=\frac{-i\,c_{0}}{2}\,\left\\{\begin{array}[]{cccc}\xi\,\left[V_{+}+V_{-}\right]&0&-V_{1}&-V_{2}\\\
0&\xi\,\left[V_{+}+V_{-}\right]&-V_{2}&-V_{1}\\\
V_{1}&-V_{2}&\xi\,\left[V_{-}-V_{+}\right]&0\\\
-V_{2}&V_{1}&0&\xi\,\left[V_{-}-V_{+}\right]\end{array}\right\\}\,,$ (44)

Now we can easily obtain the 1st order correction
$\frac{1}{\hbar\omega}\,\left[\,\hat{\mathbb{O}}_{1}(c_{0},\tau),\,\hat{\mathbb{O}}_{1}^{\dagger}(c_{0},\tau)\,\right]$
of our series expansion

$\frac{1}{\hbar\omega}\,\left[\,\hat{\mathbb{O}}_{1}(c_{0},\tau),\,\hat{\mathbb{O}}_{1}^{\dagger}(c_{0},\tau)\,\right]=\frac{c_{0}^{2}\,V_{1}}{4\,\hbar\omega}\,\left\\{\begin{array}[]{cccc}0&V_{2}&\xi\,V_{+}&0\\\
V_{2}&0&0&\xi\,V_{+}\\\ \xi\,V_{+}&0&0&-V_{2}\\\
0&\xi\,V_{+}&-V_{2}&0\end{array}\right\\}\,.$ (45)

The resulting energy dispersions are obtained as

$\displaystyle\epsilon_{\tau=\pm 1}({\bf
k}\,|\,\xi,s)=-\xi\,V_{-}k_{y}+\tau\times$ (46)
$\displaystyle\times\sqrt{\left[\,\Delta_{\xi,s}(\lambda_{0})^{2}+V_{2}k_{y}\right]^{2}+(V_{+}k_{y})^{2}+(V_{1}k_{x})^{2}}\,,$

For $k_{x}=k_{y}=0$, the energy dispersions reveal spin- and valley-dependent
bandgaps

$\Delta_{\xi,s}(\lambda_{0})=\pm\left\\{[(\xi+s\,r_{E})\Delta_{0}]^{2}+2\lambda_{0}^{2}V_{1}V_{+}\Delta_{0}(1+s\,\xi\,r_{E})+\lambda_{0}^{4}V_{1}^{2}(V_{2}^{2}+V_{+}^{2})\,\right\\}^{1/2}\,.$
(47)

The finite-$k$ terms of dispersions (46) are not affected by the circularly
polarized irradiation since only first-term correction term was taken into
account to obtain this simplified analytical expressions.

Finally, we consider the modification of the energy dispersions in the
presence of circularly polarized light taking into account the first- and the
second-order correction terms and present our numerical results in Figs. 6 \-
9.

The most interesting observation is definitely that we see a tremendous
difference between the two different valleys $\xi=1$ and $\xi=-1$, shown
separately in each of the figures. For $\xi=1$, the smaller bandgaps remain
nearly unaffected, while the larger ones demonstrate a significant increase.
Surprisingly, for $\xi=-1$ both subbands and both bandgaps show almost the
same increase.

The circularly polarized irradiation normally leads to opening a bandgap in
most known Dirac materials. If non-irradiated 1T′-MoS2 lattice has zero gap,
it becomes finite when a dressing field with circular polarization is applied.

Not only the states with no dressing field, but the irradiation-induced energy
gaps demonstrate a strong dependence on the spin and valley indices and,
therefore, are spin- and valley-polarized. Another interesting feature of the
dressed states under circular polarized radiation is that the initially bigger
bandgaps are affected stronger and become even larger in the presence of a
dressing field. This property shows a stark difference with graphene and
silicene in which the gapless states are mostly modified by circularly
polarized irradiation. In 1T′-MoS2, the slopes or Fermi velocities are also
strongly affected by this type of dressing field.

##  IV Summary and remarks

Our present work represents an attempt to construct a new type of low-energy
dispersions and electronic states in recently discovered distorted tetragonal
molybdenum disulfide 1T′-MoS2 by applying Floquet engineering.

Floquet theory and, broadly speaking, Floquet engineering represent one of the
most promising and currently studied directions in the present-day material
science and quantum optics. A fundamental interest to this topic which
connects various aspects of condensed-matter physics stems from the fact that
one can foresee some crucial and unusual electronic properties of a novel
material before it is actually tested in experiment.

The resulting electron dressed state represents a single quasiparticle which
combines the crucial features of both irradiated matter and the applied
dressing field. The properties of these states are mainly determined by the
polarization of the applied irradiation.

We have performed a rigorous theoretical and numerical investigation into the
electron dressed states in 1T′-MoS2 in the presence of external radiation with
various polarizations, which includes the derivation of closed-form analytical
expressions for the electron energy band structure driven by external
irradiation in the terahertz regime.

We have found that there is no clear distinction between the effects of
linearly and circularly polarized light for 1T′-MoS2, the obtained physical
picture is a lot richer and more involved than what we have seen for any
previously considered Dirac materials. Specifically, linearly polarized light
which has been mostly noticed to create or modify the in-plane anisotropy of a
two-dimensional lattice, in our case also substantially increases the larger
bandgap which has not been previously observed.

One of our most important findings is that not only the initial states in
1T′-MoS2 but also their modification by the dressing field substantially
depends on the valley index $\xi=\pm 1$. This includes changing both direct
and indirect bandgaps, tilting and Fermi velocities in different directions.
Therefore, the obtained dressed states also show significant degree of valley
polarization.

The energy dispersion have been also obtained analytically for both linear and
circular polarizations, but in general their expressions are so complicated
that we had to choose only very specific cases - fixed valley index $\xi=1$
and zero band gap for linearly polarized light and only first-order
approximation for circular polarized irradiation - to provide our analytical
results. However, these approximated equations give us a pretty good
understanding of how the electron states are modified by the irradiation. The
linearly polarized dressing field modified some of the Fermi velocities and,
therefore, leads to a change of the existing anisotropy and tilting, while
applying circularly polarized irradiation causes mostly the renormalization of
the existing bandgaps.

Our results for the modified electron dispersions could be later tested by
transmission electron microscopy, scanning tunneling spectroscopy, angle-
resolved photoemission spectroscopy (ARPES). This work is expected to become a
noticeable contribution to an unprecedented research effort on Floquet
engineering of all innovative two-dimensional Dirac materials.

## References

  * Kibis et al. (2020) O. Kibis, M. Boev, V. Kovalev, and I. Shelykh, Physical Review B 102, 035301 (2020).
  * Kibis et al. (2018) O. Kibis, K. Dini, I. Iorsh, and I. Shelykh, Semiconductors 52, 523 (2018).
  * Cheng et al. (2019) Q. Cheng, Y. Pan, H. Wang, C. Zhang, D. Yu, A. Gover, H. Zhang, T. Li, L. Zhou, and S. Zhu, Physical review letters 122, 173901 (2019).
  * Weitenberg and Simonet (2021) C. Weitenberg and J. Simonet, Nature Physics 17, 1342 (2021).
  * Wang et al. (2018) B. Wang, F. N. Ünal, and A. Eckardt, Physical Review Letters 120, 243602 (2018).
  * Dehghani et al. (2014) H. Dehghani, T. Oka, and A. Mitra, Physical Review B 90, 195429 (2014).
  * Nakagawa et al. (2020) M. Nakagawa, R.-J. Slager, S. Higashikawa, and T. Oka, Physical Review B 101, 075108 (2020).
  * Iurov et al. (2020a) A. Iurov, L. Zhemchuzhna, P. Fekete, G. Gumbs, and D. Huang, Physical Review Research 2, 043245 (2020a).
  * Pereira et al. (2006) V. M. Pereira, F. Guinea, J. L. Dos Santos, N. Peres, and A. C. Neto, Physical review letters 96, 036801 (2006).
  * Gumbs et al. (2014a) G. Gumbs, A. Balassis, A. Iurov, and P. Fekete, The Scientific World Journal 2014 (2014a).
  * Schnez et al. (2010) S. Schnez, J. Güttinger, M. Huefner, C. Stampfer, K. Ensslin, and T. Ihn, Physical Review B 82, 165445 (2010).
  * Castro et al. (2008) E. V. Castro, N. Peres, J. L. Dos Santos, A. C. Neto, and F. Guinea, Physical review letters 100, 026802 (2008).
  * Gumbs et al. (2014b) G. Gumbs, A. Iurov, D. Huang, and L. Zhemchuzhna, Physical Review B 89, 241407 (2014b).
  * Kibis (2010) O. Kibis, Physical Review B 81, 165433 (2010).
  * Iurov et al. (2011) A. Iurov, G. Gumbs, O. Roslyak, and D. Huang, Journal of Physics: Condensed Matter 24, 015303 (2011).
  * Ibarra-Sierra et al. (2022) V. Ibarra-Sierra, J. Sandoval-Santana, A. Kunold, S. A. Herrera, and G. G. Naumis, Journal of Physics: Materials 5, 014002 (2022).
  * Iurov et al. (2022) A. Iurov, L. Zhemchuzhna, G. Gumbs, D. Huang, and P. Fekete, Physical Review B 105, 115309 (2022).
  * Castro et al. (2022) A. Castro, U. De Giovannini, S. A. Sato, H. Hübener, and A. Rubio, arXiv preprint arXiv:2203.03387 (2022).
  * Oka and Kitamura (2018) T. Oka and S. Kitamura, arXiv preprint arXiv:1804.03212 (2018).
  * Kibis et al. (2021) O. Kibis, M. Boev, and V. Kovalev, in _Journal of Physics: Conference Series_ (IOP Publishing, 2021), vol. 2015, p. 012063.
  * Hsu and Reichl (2006) H. Hsu and L. Reichl, Physical Review B 74, 115406 (2006).
  * Sentef et al. (2015) M. Sentef, M. Claassen, A. Kemper, B. Moritz, T. Oka, J. Freericks, and T. Devereaux, Nature communications 6, 1 (2015).
  * Calvo et al. (2011) H. L. Calvo, H. M. Pastawski, S. Roche, and L. E. F. Torres, Applied Physics Letters 98, 232103 (2011).
  * Vogl et al. (2021) M. Vogl, M. Rodriguez-Vega, B. Flebus, A. H. MacDonald, and G. A. Fiete, Physical Review B 103, 014310 (2021).
  * Dal Lago et al. (2017) V. Dal Lago, E. S. Morell, and L. F. Torres, Physical Review B 96, 235409 (2017).
  * Iurov et al. (2017a) A. Iurov, G. Gumbs, and D. Huang, Journal of Modern Optics 64, 913 (2017a).
  * Kibis et al. (2017) O. Kibis, K. Dini, I. Iorsh, and I. Shelykh, Physical Review B 95, 125401 (2017).
  * Iurov et al. (2019) A. Iurov, G. Gumbs, and D. Huang, Physical Review B 99, 205135 (2019).
  * Dey and Ghosh (2019) B. Dey and T. K. Ghosh, Physical Review B 99, 205429 (2019).
  * Dey and Ghosh (2018) B. Dey and T. K. Ghosh, Physical Review B 98, 075422 (2018).
  * Weekes et al. (2021) N. Weekes, A. Iurov, L. Zhemchuzhna, G. Gumbs, and D. Huang, Physical Review B 103, 165429 (2021).
  * Tamang et al. (2021) L. Tamang, T. Nag, and T. Biswas, Physical Review B 104, 174308 (2021).
  * Wang et al. (2021) W. Wang, X. Lüu, and H. Xie, Chinese Physics B 30, 066701 (2021).
  * Tahir et al. (2016) M. Tahir, Q. Zhang, and U. Schwingenschlögl, Scientific reports 6, 1 (2016).
  * Iurov et al. (2021) A. Iurov, L. Zhemchuzhna, G. Gumbs, D. Huang, P. Fekete, F. Anwar, D. Dahal, and N. Weekes, Scientific reports 11, 1 (2021).
  * Iurov et al. (2017b) A. Iurov, L. Zhemchuzhna, G. Gumbs, and D. Huang, Journal of Applied Physics 122, 124301 (2017b).
  * Calvo et al. (2020) H. L. Calvo, J. E. B. Vargas, and L. E. F. Torres, Physical Review B 101, 075424 (2020).
  * Calvo et al. (2015) H. L. Calvo, L. F. Torres, P. M. Perez-Piskunow, C. A. Balseiro, and G. Usaj, Physical Review B 91, 241404 (2015).
  * Iurov et al. (2013) A. Iurov, G. Gumbs, O. Roslyak, and D. Huang, Journal of Physics: Condensed Matter 25, 135502 (2013).
  * Mojarro et al. (2020) M. Mojarro, V. Ibarra-Sierra, J. Sandoval-Santana, R. Carrillo-Bastos, and G. G. Naumis, Physical Review B 102, 165301 (2020).
  * Kristinsson et al. (2016) K. Kristinsson, O. V. Kibis, S. Morina, and I. A. Shelykh, Scientific reports 6, 1 (2016).
  * Iurov et al. (2020b) A. Iurov, L. Zhemchuzhna, D. Dahal, G. Gumbs, and D. Huang, Physical Review B 101, 035129 (2020b).
  * Iorsh et al. (2022) I. Iorsh, D. Zezyulin, S. Kolodny, R. Sinitskiy, and O. Kibis, Physical Review B 105, 165414 (2022).
  * Dey et al. (2020) B. Dey, P. Kapri, O. Pal, and T. K. Ghosh, arXiv preprint arXiv:2003.07143 (2020).
  * Dal Lago et al. (2015) V. Dal Lago, M. Atala, and L. F. Torres, Physical Review A 92, 023624 (2015).
  * Iorsh et al. (2017) I. Iorsh, K. Dini, O. Kibis, and I. Shelykh, Physical Review B 96, 155432 (2017).
  * Tang et al. (2017) S. Tang, C. Zhang, D. Wong, Z. Pedramrazi, H.-Z. Tsai, C. Jia, B. Moritz, M. Claassen, H. Ryu, S. Kahn, et al., Nature Physics 13, 683 (2017).
  * Zheng et al. (2016) F. Zheng, C. Cai, S. Ge, X. Zhang, X. Liu, H. Lu, Y. Zhang, J. Qiu, T. Taniguchi, K. Watanabe, et al., Advanced Materials 28, 4845 (2016).
  * Fei et al. (2017) Z. Fei, T. Palomaki, S. Wu, W. Zhao, X. Cai, B. Sun, P. Nguyen, J. Finney, X. Xu, and D. H. Cobden, Nature Physics 13, 677 (2017).
  * Liu et al. (2018) L. Liu, J. Wu, L. Wu, M. Ye, X. Liu, Q. Wang, S. Hou, P. Lu, L. Sun, J. Zheng, et al., Nature materials 17, 1108 (2018).
  * Mehmood et al. (2021) F. Mehmood, R. Pachter, T. C. Back, J. J. Boeckl, R. T. Busch, and P. R. Stevenson, Applied Optics 60, G232 (2021).
  * Gomes and Ramos (2021) Y. Gomes and R. O. Ramos, Physical Review B 104, 245111 (2021).
  * Tan et al. (2021a) C.-Y. Tan, C.-X. Yan, Y.-H. Zhao, H. Guo, H.-R. Chang, et al., Physical Review B 103, 125425 (2021a).
  * Tabert and Nicol (2013) C. J. Tabert and E. J. Nicol, Physical Review Letters 110, 197402 (2013).
  * Kara et al. (2012) A. Kara, H. Enriquez, A. P. Seitsonen, L. L. Y. Voon, S. Vizzini, B. Aufray, and H. Oughaddou, Surface science reports 67, 1 (2012).
  * Zhao et al. (2016) J. Zhao, H. Liu, Z. Yu, R. Quhe, S. Zhou, Y. Wang, C. C. Liu, H. Zhong, N. Han, J. Lu, et al., Progress in Materials Science 83, 24 (2016).
  * Drummond et al. (2012) N. Drummond, V. Zolyomi, and V. Fal’Ko, Physical Review B 85, 075423 (2012).
  * Ezawa (2012) M. Ezawa, Physical review letters 109, 055502 (2012).
  * Mannix et al. (2015) A. J. Mannix, X.-F. Zhou, B. Kiraly, J. D. Wood, D. Alducin, B. D. Myers, X. Liu, B. L. Fisher, U. Santiago, J. R. Guest, et al., Science 350, 1513 (2015).
  * Lopez-Bezanilla and Littlewood (2016) A. Lopez-Bezanilla and P. B. Littlewood, Physical Review B 93, 241405 (2016).
  * Islam (2018) S. F. Islam, Journal of Physics: Condensed Matter 30, 275301 (2018).
  * Paul et al. (2019) G. C. Paul, S. F. Islam, and A. Saha, Physical Review B 99, 155418 (2019).
  * Tan et al. (2021b) C.-Y. Tan, J.-T. Hou, C.-X. Yan, H. Guo, and H.-R. Chang, arXiv preprint arXiv:2112.09392 (2021b).
  * Park et al. (2022) S. H. Park, M. Sammon, E. Mele, and T. Low, arXiv preprint arXiv:2206.04222 (2022).
  * Mojarro et al. (2021) M. Mojarro, R. Carrillo-Bastos, and J. A. Maytorena, Physical Review B 103, 165415 (2021).
  * Goldman and Dalibard (2014) N. Goldman and J. Dalibard, Physical review X 4, 031027 (2014).