††thanks: Both authors contributed equally to this work.††thanks: Both authors
contributed equally to this work.

# Quasiclassical boundary conditions for spin-orbit coupled interfaces with
spin-charge conversion

Jacob Linder  Center for Quantum Spintronics, Department of Physics, Norwegian
University of Science and Technology, NO-7491 Trondheim, Norway  Morten
Amundsen  Nordita, KTH Royal Institute of Technology and Stockholm University,
Hannes Alfvens vag 12, SE-106 91 Stockholm, Sweden

###### Abstract

The quasiclassical theory of superconductivity provides a methodology to study
emergent phenomena in hybrid structures comprised of superconductors
interfaced with other materials. A key component in this theory is the
boundary condition that the Green functions describing the materials must
satisfy. Recently, progress has been made toward formulating such a boundary
condition for interfaces with spin-orbit coupling, the latter playing an
important role for several phenomena in spintronics. Here, we derive a
boundary condition for spin-orbit coupled interfaces that includes gradient
terms which enables the description of spin-Hall like effects with
superconductors due to such interfaces. As an example, we show that the
boundary conditions predict that a supercurrent flowing through a
superconductor that is coupled to a normal metal via a spin-orbit interface
can induce a non-local magnetization in the normal metal.

##  I Introduction

Superconductors [1] carry charge currents without resistance in their natural
state, but can also be made to transport spin currents without resistance.
This happens when the spinless Cooper pairs in a superconductor are
transformed to so-called triplet Cooper pairs [2] that carry spin. Such a
transformation has been predicted [3, 4] to take place in two ways: either by
placing the superconductor in contact with a ferromagnetic material or by
placing it in contact with spin-orbit coupled materials or impurities [5].
Besides making the superconductivity spin-polarized [6, 7], such hybrid
structures also give rise to odd-frequency superconductivity [8, 9] where
electrons are non-locally correlated in time. Although transport of charge is
already dissipationless in superconductors, it is both of fundamental and
practical interest to study dissipationless transport of spin via
superconductors [10]. Indeed, lossless spin transport can offer new types of
functionality in solid-state devices when combining superconductors with
magnetic materials that charge-based transport is not capable of. Examples of
this includes magnetization switching and domain wall/skyrmion motion with
very low dissipation of energy [11]. In addition, triplet superconductivity
can provide novel means of probing quantum materials [12], as well as giving
rise to functional material properties that are unique to the superconducting
state, without any counterpart in conventional spintronics—such as quantum
phase batteries [13].

Several experiments have confirmed theoretical predictions in the case of
hybrid structures comprised of superconductors and ferromagnets [14, 15, 16,
17, 18, 19, 20]. On the other hand, experimental studies of triplet
superconductivity in hybrid structures of superconductors and spin-orbit
coupled materials have started to appear mostly in recent years [26, 27, 28,
29, 30]. There exists several predictions [21, 22, 23, 24, 25] on emergent
phenomena in structures combining superconducting and spin-orbit coupled
materials which remain to be tested experimentally, where the main theme is
that spin-orbit interactions exists in a material due to the lack of inversion
symmetry or due to impurities.

Spin-orbit coupling can be amplified at interfaces involving heavy metals such
as Pt or W. Such interfaces are known to play an important role in
conventional spintronics [31]. Generally speaking, the Rashba-effect at
interfaces, even without inserted heavy metals, is expected to be present
simply because structural inversion symmetry is broken [32]. This makes the
description of spin-orbit coupled interfaces of general interest in any type
of hybrid structures comprised of two or more materials.

When it comes to describing quantum phenomena emerging at superconducting
interfaces, the quasiclassical theory of superconductivity [33, 34, 35] has
proven useful and to a large degree consistent with experimental results. This
theory is capable of treating non-equilibrium phenomena both in the ballistic
and diffusive limit by obtaining the quasiclassical Green function $\check{g}$
in Keldysh-space. To achieve this aim, it is necessary to supplement the
theory with boundary conditions at interfaces. For the case of interfaces with
significant spin-orbit coupling, we have recently derived such boundary
conditions in the diffusive limit [36]. These boundary conditions for
$\check{g}$ correctly predict the appearance of triplet superconductivity in
superconductors with spin-orbit coupled interfaces, but fail to capture the
effect of motion parallel to the interface. While these contributions often
cancel in the diffusive limit due to the momentum averaging brought on by
frequent impurity scattering, an important special case when it does not is
when a supercurrent is flowing parallel to the interface. In this case,
gradient terms in the boundary conditions are required in order to capture the
possible existence of spin-Hall like phenomena involving supercurrent flow and
spin accumulation. Whereas such effects are of high interest experimentally,
an accompanying theoretical description requires a derivation of boundary
conditions that include such gradient terms.

Here, we derive a boundary condition for spin-orbit coupled interfaces in the
quasiclassical theory of superconductivity that includes precisely such
gradient terms. This enables the description of spin-Hall like effects with
superconductors due to interfaces with e.g. Rashba-like spin-orbit coupling.
After deriving these boundary conditions, we apply them to predict that a
supercurrent flowing through a superconductor that is coupled to a normal
metal via a spin-orbit interface can induce a non-local magnetization in the
normal metal. This constitutes a superconducting equivalent of the
conventional, resistive spin-Hall effect where a charge current induces a
transverse spin accumulation. The boundary conditions derived here enables the
exploration of quantum phenomena that emerge in superconductors due to spin-
orbit coupled interfaces, both in and out of equilibrium.

##  II Derivation of boundary conditions

We split the derivation of the boundary conditions into two parts: one part
that takes into account the effect of tunneling through the interface and one
part that takes into account reflection at the interface.

###  II.1 Terms related to tunneling through the interface

We consider a Hamiltonian of the form

$\displaystyle H=H_{L}+H_{R}+H_{T},$ (1)

with

$\displaystyle H_{L}=$
$\displaystyle\sum_{\bm{p}\bm{p}^{\prime}ss^{\prime}}H^{L}_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})a^{\dagger}_{\bm{p}s}a_{\bm{p}^{\prime}s^{\prime}},$
(2) $\displaystyle H_{R}=$
$\displaystyle\sum_{\bm{p}\bm{p}^{\prime}ss^{\prime}}H^{R}_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})b^{\dagger}_{\bm{p}s}b_{\bm{p}^{\prime}s^{\prime}},$
(3) $\displaystyle H_{T}=$
$\displaystyle\sum_{\bm{p}\bm{p}^{\prime}ss^{\prime}}[T_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})a^{\dagger}_{\bm{p}s}b_{\bm{p}^{\prime}s^{\prime}}+T^{*}_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})b^{\dagger}_{\bm{p}^{\prime}s^{\prime}}a_{\bm{p}s}].$
(4)

From Heisenbergs equation, one obtains for the operators

$\displaystyle i\hbar\partial_{t}a_{\bm{p}s}=$
$\displaystyle+\sum_{\bm{p}^{\prime}s^{\prime}}[H^{L}_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})a_{\bm{p}^{\prime}s^{\prime}}+T_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})b_{\bm{p}^{\prime}s^{\prime}}],$
$\displaystyle i\hbar\partial_{t}a^{\dagger}_{-\bm{p}s}=$
$\displaystyle-\sum_{\bm{p}^{\prime}s^{\prime}}[[H^{L}_{ss^{\prime}}(-\bm{p},-\bm{p}^{\prime})]^{*}a^{\dagger}_{-\bm{p}^{\prime}s^{\prime}}+T^{*}_{ss^{\prime}}(-\bm{p},-\bm{p}^{\prime})b^{\dagger}_{-\bm{p}^{\prime}s^{\prime}}],$
$\displaystyle i\hbar\partial_{t}b_{\bm{p}s}=$
$\displaystyle+\sum_{\bm{p}^{\prime}s^{\prime}}[H^{R}_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})b_{\bm{p}^{\prime}s^{\prime}}+T^{*}_{s^{\prime}s}(\bm{p}^{\prime},\bm{p})a_{\bm{p}^{\prime}s^{\prime}}],$
$\displaystyle i\hbar\partial_{t}b^{\dagger}_{-\bm{p}s}=$
$\displaystyle-\sum_{\bm{p}^{\prime}s^{\prime}}[[H^{R}_{ss^{\prime}}(-\bm{p},-\bm{p}^{\prime})]^{*}b^{\dagger}_{-\bm{p}^{\prime}s^{\prime}}+T_{s^{\prime}s}(-\bm{p}^{\prime},-\bm{p})b^{\dagger}_{-\bm{p}^{\prime}s^{\prime}}].$
(5)

Using Nambu vectors,
$A_{\bm{p}}=\begin{pmatrix}a_{\bm{p}\uparrow}&a_{\bm{p}\downarrow}&a^{\dagger}_{-\bm{p}\uparrow}&a^{\dagger}_{-\bm{p}\downarrow}\end{pmatrix}$
and similarly for $B_{\bm{p}}$, this can be written as

$\displaystyle i\hbar\hat{\rho}_{3}\partial_{t}A_{\bm{p}}=$
$\displaystyle\sum_{\bm{p}^{\prime}}[\hat{H}^{L}(\bm{p},\bm{p}^{\prime})A_{\bm{p}^{\prime}}+\hat{T}(\bm{p},\bm{p}^{\prime})B_{\bm{p}^{\prime}}],$
(6) $\displaystyle i\hbar\hat{\rho}_{3}\partial_{t}B_{\bm{p}}=$
$\displaystyle\sum_{\bm{p}^{\prime}}[\hat{H}^{R}(\bm{p},\bm{p}^{\prime})B_{\bm{p}^{\prime}}+\hat{\tilde{T}}(\bm{p},\bm{p}^{\prime})A_{\bm{p}^{\prime}}],$
(7)

where $\hat{\rho}_{3}=\text{diag}(1,1,-1,-1)$ and

$\displaystyle\hat{H}^{L/R}(\bm{p},\bm{p}^{\prime})=$
$\displaystyle\begin{pmatrix}H^{L/R}(\bm{p},\bm{p}^{\prime})&0\\\
0&{H^{L/R}}^{*}(-\bm{p},-\bm{p}^{\prime})\end{pmatrix},$ (8)
$\displaystyle\hat{T}(\bm{p},\bm{p}^{\prime})=$
$\displaystyle\begin{pmatrix}T(\bm{p},\bm{p}^{\prime})&0\\\
0&T^{*}(-\bm{p},-\bm{p}^{\prime})\end{pmatrix},$ (9)
$\displaystyle\hat{\tilde{T}}(\bm{p},\bm{p}^{\prime})$
$\displaystyle=\begin{pmatrix}\tilde{T}(\bm{p},\bm{p}^{\prime})&0\\\
0&\tilde{T}^{*}(-\bm{p},-\bm{p}^{\prime})\end{pmatrix},$ (10)

where $T$ and $\tilde{T}$ are $2\times 2$ matrices in spin space and we have
defined the elements $\tilde{T}_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})\equiv
T^{*}_{s^{\prime}s}(\bm{p}^{\prime},\bm{p})$. The tunneling current from the
left to the right side is given as $J=e\langle\partial_{t}N\rangle$, with
$N=\sum_{\bm{p}s}a^{\dagger}_{\bm{p}s}a_{\bm{p}s}$. From the Heisenberg
equation we find

$\displaystyle i\hbar\partial_{t}N$ $\displaystyle=[N,\;H]=[N,\;H_{T}]$
$\displaystyle=\sum_{\bm{p}\bm{p}^{\prime}ss^{\prime}}\Big{[}T_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})a^{\dagger}_{\bm{p}s}b_{\bm{p}^{\prime}s^{\prime}}-T^{*}_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})b^{\dagger}_{\bm{p}^{\prime}s^{\prime}}a_{\bm{p}s}\Big{]}.$
(11)

We note that while we have here written $H_{L,R}$ without superconducting
terms, for notational simplicity, the approach generalizes straightforwardly
to superconducting systems by adding off-diagonal terms to Eq. 8. The
following derivation of the boundary conditions for the Green function is
valid irrespective of whether superconducting terms are included in
$\hat{H}^{L/R}$ or not. Next, we define the Green function

$\displaystyle
C^{<}_{ss^{\prime}}(\bm{p},t;\bm{p}^{\prime},t^{\prime})=i\langle
a^{\dagger}_{\bm{p}^{\prime}s^{\prime}}(t^{\prime})b_{\bm{p}s}(t)\rangle,$
(12)

which gives

$\displaystyle\langle\partial_{t}N\rangle=-\frac{1}{\hbar}\sum_{\bm{p}\bm{p}^{\prime}ss^{\prime}}$
$\displaystyle\Big{[}T_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})C^{<}_{s^{\prime}s}(\bm{p}^{\prime},t;\bm{p},t)$
$\displaystyle+T^{*}_{ss^{\prime}}(-\bm{p},-\bm{p}^{\prime})i\langle
a_{-\bm{p}s}b^{\dagger}_{-\bm{p}^{\prime}s^{\prime}}\rangle\Big{]}.$ (13)

Notice now that if we define a Nambu Green function matrix as

$\displaystyle\hat{C}^{<}(\bm{p},t;\bm{p}^{\prime},t^{\prime})=i\langle[[A^{\dagger}_{\bm{p}^{\prime}}(t^{\prime})]^{T}[B_{\bm{p}}(t)]^{T}]^{T}\rangle,$
(14)

the angular brackets in the second term of Section II.1 become the hole part
of $\hat{C}^{<}$. Hence, we may write

$\displaystyle\langle\partial_{t}N\rangle=-\frac{1}{\hbar}\sum_{\bm{p}\bm{p}^{\prime}}\text{Tr}[\hat{T}(\bm{p},\bm{p}^{\prime})\hat{C}^{<}(\bm{p}^{\prime},t;\bm{p},t)].$
(15)

For later use, we also introduce the Green function $\hat{C}^{>}$ and Keldysh
Green function $\hat{C}^{K}$ as:

$\displaystyle\hat{C}^{>}(\bm{p},t;\bm{p}^{\prime},t^{\prime})$
$\displaystyle=-i\left\langle
B_{\bm{p}}(t)A^{\dagger}_{\bm{p}^{\prime}}(t^{\prime})\right\rangle,$
$\displaystyle\hat{C}^{K}(\bm{p},t;\bm{p}^{\prime},t^{\prime})$
$\displaystyle=-i\hat{\rho}_{3}\langle[B_{\bm{p}}(t),A^{\dagger}_{\bm{p}^{\prime}}(t^{\prime})]\rangle.$
(16)

To find $\hat{C}^{<}$, we first find its time-ordered equivalent,

$\displaystyle\hat{C}(\bm{p},t;\bm{p}^{\prime},t^{\prime})=$
$\displaystyle-i\hat{\rho}_{3}\langle\mathcal{T}B_{\bm{p}}(t)A^{\dagger}_{\bm{p}^{\prime}}(t^{\prime})\rangle$
$\displaystyle=$ $\displaystyle-i\hat{\rho}_{3}\theta(t-t^{\prime})\langle
B_{\bm{p}}(t)A^{\dagger}_{\bm{p}^{\prime}}(t^{\prime})\rangle$
$\displaystyle+i\hat{\rho}_{3}\theta(t^{\prime}-t)\langle[[A^{\dagger}_{\bm{p}^{\prime}}(t^{\prime})]^{T}[B_{\bm{p}}(t)]^{T}]^{T}\rangle$
(17)

where $\mathcal{T}$ is the time-ordering operator, which acts as a normal-
operator at equal times. Differentiating with respect to $t$ gives us

$\displaystyle
i\hbar\hat{\rho}_{3}\partial_{t}\hat{C}(\bm{p},t;\bm{p}^{\prime},t^{\prime})$
$\displaystyle=\sum_{\bm{q}}\Big{[}\hat{H}^{R}(\bm{p},\bm{q})\hat{C}(\bm{q},t;\bm{p}^{\prime},t^{\prime})$
$\displaystyle+\hat{\tilde{T}}(\bm{p},\bm{q})\hat{G}_{L}(\bm{q},t;\bm{p}^{\prime},t^{\prime})\Big{]},$
(18)

where

$\displaystyle\hat{G}_{L}(\bm{q},t;\bm{p}^{\prime},t^{\prime})=-i\hat{\rho}_{3}\langle\mathcal{T}A_{\bm{q}}(t)A^{\dagger}_{\bm{p}^{\prime}}(t^{\prime})\rangle.$

This may be written as

$\displaystyle\sum_{\bm{q}}(i\hbar\hat{\rho}_{3}\partial_{t}\delta_{\bm{q}\bm{p}}-$
$\displaystyle\hat{H}^{R}(\bm{p},\bm{q}))\hat{C}(\bm{q},t;\bm{p}^{\prime},t^{\prime})$
$\displaystyle=\sum_{\bm{q}}\hat{\tilde{T}}(\bm{p},\bm{q})\hat{G}_{L}(\bm{q},t;\bm{q}^{\prime},t^{\prime}).$
(19)

A solution to this equation is given by

$\displaystyle\hat{C}(\bm{p},t;\bm{p}^{\prime},t^{\prime})$
$\displaystyle=\int
dt_{1}\sum_{\bm{q}\bm{q}^{\prime}}\hat{G}_{0,R}(\bm{p},t;\bm{q},t_{1})\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})$
$\displaystyle\times\hat{G}_{L}(\bm{q}^{\prime},t_{1};\bm{p}^{\prime},t^{\prime}),$
(20)

with $\hat{G}_{0,R}$ the right-side unperturbed Green function, satisfying

$\displaystyle\sum_{\bm{q}}(i\hbar\hat{\rho}_{3}\partial_{t}\delta_{\bm{p}\bm{q}}-\hat{H}^{R}(\bm{p},\bm{q}))\hat{G}_{0,R}(\bm{q},t;\bm{p}^{\prime},t^{\prime})=\delta(t-t^{\prime})\delta_{\bm{p}\bm{p}^{\prime}},$
(21)

as can be seen by direct insertion. Next, we deform the time axis to the
Keldysh contour,

$\displaystyle\hat{C}(\bm{p},\tau;\bm{p}^{\prime},\tau^{\prime})$
$\displaystyle=\int_{C}d\tau_{1}\sum_{\bm{q}\bm{q}^{\prime}}\hat{G}_{0,R}(\bm{p},\tau;\bm{q},\tau_{1})\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})$
$\displaystyle\times\hat{G}_{L}(\bm{q}^{\prime},\tau_{1};\bm{p}^{\prime},\tau^{\prime}),$
(22)

to obtain via the Langreth rules [37]

$\displaystyle\hat{\rho}_{3}\hat{C}^{<}(\bm{p},t;\bm{p}^{\prime},t^{\prime})=$
$\displaystyle\int
dt_{1}\sum_{\bm{q}\bm{q}^{\prime}}[\hat{G}^{R}_{0,R}(\bm{p},t;\bm{q},t_{1})\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})\hat{\rho}_{3}\hat{G}^{<}_{L}(\bm{q}^{\prime},t_{1};\bm{p}^{\prime},t^{\prime})+\hat{\rho}_{3}\hat{G}^{<}_{0,R}(\bm{p},t;\bm{q},t_{1})\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})\hat{G}^{A}_{L}(\bm{q}^{\prime},t_{1};\bm{p}^{\prime},t^{\prime})],$
(23)
$\displaystyle\hat{\rho}_{3}\hat{C}^{>}(\bm{p},t;\bm{p}^{\prime},t^{\prime})=$
$\displaystyle\int
dt_{1}\sum_{\bm{q}\bm{q}^{\prime}}[\hat{G}^{R}_{0,R}(\bm{p},t;\bm{q},t_{1})\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})\hat{\rho}_{3}\hat{G}^{>}_{L}(\bm{q}^{\prime},t_{1};\bm{p}^{\prime},t^{\prime})+\hat{\rho}_{3}\hat{G}^{>}_{0,R}(\bm{p},t;\bm{q},t_{1})\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})\hat{G}^{A}_{L}(\bm{q}^{\prime},t_{1};\bm{p}^{\prime},t^{\prime})].$
(24)

Note that the lesser and greater Green functions were defined without
$\hat{\rho}_{3}$, whereas the retarded, advanced and time-ordered were defined
with. Hence, when going from time-ordered to lesser/greater, there has to be
additional $\hat{\rho}_{3}$ appearing, as shown in the above equations. The
lesser and greater Green functions are further related to the Keldysh Green
function as

$\displaystyle\hat{G}^{K}(\bm{p},t;\bm{p}^{\prime},t^{\prime})=\hat{\rho}_{3}[\hat{G}^{<}(\bm{p},t;\bm{p}^{\prime},t^{\prime})+\hat{G}^{>}(\bm{p},t;\bm{p}^{\prime},t^{\prime})].$
(25)

The same relation also holds if one replaces $\hat{G}$ with $\hat{C}$. We also
have the relationship

$\displaystyle\hat{C}^{<}(\bm{p},t;\bm{p}^{\prime},t)=\hat{C}^{>}(\bm{p},t;\bm{p}^{\prime},t),$
(26)

By adding Eqs. (23) and then multiplying the whole equation with
$\hat{\rho}_{3}/2$, one therefore obtains:

$\displaystyle\hat{C}^{<}(\bm{p},t;\bm{p}^{\prime},t)=$
$\displaystyle\frac{1}{2}\hat{\rho}_{3}\int
dt_{1}\sum_{\bm{q}\bm{q}^{\prime}}[\hat{G}^{R}_{0,R}(\bm{p},t;\bm{q},t_{1})\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})\hat{G}^{K}_{L}(\bm{q}^{\prime},t_{1};\bm{p}^{\prime},t)+\hat{G}^{K}_{0,R}(\bm{p},t;\bm{q},t_{1})\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})\hat{G}^{A}_{L}(\bm{q}^{\prime},t_{1};\bm{p}^{\prime},t)].$
(27)

Insertion into Eq. 15 then gives

$\displaystyle\langle\partial_{t}N\rangle$
$\displaystyle=-\frac{1}{2\hbar}\sum_{\bm{p}\bm{p}^{\prime}\bm{q}\bm{q}^{\prime}}\int
dt_{1}\text{Tr}[\hat{\rho}_{3}(\hat{T}(\bm{p},\bm{p}^{\prime})\check{G}_{0,R}(\bm{p}^{\prime},t;\bm{q},t_{1})$
$\displaystyle\times\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})\check{G}_{L}(\bm{q}^{\prime},t_{1};\bm{p},t))^{K}],$
(28)

where it has been used that $\hat{T}$ commutes with $\hat{\rho}_{3}$. The
superscript K means in this context that the Keldysh component of the matrix
product enclosed in the parentheses should be extracted. Repeating the process
by calculating the currents flowing into the opposite material gives

$\displaystyle\langle\partial_{t}N\rangle$
$\displaystyle=+\frac{1}{2\hbar}\sum_{\bm{p}\bm{p}^{\prime}\bm{q}\bm{q}^{\prime}}\int
dt_{1}\text{Tr}\Big{[}\hat{\rho}_{3}(\check{G}_{0,L}(\bm{p}^{\prime},\bm{q})\hat{T}(\bm{q},\bm{q}^{\prime})$
$\displaystyle\times\check{G}_{R}(\bm{q}^{\prime},t;\bm{p},t_{1})\hat{\tilde{T}}(\bm{p},\bm{p}^{\prime}))^{K}\Big{]}.$
(29)

Averaging the two and Fourier transforming then gives

$\displaystyle\langle\partial_{t}N\rangle=$
$\displaystyle+\frac{1}{4\hbar}\sum_{\bm{p}\bm{p}^{\prime}\bm{q}\bm{q}^{\prime}}\int\frac{d\varepsilon}{2\pi}\text{Tr}\Big{[}\hat{\rho}_{3}(\check{G}_{0,L}(\bm{p},\bm{q};\varepsilon)\hat{T}(\bm{q},\bm{q}^{\prime})$
$\displaystyle\times\check{G}_{R}(\bm{q}^{\prime},\bm{p}^{\prime};\varepsilon)\hat{\tilde{T}}(\bm{p}^{\prime},\bm{p}))^{K}\Big{]}$
$\displaystyle-\frac{1}{4\hbar}\sum_{\bm{p}\bm{p}^{\prime}\bm{q}\bm{q}^{\prime}}\int\frac{d\varepsilon}{2\pi}\text{Tr}\Big{[}\hat{\rho}_{3}(\hat{T}(\bm{p}^{\prime},\bm{p})\check{G}_{0,R}(\bm{p},\bm{q};\varepsilon)$
$\displaystyle\times\check{\hat{T}}(\bm{q},\bm{q}^{\prime})\check{G}_{L}(\bm{q}^{\prime},\bm{p}^{\prime};\varepsilon))^{K}\Big{]}$
(30)

where we defined the $8\times 8$ matrix $\check{G}$ in Keldysh-space:

$\displaystyle\check{G}=\begin{pmatrix}\hat{G}^{R}&\hat{G}^{K}\\\
\hat{0}&\hat{G}^{A}\\\ \end{pmatrix}.$ (31)

A product $\hat{A}\check{B}$ between a $\hat{A}$-matrix and $\check{B}$-matrix
is to be understood as $\check{A}\check{B}$ where
$\check{A}=\text{diag}(\hat{A},\hat{A})$. In quasiclassical theory, the
electric current $J$ is computed from a quantity known as the matrix current
$\check{I}$ via the expression:

$\displaystyle J\sim\int dE\;\text{Tr}\\{\hat{\rho}_{3}\check{I}^{K}\\}$ (32)

Therefore, one may identify the matrix current as

$\displaystyle\check{I}\sim\sum_{\bm{p}\bm{p}^{\prime}\bm{q}\bm{q}^{\prime}}[\check{G}_{L},\hat{T}\check{G}_{R}\hat{\tilde{T}}]+\text{traceless
terms}.$ (33)

where we has set approximated $\check{G}_{0,j}\simeq\check{G}_{j}$ for
material $j$. The designation ”traceless terms” indicates that there may exist
terms in the matrix current which do not give any contribution to the charge
current, in effect which give zero when multiplied with $\hat{\rho}_{3}$ and
taken the trace over. In the next subsection, we shall identify these terms.

To compute the boundary conditions specific for an interface with Rashba spin-
orbit coupling, we assume an interface potential of the form

$\displaystyle
T(\bm{r})=T_{0}(\bm{r}_{\parallel})\delta(r_{\perp}-R_{0})-i\\{T_{ij}(\bm{r}_{\parallel})\delta(r_{\perp}-R_{0})\sigma_{j},\partial_{i}\\},$
(34)

where $r_{\parallel}$ and $r_{\perp}$ are coordinates parallel and
perpendicular to the interface. Note that the term involving the differential
operator $\partial_{i}$ has been symmetrized in order to ensure hermiticity.
The tunneling Hamiltonian is

$\displaystyle H_{T}=\sum_{ss^{\prime}}\int
d\bm{r}\;a_{s}^{\dagger}(\bm{r})T_{ss^{\prime}}(\bm{r})b_{s^{\prime}}(\bm{r})+\text{h.c.}$
(35)

with $a(\bm{r})$ and $b(\bm{r})$ the field operators on the left and right
side of the interface, respectively. We may now define
$a_{s}(\bm{r})=\sum_{\bm{p}}a_{\bm{p}s}e^{i\bm{p}\cdot\bm{r}}$, and
$b_{s}(\bm{r})=\sum_{\bm{p}}b_{\bm{p}s}e^{i\bm{p}\cdot\bm{r}}$ which produces
precisely $H_{T}$ in Eq. (2) with

$\displaystyle T_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})$ $\displaystyle=\int
d\bm{r}_{\parallel}[T_{0}(\bm{r}_{\parallel})+T_{ij}(\bm{r}_{\parallel})(p_{i}+p_{i}^{\prime})\sigma_{ss^{\prime},j}]$
$\displaystyle e^{-i(\bm{p}-\bm{p}^{\prime})\cdot\bm{r}_{\parallel}}.$ (36)

We now approximate the tunneling matrix element as follows to describe the
effect of interfacial Rashba-coupling:

$\displaystyle T_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})\simeq
T_{0}+T_{ij}(p_{i}+p_{i}^{\prime})\sigma_{ss^{\prime},j},$ (37)

where $T_{0}$ and $T_{ij}$ are real phenomenological constants. As an aside,
this approximation is equivalent to saying that both $T_{0}$ and $T_{ij}$ have
a $\delta(\bm{r}_{\parallel})$ spatial dependency, which means that the
surface is approximated as a point impurity. Moreover, we note that from the
definition of $\tilde{T}(\bm{p},\bm{p}^{\prime})$, it then follows that

$\displaystyle\tilde{T}_{ss^{\prime}}(\bm{p},\bm{p}^{\prime})\simeq
T_{0}+T_{ij}(p_{i}^{\prime}+p_{i})\sigma_{ss^{\prime},j}=T_{ss^{\prime}}(\bm{p},\bm{p}^{\prime}),$
(38)

Next, we need to calculate the integral

$\displaystyle\check{M}$ $\displaystyle=\int
d\bm{p}d\bm{p}^{\prime}d\bm{q}d\bm{q}^{\prime}\;\hat{T}(\bm{p}^{\prime},\bm{p})\check{G}_{R}(\bm{p},\bm{q})\hat{\tilde{T}}(\bm{q},\bm{q}^{\prime})\check{G}_{L}(\bm{q}^{\prime},\bm{p}^{\prime})$
(39)

which is one of the two main terms appearing in Eq. (30) after going to the
continuum limit in momentum-space. Defining

$\displaystyle\hat{T}(\bm{p},\bm{p}^{\prime})$
$\displaystyle=T_{0}+(p_{i}+p_{i}^{\prime})\hat{t}_{i},$
$\displaystyle\hat{t}_{i}$
$\displaystyle=T_{ij}\hat{\rho}_{3}\hat{\sigma}_{j},$
$\displaystyle\hat{\bm{\sigma}}$
$\displaystyle=\begin{pmatrix}\bm{\sigma}&0\\\ 0&\bm{\sigma}^{*}\\\
\end{pmatrix}.$ (40)

we get

$\displaystyle M$ $\displaystyle=\int
d\bm{p}d\bm{p}^{\prime}d\bm{q}d\bm{q}^{\prime}\;[T_{0}^{2}\check{G}_{R}(\bm{p},\bm{q})\check{G}_{L}(\bm{q}^{\prime},\bm{p}^{\prime})$
$\displaystyle+T_{0}\check{G}_{R}(\bm{p},\bm{q})\hat{t}_{k}(q_{k}+q_{k}^{\prime})\check{G}_{L}(\bm{q}^{\prime},\bm{p}^{\prime})$
$\displaystyle+T_{0}\hat{t}_{i}\check{G}_{R}(\bm{p},\bm{q})(p_{i}+p_{i}^{\prime})\check{G}_{L}(\bm{q}^{\prime},\bm{p}^{\prime})$
$\displaystyle+\hat{t}_{i}\check{G}_{R}(\bm{p},\bm{q})\hat{t}_{k}(p_{i}+p_{i}^{\prime})(q_{k}+q_{k}^{\prime})\check{G}_{L}(\bm{q}^{\prime},\bm{p}^{\prime})].$
(41)

We make use of the fact that the Green function
$\check{G}_{L/R}(\bm{p},\bm{q})$ is strongly peaked for
$\bm{p}\simeq\bm{q}\simeq{\bm{p}}_{F}$, since the aim is to construct a low-
energy theory relative the Fermi level. For such low energies, particle have
momenta close to ${\bm{p}}_{F}$. Since both $\bm{p}$ and $\bm{q}$ are close to
${\bm{p}}_{F}$, the center-of-mass momentum $({\bm{p}}+{\bm{q}})/2$ must also
be close to ${\bm{p}}_{F}$. Therefore, we may approximate the individual
momenta variables $\\{p_{i},p_{i}^{\prime},q_{k},q_{k}^{\prime}\\}$ entering
the above equation as center-of-mass variables with magnitude equal to the
Fermi momentum. Making use of the fact that the Jacobian of a transformation
from the momenta of the field operators to a center-of-mass and relative
momentum representation is unity, we find within the quasiclassical
approximation that

$\displaystyle\check{M}$ $\displaystyle\sim
T_{0}^{2}\check{g}_{R,s}\check{g}_{L,s}-m_{R}D_{R}T_{0}\check{g}_{R,s}\partial_{k}\check{g}_{R,s}\hat{t}_{k}\check{g}_{L,s}$
$\displaystyle-
m_{L}D_{L}T_{0}\check{g}_{R,s}\hat{t}_{k}\check{g}_{L,s}\partial_{k}\check{g}_{L,s}-m_{R}D_{R}T_{0}\hat{t_{i}}\check{g}_{R,s}\partial_{i}\check{g}_{R,s}\check{g}_{L,s}$
$\displaystyle-
m_{L}D_{L}T_{0}\hat{t}_{i}\check{g}_{R,s}\check{g}_{L,s}\partial_{i}\check{g}_{L,s}+\frac{2p_{F}^{2}}{3}\hat{t}_{i}\check{g}_{R,s}\hat{t}_{i}\check{g}_{L,s}.$
(42)

To obtain this result, we used that

$\displaystyle\check{g}_{R/L}(\bm{p}_{F},{\bm{R}}=0)=\frac{1}{(2\pi)^{3}}\frac{i}{\pi}\int
d\bm{q}\int d\xi_{\bm{p}}\check{G}_{R/L}(\bm{p},\bm{q})$ (43)

is the quasiclassical Green function evaluated at the interface and $\bm{q}$
in the equation above is the relative momentum-coordinate whereas $\bm{p}$ is
the center-of-mass momentum. We also inserted
$\check{g}_{R}=\check{g}_{R,s}+\hat{\bm{p}}_{F}\cdot\check{\bm{g}}_{R}$, where
$\check{\bm{g}}_{R}=-\tau v_{F}\check{g}_{R,s}\nabla\check{g}_{R,s}$, as
follows from the Usadel equation, and ${\bm{p}}$ is now the center-of-mass
momentum coordinate. For the angular averaging, we used that

$\displaystyle\int\frac{d\Omega_{p}}{4\pi}\;p_{F,i}p_{F,j}/p_{F}^{2}=\frac{1}{3}\delta_{ij},$

. Here, $m$ is the electron mass and we used that $p_{F}=mv_{F}$ and
$D=v_{F}^{2}\tau/3$. Since the parameters $T_{0}$ and $T_{ij}$ will be treated
as phenomenological interface parameters, we absorb the proportionality
constant in Eq. (II.1) into them. We also absorb a factor 1/3 into $p_{F}^{2}$
in the last term of Eq. (II.1) for brevity of notation. Finally, we neglected
terms in Eq. (II.1) which were second order in gradient terms of the Green
functions under the assumption that such terms are small.

The boundary conditions thus become

$\displaystyle\bm{n}\cdot D\check{g}_{L}\nabla\check{g}_{L}$
$\displaystyle=T_{0}^{2}[\check{g}_{L},\;\check{g}_{R}]+2T_{ij}T_{il}p_{F}^{2}[\check{g}_{L}\;,\;\hat{\rho}_{3}\hat{\sigma}_{j}\check{g}_{R}\hat{\rho}_{3}\hat{\sigma}_{l}]$
$\displaystyle-
mDT_{ij}T_{0}[\check{g}_{L},\;\\{\hat{\rho}_{3}\hat{\sigma}_{j}\;,\;\check{g}_{R}\partial_{i}\check{g}_{R}\\}]$
$\displaystyle-
mDT_{ij}T_{0}[\check{g}_{L}\partial_{i}\check{g}_{L},\;\\{\hat{\rho}_{3}\hat{\sigma}_{j}\;,\;\check{g}_{R}\\}]$
(44)

For the special case of Rashba SOC, we have
$T_{ij}=T_{1}n_{k}\varepsilon_{kji}$, where $n_{k}$ is component $k$ of the
interface normal. This gives

$\displaystyle\bm{n}\cdot D\check{g}_{L}\nabla\check{g}_{L}$
$\displaystyle=T_{0}^{2}[\check{g}_{L},\;\check{g}_{R}]+2T_{1}^{2}p_{F}^{2}[\check{g}_{L},\;\bm{\sigma}_{||}\check{g}_{R}\bm{\sigma}_{||}]$
$\displaystyle-
mDT_{1}T_{0}[\check{g}_{L},\;\\{\bm{\sigma}_{||},\;(\check{g}_{R}\nabla\check{g}_{R})\times\bm{n}\\}]$
$\displaystyle-
mDT_{1}T_{0}[(\check{g}_{L}\nabla\check{g}_{L})\times\bm{n},\;\\{\bm{\sigma}_{||},\;\check{g}_{R}\\}].$
(45)

with
$\hat{\bm{\sigma}}_{||}=\hat{\rho}_{3}[\hat{\bm{\sigma}}-\bm{n}(\hat{\bm{\sigma}}\cdot\bm{n})]$.
We see that the second term on the right-hand side of Eq. (II.1) reproduces
the tunneling term due to spin-orbit coupling derived in Ref. [36]. The third
and fourth terms on the right-hand side is new and involves the gradient of
the Green function. This allows us to describe for instance spin-Hall like
phenomena with supercurrents, as we will demonstrate below.

###  II.2 Terms related to reflection at the interface

Figure 1: (Color online) Spin-orbit coupled interface separating left (L) and
right (R) region in the model used to derive the reflection boundary condition
terms.

Here, we follow a similar procedure as [38, 41] to derive the terms related to
reflection at the interface. The Usadel equation including antisymmetric spin-
orbit coupling reads:

$\displaystyle
D\tilde{\nabla}(\check{g}\tilde{\nabla}\check{g})+i[\varepsilon\hat{\rho}_{3}+\hat{\Sigma},\check{g}]=0.$
(46)

Here, we have defined

$\displaystyle\tilde{\nabla}\cdot\ldots=\nabla\cdot\ldots-i[\hat{\bm{A}},\ldots].$
(47)

We use the convention that when multiplying a matrix $\hat{A}$ ($4\times 4$)
with a matrix $\check{B}$ ($8\times 8$), that product should be understood as
$\text{diag}(\hat{A},\hat{A})\cdot\check{B}$. The matrix $\hat{\bm{A}}$ has
vectors as elements and classifies the spin-orbit coupling according to

$\displaystyle\hat{\bm{A}}=\begin{pmatrix}{\bm{A}}&0\\\ 0&-{\bm{A}}^{*}\\\
\end{pmatrix}$ (48)

where ${\bm{A}}$ defines the spin-orbit coupling part of the Hamilton-operator
via:

$\displaystyle H_{\text{soc}}=-\frac{1}{m}{\bm{k}}\cdot{\bm{A}}.$ (49)

This is a generally valid Hamilton-operator for an antisymmetric spin-orbit
coupling that is linear in momentum. For instance, the Rashba case
$H_{\text{soc}}=-\frac{\alpha}{m}({\bm{k}}\times{\bm{\sigma}})\cdot{\bm{z}}$
for structures breaking inversion symmetry in the $z$-direction is obtained
using ${\bm{A}}=\alpha(-\sigma_{y},\sigma_{x},0)$. Here, we will derive
boundary conditions valid for any ${\bm{A}}$ that satisfies
${\bm{n}}\cdot{\bm{A}}=0$ where ${\bm{n}}$ is the interface normal vector.
This includes the Rashba interaction as a limiting case. The criterion
${\bm{n}}\cdot{\bm{A}}=0$ is chosen as to avoid the complication of
symmetrizing the spin-orbit coupling Hamiltonian, as is required for the
hermiticity of the Hamilton-operator when the spin-orbit coupling is confined
to a certain spatial region. Nevertheless, the condition
${\bm{n}}\cdot{\bm{A}}$ still allows for several interesting new phenomena
related to spin-Hall physics to be predicted.

The underlying assumption is that the region with spin-orbit coupling is
sufficiently thin that we may neglect the spatial variation of the Green
function along the interface normal ${\bm{n}}$. Let $n$ denote the coordinate
along the ${\bm{n}}$-axis (chosen as $z$-direction in Fig. 1). Integrating the
Usadel-equation from $n=-d$ to $n=d$ gives:

$\displaystyle(\check{g}\partial_{n}$
$\displaystyle\check{g})|_{d}-(\check{g}\partial_{n}\check{g})|_{-d}-id\nabla\check{g}[\hat{{\bm{A}}},\check{g}]|_{d}$
$\displaystyle-i\check{g}d[\hat{{\bm{A}}},\nabla{\check{g}}]-id[\hat{{\bm{A}}},\check{g}\nabla\check{g}]|_{d}-d[\hat{{\bm{A}}},\check{g}[\hat{{\bm{A}}},\check{g}]]|_{d}=0.$
(50)

Here, we made use of the fact that ${\bm{n}}\cdot{\bm{A}}=0$ and approximated
the value of all terms not involving a derivative in the $n$-direction with
their value at $n=d$. This is consistent with the assumption that $d$ is
smaller than the length scale over which the Green function varies in the
$n$-direction. Now, if there was no spin-orbit coupling present, the tunneling
current alone derived in the previous section (with scalar tunneling
amplitudes due to the absence of SOC) provides the value of
$\check{g}\partial_{n}\check{g}|_{d}$, constituting the effect boundary
condition for the Green function $\check{g}$. Therefore, we may identify:

$\displaystyle(\check{g}\partial_{n}\check{g})|_{-d}=\text{tunneling current
contribution}\equiv I_{T},$ (51)

and consequently write Eq. (II.2) as:

$\displaystyle\check{g}\partial_{n}$
$\displaystyle\check{g}|_{d}=I_{T}+id\nabla\check{g}[\hat{{\bm{A}}},\check{g}]_{d}+id\check{g}[\hat{{\bm{A}}},\nabla\check{g}]_{d}$
$\displaystyle+id[\hat{{\bm{A}}},\check{g}\nabla\check{g}]_{d}+d[\hat{{\bm{A}}},\check{g}[\hat{{\bm{A}}},\check{g}]]_{d}.$
(52)

After cancelling some terms, we end up with the final form which is generally
valid for any SOC at the interface so long that ${\bm{n}}\cdot{\bm{A}}=0$:

$\displaystyle\bm{n}\cdot\check{g}\nabla$
$\displaystyle\check{g}|_{d}=I_{T}+id\nabla\check{g}[\hat{{\bm{A}}},\check{g}]|_{d}-2id\check{g}\nabla\check{g}\hat{{\bm{A}}}_{d}+id\check{g}\hat{{\bm{A}}}\nabla\check{g}|_{d}$
$\displaystyle+id\hat{{\bm{A}}}\check{g}\nabla\check{g}_{d}+d[\hat{{\bm{A}}},\check{g}\hat{{\bm{A}}}\check{g}]_{d}.$
(53)

This constitutes the boundary condition for the Green function on the right
side of the interface $\check{g}_{R}=\check{g}|_{d}$. The boundary condition
for the Green function on the left side of the interface
$\check{g}_{L}=\check{g}|_{-d}$ is obtained from Eq. (II.2) by letting
$d\to-d$. This is derived in exactly the same way we did for $\check{g}_{R}$
except that we integrate $\int_{d}^{-d}$ instead of $\int_{-d}^{d}$.

As an application of these boundary conditions, consider now a special case.
We examine a 2D structure ($yz$-plane) with a Rashba spin-orbit coupled
interface. The $z$-direction is perpendicular to the interface plane, so that
$\bm{n}=\hat{\bm{z}}$. We permit not only $\partial_{z}\check{g}_{R}\neq 0$,
but also $\partial_{y}\check{g}_{R}\neq 0$ in order to capture spin Hall like
phenomena involving motion parallell with the interface. In the Rashba case,
we have as previously mentioned ${\bm{A}}=\alpha(-\sigma_{y},\sigma_{x},0)$,
which gives

$\displaystyle\check{g}_{R}\partial_{z}$
$\displaystyle\check{g}_{R}=\check{I}_{T}+id\partial_{y}\check{g}_{R}[\hat{A}_{y},\check{g}_{R}]-2id\check{g}_{R}\partial_{y}\check{g}_{R}\hat{A}_{y}+id\check{g}_{R}\hat{A}_{y}\partial_{y}\check{g}_{R}$
$\displaystyle+id\hat{A}_{y}\check{g}_{R}\partial_{y}\check{g}_{R}+d[\hat{{\bm{A}}},\check{g}_{R}\hat{{\bm{A}}}\check{g}_{R}].$
(54)

We define

$\displaystyle\hat{\rho}_{x}=\begin{pmatrix}\sigma_{x}&0\\\
0&-\sigma_{x}^{*}\\\
\end{pmatrix},\;\hat{\rho}_{y}=\begin{pmatrix}\sigma_{y}&0\\\
0&-\sigma_{y}^{*}\\\ \end{pmatrix},$ (55)

so that $\hat{A}_{x}=-\alpha\hat{\rho}_{y}$ and
$\hat{A}_{y}=\alpha\hat{\rho}_{x}$. We then get:

$\displaystyle\check{g}_{R}\partial_{z}\check{g}_{R}$
$\displaystyle=\check{I}_{T}+id\alpha\partial_{y}\check{g}_{R}[\hat{\rho}_{x},\check{g}_{R}]-2id\alpha\check{g}_{R}\partial_{y}\check{g}_{R}\hat{\rho}_{x}+id\alpha\check{g}_{R}\hat{\rho}_{x}\partial_{y}\check{g}_{R}+id\alpha\hat{\rho}_{x}\check{g}_{R}\partial_{y}\check{g}_{R}$
$\displaystyle+d\alpha^{2}[\hat{\rho}_{x},\check{g}_{R}\hat{\rho}_{x}\check{g}_{R}]+d\alpha^{2}[\hat{\rho}_{y},\check{g}_{R}\hat{\rho}_{y}\check{g}_{R}].$
(56)

It is clear from the boundary conditions that any variation in the Green
function $\check{g}_{R}$ in the direction parallell to the interface
($y$-direction) will now couple directly to variations of the Green function
in direction normal to the interface ($z$-direction) due to the Rashba spin-
orbit coupling at the interface.

###  II.3 Final form of complete boundary condition for Rashba-coupling

In the preceeding sections, we have derived the contribution to the boundary
condition for the quasiclassical Green function coming both from tunneling
terms and reflection terms at an interface with spin-orbit coupling. Here, we
give the final complete form of the boundary conditions for the Green function
on both the left and right side of the interface. Let $\bm{n}$ point from left
to right.

For the quasiclassical Green function matrix $\check{g}_{L}$ on the left side
of the interface, the boundary condition generally reads

$\displaystyle\bm{n}\cdot D\check{g}_{L}\nabla\check{g}_{L}$
$\displaystyle=T_{0}^{2}[\check{g}_{L},\;\check{g}_{R}]-mDT_{ij}T_{0}[\check{g}_{L}\;,\;\\{\hat{\rho}_{3}\hat{\sigma}_{j}\;,\;\check{g}_{R}\partial_{i}\check{g}_{R}\\}]\;-mDT_{ij}T_{0}[\check{g}_{L}\partial_{i}\check{g}_{L}\;,\;\\{\hat{\rho}_{3}\hat{\sigma}_{j}\;,\;\check{g}_{R}\\}]+T_{ij}T_{il}p_{F}^{2}[\check{g}_{L}\;,\;\hat{\rho}_{3}\hat{\sigma}_{j}\check{g}_{R}\hat{\rho}_{3}\hat{\sigma}_{l}]$
$\displaystyle-
iDd\nabla\check{g}_{L}[\hat{{\bm{A}}},\check{g}_{L}]+2iDd\check{g}_{L}\nabla\check{g}_{L}\hat{{\bm{A}}}-iDd\check{g}_{L}\hat{{\bm{A}}}\nabla\check{g}_{L}-iDd\hat{{\bm{A}}}\check{g}_{L}\nabla\check{g}_{L}-Dd[\hat{{\bm{A}}},\check{g}_{L}\hat{{\bm{A}}}\check{g}_{L}].$
(57)

while for the Green function on the right side:

$\displaystyle\bm{n}\cdot D\check{g}_{R}\nabla\check{g}_{R}$
$\displaystyle=T_{0}^{2}[\check{g}_{L},\;\check{g}_{R}]+mDT_{ij}T_{0}[\check{g}_{R}\;,\;\\{\hat{\rho}_{3}\hat{\sigma}_{j}\;,\;\check{g}_{L}\partial_{i}\check{g}_{L}\\}]\;+mDT_{ij}T_{0}[\check{g}_{R}\partial_{i}\check{g}_{R}\;,\;\\{\hat{\rho}_{3}\hat{\sigma}_{j}\;,\;\check{g}_{L}\\}]-T_{ij}T_{il}p_{F}^{2}[\check{g}_{R}\;,\;\hat{\rho}_{3}\hat{\sigma}_{j}\check{g}_{L}\hat{\rho}_{3}\hat{\sigma}_{l}]$
$\displaystyle+iDd\nabla\check{g}_{R}[\hat{{\bm{A}}},\check{g}_{R}]-2iDd\check{g}_{R}\nabla\check{g}_{R}\hat{{\bm{A}}}+iDd\check{g}_{R}\hat{{\bm{A}}}\nabla\check{g}_{R}+iDd\hat{{\bm{A}}}\check{g}_{R}\nabla\check{g}_{R}+Dd[\hat{{\bm{A}}},\check{g}_{R}\hat{{\bm{A}}}\check{g}_{R}].$
(58)

For the special case of Rashba spin-orbit coupling and an interface normal
$\bm{n}$ along the $z$-direction with $\check{g}_{L}=\check{g}_{L}(y)$ due to
supercurrent flow in the $y$-direction and a thin region to the right of the
interface so that $\check{g}_{R}=\check{g}_{R}(z)$, as shown in Fig. 2, one
obtains on the left side:

$\displaystyle D\check{g}_{L}\partial_{z}\check{g}_{L}$
$\displaystyle=T_{0}^{2}[\check{g}_{L},\;\check{g}_{R}]+T_{1}^{2}p_{F}^{2}[\check{g}_{L}\;,\;\hat{\bm{\sigma}}_{||}\check{g}_{R}\hat{\bm{\sigma}}_{||}]-mDT_{1}T_{0}[\check{g}_{L}\;,\;\\{\hat{\sigma}_{||,x}\;,\;\check{g}_{R}\partial_{y}\check{g}_{R}\\}]-mDT_{1}T_{0}[\check{g}_{L}\partial_{y}\check{g}_{L}\;,\;\\{\hat{\sigma}_{||,x}\;,\;\check{g}_{R}\\}]$
$\displaystyle-
Dd\alpha^{2}[\hat{\rho}_{x},\check{g}_{L}\hat{\rho}_{x}\check{g}_{L}]-Dd\alpha^{2}[\hat{\rho}_{y},\check{g}_{L}\hat{\rho}_{y}\check{g}_{L}],$
(59)

while on the right side:

$\displaystyle D\check{g}_{R}\partial_{z}\check{g}_{R}$
$\displaystyle=T_{0}^{2}[\check{g}_{L},\;\check{g}_{R}]-T_{1}^{2}p_{F}^{2}[\check{g}_{R}\;,\;\hat{\bm{\sigma}}_{||}g_{L}\hat{\bm{\sigma}}_{||}]+mDT_{1}T_{0}[\check{g}_{R}\;,\;\\{\hat{\sigma}_{||,x}\;,\;\check{g}_{L}\partial_{y}\check{g}_{L}\\}]+mDT_{1}T_{0}[\check{g}_{R}\partial_{y}\check{g}_{R}\;,\;\\{\hat{\sigma}_{||,x}\;,\;\check{g}_{L}\\}]$
$\displaystyle+Dd\alpha^{2}[\hat{\rho}_{x},\check{g}_{R}\hat{\rho}_{x}\check{g}_{R}]+Dd\alpha^{2}[\hat{\rho}_{y},\check{g}_{R}\hat{\rho}_{y}\check{g}_{R}].$
(60)

In this particular case, one has
$\bm{\hat{\sigma}}_{||}=(\hat{\rho}_{x},\hat{\rho}_{y},\hat{0})$. The
derivation of the boundary conditions provided here treats the interface in a
phenomenological manner, in particular when it comes to the reflection terms.
A stricter microscopic derivation makes use of concepts like isotropization
zone and ballistic zones near an interface [39]. Nevertheless, the methodology
used in this paper to derive the boundary conditions is known to give the same
results as the microscopic derivation in the case of magnetic boundary
conditions [38], a fact which motivates the present approach.

##  III Application: supercurrent-induced magnetization via Rashba interface

We here show how the boundary conditions can be used on a concrete system and
that they predict that a supercurrent can induce a non-local magnetization via
structural inversion symmetry breaking. The polarization of the induced
magnetization is perpendicular to both the current flow and the direction of
inversion symmetry breaking and can be viewed as a type of superconducting
spin Hall effect. We consider a superconductor in contact with a normal metal
through a thin heavy metal interface. Two possible experimental realizations
of this is shown in Fig. 1. When a supercurrent runs through the
superconductor, a non-local magnetization appears in the normal metal despite
the absence of any magnetic elements in the system. To show this analytically,
we consider the weak proximity effect regime where the Usadel equation is
linearized in the anomalous Green function.

Figure 2: (Color online) Two possible experimental setups for observation of
supercurrent-induced magnetization via a Rashba interface. In both cases,
current-biasing the superconductor induces a non-local spin magnetization in a
normal metal without any magnetic elements in the system. The magnetization
points in the direction orthogonal to both the supercurrent flow and the
broken inversion symmetry axis, similarly to the conventional spin Hall
effect.

The quasiclassical retarded Green function in a supercurrent-carrying
superconductor is given by

$\displaystyle\hat{g}_{S}=\begin{pmatrix}c&s\mathrm{e}^{iJy}i\sigma_{y}\\\
s\mathrm{e}^{-iJy}i\sigma_{y}&-c\\\ \end{pmatrix}$ (61)

for a current density $J$ well below the critical supercurrent density. Here,
$c=\cosh(\theta)$, $s=\sinh(\theta)$, where
$\theta=\text{atanh}(\Delta/\varepsilon)$ where $\Delta$ is the
superconducting gap. In the normal metal, the retarded Green function matrix
reads

$\displaystyle\hat{g}_{N}=\begin{pmatrix}1&f\\\ -\tilde{f}&-1\\\
\end{pmatrix}$ (62)

in the case of a weak proximity effect where the anomalous Green function
matrix has the form

$\displaystyle f=\begin{pmatrix}f_{\uparrow}&f_{t}+f_{s}\\\
f_{t}-f_{s}&f_{\downarrow}\\\ \end{pmatrix}.$ (63)

The $\tilde{\ldots}$ operator denotes reversal of energy
($\varepsilon\to-\varepsilon$) and complex conjugation. In the normal metal,
the Usadel equation then takes the form

$\displaystyle D\partial_{z}^{2}f_{s,t,\sigma}+2i\varepsilon f_{s,t,\sigma}=0$
(64)

and has a general solution

$\displaystyle
f_{s,t,\sigma}=A_{s,t,\sigma}\mathrm{e}^{ikz}+B_{s,t,\sigma}\mathrm{e}^{-ikz}$
(65)

where $k=\sqrt{2i\varepsilon/D}$. The boundary condition at the vacuum edge
$z=L$ of the normal metal is

$\displaystyle\partial_{z}f_{s,t,\sigma}=0$ (66)

in order to ensure that no current flows into vacuum. The boundary condition
between the normal metal and the superconductor at $z=0$, which are connected
via an atomically thin Pt layer, is the key ingredient for achieving the
superspin Hall effect. In the weak proximity effect regime, we obtain

$\displaystyle D\partial_{z}f_{\sigma}+4i\sigma
DJT_{0}T_{1}csm\mathrm{e}^{iJy_{0}}-2T_{0}^{2}cf_{\sigma}$
$\displaystyle\;\;-4i\sigma
DmT_{1}Js^{2}T_{0}f_{s}-4(Dd\alpha^{2}+2cT_{1}^{2}p_{F}^{2})f_{\sigma}=0,$
$\displaystyle
D\partial_{z}f_{t}-8\alpha^{2}Ddf_{t}-c(T_{0}^{2}+4T_{1}^{2}p_{F}^{2})f_{t}=0,$
$\displaystyle
D\partial_{z}f_{s}-8cT_{1}^{2}p_{F}^{2}f_{s}-4iDT_{1}T_{0}s^{2}Jm(f_{\uparrow}-f_{\downarrow})$
$\displaystyle\;\;-4cT_{0}^{2}f_{s}+4s\mathrm{e}^{iJy_{0}}(T_{0}^{2}+4T_{1}^{2}p_{F}^{2})=0.$
(67)

Here, $y_{0}$ is the position along the $y$-axis where the normal metal is
coupled to the superconductor.

From these boundary conditions, one obtains the solution for the unknown
coefficients $A_{s,t,\sigma}$ and $B_{s,t,\sigma}$. With the anomalous Green
function in hand, one may compute the spin magnetization induced in the normal
metal:

$\displaystyle{\bm{M}}=M_{0}\int^{\infty}_{-\infty}\frac{d\varepsilon}{\varepsilon_{\text{Th}}}\text{Re}\text{Tr}\\{{\bm{\sigma}}\hat{g}^{K}\\}$
(68)

where $\varepsilon_{\text{Th}}=D/L^{2}$ is the Thouless energy and $M_{0}$ is
a normalization constant of the same dimension as ${\bm{M}}$.

The solution for the coefficients $A_{s,t,\sigma}$ and $B_{s,t,\sigma}$ are
given in the appendix, and show that the anomalous Green function $f_{t}$ is
zero whereas $f_{\uparrow}=-f_{\downarrow}$. Therefore, the only component of
the induced magnetization that exists is $M_{x}$ which is determined by the
product of the singlet anomalous Green function $f_{s}$ and the $x$-component
of the triplet $d$-vector [2], $(f_{\downarrow}-f_{\uparrow})/2$. From the
expressions in the appendix, it then follows that

$\displaystyle M_{x}\sim JT_{0}T_{1}\Delta^{2}.$ (69)

Therefore, the induced magnetization only exists in the presence of a
supercurrent $(J)$ and in the presence of both normal tunneling $(T_{0})$ and
a Rashba-like tunneling $(T_{1})$. We note that spin-orbit impurity scattering
[40] can cause a similar effect due to supercurrent flow. Moreover, we note
that the effect predicted here is different from Ref. [41] since no exchange
field is required anywhere in the junction.

##  IV Conclusion

In summary, we have derived a set of quasiclassical boundary conditions for
spin-orbit coupled interfaces. The boundary conditions contain both terms
related to transmission and reflection of particles and are valid both in and
out of equilibrium. As an application, we have shown that a supercurrent
flowing in a superconductor that is coupled to a normal metal through a
Rashba-interface will induce an non-local magnetization in the normal metal.
The magnetization is polarized in the direction perpendicular to both the
supercurrent flow and the interface normal. These boundary conditions may
prove useful to predict new spin-Hall like phenomena in superconducting hybrid
structures and to model experimental data.

###### Acknowledgements.

E. H. Fyhn and S. Aunsmo is thanked for useful discussions. This work was
supported by the Research Council of Norway through grant 240806, and its
Centres of Excellence funding scheme grant 262633 “ _QuSpin_ ”.

##  V Appendix

The coefficients determining the anomalous Green functions are obtained as

$\displaystyle B_{s,t,\sigma}=A_{s,t,\sigma}\mathrm{e}^{2ikL}$ (70)

and $A_{\downarrow}=-A_{\uparrow}$ where $e_{\pm}\equiv\mathrm{e}^{\pm ikL}$
and

$\displaystyle A_{\uparrow}$
$\displaystyle=\frac{4iT_{1}DT_{0}mJ\mathrm{e}^{iJy_{0}}s\Big{[}e_{+}\Big{(}(-4T_{0}^{2}-16T_{1}^{2}p_{F}^{2})c^{2}+4s^{2}(4T_{1}^{2}p_{F}^{2}+T_{0}^{2})\Big{)}+iDce_{-}k\Big{]}}{\Big{(}-(T_{0}^{2}+4T_{1}^{2}p_{F}^{2})^{2}c^{2}-2D\alpha^{2}d(T_{0}^{2}+4T_{1}^{2}p_{F}^{2})c-8D^{2}m^{2}T_{0}^{2}T_{1}^{2}p_{F}^{2}J^{2}s^{4}\Big{)}8e_{+}^{2}+4i\Big{(}(3T_{0}^{2}/2+6T_{1}^{2}p_{F}^{2})c+D\alpha^{2}d\Big{)}ke_{-}De_{+}+D^{2}e_{-}^{2}k^{2}}$
$\displaystyle A_{s}$
$\displaystyle=\frac{4\mathrm{e}^{iJy_{0}}s\Big{[}\Big{(}-16D^{2}J^{2}T_{0}^{2}T_{1}^{2}cs^{2}m^{2}-4\alpha^{2}d(4T_{1}^{2}p_{F}^{2}+T_{0}^{2})D-2c(T_{0}^{2}+4T_{1}^{2}p_{F}^{2})(4T_{1}^{2}p_{F}^{2}+T_{0})\Big{)}e_{+}+i(4T_{1}^{2}p_{F}^{2}+T_{0}^{2})kDe_{-}\Big{]}}{e_{+}^{2}\Big{(}-8(T_{0}^{2}+4T_{1}^{2}p_{F}^{2})^{2}c^{2}-16D\alpha^{2}d(T_{0}^{2}+4T_{1}^{2}p_{F}^{2})c-64D^{2}m^{2}T_{0}^{2}T_{1}^{2}J^{2}s^{4}\Big{)}+4ikDe_{+}e_{-}\Big{(}D\alpha^{2}d+3c(T_{0}^{2}+4T_{1}^{2}p_{F}^{2})/2\Big{)}+D^{2}e_{-}^{2}k^{2}}$
(71)

Note that the value of the interface location along the $y$-axis (the
coordinate $y_{0}$) does not affect the value of physical observables such as
the induced non-local magnetization.

## References

  * [1] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175–1204 (1957).
  * [2] A. P. Mackenzie and Y. Maeno, Rev. Mod. Phys. 75, 657 (2003).
  * [3] F. S. Bergeret, A. F. Volkov, and K. B. Efetov, Rev. Mod. Phys. 77, 1321 (2005).
  * [4] A. I. Buzdin, Rev. Mod. Phys. 77, 935 (2005)
  * [5] F. S. Bergeret and I. V. Tokatly, Phys. Rev. B 89, 134517 (2014).
  * [6] J. Linder and J. W. A. Robinson, Nat. Phys. 11, 307–315 (2015).
  * [7] M. Eschrig, Rep. Prog. Phys. 78, 104501 (2015).
  * [8] V. L. Berezinskii, Pis’ma Zh. Eksp. Teor. Fiz. 20, 628–631 (1974).
  * [9] J. Linder and A. V. Balatsky, Rev. Mod. Phys. 91, 045005 (2019).
  * [10] K. Ohnishi, S. Komori, G. Yang, K.-R. Jeon, L. A. B. Olde Olthof, X. Montiel, M. G. Blamire, and J. W. A. Robinson, Appl. Phys. Lett. 116, 130501 (2020).
  * [11] A. Brataas, A. D. Kent, and H. Ohno, Nature Materials 11, 372–381 (2012).
  * [12] W. Han, S. Maekawa, and X.-C. Xie, Nat. Mat. 19, 139 (2020).
  * [13] D. B. Szombati et al. , Nat. Phys. 12, 568 (2016).
  * [14] V. V. Ryazanov, V. A. Oboznov, A. Yu. Rusanov, A. V. Veretennikov, A. A. Golubov, and J. Aarts, Phys. Rev. Lett. 86, 2427 (2001).
  * [15] T. Kontos, M. Aprili, J. Lesueur, and X. Grison, Phys. Rev. Lett. 86, 304 (2001).
  * [16] R. S. Keizer, S. T. B. Goennenwein, T. M. Klapwijk, G. Miao, G. Xiao, A. Gupta, Nature 439, 825 (2006).
  * [17] T. S. Khaire, M. A. Khasawneh, W. P. Pratt, Jr., and N. O. Birge, Phys. Rev. Lett. 104, 137002 (2010).
  * [18] J. W. A. Robinson, J. D. S. Witt, and M. G. Blamire, Science 329, 59-61 (2010).
  * [19] S. Kolenda, M. J. Wolf, and D. Beckmann, Phys. Rev. Lett. 116, 097001 (2016).
  * [20] E. C. Gingrich, B. M. Niedzielski, J. A. Glick, Y. Wang, D. L. Miller, R. Loloee, W. P. Pratt Jr, and N. O. Birge, Nature Physics 12, 564–567 (2016).
  * [21] C. R. Reeg and D. L. Maslov, Phys. Rev. B 92, 134512 (2015)
  * [22] P. Högl, A. Matos-Abiague, I. Žutić, and J. Fabian, Phys. Rev. Lett. 115, 116601 (2015),
  * [23] I. V. Bobkova and A. M. Bobkov, Phys. Rev. B 95, 184518 (2017)
  * [24] V. Mishra, Y. Li, F.-C. Zhang, and S. Kirchner, Phys. Rev. B 103, 184505 (2021)
  * [25] A. Costa and J. Fabian, Phys. Rev. B 104, 174504 (2021)
  * [26] N. Banerjee, J. A. Ouassou, Y. Zhu, N. A. Stelmashenko, J. Linder, and M. G. Blamire, Phys. Rev. B 97, 184521 (2018)
  * [27] N. Satchell and N. O. Birge, Phys. Rev. B 97, 214509 (2018)
  * [28] N. Satchell, R. Loloee, and N. O. Birge, Phys. Rev. B 99, 174519 (2019)
  * [29] I. Martínez, P. Högl, C. González-Ruano, J. P. Cascales, C. Tiusan, Y. Lu, M. Hehn, A. Matos-Abiague, J. Fabian, I. Žutić, and F. G. Aliev, Phys. Rev. Applied 13, 014030 (2020)
  * [30] C. González-Ruano, L. G. Johnsen, D. Caso, C. Tiusan, M. Hehn, N. Banerjee, J. Linder, and F. G. Aliev, Phys. Rev. B 102, 020405(R) (2020)
  * [31] V. P. Amin, P. M. Haney, and M. D. Stiles, Journal of Applied Physics 128, 151101 (2020)
  * [32] J. Linder and T. Yokoyama, Phys. Rev. Lett. 106, 237201 (2011)
  * [33] J. W. Serene, and D. Rainer, Phys. Rep. 101, 221 (1983).
  * [34] W. Belzig, F. K. Wilhelm, C. Bruder, G. Schön, and A. D. Zaikin, Superlatt. Microstruct. 25, 1251 (1999).
  * [35] J. Rammer, and H. Smith, Rev. Mod. Phys. 58, 323 (1986).
  * [36] M. Amundsen and J. Linder, Phys. Rev. B 100, 064502 (2019) [Erratum: Phys. Rev. B 103, 059902 (2021)].
  * [37] H. Haug and A. P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors, Springer, Berlin, 1996, ISBN: 3540616020.
  * [38] F. S. Bergeret, A. Verso, and A. F. Volkov, Phys. Rev. B 86, 214516 (2012).
  * [39] A. Cottet, D. Huertas-Hernando, W. Belzig, and Y. V. Nazarov, Phys. Rev. B 80, 184511 (2009).
  * [40] F. Sebastian Bergeret and Ilya V. Tokatly, Phys. Rev. B 94, 180502(R) (2016).
  * [41] M. A. Silaev, I. V. Bobkova, and A. M. Bobkov, Phys. Rev. B 102, 100507(R)


# Ultrafast dynamics of a fermion chain in a terahertz field-driven optical
cavity

Mohsen Yarmohammadi<EMAIL_ADDRESS>Department of Physics,
The University of Texas at Dallas, Richardson, Texas 75080, USA  John Sous
<EMAIL_ADDRESS>Department of Chemistry and Biochemistry, University of
California, San Diego, La Jolla, California 92093, USA  Marin Bukov
<EMAIL_ADDRESS>Max Planck Institute for the Physics of Complex Systems,
Nöthnitzer Str. 38, 01187 Dresden, Germany  Michael H. Kolodrubetz
<EMAIL_ADDRESS>Department of Physics, The University of Texas at
Dallas, Richardson, Texas 75080, USA

###### Abstract

We study the effect of a terahertz field-driven single cavity mode for
ultrafast control of a fermion chain with dissipation-induced nonlinearity and
quadratic coupling to an infrared-active phonon mode. Without photon loss from
the cavity, we uncover a first-order phase transition in the nonequilibrium
steady state only for the lower phonon-polariton, accompanied by polaritons
whose frequency response is asymmetric with respect to the photon frequency
due to the direct laser-induced dressing effect on the photon. A weak laser
field fails to induce the phase transition but renders the polaritons
symmetrical. Finally, we show that sufficiently strong photon loss from the
cavity eliminates the polaritons and the associated phase transition. The
experimental feasibility of these phenomena is also proposed.

Introduction.—Laser fields [1, 2, 3, 4, 5, 6] may provide control over quantum
materials, possibly revealing phenomena such as light-induced topology [7, 8],
laser-induced magnetism [9, 10], ferroelectric phase-transitions [11], and
superconductivity [12]. Classical light fields often introduce challenges
associated with heating which may be circumvented via the use of optical
cavities [13, 14, 15, 16] to extend the lifetime of laser-induced states.
Placing laser-driven systems in optical cavities not only establishes steady
states but also extends their preservation by suppressing disruptive heating
effects [13, 14, 15, 16, 17, 18]. The effectiveness of light-matter coupling
in terahertz cavities, where quantum fluctuations play a substantial role in
influencing the matter degrees of freedom, is of great importance for tuning
the lifetime of nonequilibrium states [19].

Understanding steady-state behavior requires systematic consideration of
dissipation [20, 21, 10, 22], affecting the balance between injected laser
energy and energy loss to the environment. This interplay facilitates the
emergence of intriguing nonthermal states of matter. In a recent study [21],
we uncovered a novel first-order phase transition in the steady state of a
fermion chain coupled to laser-driven infrared (IR)-active phonons. This
transition arises from the dressed dispersion of electrons in dissipation
processes. For practical applications in laser-driven switching and control
processes [23, 24, 25], tuning the timescale of the transient process for
achieving nonequilibrium steady states (NESS) becomes crucial. In light of the
significant role of the cavity in adjusting the timescale of steady states in
a driven-dissipative material, placing the above fermion chain into an optical
cavity [26] emerges as a promising approach to engineering the system’s
dynamic responses.

This Letter addresses three questions. First, with the inclusion of laser-
induced transient excitations into the dissipation process in a fermion chain
[21], how can we manipulate the NESS by coupling electrons and phonons to a
single cavity mode? Second, does the observed first-order phase transition in
the cavity-free model persist in the coupled cavity-fermion chain, or do we
observe different phase transitions? Lastly, what effect does photon loss from
the cavity have on the dynamic characteristics of the system? In answering
these questions, we discover a first-order phase transition only specific to
lower phonon-polaritons, linked to photon’s asymmetric responses, a
consequence of dressing effects induced by the laser. Our findings reveal that
a strong photon loss from the cavity destroys both the phase transition and
polaritons.

Model.—The model, Fig. 1, comprises of seven terms, including spinless
electrons (assuming spin-polarized bands for simplicity), phonons, a single
photon mode, electron-phonon coupling, phonon-photon coupling, electron-photon
coupling, and the interaction between the laser and the single photon mode
[18, 27, 22, 21, 28]. We employ the dipole approximation for long-wavelength
optical cavities, neglecting the spatial dependence of the electromagnetic
vector potential [18, 27].

The total time-dependent Hamiltonian (we set $\hbar=1$) is:

$\displaystyle\mathcal{H}(t)=-t_{0}\sum_{\ell}\big{(}c^{\dagger}_{\ell}c_{\ell+1}+{\rm
H.c.}\big{)}+\omega_{0}\sum_{\ell}a^{\dagger}_{\ell}a_{\ell}$ (1)
$\displaystyle+\eta b^{\dagger}b+g_{\rm
q}\sum_{\ell}\big{(}a^{\dagger}_{\ell}+a_{\ell}\big{)}^{2}\Big{(}c^{\dagger}_{\ell}c_{\ell}-\frac{1}{2}\Big{)}$
$\displaystyle-\frac{i\Delta}{\sqrt{L}}(b^{\dagger}+b)\sum_{\ell}(a^{\dagger}_{\ell}-a_{\ell})+\frac{2\Delta^{2}}{L\omega_{0}}(b^{\dagger}+b)^{2}$
$\displaystyle-
t_{0}\sum_{\ell}\big{(}e^{-i\frac{\lambda}{\sqrt{L}}(b^{\dagger}+b)}c^{\dagger}_{\ell}c_{\ell+1}+{\rm
H.c.}\big{)}+\frac{\mathcal{E}(t)}{\sqrt{L}}\big{(}b^{\dagger}+b\big{)}\,,$

which incorporates the electronic nearest-neighbor hopping amplitude $t_{0}$,
optical phonon ($\omega_{0}$), and photon ($\eta$) frequencies. Operators
$c^{\dagger}_{\ell}~{}(c_{\ell})$, $a^{\dagger}_{\ell}~{}(a_{\ell})$, and
$b^{\dagger}~{}(b)$ create (annihilate) electrons, phonons, and photon at site
$\ell\in\\{1,\dots,L\\}$ in a chain of length $L$. Quadratic electron-phonon
coupling is denoted by $g_{\rm q}$, and $\lambda$ represents the dimensionless
electron-photon coupling strength.
$\Delta=\alpha\,\eta\sqrt{\omega_{0}/8\eta}$ characterizes the phonon-photon
coupling with $\alpha\approx 0.02$ to 0.2 being the polariton-to-photon
frequency ratio in our simulations [18, 27, 28]. To coherently drive the
cavity, we use a steady laser field
$\mathcal{E}(t)=\mathcal{A}_{0}\cos(\omega\,t)$ with frequency $\omega$ and
amplitude $\mathcal{A}_{0}$.

Figure 1: A fermion chain in a laser-driven optical cavity. The model depicts
a driven cavity-fermion chain with a laser described as a classical field
$\mathcal{E}(t)$, driving a resonant single photon mode $\eta$ in the cavity.
The chain includes electrons that hop with energy $t_{0}$ and optical IR-
active phonons $\omega_{0}$, where $g_{\rm q}$ denotes quadratic electron-
phonon coupling. The electron-photon (phonon-photon) coupling strength is
denoted by the parameter $\lambda$ ($\Delta$). Damping of the driven photon,
phonon, and electron to the internal bath (a combination of other phonons and
an electromagnetic field of the non-driven photons) is represented by the
rates $\gamma_{\rm pt}$, $\gamma_{\rm pn}$, and $\gamma_{\rm e}$,
respectively. Additionally, the damping of the photon to the external bath
(vacuum) from the cavity is defined by $\kappa$.

Our model serves as a minimal descriptor of the ultrafast dynamics induced by
dissipation and nonlinear electron-phonon coupling in a terahertz field-driven
cavity and may be relevant to materials with symmetry constraints on their
electron-phonon coupling if placed in a cavity, such as ET-F2TCNQ,
(Pb,Bi)2Sr2CaCu2O8, YBa2Cu3O6+x, K3C60, and $\kappa$-salt [29, 30, 31, 5, 6,
32, 33, 34]. Application of our model to these materials requires the
inclusion of other material-specific details, which we leave to future work.

We study the dynamics for small values of the electron-photon coupling [19] in
compliance with the dipole approximation, first, and second, to have a
negligible energy density effect on the electrons. We focus only on the zero-
momentum mode of the dispersionless phonon $a_{0}$, as it is the only mode
that couples to electrons or photons in our translation-invariant model under
the dipole approximation. This comes from the fact that phonons are generated
by atoms spaced by order of an Å, while the laser wavelengths are on the order
of nanometers. Since the relative quantum fluctuations in phonon occupation
approach zero in the thermodynamic limit, the mean-field decoupling
approximation is then employed for describing the interaction between
particles [22].

To consider dissipation across all particles, we introduce an internal bath (a
combination of other phonons of the chain and the electromagnetic field inside
the cavity from non-driven cavity photon modes) and an external vacuum bath
from the cavity to simulate photon loss. This is implemented using the
Lindblad quantum master equation framework, incorporating phenomenological
damping rates [35, 36]. The equation governing coherent evolution and the
dissipator for any arbitrary observable $O$ is then expressed as
$\langle\dot{O}\rangle(t)=i\langle[\mathcal{H},O(t)]\rangle+\sum_{j}\Gamma_{j}\Big{<}\mathcal{L}_{j}^{\dagger}O(t)\mathcal{L}_{j}-\frac{1}{2}\big{\\{}\mathcal{L}_{j}^{\dagger}\mathcal{L}_{j},O(t)\big{\\}}\Big{>}$.
The Lindblad jump operators $\mathcal{L}_{j}=O$ and $O^{\dagger}$, with $O$
being any of $\\{a,b,c\\}$ operators [21], have corresponding damping
parameters $\Gamma_{j}=\gamma_{\rm pn}$, $\gamma_{\rm pt}$, and $\gamma_{\rm
e}$ for the phonon, photon, and electron subsystems inside the cavity.
Additionally, $\Gamma_{j}=\kappa$ represents the damping parameter for cavity
photon loss.

By incorporating laser-induced transient excitations into the dissipation
processes, the ratio between the damping of two dissipators in the electronic
subsystem results in a modulated dispersion $\widetilde{\varepsilon}_{k}(t)$
[35, 36]. Modes $k$ in the vicinity of the Fermi level $k_{\rm F}$ play a key
role in the laser modulation of electron occupation, influenced by electron-
photon and electron-phonon couplings. Thus, by expanding the electronic
dispersion around the Fermi level, we obtain

$\widetilde{\varepsilon}_{k}(t)\approx-2\,t_{0}(k-k_{\rm F})+g_{\rm
q}q^{2}_{\rm pn}(t)+2\,t_{0}\,\lambda\,q_{\rm pt}(t)\,,$ (2)

where $q_{\rm pt}(t)=L^{-1/2}\langle b^{\dagger}+b\rangle(t)$ and $q_{\rm
pn}(t)=L^{-1/2}\langle a^{\dagger}_{0}+a_{0}\rangle(t)$ refer to the photon
and phonon displacements, stemming from Eq. (1), while keeping the linear
order of $\lambda$.

Results.—Among various observables in the system, we focus on the occupations
of subsystems; $n_{\rm pt}(t)=L^{-1}\langle b^{\dagger}b\rangle(t)$, $n_{\rm
pn}(t)=L^{-1}\langle a^{\dagger}_{0}a_{0}\rangle(t)$, and $n_{\rm
e}(t)=L^{-1}\sum_{k}\langle c^{\dagger}_{k}c_{k}\rangle(t)$ and derive
corresponding equations of motion. To achieve sufficient accuracy of time
propagation for a true NESS, we numerically solve the equations of motion up
to a sufficiently long time, see Sec. S1 of the Supplemental Material (SM)
[37].

We consider a resonantly excited photon mode through laser interaction,
setting $\omega=\eta$. Additionally, we resonantly couple phonons to the
excited photon mode, ensuring that the laser field indirectly influences the
phonons (i.e., we set $\omega_{0}=\eta$). Previous studies on $\kappa$-salt
[34, 18, 27] demonstrated that the effective phonon frequency is
$\omega_{0}\approx 2\,t_{0}$, positioned near the electronic band edge. As
previously demonstrated [21], our mean-field-type approximation results in a
minor deviation of the electron occupation from its equilibrium value of 1/2
in the cavity-free model, maintaining the model stability. This allows for the
arbitrary setting of strengths for $g_{\rm q}$. In our system, the most
important parameter is the effective phonon-photon coupling $\Delta$, which is
responsible for transferring laser energy to the IR-active phonon, first, and
second, to the electron. For the majority of parameters, we reach a NESS, as
detailed in Sec. S2 of the SM.

Figure 2: Asymmetrical nonequilibrium first-order phase transition. The
dynamics of occupied steady states for (a) photon, (b) phonon, and (c)
electron are depicted across various phonon-photon couplings ($\Delta$) and
driving frequencies ($\omega$). The laser amplitude is set to
$\mathcal{A}_{0}/\gamma_{\rm pt}=0.25$. Due to dissipation-induced
nonlinearity in the relaxation process of electrons [21], leading to a
reduction in the number of electrons in (c), shift and splitting of resonances
in photon and phonon in (a) and (b) take place. This process facilitates the
first-order phase transition, as highlighted in inset panels, and gives rise
to the creation of hybridized states – phonon-polaritons – particularly for
$\Delta>0.05\,\eta$ values. The axes in (c) are swapped to increase the
visibility of jumps at the phase transitions. Panel (d) shows the time
evolution of dressed phonon occupation away from resonance for different
driving frequencies at $\Delta/\eta=0.088$, confirming the transition.

In addition to the provided data, Sec. S3 of the SM offers supplementary
information on various $g_{\rm q}$ and $\lambda$, confirming that the chosen
values $(g_{\rm q}/\omega_{0},\lambda)=(2,0.1)$ are proper choices following
Ref. [21]. Further, we proceed with weak damping rates for all degrees of
freedom; throughout the paper, we fix $\gamma_{\rm pn}/\omega_{0}=0.05$,
$\gamma_{\rm pt}/\eta=0.05$, and $\gamma_{\rm e}/t_{0}=0.01$. While finite
electron damping $\gamma_{\rm e}\neq 0$ is necessary to reflect the transient
feedback effects between phonons/photons and electrons in the relaxation
processes, its strength is not important. For further details on the effect of
phenomenological photon and phonon damping rates, see Sec. S4 of the SM. As
for the laser amplitude, we set $\mathcal{A}_{0}/\gamma_{\rm pt}\ll 1$ to be
weak enough to avoid system heating; we normalize it to the photon damping
because a driven damped harmonic oscillator has an occupation of
$(\mathcal{A}_{0}/\gamma_{\rm pt})^{2}$ in the decoupled phase, which we will
use as a benchmark for the coupled phase [22]. Finally, to eliminate finite-
size effects, we fix $L=2001$ for the chain length.

In what follows, we consider two scenarios: one without photon loss
($\kappa=0$) from the cavity and one with photon loss ($\kappa\neq 0$). To
investigate the properties of the NESS, we compute the average of temporal
coherent signals over one period.

We start with the effect of $\Delta$ on occupations in the NESS in Fig.
2(a-c). Panel (d) specifically considers $\Delta/\eta=0.088$. When there is no
phonon-photon coupling ($\Delta=0$), all entities peak at the resonance
$\omega=\omega_{0}=\eta$. However, once the photon interacts with the phonon
and electron, resonances experience shifts and splittings, attributed to the
formation of a “phonon-polariton”. This observation aligns well with findings
in Refs. [27, 38, 18]. Nevertheless, the behavior of the lower polariton in
photon occupation differs from that in phonon and electron occupations,
resulting in an asymmetrical presence in the photon occupation. The primary
reason for this distinction lies in the dressing effect in the model, where
the electric field directly influences the photon occupation through
$\widetilde{\mathcal{E}}(t)\approx{}\mathcal{E}(t)-\Delta\,p_{\rm pn}(t)$,
where $p_{\rm pn}(t)=i\,L^{-1/2}\langle a^{\dagger}_{0}-a_{0}\rangle(t)$
refers to the phonon momentum. In contrast, the strength of this effect is
comparatively weaker in the electron and photon occupations, where they are
indirectly influenced by the laser field.

Remarkably, simultaneous with the formation of asymmetrical and symmetrical
branches in photon and phonon/electron occupations during phonon-polariton
development, a notable sharp jump is observed only for the lower polariton at
various critical drive frequencies and $\Delta>0.05\,\eta$, as illustrated in
the inset panels. The time evolution of the photon occupation in Fig. 2(d)
further supports the occurrence of a phase transition at $\Delta/\eta=0.088$
and $\omega/\eta\approx 0.867$. Indeed, at the critical drive frequency, the
phonon occupation initially surges, reaching a quasi-divergence, but
immediately descends to another NESS plateau (the highest value). The
emergence of quasi-divergence signifies the first-order phase transition. This
represents a novel theoretical finding regarding phonon-polaritons out of
equilibrium that has not been uncovered previously.

By leveraging the dominant harmonics of oscillations in the NESS for all
entities, we analytically confirm that the singularities in phonon
displacement emerge at drive frequencies smaller than single-photon frequency,
$\omega<\eta$, see Sec. S5 of the SM. This arises from the solutions of a
polynomial equation for $|q^{\rm pn}_{1}|$ (the first harmonic of $q_{\rm
pn}(t)$), wherein its minimum reaches zero. From a physical standpoint, a zero
minimum of phonon and photon displacements means that the objects are
instantaneously frozen at the critical drive frequencies for
$\Delta>0.05\,\eta$. The first-order phase transition is a direct consequence
of such freezing-like effects.

Figure 3: Elimination of nonequilibrium first-order phase transition with weak
laser field. The role of laser amplitude on the first-order phase transition
and the symmetrical responses of lower and upper polaritons in the NESS is
illustrated for (a) photon, (b) phonon, and (c) electron at
$\Delta/\eta\approx 0.1$. As $\mathcal{A}_{0}$ decreases, the lower and upper
branches of the polariton display symmetry in the photon sector, and the
first-order phase transition disappears in the entire system. This happens
because of the reduced dressing effect caused by the laser field.

In the following analysis, we investigate the robustness of the phase
transition and symmetrical responses of lower and upper polaritons, depicted
in Fig. 3, under various laser amplitudes $\mathcal{A}_{0}$. In experimental
settings, laser amplitude plays a crucial role in unveiling unique behaviors
and potentially tracking phase transitions. As the laser amplitude decreases,
the lower and upper branches of the phonon-polaritons symmetrically align
around the resonance peak at $\omega/\eta=1$, as seen in the inset panel of
Fig. 3(a). The rationale behind this remains consistent with previous
explanations; the weakening of the dressing effect with smaller
$\mathcal{A}_{0}$ values causes the photon occupation to behave similarly to
the phonon and electron. Moreover, we observe that such a weak drive field
leads to the revocation of the first-order phase transition. This implies that
the minimum values of phonon and photon displacements persist as finite
without vanishing if the laser amplitude is sufficiently weak. This suggests
that under such conditions, these objects never experience instantaneous
freezing in the system, in contrast to the regime with strong
$\mathcal{A}_{0}$.

Figure 4: Elimination of nonequilibrium first-order phase transition and
phonon-polariton with strong photon loss. The influence of photon loss
$\kappa$ on the dynamical responses of phonon-polaritons in the NESS of the
(a) photon, (b) phonon, and (c) electron at $\mathcal{A}_{0}/\gamma_{\rm
pt}=0.25$ and $\Delta/\eta\approx 0.1$. When photon loss $\kappa/\eta$
surpasses 0.03 and 0.2, the first-order phase transition and polaritons
vanish, respectively, due to the blocking of input power from entering the
cavity, as it is dissipated to the vacuum through photon loss.

Up to this point, we have considered a scenario where the system lacks photon
loss from the cavity, leading to the formation of phonon-polaritons and the
occurrence of the first-order phase transition under specific parameter
conditions. For the sake of completeness, in Fig. 4, our focus turns on the
effect of photon loss, addressing how the aforementioned dynamic features
respond to changes in $\kappa$. In this scenario, we observe that for
$\kappa/\eta>0.03$, the first-order phase transition does not occur for our
choice of remaining parameters. Additionally, the presence of the phonon-
polariton is noticeable up to $\kappa/\eta\approx 0.2$, but beyond this
threshold, polaritons fail to form. This can be understood by noting that with
the inclusion of photon loss $\kappa$, a significant portion of the input
laser power is dissipated to the vacuum before transferring energy to the
single cavity mode and other components. Essentially, the resonance response
of the system remains robust against couplings when exceeding
$\kappa/\eta=0.2$.

Experimental perspective.—While the details of the couplings strongly depend
on the actual material at hand, particularly the electron-phonon coupling, we
can propose a physical realization of the first-order phase transition. Our
findings indicate that the interaction between phonons and photons plays a key
role in determining novel phenomena, motivating us to focus on the phonon as
an easily accessible lattice degree of freedom in reality. By driving phonons
through a terahertz field-driven optical photon, a portion of the absorbed
laser power, transitioning from the excited phonon to the environment inside
the cavity, can be described by $I(t)=\mathcal{M}\omega_{0}\gamma_{\rm
ph}n_{\rm ph}(t)$, where $\mathcal{M}$ represents the mole amount of thin film
of our sample determined by its volume and molar density. To mitigate thermal
effects, the sample can be placed in contact with a cold sink (characterized
by temperature $T$, mass $m$, and a small heat capacity $C(T)\propto T$). With
this, we can estimate the timescale at which the system maintains the NESS to
observe phase transition, given by $t_{\rm
exp}=m\int_{0}^{T}C(T)\,d\,T/\overline{I}$. For instance, consider a chain
with hopping energy of $t_{0}\approx 1$ THz, where the phonon frequency
$\omega_{0}\approx 2\,t_{0}$ in our simulations lies near the electronic band
edge. For a chain with approximately $1\%$ phonon occupation per site (see
panel (b) of above figures), a cold sink possessing a mass of 2 g and a
temperature of 2 K, the NESS can be maintained for up to $t_{\rm exp}\approx
500$ ns, providing sufficient time to observe the phase transition in the
NESS. This is reachable using existing terahertz laser sources [39, 40, 41,
42].

Summary and outlook.—The occurrence of a nonequilibrium first-order phase
transition stems from nonlinear dissipative dynamics in a driven fermion
chain, where the electron density is coupled quadratically to an infrared-
active phonon. Our current investigation unveiled that a driven optical cavity
can finely adjust the ultrafast dynamics of particles in the non-equilibrium
steady states (NESS) of such a chain. To analyze this system, we employed the
Lindblad formalism, dipole approximation, and mean-field theory. Laser
excitation of a single photon mode gives rise to hybridized states,
particularly phonon-polaritons in the NESS primarily formed by the coupling
between the driven photon and phonon.

We determine that, in the absence of photon loss from the cavity, only the
lower phonon-polariton experiences the first-order phase transition in the
NESS, accompanied by the signature of the terahertz field-induced dressing
effect in the (a)symmetrical responses of polaritons. Our findings suggest
that a laser field lacking sufficient strength is incapable of triggering the
first-order phase transition, but it does restore the symmetry of the lower
and upper polaritons. Additionally, photon loss disrupts the phonon-polariton
and its distinctive phase transition characteristic. These results bear
significance for laser-driven switching and control processes in driven open
quantum systems.

We briefly comment that a chirp drive field can effectively modulate all these
dynamic features via hysteresis processes, enabling the manifestation of a
phase transition for an alternative set of parameters as well [21, 10]. For
example, it remains interesting to explore whether the chirp can prompt a
first-order phase transition to the upper polariton of our current model.
Moreover, the induction of a first-order phase transition in the presence of
dissipation and the nonlinear coupling of electrons to photon-driven phonons
in a cavity raises important questions. Specifically, it prompts inquiries
about the cavity-mediated induction of phonon-induced interactions between
electrons with spin degrees of freedom for real materials in the cavity. For
example, it will be interesting to study the modalities afforded by a cavity
field in the presence of dissipation for the induction of coherent phase
transitions of the electrons such as cavity-mediated superconductivity and
density waves [43, 44, 45]. We will address these questions in future work.

Acknowledgments.—We greatly thank Edoardo Baldini for helpful discussions.
This work was performed with support from the National Science Foundation
(NSF) through award numbers MPS-2228725 and DMR-1945529. Part of this work was
performed at the Aspen Center for Physics, which is supported by NSF grant No.
PHY-1607611. Funded by the European Union (ERC, QuSimCtrl, 101113633). Views
and opinions expressed are however those of the authors only and do not
necessarily reflect those of the European Union or the European Research
Council Executive Agency. Neither the European Union nor the granting
authority can be held responsible for them.

## References

  * Först  _et al._ [2011] M. Först, C. Manzoni, S. Kaiser, Y. Tomioka, Y. Tokura, R. Merlin, and A. Cavalleri, Nonlinear phononics as an ultrafast route to lattice control, Nature Physics 7, 854 (2011).
  * Mitrano  _et al._ [2016] M. Mitrano, A. Cantaluppi, D. Nicoletti, S. Kaiser, A. Perucchi, S. Lupi, P. Di Pietro, D. Pontiroli, M. Riccò, S. R. Clark, D. Jaksch, and A. Cavalleri, Possible light-induced superconductivity in K3C60 at high temperature, Nature 530, 461 (2016).
  * Knap  _et al._ [2016] M. Knap, M. Babadi, G. Refael, I. Martin, and E. Demler, Dynamical cooper pairing in nonequilibrium electron-phonon systems, Phys. Rev. B 94, 214504 (2016).
  * Sentef  _et al._ [2017] M. A. Sentef, A. Tokuno, A. Georges, and C. Kollath, Theory of laser-controlled competing superconducting and charge orders, Phys. Rev. Lett. 118, 087002 (2017).
  * Kennes  _et al._ [2017] D. M. Kennes, E. Y. Wilner, D. R. Reichman, and A. J. Millis, Transient superconductivity from electronic squeezing of optically pumped phonons, Nature Physics 13, 479 (2017).
  * Sentef [2017] M. A. Sentef, Light-enhanced electron-phonon coupling from nonlinear electron-phonon coupling, Phys. Rev. B 95, 205111 (2017).
  * Claassen  _et al._ [2016] M. Claassen, C. Jia, B. Moritz, and T. P. Devereaux, All-optical materials design of chiral edge modes in transition-metal dichalcogenides, Nature Communications 7, 13074 (2016).
  * Hübener  _et al._ [2017] H. Hübener, M. A. Sentef, U. De Giovannini, A. F. Kemper, and A. Rubio, Creating stable Floquet–Weyl semimetals by laser-driving of 3d Dirac materials, Nature Communications 8, 13940 (2017).
  * Shin  _et al._ [2018] D. Shin, H. Hübener, U. De Giovannini, H. Jin, A. Rubio, and N. Park, Phonon-driven spin-Floquet magneto-valleytronics in MoS2, Nature Communications 9, 638 (2018).
  * Yarmohammadi  _et al._ [2023a] M. Yarmohammadi, M. Bukov, V. Oganesyan, and M. H. Kolodrubetz, Laser-enhanced magnetism in SmFeO3 arXiv:2308.14939 (2023a).
  * Ashida  _et al._ [2020] Y. Ashida, A. m. c. İmamoğlu, J. Faist, D. Jaksch, A. Cavalleri, and E. Demler, Quantum electrodynamic control of matter: Cavity-enhanced ferroelectric phase transition, Phys. Rev. X 10, 041027 (2020).
  * Thomas  _et al._ [2019] A. Thomas, E. Devaux, K. Nagarajan, T. Chervy, M. Seidel, D. Hagenmüller, S. Schütz, J. Schachenmayer, C. Genet, G. Pupillo, and T. W. Ebbesen, Exploring superconductivity under strong coupling with the vacuum electromagnetic field arXiv:1911.01459 (2019).
  * Ruggenthaler  _et al._ [2018] M. Ruggenthaler, N. Tancogne-Dejean, J. Flick, H. Appel, and A. Rubio, From a quantum-electrodynamical light–matter description to novel spectroscopies, Nature Reviews Chemistry 2, 0118 (2018).
  * Frisk Kockum  _et al._ [2019] A. Frisk Kockum, A. Miranowicz, S. De Liberato, S. Savasta, and F. Nori, Ultrastrong coupling between light and matter, Nature Reviews Physics 1, 19 (2019).
  * Hübener  _et al._ [2021] H. Hübener, U. De Giovannini, C. Schäfer, J. Andberger, M. Ruggenthaler, J. Faist, and A. Rubio, Engineering quantum materials with chiral optical cavities, Nature Materials 20, 438 (2021).
  * Genet  _et al._ [2021] C. Genet, J. Faist, and T. W. Ebbesen, Inducing new material properties with hybrid light–matter states, Physics Today 74, 42 (2021).
  * Schlawin  _et al._ [2022] F. Schlawin, D. M. Kennes, and M. A. Sentef, Cavity quantum materials, Applied Physics Reviews 9, 011312 (2022).
  * Dé  _et al._ [2022] B. L. Dé, C. J. Eckhardt, D. M. Kennes, and M. A. Sentef, Cavity engineering of Hubbard u via phonon polaritons, Journal of Physics: Materials 5, 024006 (2022).
  * Li  _et al._ [2018] X. Li, M. Bamba, Q. Zhang, S. Fallahi, G. C. Gardner, W. Gao, M. Lou, K. Yoshioka, M. J. Manfra, and J. Kono, Vacuum Bloch–Siegert shift in Landau polaritons with ultra-high cooperativity, Nature Photonics 12, 324 (2018).
  * Yarmohammadi  _et al._ [2023b] M. Yarmohammadi, M. Krebs, G. S. Uhrig, and B. Normand, Strong-coupling magnetophononics: Self-blocking, phonon-bitriplons, and spin-band engineering, Phys. Rev. B 107, 174415 (2023b).
  * Yarmohammadi  _et al._ [2023c] M. Yarmohammadi, M. Bukov, and M. H. Kolodrubetz, Nonequilibrium phononic first-order phase transition in a driven fermion chain, Phys. Rev. B 108, L140305 (2023c).
  * Yarmohammadi  _et al._ [2021] M. Yarmohammadi, C. Meyer, B. Fauseweh, B. Normand, and G. S. Uhrig, Dynamical properties of a driven dissipative dimerized ${S}=\frac{1}{2}$ chain, Phys. Rev. B 103, 045132 (2021).
  * Parag  _et al._ [2015] A. K. Parag, L. I. Giri, and R. Basu, Design and construction of a laser switch based electrical control circuit, in  _2015 6th International Conference on Computers and Devices for Communication (CODEC)_ (2015) pp. 1–4.
  * Kuzikova  _et al._ [2023] A. V. Kuzikova, L. A. Shelukhin, F. M. Maksimov, A. I. Chernov, R. V. Pisarev, and A. M. Kalashnikova, Laser-driven first-order spin reorientation and verwey phase transitions in magnetite Fe3O4 beyond the range of thermodynamic equilibrium, Phys. Rev. B 107, 024413 (2023).
  * Boschker  _et al._ [2015] J. E. Boschker, R. N. Wang, V. Bragaglia, P. Fons, A. Giussani, L. L. Guyader, M. Beye, I. Radu, A. V. Kolobov, K. Holldack, and R. Calarco, Laser-driven switching dynamics in phase change materials investigated by time-resolved x-ray absorption spectroscopy, Phase Transitions 88, 82 (2015).
  * Mazza and Georges [2019] G. Mazza and A. Georges, Superradiant quantum materials, Phys. Rev. Lett. 122, 017401 (2019).
  * Eckhardt  _et al._ [2022] C. J. Eckhardt, G. Passetti, M. Othman, C. Karrasch, F. Cavaliere, M. A. Sentef, and D. M. Kennes, Quantum Floquet engineering with an exactly solvable tight-binding chain in a cavity, Communications Physics 5, 122 (2022).
  * Mahan [2013] G. Mahan,  _Many-Particle Physics_, Physics of Solids and Liquids (Springer US, 2013).
  * Kaiser  _et al._ [2014] S. Kaiser, S. R. Clark, D. Nicoletti, G. Cotugno, R. I. Tobey, N. Dean, S. Lupi, H. Okamoto, T. Hasegawa, D. Jaksch, and A. Cavalleri, Optical properties of a vibrationally modulated solid state Mott insulator, Scientific Reports 4, 3823 (2014).
  * Singla  _et al._ [2015] R. Singla, G. Cotugno, S. Kaiser, M. Först, M. Mitrano, H. Y. Liu, A. Cartella, C. Manzoni, H. Okamoto, T. Hasegawa, S. R. Clark, D. Jaksch, and A. Cavalleri, THz-frequency modulation of the Hubbard ${U}$ in an organic Mott insulator, Phys. Rev. Lett. 115, 187401 (2015).
  * Kim  _et al._ [2003] T. K. Kim, A. A. Kordyuk, S. V. Borisenko, A. Koitzsch, M. Knupfer, H. Berger, and J. Fink, Doping dependence of the mass enhancement in (Pb,Bi)2Sr2CaCu2O8 at the antinodal point in the superconducting and normal states, Phys. Rev. Lett. 91, 167002 (2003).
  * Puviani and Sentef [2018] M. Puviani and M. A. Sentef, Quantum nonlinear phononics route towards nonequilibrium materials engineering: Melting dynamics of a ferrielectric charge density wave, Phys. Rev. B 98, 165138 (2018).
  * Sous  _et al._ [2021] J. Sous, B. Kloss, D. M. Kennes, D. R. Reichman, and A. J. Millis, Phonon-induced disorder in dynamics of optically pumped metals from nonlinear electron-phonon coupling, Nature Communications 12, 5803 (2021).
  * Buzzi  _et al._ [2020] M. Buzzi, D. Nicoletti, M. Fechner, N. Tancogne-Dejean, M. A. Sentef, A. Georges, T. Biesner, E. Uykur, M. Dressel, A. Henderson, T. Siegrist, J. A. Schlueter, K. Miyagawa, K. Kanoda, M.-S. Nam, A. Ardavan, J. Coulthard, J. Tindall, F. Schlawin, D. Jaksch, and A. Cavalleri, Photomolecular high-temperature superconductivity, Phys. Rev. X 10, 031028 (2020).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Comm. Math. Phys. 48, 119 (1976).
  * Breuer and Petruccione [2007] H. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (OUP Oxford, 2007).
  * [37] See the Supplemental Materials at http://xxx.yyy.zzz for further details of the equations of motion and photonic, phononic, and electronic responses out of equilibrium .
  * Sentef  _et al._ [2018] M. A. Sentef, M. Ruggenthaler, and A. Rubio, Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity, Science Advances 4, 10.1126/sciadv.aau6969 (2018).
  * Biasco  _et al._ [2018] S. Biasco, K. Garrasi, F. Castellano, L. Li, H. E. Beere, D. A. Ritchie, E. H. Linfield, A. G. Davies, and M. S. Vitiello, Continuous-wave highly-efficient low-divergence terahertz wire lasers, Nature Communications 9, 1122 (2018).
  * Khalatpour  _et al._ [2021] A. Khalatpour, A. K. Paulsen, C. Deimert, Z. R. Wasilewski, and Q. Hu, High-power portable terahertz laser systems, Nature Photonics 15, 16 (2021).
  * Liebermeister  _et al._ [2021] L. Liebermeister, S. Nellen, R. B. Kohlhaas, S. Lauck, M. Deumer, S. Breuer, M. Schell, and B. Globisch, Optoelectronic frequency-modulated continuous-wave terahertz spectroscopy with 4 THz bandwidth, Nature Communications 12, 1071 (2021).
  * Yang  _et al._ [2023] C.-J. Yang, J. Li, M. Fiebig, and S. Pal, Terahertz control of many-body dynamics in quantum materials, Nature Reviews Materials 8, 518 (2023).
  * Allocca  _et al._ [2019] A. A. Allocca, Z. M. Raines, J. B. Curtis, and V. M. Galitski, Cavity superconductor-polaritons, Phys. Rev. B 99, 020504 (2019).
  * Curtis  _et al._ [2019] J. B. Curtis, Z. M. Raines, A. A. Allocca, M. Hafezi, and V. M. Galitski, Cavity quantum eliashberg enhancement of superconductivity, Phys. Rev. Lett. 122, 167002 (2019).
  * Ruggenthaler  _et al._ [2023] M. Ruggenthaler, D. Sidler, and A. Rubio, Understanding polaritonic chemistry from ab initio quantum electrodynamics, Chemical Reviews 123, 11191 (2023).


* Huang  _et al._ [2017] B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero, and X. Xu, Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit, Nature (London) 546, 270 (2017).
  * Fei  _et al._ [2018] Z. Fei, B. Huang, P. Malinowski, W. Wang, T. Song, J. Sanchez, W. Yao, D. Xiao, X. Zhu, A. F. May, W. Wu, D. H. Cobden, J.-H. Chu, and X. Xu, Two-dimensional itinerant ferromagnetism in atomically thin Fe3GeTe2, Nat. Mater. 17, 778 (2018).
  * Deng  _et al._ [2018] Y. Deng, Y. Yu, Y. Song, J. Zhang, N. Z. Wang, Z. Sun, Y. Yi, Y. Z. Wu, S. Wu, J. Zhu, J. Wang, X. H. Chen, and Y. Zhang, Gate-tunable room-temperature ferromagnetism in two-dimensional Fe3GeTe2, Nature (London) 563, 94 (2018).
  * Xu  _et al._ [2018] C. Xu, J. Feng, H. Xiang, and L. Bellaiche, Interplay between Kitaev interaction and single ion anisotropy in ferromagnetic CrI3 and CrGeTe3 monolayers, npj Computational Materials 4, 57 (2018).
  * Mackintosh and Andersen [1980] A. R. Mackintosh and O. K. Andersen, The electronic structure of transition metals, in  _Electrons at the Fermi Surface_, edited by M. Springford (Cambridge University Press, Cambridge, 1980) pp. 149–224.
  * Heine [1980] V. Heine, Electronic structure from the point of view of the local atomic environment, in  _Solid State Physics_, Vol. 35, edited by H. Ehrenreich, F. Seitz, and D. Turnbull (Academic Press, New York, 1980) pp. 1–123.
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Daalderop  _et al._ [1990] G. H. O. Daalderop, P. J. Kelly, and M. F. H. Schuurmans, First-principles calculation of the magnetocrystalline anisotropy energy of iron, cobalt and nickel, Phys. Rev. B 41, 11919 (1990).
  * Daalderop  _et al._ [1992] G. H. O. Daalderop, P. J. Kelly, and F. J. A. den Broeder, Prediction and confirmation of perpendicular magnetic anisotropy in Co/Ni multilayers, Phys. Rev. Lett. 68, 682 (1992).
  * Daalderop  _et al._ [1994] G. H. O. Daalderop, P. J. Kelly, and M. F. H. Schuurmans, Magnetic anisotropy of a free-standing Co monolayer and of multilayers which contain Co monolayers, Phys. Rev. B 50, 9989 (1994).
  * Hong  _et al._ [2015] J. Hong, Z. Hu, M. Probert, K. Li, D. Lv, X. Yang, L. Gu, N. Mao, Q. Feng, L. Xie, J. Zhang, D. Wu, Z. Zhang, C. Jin, W. Ji, X. Zhang, J. Yuan, and Z. Zhang, Exploring atomic defects in molybdenum disulphide monolayers, Nat. Commun. 6, 6293 (2015).
  * Wang  _et al._ [2018] S. Wang, A. Robertson, and J. H. Warner, Atomic structure of defects and dopants in 2D layered transition metal dichalcogenides, Chem. Soc. Rev. 47, 6764 (2018).
  * Rasool  _et al._ [2015] H. I. Rasool, C. Ophus, and A. Zettl, Atomic defects in two dimensional materials, Adv. Mater. 27, 5771 (2015).
  * Komsa and Krasheninnikov [2015] H. Komsa and A. V. Krasheninnikov, Native defects in bulk and monolayer mos2 from first principles, Phys. Rev. B 91, 125304 (2015).
  * Li  _et al._ [2016a] W.-F. Li, C. Fang, and M. A. van Huis, Strong spin-orbit splitting and magnetism of point defect states in monolayer WS2, Phys. Rev. B 94, 195425 (2016a).
  * Haldar  _et al._ [2015] S. Haldar, H. Vovusha, M. K. Yadav, O. Eriksson, and B. Sanyal, A systematic study of structural, electronic and optical properties of atomic scale defects in 2D transition metal dichalcogenides MX2 (M = Mo,W; X = S, Se, Te), Phys. Rev. B 92, 235408 (2015).
  * Yang  _et al._ [2019b] S.-Z. Yang, W. Sun, Y.-Y. Zhang, Y. Gong, M. P. Oxley, A. R. Lupini, P. M. Ajayan, M. F. Chisholm, S. T. Pantelides, and W. Zhou, Direct Cation Exchange in Monolayer MoS2 via Recombination-Enhanced Migration, Phys. Rev. Lett. 122, 106101 (2019b).
  * Tanabe  _et al._ [2016] I. Tanabe, M. Gomez, W. C. Coley, D. Le, E. M. Echeverria, G. Stecklein, V. Kandyba, S. K. Balijepalli, V. Klee, A. E. Nguyen, E. Preciado, I.-H. Lu, S. Bobek, D. Barroso, D. Martinez-Ta, A. Barinov, T. S. Rahman, P. A. Dowben, P. A. Crowell, and L. Bartels, Band structure characterization of WS2 grown by chemical vapor deposition, Appl. Phys. Lett. 108, 252103 (2016).
  * Yeh  _et al._ [2015] P.-C. Yeh, W. Jin, N. Zaki, D. Zhang, J. T. Liou, J. T. Sadowski, A. Al-Mahboob, J. I. Dadap, I. P. Herman, P. Sutter, and R. M. Osgood, Jr., Layer-dependent electronic structure of an atomically heavy two-dimensional dichalcogenide, Phys. Rev. B 91, 041407 (2015).
  * Kübler [2000] J. Kübler,  _Theory of Itinerant Electron Magnetism_ (Oxford University Press, Oxford, 2000).
  * Mohn [2003] P. Mohn,  _Magnetism in the Solid State: An Introduction_ , edited by M. Cardona, P. Fulde, K. von Klitzing, and H.-J. Queisser, Springer Series in Solid-State Sciences, Vol. 134 (Springer, Berlin, 2003).
  * Lejaeghere  _et al._ [2016] K. Lejaeghere, G. Bihlmayer, T. Björkman, P. Blaha, S. Blügel, V. Blum, D. Caliste, I. E. Castelli, S. J. Clark, A. D. Corso, S. de Gironcoli, T. Deutsch, J. K. Dewhurst, I. Di Marco, C. Draxl, M. Dułak, O. Eriksson, J. A. Flores-Livas, K. F. Garrity, L. Genovese, P. Giannozzi, M. Giantomassi, S. Goedecker, X. Gonze, O. Grånäs, E. K. U. Gross, A. Gulans, F. Gygi, D. R. Hamann, P. J. Hasnip, N. A. W. Holzwarth, D. Iuşan, D. B. Jochym, F. Jollet, D. Jones, G. Kresse, K. Koepernik, E. Küçükbenli, Y. O. Kvashnin, I. L. M. Locht, S. Lubeck, M. Marsman, N. Marzari, U. Nitzsche, L. Nordström, T. Ozaki, L. Paulatto, C. J. Pickard, W. Poelmans, M. I. J. Probert, K. Refson, M. Richter, G.-M. Rignanese, S. Saha, M. Scheffler, M. Schlipf, K. Schwarz, S. Sharma, F. Tavazza, P. Thunström, A. Tkatchenko, M. Torrent, D. Vanderbilt, M. J. van Setten, V. Van Speybroeck, J. M. Wills, J. R. Yates, G.-X. Zhang, and S. Cottenier, Reproducibility in density functional theory calculations of solids, Science 351, aad3000 (2016).
  * Gunnarsson  _et al._ [1974] O. Gunnarsson, B. I. Lundqvist, and J. W. Wilkins, Contribution to the cohesive energy of simple metals: Spin-dependent effect, Phys. Rev. B 10, 1319 (1974).
  * Anisimov  _et al._ [1991] V. I. Anisimov, J. Zaanen, and O. K. Andersen, Band theory and mott insulators: Hubbard $u$ instead of stoner $i$, Phys. Rev. B 44, 943 (1991).
  * Solovyev  _et al._ [1994] I. V. Solovyev, P. H. Dederichs, and V. I. Anisimov, Corrected atomic limit in the local-density approximation and the electronic structure of $d$ impurities in Rb, Phys. Rev. B 50, 16861 (1994).
  * Liechtenstein  _et al._ [1995] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Density-functional theory and strong-interactions - orbital ordering in mott-hubbard insulators, Phys. Rev. B 52, 5467 (1995).
  * Bethe and Jackiw [1986] H. A. Bethe and R. Jackiw,  _Intermediate Quantum Mechanics_, 3rd ed., edited by D. Pines (Benjamin Cummings, Menlo Park, CA, USA, 1986).
  * Dolui  _et al._ [2013] K. Dolui, I. Rungger, C. Pemmaraju, and S. Sanvito, Possible doping strategies for MoS2 monolayers: An _a_ b initio study, Phys. Rev. B 88, 075420 (2013).
  * Andriotis and Menon [2014] A. N. Andriotis and M. Menon, Tunable magnetic properties of transition metal doped MoS2, Phys. Rev. B 90, 125304 (2014).
  * Yue  _et al._ [2013] Q. Yue, S. Chang, S. Qin, and J. Li, Functionalization of monolayer MoS2 by substitutional doping: A first-principles study, Phys. Lett. A 377, 1362 (2013).
  * Lu and Leburton [2014] S.-C. Lu and J.-P. Leburton, Electronic structures of defects and magnetic impurities in MoS2 monolayers, Nanoscale Research Letters 9, 676 (2014).
  * Singh and Schwingenschlögl [2016] N. Singh and U. Schwingenschlögl, Extended moment formation in monolayer WS2 doped with 3d transition-metals, ACS Applied Materials & Interfaces 8, 23886 (2016).
  * Li  _et al._ [2016b] H. Li, S. Liu, S. Huang, D. Yin, C. Li, and Z. Wang, Impurity-induced ferromagnetism and metallicity of WS2 monolayer, Ceramics International 42, 2364 (2016b).
  * Wu  _et al._ [2018] M. Wu, X. Yao, Y. Hao, H. Dong, Y. Cheng, H. Liu, F. Lu, W. Wang, K. Cho, and W.-H. Wang, Electronic structures, magnetic properties and band alignments of 3d transition metal atoms doped monolayer MoS2, Phys. Lett. A 382, 111 (2018).
  * Duong  _et al._ [2019] D. L. Duong, S. J. Yun, Y. Kim, S.-G. Kim, and Y. H. Lee, Long-range ferromagnetic ordering in vanadium-doped WSe2 semiconductor, Appl. Phys. Lett. 115, 242406 (2019).
  * Pan  _et al._ [2020] M. Pan, J. T. Mullen, and K. W. Kim, First-principles analysis of magnetically doped transition-metal dichalcogenides, J. Phys. D: Appl. Phys. 54, 025002 (2020).
  * Tiwari  _et al._ [2021] S. Tiwari, M. L. Van de Put, B. Sorée, and W. G. Vandenberghe, Magnetic order and critical temperature of substitutionally doped transition metal dichalcogenide monolayers, npj 2D Materials and Applications 5, 54 (2021).
  * Zhao  _et al._ [2016] X. Zhao, C. Xia, T. Wang, and X. Dai, Electronic and magnetic properties of X-doped (X = Ti, Zr, Hf) tungsten disulphide monolayer, J. Alloys Compd. 654, 574 (2016).
  * Wiesendanger [2009] R. Wiesendanger, Spin mapping at the nanoscale and atomic scale, Rev. Mod. Phys. 81, 1495 (2009).
  * Gil  _et al._ [2014] C. J. Gil, A. Pham, A. Yu, and S. Li, An _a_ b initio study of transition metals doped with WSe2 for long-range room temperature ferromagnetism in two-dimensional transition metal dichalcogenide, J. Phys.: Condens. Matter 26, 306004 (2014).
  * Mallet  _et al._ [2020] P. Mallet, F. Chiapello, H. Okuno, H. Boukari, M. Jamet, and J.-Y. Veuillen, Bound hole states associated to individual vanadium atoms incorporated into monolayer WSe2, Phys. Rev. Lett. 125, 036802 (2020).
  * Edwards and Katsnelson [2006] D. M. Edwards and M. I. Katsnelson, High-temperature ferromagnetism of sp electrons in narrow impurity bands: application to CaB6, J. Phys.: Condens. Matter 18, 7209 (2006).
  * Su  _et al._ [2019] J. Su, K. Liu, F. Wang, B. Jin, Y. Guo, G. Liu, H. Li, and T. Zhai, Van der Waals 2D transition metal tellurides, Advanced Materials Interfaces 6, 1900741 (2019).
  * Ataca  _et al._ [2012] C. Ataca, H. Şahin, and S. Ciraci, Stable, Single-Layer MX2 Transition-Metal Oxides and Dichalcogenides in a Honeycomb-Like Structure, J. Phys. Chem. C 116, 8983 (2012).
  * Janak [1978] J. F. Janak, Proof that $\partial e/\partial n_{i}=\epsilon_{i}$ in density-functional theory, Phys. Rev. B 18, 7165 (1978).
  * Jones and Gunnarsson [1989] R. O. Jones and O. Gunnarsson, The density functional formalism, its applications and prospects, Rev. Mod. Phys. 61, 689 (1989).