# Exceptional points and phase transitions in non-Hermitian binary systems

Amir Rahmani  Institute of Physics Polish Academy of Sciences, Al. Lotników
32/46, 02-668 Warsaw, Poland  Andrzej Opala  Institute of Physics Polish
Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland  Institute of
Experimental Physics, Faculty of Physics, University of Warsaw, ul. Pasteura
5, PL-02-093 Warsaw, Poland  Michał Matuszewski  Institute of Physics Polish
Academy of Sciences, Al. Lotników 32/46, 02-668 Warsaw, Poland

###### Abstract

Recent study demonstrated that steady states of a polariton system may
demonstrate a first-order dissipative phase transition with an exceptional
point that appears as an endpoint of the phase boundary [R. Hanai et al.,
Phys. Rev. Lett. 122, 185301 (2019)]. Here, we show that this phase transition
is strictly related to the stability of solutions. In general, the exceptional
point does not correspond to the endpoint of a phase transition, but rather it
is the point where stable and unstable solutions coalesce. Moreover, we show
that the transition may occur also in the weak coupling regime, which was
excluded previously. In a certain range of parameters, we demonstrate
permanent Rabi-like oscillations between light and matter fields. Our results
contribute to the understanding of nonequilibrium light-matter systems, but
can be generalized to any two-component oscillatory systems with gain and
loss.

Phase transitions correspond to significant alterations of the properties of a
system caused by the modification of physical parameters. Examples include the
ferromagnetic-paramagnetic phase transition Jiles (2015), gas-liquid-solid
transition Stishov (2018), Bose-Einstein condensation in bosonic and fermionic
systems Proukakis  _et al._ (2017), metal–insulator transition in solid state
Imada  _et al._ (1998), and topological phase transitions Brink  _et al._
(2018). Phase transitions may also occur in non-Hermitian systems, which are
systems that do not satisfy the condition of Hermiticity, which is embedded in
quantum mechanics Ashida  _et al._ (2020). Here the non-Hermitian
contributions may stem from dissipation Miri and Alù (2019) or asymmetric
coupling Gong  _et al._ (2018) and lead to a number of unique properties such
as non-reciprocity Peng  _et al._ (2022), mutually interlinked non-Hermitian
phase transitions Weidemann  _et al._ (2022) and the non-Hermitian skin effect
Okuma and Sato (2023).

A striking example of non-Hermitian physics that deviates significantly from
the Hermitian case is the coalescence of eigenstates and energy eigenvalues at
so-called exceptional points (EPs). These spectral singularities may be
accompanied by a non-Hermitian phase transition Öztürk  _et al._ (2021).
Standard procedure to investigate these phase transitions is through the study
of the spectrum of the system as some controllable parameters are changed Miri
and Alù (2019). Typically, the process involves meticulous adjustment of loss
and gain in order to achieve the desired outcome. In general, in a linear
system the presence of EPs is independent of the stability of the stationary
state that the system evolves to Khurgin (2020). However, in a nonlinear
system, more than one solution may be stable, which gives rise to the
phenomena of bistability and multistability Boyd (2020); Paraïso  _et al._
(2010); Gippius  _et al._ (2007); Cancellieri  _et al._ (2011). The existence
of nonlinear features may affect the non-Hermitian effects realized in linear
cases or give rise to entirely new phenomena Yu  _et al._ (2021); Wingenbach
_et al._ (2023); Ramezanpour and Bogdanov (2021); Xia  _et al._ (2021); Hassan
_et al._ (2015); Ramezani  _et al._ (2010); El-Ganainy  _et al._ (2018);
Wimmer  _et al._ (2015).

In order to examine the relationship between nonlinearity and non-Hermitian
physics, it is necessary to study systems that possess variable nonlinearity
and controllable gain and loss. Particularly suitable systems for this study
are those where matter couples with light, as they allow to take advantage of
the difference in physical properties of these components. For example, it was
demonstrated that exceptional points appear naturally in light-matter systems
of exciton-polaritons and subtreshold Fabry-Perot lasers Khurgin (2020); Hanai
_et al._ (2019). Moreover, it is possible to induce exceptional points by
manipulating spatial and spin degrees of freedom of exciton-polaritons in
various configurations Wingenbach  _et al._ (2023); Gao  _et al._ (2018a); Li
_et al._ (2022); Król  _et al._ (2022); Gao  _et al._ (2015); Su  _et al._
(2021); Song  _et al._ (2021); Hanai and Littlewood (2020); Rahmani  _et al._
(2023); Opala  _et al._ (2023); Gao  _et al._ (2018b); Liao  _et al._ (2021).
In the case of bosonic condensates of exciton-polaritons, it was predicted
that a dissipative first-order phase transition line exists in the phase
diagram Hanai  _et al._ (2019), similar to a critical point in a liquid-gas
phase transition. According to this study, this phase transition line exists
in the regime of strong light-matter coupling and has an endpoint which
corresponds to an exceptional point Hanai  _et al._ (2019).

In this letter, we investigate a non-Hermitian model describing interaction
between two oscillating modes. We use it to examine the significance of
nonlinearity in a non-Hermitian phase transition. This model can describe
light and matter modes in exciton-polariton condensation and lasing, as
investigated in Ref. Hanai  _et al._ (2019). We find that the model is
incomplete unless nonlinear saturation of gain is taken into account.
Importantly, saturation increases the complexity of the phase diagram and
leads to the appearance of bistability. It has also profound consequences on
the physics of the system. We find that while the first-order phase transition
line with an endpoint is present, the equivalence of the endpoint to an
exceptional point as found in Hanai  _et al._ (2019) is no longer valid in the
general case. The phase diagram of Ref. Hanai  _et al._ (2019) can be restored
in the limit of strong saturation. In contrast to the results of Ref. Hanai
_et al._ (2019), the transition between solutions can occur also in the weak
coupling regime. This suggests that the second threshold from polariton to
photon lasing, observed in experiments Hu  _et al._ (2021); Pieczarka  _et
al._ (2022); Tempel  _et al._ (2012), may be related to a dissipative phase
transition in the weak coupling regime. Moreover, we find a regime of
permanent Rabi-like oscillations between two stable solutions. This regime
corresponds to a line in the phase diagram that ends with an exceptional
point.

_Model and Analytical Solutions._ We consider a system of two coupled
oscillators described by a non-Hermitian Hamiltonian with gain and loss. The
imbalance between gain and loss in a linear system leads in general to
solutions exponentially growing or decaying in time. To obtain non-trivial
stationary solutions it is necessary to include nonlinearity. Here we adopt
cubic nonlinearity that appears naturally in symmetric systems with no
dependence on the complex phase. Such a model can be realized, among many
other physical systems, in the case of cavity photons coupled to excitons,
where the nonlinearity occurs only in the matter (exciton) component Kavokin
_et al._ (2016). The system is described by complex functions
$\psi_{C}=n_{C}e^{i\varphi_{C}}$ and $\psi_{X}=n_{X}e^{i\varphi_{X}}$,
corresponding to amplitudes of cavity photons and excitons, respectively. The
dynamics is governed by equations $i\hbar\partial\psi/\partial t=$
$i\hbar\partial_{t}|\Psi\rangle=H|\Psi\rangle$ with
$|\Psi\rangle=(\psi_{C},\psi_{X})^{\mathrm{T}}$, where non-Hermitian
Hamiltonian $H$ is given by Hanai  _et al._ (2019)

$\displaystyle H=\left(\begin{matrix}E_{C}-i\hbar\gamma_{C}&\hbar\Omega_{R}\\\
\hbar\Omega_{R}&E_{X}+g|\psi_{X}|^{2}+ip\end{matrix}\right)\,.$ (1)

Here $\hbar\Omega_{R}$ is the coupling strength, $\gamma_{C}$ is the decay
rate of the photon field, and $p$ represents the gain to the exciton field.
This gain can be realized in practice by nonresonant optical or electrical
pumping. We define the complex nonlinear coefficient as $g=g_{1}-ig_{2}$,
where $g_{1}$ is the strength of two body interactions (Kerr-like
nonlinearity) and $g_{2}|\psi_{X}|^{2}$ is the saturation term that allows to
avoid instability. Spectrum of Hamiltonian (1) can be found analytically

$\displaystyle E=$
$\displaystyle\frac{1}{2}\big{[}E_{c}+\mathcal{E}+i(\mathcal{P}-\hbar\gamma_{c})$
$\displaystyle\pm\sqrt{4\hbar^{2}\Omega_{R}^{2}+[\mathcal{E}-E_{c}+i(\mathcal{P}+\hbar\gamma_{c})]^{2}}\big{]}\,,$
(2)

where $\mathcal{P}=p-g_{2}(n_{X}^{\mathrm{SS}})^{2}$ and
$\mathcal{E}=E_{x}+g_{1}(n_{X}^{\mathrm{SS}})^{2}$. For convenience, we denote
the solution associated with plus (minus) by $U(L)$. The respective steady
state analytical solutions $|\Psi\rangle=|\Psi_{0}\rangle e^{-iEt}$ can be
found from the condition $\mathrm{Im}[E]=0$, that is, the imaginary part of
the eigenvalue of (1) must be zero. In Hanai  _et al._ (2019), it was argued
that one or two real energy solutions exist in certain regions in parameter
space. However, it can be seen from (Exceptional points and phase transitions
in non-Hermitian binary systems) that except from special values of
parameters, real energy solutions can exist only when saturation represented
by $g_{2}$ is taken into account. We will show below that accounting for the
nonlinear $g_{2}$ term does in fact lead to the appearance of up to three
real-energy solutions, each of them of the form (Exceptional points and phase
transitions in non-Hermitian binary systems).

The condition $\mathrm{Im}[E]=0$ allows one to find analytical expression for
$n_{X}^{\text{SS}}$

$\displaystyle(n_{X}^{\mathrm{SS}})^{2}=\frac{1}{g}\big{(}\mathrm{Re}[E]-E_{X}-iP-\frac{(\hbar\Omega_{R})^{2}}{\mathrm{Re}[E]-E_{C}+i\hbar\gamma_{C}}\big{)}.$
(3)

The resulting explicit formula for $n_{X}^{\text{SS}}$ is tedious, but for a
given $n_{X}^{\text{SS}}$, one can find closed forms of steady state
$n_{C}^{\text{SS}}$ and $\varphi_{CX}=\varphi_{C}-\varphi_{X}$

$\displaystyle n^{\mathrm{SS}}_{C}=$ $\displaystyle
n^{\mathrm{SS}}_{X}\sqrt{\frac{p}{\hbar\gamma_{C}}-\frac{\left(n_{X}^{\mathrm{SS}}\right)^{2}g_{2}}{\hbar\gamma_{C}}}\,,$
(4) $\displaystyle\varphi_{CX}^{SS}=$ $\displaystyle\arg\left(\frac{\delta-
g_{1}(n_{X}^{\mathrm{SS}})^{2}}{\hbar\Omega_{R}\left(n^{\mathrm{SS}}_{C}/n_{X}^{\mathrm{SS}}-n_{X}^{\mathrm{SS}}/n^{\mathrm{SS}}_{C}\right)}-i\frac{\gamma_{C}n^{\mathrm{SS}}_{C}}{\Omega_{R}n_{X}^{\mathrm{SS}}}\right)\,,$
(5)

where we introduced photon-exciton energy detuning $\delta=E_{C}-E_{X}$.

Figure 1:  (a-b): Phase diagrams of generic nonresonantly driven binary system
(1). In (a) the number of stationary states is marked with colors in the
function of normalized photon decay rate ($\gamma_{c}$) and pumping strength
($p$). In (b) only the lowest-energy stable state is shown. Here colors
indicate the real part of the energy of the corresponding solution. In (a) and
(b) the exceptional point (EP, marked with a red star) appears both as a point
on the phase boundary and as an endpoint of the R-line. The R-Line corresponds
to non-decaying Rabi-like oscillations between two modes. The endpoint of the
phase transition (ET) is marked with a black star. Cross-sections of constant
$\gamma_{c}$ with different numbers of thresholds (th) are marked with
horizontal lines. Panels (c-f) show bistability and phase transitions in more
detail. In (c) we show the case $\gamma_{C}=\Omega_{R}$, for which the energy
eigenvalues coalesce at the EP, which is also a turning point of a bistability
curve. Stable solutions are marked with S and black lines, while unstable
solutions are marked with US and orange lines. Panel (d) shows real part of
energies for different pumping and decay rates. The ET point corresponds to
the transition to bistability at $\gamma_{C}>\Omega_{R}$. This cross-section
is depicted in panel (e), while in panel (f) we show the case
$\gamma_{C}<\Omega_{R}$, where the unstable solution is split into two
branches, and the lowest-energy solution becomes unstable in certain regions
of pumping. Other parameters are $\delta=0.2~{}\hbar\Omega_{R}$,
$g_{1}=0.1~{}\hbar\Omega_{R}$, $E_{X}=0$, $E_{C}=0.2~{}\hbar\Omega$ and
$g_{2}=0.3~{}g_{1}$.

_Non-Hermitian Phase Transitions._ We use the analytical solutions from the
previous section to determine the phase diagram of the system, looking at it
from two perspectives. We analyze the steady state solutions and their
multiplicity, as in Fig. 1(a). On the other hand, we consider the lowest-
energy state among the dynamically stable ones and investigate its properties
and possible transitions, see Fig. 1(b). The latter approach is equivalent to
analyzing a system that is weakly coupled to an energy sink, which does not
perturb the spectrum, but picks the lowest-energy stable solution after a
sufficiently long evolution due to its energetic stability.

In the case when the conservative nonlinearity $g_{1}$ is stronger than the
dissipative nonlinearity $g_{2}$, representative phase diagrams are shown in
Fig. 1. We focus on the blue-detuned case ($\delta>0$), which is much richer
that the red-detuned case. In Fig. 1(a) the number of steady state solutions
is shown. Up to three non-zero solutions, corresponding to both upper and
lower branches of Eq. (Exceptional points and phase transitions in non-
Hermitian binary systems) can exist, which results from the nonlinearity of
the system. The region of zero solutions corresponds to the situation where
pumping cannot overcome losses and no lasing nor polariton condensation
occurs. For given $\Omega$ and $\gamma_{C}$, increasing pumping $p$ can lead
to one or several thresholds, as indicated with horizontal lines.

Special points in the phase diagram (marked by stars in Fig. 1) include the
exceptional point (EP) and the endpoint of the first-order phase transition
(ET). In contrast to Hanai  _et al._ (2019), we find that in general they do
not coincide. To determine the position of the EP, one can find the following
conditions for which the real and imaginary parts of eigenvalues are zero in
Eq. (Exceptional points and phase transitions in non-Hermitian binary systems)

$p^{\mathrm{EP}}=\hbar\Omega_{R}+\frac{g_{2}\delta}{g_{1}}\,,~{}~{}\gamma_{C}=\Omega_{R}\,.$
(6)

This can occur when $n_{X}^{\mathrm{SS}}=\delta/g_{1}$, that is, whenever the
system is blue-detuned ($\delta>0$). On the other hand, the ET point is
clearly visualised in the phase diagram that takes into account the energetic
instability in panel Fig. 1(b). The first-order phase transition line begins
at the ET point in the weak coupling regime ($\gamma_{C}>\Omega_{R}$) and
follows the arc represented by the ET-EP line towards the EP point. Below the
EP, the phase transition line follows into the strong coupling regime. We
conclude that, contrary to the results of Hanai  _et al._ (2019), the first-
order phase transition can occur also in the weak coupling regime. This can be
explained by a simple physical argument. Since the pumping influences the
effective photon-exciton detuning
$\tilde{\delta}=E_{C}-(E_{X}+g(n^{\mathrm{SS}}_{X})^{2})$, the increase of
pumping can change of the sign of $\tilde{\delta}$, leading to an abrupt
change of the lowest-energy state in the weak-coupling regime.

Figure 1(d) shows the dependence of the real part of the energy of solutions
shown in Figs. 1(a,b), in the vicinity of the ET-EP line. As can be seen, the
ET point is the point of the transition to bistability. On the other hand, the
EP point corresponds to a turning point in the bistability curve. The cross-
section including the EP point ($\gamma_{C}=\Omega$) is depicted in more
detail in Figure 1(c), which shows the occurrence of two stable branches from
the upper and lower branches of Eq. (Exceptional points and phase transitions
in non-Hermitian binary systems) and one unstable branch. At the EP, the
unstable upper branch coalesces with the lower stable branch, leading to the
first-order phase transition. The cross-section with the ET point
($\gamma_{C}>\Omega_{R}$) is shown in Fig. 1(e), where the bistability curve
closes, and the transition from the upper to lower branch becomes smooth. This
leads to the possibility to encircle the exceptional point as indicated with
arrows in Fig. 1(d).

Interestingly, additional features that have an influence on the physics of
the system can occur in the strong coupling case ($\gamma_{C}<\Omega_{R}$),
see Fig. 1(f). These include the disappearance of one of the solutions in a
certain parameter range and the dynamical instability of the lowest-energy
branch (marked with orange line). Consequently, the upper, higher-energy
solution may become the only viable solution despite the existence of lower-
energy solutions.

Figure 2: Phase diagrams in the case when dissipative nonlinearity $g_{2}$
dominates over the conservative nonlinearity $g_{1}$. In this case, endpoint
of the phase transition (ET) and exceptional point (EP) correspond to the same
point in parameter space, recovering the results of Hanai  _et al._ (2019).
Parameters are the same as in Fig. 1, except for $g_{2}=4.5~{}g_{1}$.

In the opposite case when the dissipative nonlinearity dominates over the
conservative one, we find that the phase diagram of energetically stable
solutions recovers the results of Hanai  _et al._ (2019), see Fig. 2. As the
dissipative nonlinearity is increased, the length of the ET-EP arc decreases,
and finally the two points coalesce. In this specific case, the exceptional
point is characterized by a jagged crest in the phase diagram, embodying a
third-order exceptional point (see supplementary materials). This phenomenon
arises from the coalescence of two stable solutions and a single unstable
solution.

_Permanent Rabi-Like Oscillations: R-Line._ Our analysis allows to predict
that a peculiar oscillating state may form, as indicated in Fig. 1(a) by _R-
Line_. In this case, long evolution leads to permanent oscillations,
resembling Rabi oscillations in a two-level system, instead of stationary
solutions. To explain this phenomenon, we examine imaginary and real parts of
eigenvalues given in Eq. (Exceptional points and phase transitions in non-
Hermitian binary systems). An example is shown in Figs. 3(a) and 3(b). In
general, two kinds of stationary solutions corresponding to
$\mathrm{Im}[E(n_{X})]=0$ may exist. As shown in Fig. 3(a), in this particular
case there are two solutions from the upper branch and one solution from the
lower branch (the black dashed vertical lines denote the emergent solutions).
Our interest is in solutions from upper and lower branches that occur at the
same $n_{X}$, while there is a gap in respective real parts, see Fig. 3(b).
Such solutions occur when $p=(g_{2}/g_{1})\delta+\hbar\gamma_{C}$, which
corresponds to a straight line (marked by R-line) in the phase diagram of Fig.
1(c).

An example of such permanent oscillations is shown in Fig. 3(c). After initial
transient time, the oscillations stabilize at a cetain amplitude. When
different initial conditions are used, the system may end up in one of the
steady state solutions, as shown in Fig. 3(d). The frequency of oscillations
is given by the gap, $\Omega=2\sqrt{\Omega_{R}^{2}-\gamma_{C}^{2}}$. When the
parameters of the system approach the exceptional point along the R-line, the
gap decreases and the period of oscillations increases. At the exceptional
point ($\Omega_{R}=\gamma_{C}$), the solutions coalesce and the period becomes
infinite. Therefore, the exceptional point is the endpoint of the R-line.

Figure 3: Permanent oscillations. (a) Imaginary part of energy eigenvalues
versus exciton density $n_{X}$. The condition $\mathrm{Im}[E(n_{X})]=0$
provides the possible stationary solutions. Here the solution from the lower
energy branch (denoted by L) coincides with one of the solutions from the
upper branch (denoted by U). (b) Real part of eigenvalues versus $n_{X}$. The
nonzero energy gap between the solutions results in Rabi oscillations with
frequency $\Omega_{R}$. Examples of simulations are shown in (c) and (d). In
(c) oscillations stabilize after initial transient. In (d), different initial
condition for the same parameters leads to a steady state solution from the
upper branch. Using parameters: $p=0.82~{}\hbar\Omega_{R}$,
$\gamma_{C}=0.75~{}\Omega_{R}$, $g_{1}=0.1~{}\hbar\Omega_{R}$,
$g_{2}=0.3~{}g_{1}$ and $\delta=0.2~{}\hbar\Omega_{R}$.

Discussion. We showed that, contrary to previous understanding, non-Hermitian
polariton systems exhibit first-order phase transition with an endpoint that
in general does not coincide with the exceptional point. Explanation of this
phenomenon requires taking into account the nonlinear gain saturation and the
consideration of the bistability curve. While the endpoint of the phase
transition is where the bistability appears, the exceptional point is where
the stable and unstable solutions coalesce. In addition, we demonstrated that
first-order phase transition may occur in the weak coupling regime, and that
for certain values of parameters one can predict permanent oscillations, whose
frequency vanishes at the exceptional point.

The predicted results contribute to the ongoing debate surrounding
polariton/photon lasing. The presence of an exceptional point has been
identified as the possible underlying factor for the observed second threshold
Hanai  _et al._ (2019). Here, we provide further insights by identifying
several other thresholds in phase diagrams and pointing out that multiplicity
and stability of solutions are also crucial factors, so far overlooked.

The presented results may be applied to much broader class of systems. The
non-Hermitian Hamiltonian represented by the $2\times 2$ matrix in Eq. (1)
describes in general an arbitrary two-mode oscillatory system with gain and
loss in the two modes, and the cubic nonlinearity in one of them. This term
appears naturally in any oscillatory system in the first order as long as the
nonlinearity respects the global $U(1)$ symmetry of the oscillations. Examples
include not only all quantum mechanical systems such as Bose-Einstein
condensates, but also high-frequency coupled classical oscillators, where
phase of oscillations is irrelevant on the time scale of a slowly varying
envelope. The results presented here should be applicable to any such system
that exhibits exceptional points and nonlinearity.

A.R. and M.M. acknowledge support from National Science Center, Poland (PL),
Grant No. 2016/22/E/ST3/00045. A.O. acknowledges support from Grant No.
2019/35/N/ST3/01379.

## References

  * Jiles (2015) D. Jiles,  _Introduction to Magnetism and Magnetic Materials_ (CRC Press, 2015).
  * Stishov (2018) S. M. Stishov,  _Phase Transitions for Beginners_ (WORLD SCIENTIFIC, 2018) https://www.worldscientific.com/doi/pdf/10.1142/11096 .
  * Proukakis  _et al._ (2017) N. Proukakis, D. Snoke, and P. Littlewood,  _Universal Themes of Bose-Einstein Condensation_ (Cambridge University Press, 2017).
  * Imada  _et al._ (1998) M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * Brink  _et al._ (2018) L. Brink, M. Gunn, J. Jose, J. Kosterlitz, and K. Phua,  _Topological Phase Transitions And New Developments_ (World Scientific Publishing Company, 2018).
  * Ashida  _et al._ (2020) Y. Ashida, Z. Gong, and M. Ueda, Advances in Physics 69, 249 (2020), https://doi.org/10.1080/00018732.2021.1876991 .
  * Miri and Alù (2019) M.-A. Miri and A. Alù, Science 363, eaar7709 (2019), https://www.science.org/doi/pdf/10.1126/science.aar7709 .
  * Gong  _et al._ (2018) Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Phys. Rev. X 8, 031079 (2018).
  * Peng  _et al._ (2022) B. Peng, S. K. Ozdemir, F. Lei, F. Monifi, M. Gianfreda, G. L. Long, S. Fan, F. Nori, C. M. Bender, and L. Yang, Nature Physics 30, 394 (2022).
  * Weidemann  _et al._ (2022) S. Weidemann, M. Kremer, S. Longhi, and A. Szameit, Nature 601, 354 (2022).
  * Okuma and Sato (2023) N. Okuma and M. Sato, Annual Review of Condensed Matter Physics 14, 83 (2023), https://doi.org/10.1146/annurev-conmatphys-040521-033133 .
  * Öztürk  _et al._ (2021) F. E. Öztürk, T. Lappe, G. Hellmann, J. Schmitt, J. Klaers, F. Vewinger, J. Kroha, and M. Weitz, Science 372, 88 (2021), https://www.science.org/doi/pdf/10.1126/science.abe9869 .
  * Khurgin (2020) J. B. Khurgin, Optica 7, 1015 (2020).
  * Boyd (2020) R. W. Boyd,  _Nonlinear optics_ (Academic press, 2020).
  * Paraïso  _et al._ (2010) T. Paraïso, M. Wouters, Y. Léger, F. Morier-Genoud, and B. Deveaud-Plédran, Nature materials 9, 655 (2010).
  * Gippius  _et al._ (2007) N. Gippius, I. Shelykh, D. Solnyshkov, S. Gavrilov, Y. G. Rubo, A. Kavokin, S. Tikhodeev, and G. Malpuech, Physical review letters 98, 236401 (2007).
  * Cancellieri  _et al._ (2011) E. Cancellieri, F. Marchetti, M. Szymańska, and C. Tejedor, Physical Review B 83, 214507 (2011).
  * Yu  _et al._ (2021) Z.-F. Yu, J.-K. Xue, L. Zhuang, J. Zhao, and W.-M. Liu, Phys. Rev. B 104, 235408 (2021).
  * Wingenbach  _et al._ (2023) J. Wingenbach, S. Schumacher, and X. Ma, preprint at https://arxiv.org/abs/2305.04855 (2023), https://doi.org/10.48550/arXiv.2305.04855.
  * Ramezanpour and Bogdanov (2021) S. Ramezanpour and A. Bogdanov, Phys. Rev. A 103, 043510 (2021).
  * Xia  _et al._ (2021) S. Xia, D. Kaltsas, D. Song, I. Komis, J. Xu, A. Szameit, H. Buljan, K. G. Makris, and Z. Chen, Science 372, 72 (2021).
  * Hassan  _et al._ (2015) A. U. Hassan, H. Hodaei, M.-A. Miri, M. Khajavikhan, and D. N. Christodoulides, Phys. Rev. A 92, 063807 (2015).
  * Ramezani  _et al._ (2010) H. Ramezani, T. Kottos, R. El-Ganainy, and D. N. Christodoulides, Phys. Rev. A 82, 043803 (2010).
  * El-Ganainy  _et al._ (2018) R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nature Physics 14, 11 (2018).
  * Wimmer  _et al._ (2015) M. Wimmer, A. Regensburger, M.-A. Miri, C. Bersch, D. N. Christodoulides, and U. Peschel, Nature communications 6, 7782 (2015).
  * Hanai  _et al._ (2019) R. Hanai, A. Edelman, Y. Ohashi, and P. B. Littlewood, Phys. Rev. Lett. 122, 185301 (2019).
  * Gao  _et al._ (2018a) T. Gao, G. Li, E. Estrecho, T. C. H. Liew, D. Comber-Todd, A. Nalitov, M. Steger, K. West, L. Pfeiffer, D. Snoke, _et al._ , Physical review letters 120, 065301 (2018a).
  * Li  _et al._ (2022) Y. Li, X. Ma, Z. Hatzopoulos, P. G. Savvidis, S. Schumacher, and T. Gao, ACS Photonics 9, 2079 (2022), https://doi.org/10.1021/acsphotonics.2c00288 .
  * Król  _et al._ (2022) M. Król, I. Septembre, P. Oliwa, M. Kedziora, K. Łempicka-Mirek, M. Muszyński, R. Mazur, P. Morawiak, W. Piecek, P. Kula, _et al._ , Nature Communications 13, 5340 (2022).
  * Gao  _et al._ (2015) T. Gao, E. Estrecho, K. Bliokh, T. Liew, M. Fraser, S. Brodbeck, M. Kamp, C. Schneider, S. Höfling, Y. Yamamoto, _et al._ , Nature 526, 554 (2015).
  * Su  _et al._ (2021) R. Su, E. Estrecho, D. Biegańska, Y. Huang, M. Wurdack, M. Pieczarka, A. G. Truscott, T. C. Liew, E. A. Ostrovskaya, and Q. Xiong, Science Advances 7, eabj8905 (2021).
  * Song  _et al._ (2021) H. G. Song, M. Choi, K. Y. Woo, C. H. Park, and Y.-H. Cho, Nature Photonics 15, 582 (2021).
  * Hanai and Littlewood (2020) R. Hanai and P. B. Littlewood, Physical Review Research 2, 033018 (2020).
  * Rahmani  _et al._ (2023) A. Rahmani, M. Kedziora, A. Opala, and M. Matuszewski, Physical Review B 107, 165309 (2023).
  * Opala  _et al._ (2023) A. Opala, M. Furman, M. Król, R. Mirek, K. Tyszka, B. Seredyński, W. Pacuski, J. Szczytko, M. Matuszewski, and B. Pietka, arXiv preprint arXiv:2306.01366 (2023).
  * Gao  _et al._ (2018b) W. Gao, X. Li, M. Bamba, and J. Kono, Nature Photonics 12, 362 (2018b).
  * Liao  _et al._ (2021) Q. Liao, C. Leblanc, J. Ren, F. Li, Y. Li, D. Solnyshkov, G. Malpuech, J. Yao, and H. Fu, Physical Review Letters 127, 107402 (2021).
  * Hu  _et al._ (2021) J. Hu, Z. Wang, S. Kim, H. Deng, S. Brodbeck, C. Schneider, S. Höfling, N. H. Kwong, and R. Binder, Phys. Rev. X 11, 011018 (2021).
  * Pieczarka  _et al._ (2022) M. Pieczarka, D. Biegańska, C. Schneider, S. Höfling, S. Klembt, G. Sek, and M. Syperek, Opt. Express 30, 17070 (2022).
  * Tempel  _et al._ (2012) J.-S. Tempel, F. Veit, M. Aßmann, L. E. Kreilkamp, A. Rahimi-Iman, A. Löffler, S. Höfling, S. Reitzenstein, L. Worschech, A. Forchel, and M. Bayer, Phys. Rev. B 85, 075318 (2012).
  * Kavokin  _et al._ (2016) A. V. Kavokin, J. J. Baumberg, G. Malpuech, and F. P. Laussy,  _Microcavities_ (Oxford University Press, Oxford, 2016).


Floquet non-equilibrium Green's function and Floquet quantum master equation for electronic transport: The role of electron-electron interactions

Vahid Mosallanejad
Department of Chemistry, Westlake University, Hangzhou, Zhejiang 310024, China $\&$ Institute of Natural Sciences, Westlake Institute for Advanced Study, Hangzhou, Zhejiang 310024, China

Yu Wang
Department of Chemistry, Westlake University, Hangzhou, Zhejiang 310024, China $\&$ Institute of Natural Sciences, Westlake Institute for Advanced Study, Hangzhou, Zhejiang 310024, China

Wenjie Dou
Department of Chemistry, Westlake University, Hangzhou, Zhejiang 310024, China $\&$ Institute of Natural Sciences, Westlake Institute for Advanced Study, Hangzhou, Zhejiang 310024, China

Non-equilibrium Green's function (NEGF) and quantum master equation (QME) are two main classes of approaches for electronic transport. We discuss various Floquet variances of these formalisms for transport properties of a quantum dot driven via interaction with an external periodic field. We first derived two versions of the Floquet NEGF. We also explore an ansatz of the Floquet NEGF formalism for the interacting systems. In addition, we derived two versions of Floquet QME in the weak interaction regime. With each method, we elaborate on the evaluation of the expectation values of the number and current operators. We examined these methods for transport through a two-level system that is subject to periodic driving. The results of all four methods show great consistency. We expect these Floquet quantum transport methods to be useful in studying molecular junctions exposed to light.

§ INTRODUCTION

Light-matter interactions on quantum scales have long been an attractive research topic<cit.>. In particular, strong electron-photon interaction opens a new avenue toward engineering material properties <cit.>. In the interpretation of phenomena related to light-matter interaction, the strength of light plays an essential role such that a perturbative treatment of light-matter interaction becomes insufficient <cit.>.

Technological advancements, for example, in time-resolved photo spectroscopy <cit.>, have allowed us to experimentally demonstrate exotic states such as Floquet-Bloch states on topological insulators <cit.>. Electronic states in such a system can be understood via Floquet theory as if a classical monochromatic light is coupled to a closed fermionic system <cit.>. Floquet theory is a non-perturbative mathematical technique that essentially turns the time-dependent schrödinger equation with a time-periodic Hamiltonian to a time-independent equation in the price of increasing the dimensionality of the basis set <cit.>.

Floquet theory has been used in the interpretation of Kondo effects, and noise effects <cit.>, spin currents in DNA <cit.>, electron transfer in donor-bridge-acceptor systems <cit.>, driven quantum dots <cit.>, laser absorption properties<cit.>, nonadiabatic dynamics <cit.> and many other applications.

To further enhance our understanding of light-matter interactions, it seems essential to develop approaches that describe how light interacts with open quantum systems. This requirement can be met by properly combining the Floquet theory with quantum transport theories.

In other words, two classes of quantum transport approaches, namely the non-equilibrium Green's function (NEGF) and the quantum master equation (QME), can be complemented with the Floquet theory.

There have been several works on extending the Floquet theory to open quantum systems. Floquet Green’s function is defined in a number of different ways <cit.>. While the Floquet QME was less well-established in the past <cit.>, recent demands for dynamically controlling the material properties by the Floquet engineering stimulated the development of new methods <cit.>.

NEGF approaches are often employed when a one-body Hamiltonian (e.g., the tight-binding Hamiltonian) can sufficiently describe a system of interest (the dot) whereas QME is suitable when many-body interactions (such as electron-electron Coulomb repulsion) are important <cit.>. In addition, QME is known to work well when the coupling between the System and environment is weak (weak-coupling limit) whereas NEGF is insensitive to the System-environment coupling strength.

Historically, two main flavors for Floquet theory exist in the context of quantum mechanics/solid state physics. The first one that focused on the evolution operator was developed by Shirley in the mid-sixties <cit.> while the second one that focused on the discreet expansion of wavefunction was developed by Sambe <cit.>.

In this work, we intend to discuss that there are multiple ways in which both quantum transport approaches can benefit from the Floquet theory.

We will offer compact derivations for Floquet non-equilibrium Green's functions (which rely on the expansion of Green's function pioneered by Sambe) and we will also microscopically derive two versions of Floquet quantum master equations (one based on the Floquet-based evolution operator pioneered by Shirley and the other relies on the concept of Floquet density matrix). We will also examine presented methods for a model system.

To the best of our knowledge, there is no work dedicated to the consistency check between different Floquet quantum transport theories.
Here, we sketch general Floquet theories which are neither based on an average Hamiltonian approach (Floquet-Magnus) nor based on stroboscopic and micromotion evolution (Van Vleck approximation) <cit.>. Hence, the approaches presented here are not sensitive to the value of driving frequency, in principle.

§ THEORY

§.§ Model Hamiltonian

With this work, we are concerned with the electronic transport through a quantum dot (the System) that is weakly connected to the thermal fermionic bath (lead) while it strongly interacts with an external monochromatic light with the frequency $\omega$. In such scenarios, often a time-periodic off-diagonal term (that represents the dipole approximation) should be added to the unperturbed System Hamiltonian ($\hat{H}_S$) such that the total System Hamiltonian will become time-periodic $\hat{H}_S(t)=\hat{H}_S(t+T)$ [$T=2\pi/\omega$]. The bath Hamiltonian, $\hat{H}_B$, and the system-bath interaction, $\hat{H}_{SB}$, will be remained time-independent.

The spinless model Hamiltonian can be given by

\begin{eqnarray}
\label{eq:1}
\hat{H}_{tot}(t) &=& \hat{H}_S(t)  + \hat{H}_B + \hat{H}_{SB}
\\
\label{eq:2}
\hat{H}_S(t) &=& \sum_{ij} h_{ij} (t) \hat{d}_i^\dagger \hat{d}_j \\
\label{eq:3}
\hat{H}_B &=& \sum_{l k} \epsilon_{l k} \hat{c}_{l k}^\dagger \hat{c}_{l k}
\\
\label{eq:4}
\hat{H}_{SB}  &=&  \sum_{l k,i}
V_{l k, i}
\hat{c}_{l k}^\dagger \hat{d}_i +\mathrm{H.c.}
\end{eqnarray}

Here, $\hat{d}_i$ ($\hat{d}_i^{\,\dagger}$) is the System's electronic annihilation (creation) operator in the many-body space and $h_{ij}(t)=h_{ij}(t+T)$ represents a periodic one-body Hamiltonian. Likewise, $\hat{c}_{l k}$ ($\hat{c}_{l k}^\dagger$) is the annihilation (creation) operator for the $k$th electronic orbital in the bath $l$.

The quantity $V_{l k, i}$ is the coupling strength between the $k$th orbital of the bath (lead) $l$, and the System's orbital $\hat{d}_i$.
In addition the number operator is $\hat{n}=\sum_{i} \hat{d}_i^\dagger \hat{d}_i$.
We can further simplify the notation of system-bath interaction Hamiltonian as
$\hat{H}_{SB} = \sum_{i}
\hat{C}_{i}^\dagger \hat{d}_i +\mathrm{H.c}$.

We shall also associate the bath $l$ with the electrochemical potential $\mu_{l}$. Note that the System Hamiltonian is quadratic whereas the bath Hamiltonian is non-interacting.

§.§ Floquet Green's function

In the following we explore two types of Floquet Green's function whereby we named 1) Vector-like NEGF, and 2) Matrix-like NEGF. Our starting point for deriving these two flavors of Floquet Green's function is the two-time Kadanoff-Baym (KB) equation <cit.> for the retarded Green's function, $G^r(t, t^{\prime})$, as

\begin{eqnarray}
\label{eq:5}
\begin{aligned}
&\Big(i \frac{d}{d t}-h(t)\Big)G^r(t, t^{\prime})-\int_{-\infty}^{+\infty} d t_1 \Sigma^r(t, t_1) \times
\\&
G^r(t_1, t^{\prime})=I \delta(t-t^{\prime}).
\end{aligned}
\end{eqnarray}

Hereafter, we will set $\hbar=1$.

$\Sigma^r=\sum_{l} \Sigma^r_{l}$ refers to the total retarded self-energy. $\Sigma^r_l$ refers to the retarded self-energy of the non-interacting bath $l$, which only depends on the time difference, such that $\Sigma^r(t, t^{\prime})=\Sigma^r(t-t^{\prime})$.

§.§.§ Vector-like (V-like) Floquet NEGF

One way to derive the EOM within the V-like Floquet NEGF is to split the derivation into two parts. In the first part, we define following time-energy retarded and advanced Green's functions

\begin{eqnarray}
\label{eq:6}
\begin{aligned}
&\int_{-\infty}^{\infty} d t^{\prime} G^{r/a}(t, t^{\prime}) e^{i \mathcal{E}(t-t^{\prime})} = G^{r/a}(t, \mathcal{E}),
\end{aligned}
\end{eqnarray}

where $\mathcal{E}$ is the quasi-energy variable. Then, after performing the continuous Fourier transformation, the KB equation for retarded Green's function simplifies as

\begin{eqnarray}
\label{eq:7}
\begin{aligned}
&\Big(\mathcal{E}\,I+i \frac{d}{d t} -h(t)\Big) G^r(t, \mathcal{E})-\int_0^{+\infty} d \tau \Sigma^r(\tau) \times
\\&
e^{\tfrac{i \mathcal{E}}{\hbar} \tau} G^r(t-\tau, \mathcal{E})= I
\end{aligned}
\end{eqnarray}

Note that $\tau=t-t_1$, and we have employed convolution identity in the self-energy term. In the second part, we expand the $G^r(t, \mathcal{E})$ with the following complex Fourier series

\begin{eqnarray}
\label{eq:8}
\begin{aligned}
G^r(t, \mathcal{E})=\sum_{n=-N}^N e^{-i n \omega t} g^r_n(\mathcal{E}).
\end{aligned}
\end{eqnarray}

$N$ is an integer that determines the truncation of the Fourier space basis set. After substitution Eq. (<ref>) into Eq. (<ref>), we arrive at

\begin{eqnarray}
\label{eq:9}
\begin{aligned}
\left(e^{-i n \omega t} \mathcal{E}~I+e^{-i n \omega t} n
\omega~I-h(t) e^{-i n \omega t}
\right) g^r_n(\mathcal{E})-
\\&
\sum_{n} e^{-i n \omega(t)} \int_0^{+\infty} d \tau \Sigma^r(\tau) e^{i(\mathcal{E}+n \omega) \tau} g^r_n(\mathcal{E})= I.
\end{aligned}
\end{eqnarray}

We then desire to find an EOM for $g^r_n(\mathcal{E})$ which requires multiplying both sides to $e^{i m \omega t}$ and taking the time average over one period. EOM for $g^r_n(\mathcal{E})$ reads as

\begin{eqnarray}
\label{eq:10}
\begin{aligned}
\big(
(\mathcal{E}+m \omega) I
-\bar{\Sigma}^{r}(\mathcal{E}+m \omega)
\big)\delta_{mn}-
\sum_{n} h_{m n}\Big)\times
\\&
g^r_n(\mathcal{E})=\delta_{m 0} I,
\end{aligned}
\end{eqnarray}

where $\bar{\Sigma}^r$ and $[h]_{m n}$ are given by

\begin{eqnarray}
\label{eq:11}
\bar{\Sigma}^{r}(\mathcal{E}+m \hbar \omega) &=&\int_0^{+\infty} d \tau \Sigma^r(\tau) e^{i(\mathcal{E}+m \omega) \tau}, \\
\label{eq:12}
h_{m n} &=& \frac{1}{T} \int_0^T d t e^{i m \omega t}h(t) e^{-i n \omega t}.
\end{eqnarray}

By running Eq. (<ref>) over indices $n$ and $m$, the compact form of V-like Floquet NEGF appears as

\begin{eqnarray}
\label{eq:13}
\begin{aligned}
\mathcal{E}[I^F]
-[h^F]-[\bar{\Sigma}^{r\,F}] \Big)
\mathbf{g}^{r\,F}(\mathcal{E})=\mathbf{I}_o^{F},
\end{aligned}
\end{eqnarray}

where the $[h^F]=\sum_{n,m}h_{m n}-I (m \omega) \delta_{mn}$ is a standard definition for the Floquet Hamiltonian. $[I^F]$ and $[\bar{\Sigma}^{r/a\,F}]$ are block diagonal matrices made of $I$ and $\bar{\Sigma}^{r/a}$.

In the above block matrix EOM, the $\mathbf{g}^{r\,F}$ refers to a stack of $g^r_n$ placed over each other (V-like), and $\mathbf{I}_o^{F}$ refers to a stacking zero matrices symmetrically placed around a central identity matrix. Note that with the energy-independent wide-band approximation, $\bar{\Sigma}^{r}(\mathcal{E}+m \omega)=\bar{\Sigma}^{r}$, there would be no implementation complexity on the Floquet self-energy term. In addition, $\mathcal{E}$ is unbounded. The most important point about Eq. (<ref>) is that it is a time-independent EOM.

Such a two step procedure was first employed by Sambe to directly solve Schrödinger equation with the time-periodic Hamiltonian.

Also, one can combine the inverse of Eq. (<ref>) with Eq. (<ref>) to expand the two-time quantity $Q(t, t^{\prime})$ in the V-like form as

\begin{eqnarray}
\label{eq:14}
\begin{aligned}
&Q(t, t^{\prime})=\sum_{n} \int_{-\infty}^{\infty}
\frac{d \mathcal{E}}{2\pi}
e^{-i(\mathcal{E}+n\omega) t} e^{i \mathcal{E} t^{\prime}} q_n(\mathcal{E}).
\end{aligned}
\end{eqnarray}

From the definition given in Eq. (<ref>), the requirement $G^{a}( t^{\prime},t)=G^r(t, t^{\prime})^{\, \dagger}$ indicates $g^a_{n}(\mathcal{E})=g^r_{-n}(\mathcal{E}+n \omega) ^{\,\dagger}$.

§.§.§ Matrix-like (M-like) Floquet NEGF

We first define the following one-step transformation for a two-time function

\begin{eqnarray}
\label{eq:15}
\begin{aligned}
&Q(t, t^{\prime})=\sum_{n,m} \int_0^{ \omega}
\frac{d \mathcal{E}}{2\pi}
e^{-i(\mathcal{E}+n \omega) t} e^{i(\mathcal{E}+m \omega) t^{\prime}} q_{n m}(\mathcal{E}).~~
\end{aligned}
\end{eqnarray}

The inverse transformation, for the matrix coefficients $q_{n m}(\mathcal{E})$, is given by

\begin{eqnarray}
\label{eq:16}
\begin{aligned}
q_{n m}(\mathcal{E})=
%d t^{\prime}
\frac{d t^{\prime}}{T} \int_0^T d t e^{i(\mathcal{E}+n \omega) t}
e^{-i(\mathcal{E}+m \omega) t^{\prime}}
Q(t, t^{\prime}).~~
\end{aligned}
\end{eqnarray}

Contrary to V-like Floquet NEGF, $\mathcal{E}$ is a bounded (independent) variable on the range [$0$ $\omega$]. We proceed by expanding $G^{r}$ using Eq. (<ref>) and substituting the result into the KB equation, and solving for the matrix coefficient $g^r_{nm}$. In this process, one can first notice that RHS can be reduced to $I \delta_{n m}$. After some lengthy algebra, we arrive at the following EOM for the retarded Green’s

\begin{eqnarray}
\label{eq:17}
\begin{aligned}
(\mathcal{E}+n \omega) -
\bar{\Sigma}^{r}
\mathcal{E}+n\omega
\Big)
g^r_{n m}(\mathcal{E})-
\\&
\sum_{k}h_{n k}~g^r_{k m}(\mathcal{E})=I \delta_{n m}.
\end{aligned}
\end{eqnarray}

The definition of $\bar{\Sigma}^{r}$ and $h_{n k}$ is identical to what is given in Eqs. (<ref>) and (<ref>). Comparing Eq. (<ref>) with Eq. (<ref>), it is evident why this form of EOM can be referred to as the Matrix-like Floquet NEGF. The compact form of M-like Floquet Green's function reads as

\begin{eqnarray}
\label{eq:18}
\begin{aligned}
\mathcal{E}[I^F]
-[h^F]-[\bar{\Sigma}^{r\,F}] \Big)
\end{aligned}
\end{eqnarray}

In the M-like framework, the condition $G^{a}( t^{\prime},t)=G^r(t, t^{\prime})^\dagger$ imposes $g^a_{n m}=(g^r_{m n})^\dagger$ which in turn confirms that $[G^{a\,F}]=[G^{r\,F}]^{\,\dagger}$. This property makes the M-like framework more convenient compared to the V-like.

§.§ Observable in Floquet NEGF

To explore the transport properties, the expectation values of two (operators) observable are essential. 1) the number operator $\langle \hat{n} \rangle (t)$, and 2) the terminal current $\langle \hat{J} \rangle={J}$. For $\langle \hat{n} \rangle (t)$, we shall trace the two-time lesser Green's function as $\langle \hat{n} \rangle(t)=-i \operatorname{Tr}(G^{<}(t,t))$ <cit.>. The $G^{<}(t,t^{\prime})$ given by

\begin{eqnarray}
\label{eq:19}
\begin{aligned}
\int_{-\infty}^{+\infty} d t_1 \int_{-\infty}^{+\infty} d t_2
&G^r(t, t_1)
\times\\
\Sigma^{<}(t_1-t_2) G^a(t_2, t^{\prime}).
\end{aligned}
\end{eqnarray}

The total lesser self-energy is denoted by: $\Sigma^{<} = \sum_{l} \Sigma^{<}_{l}$. The two-time particle current per spin passing through the contact $l$ is given by

\begin{eqnarray}
\label{eq:20}
\begin{aligned}
{J}_l(t, t^{\prime})=
\int_{-\infty}^{+\infty} d t_1
\Big(
\Sigma_{l}^{<} (t_1,t^{\prime})+
\\&
G^{<}(t, t_1)
\Sigma_{l}^a (t_1,t^{\prime})
\Big).
\end{aligned}
\end{eqnarray}

In the following (without details of derivation), we summarize general expressions for Floquet components of these two observable within the V-like and M-like Floquet NEGF frameworks.

Observable in the V-like Floquet NEGF.–

The expansion given in Eq. (<ref>) will be used to expand $G^{</r/a}$ (in terms of $g_n^{</r/a}$) and resulting expressions will be substituted in Eq. (<ref>). Then, we proceed to derive the matrix coefficients $g_n^{<}(\mathcal{E})$. After doing some algebra, it reads as

\begin{eqnarray}
\label{eq:21}
\begin{aligned}
\Sigma^{<}(\mathcal{E}+m\omega)
\\
\sum_{m}
\Sigma^{<}(\mathcal{E}+m\omega)
\end{aligned}
\end{eqnarray}

In the energy domain, the lesser self-energy is: $\bar{\Sigma}^{<}_{l} (\mathcal{E})=i
\Gamma^r_{l}(\mathcal{E}) f_l(\mathcal{E})$, $\Gamma^r_{l}=2 \operatorname{Im}(\bar{\Sigma}^r_{l})$. Here, $ f_l(\mathcal{E})\equiv f(\mathcal{E},\mu_{l})$ and $\mu_l$ are the Fermi function and the electrochemical potential associated with the terminal $l$.

For the current, we will follow the same procedure with expanding the ${J}_l(t, t^{\prime})$ in terms of the matrix coefficient $J_{l,n}(\mathcal{E})$. Solving for $J_{l,n}(\mathcal{E})$, we arrive at

\begin{eqnarray}
\label{eq:22}
\begin{aligned}
%\frac{2 e}{\hbar}
\Big(
\bar{\Sigma}_{l}^{<} (\mathcal{E})+
g_n^{<}(\mathcal{E})\bar{\Sigma}_{l}^a (\mathcal{E})
\Big).
\end{aligned}
\end{eqnarray}

Observable in the M-like Floquet NEGF.– Similar to above, one can employ Eq. (<ref>) to expand $G^{</a/r}$ and substitute the resulting expressions in Eq. (<ref>). After doing some algebra, we arrive at the following relation for $g^<_{nm}(\mathcal{E})$

\begin{eqnarray}
\label{eq:23}
\begin{aligned}
\sum_{k}
\bar{\Sigma}^{<}(\mathcal{E}+m\omega)
\end{aligned}
\end{eqnarray}

For the current, ${J}_{l}(t,t^{\prime})$ has to be expanded in terms of the coefficient matrix $J_{l,nm}(\mathcal{E})$ using Eq. (<ref>). Substituting expanded forms of all components into Eq. (<ref>) and solving for $J_{l,nm}(\mathcal{E})$ results in

\begin{eqnarray}
\label{eq:24}
\begin{aligned}
\Big(
\bar{\Sigma}_{l}^{<} (\mathcal{E}+m\omega)+
\bar{\Sigma}_{l}^a (\mathcal{E}+m\omega)
\Big).
\end{aligned}
\end{eqnarray}

Eqs. (<ref>-<ref>) may look complex at first glance, however, these equations can be further simplified when we evaluate the time-averaged of an observable.

§.§ Time-average of an observable

Within the V-like Floquet NEGF framework, the time average of an observable over one period can be obtained by changing the $t^{\prime}$ to $t$ in Eq. (<ref>) and performing the time average over one period. Using Eq. (<ref>), the time average of the number operator over one period, $\overline{ n }$, in the V-like framework read as

\begin{eqnarray}
\label{eq:25}
\begin{aligned}
\overline{ n }
\frac{-i}{T} \int_0^T d t~ \operatorname{Tr}(G^{<}(t, t))=
\int_{-\infty}^{\infty}
\frac{d \mathcal{E}}{2\pi}
~ \operatorname{Tr}(g^<_0(\mathcal{E}))
\\&
\int_{-\infty}^{\infty}
\frac{d \mathcal{E}}{2\pi}
\sum_{m}
\operatorname{Tr}
\Big(
\bar{\Sigma}^{<} (\mathcal{E}+m\omega)
\Big).
\end{aligned}
\end{eqnarray}

In the above expression, we have used the property
=g_m^r(\mathcal{E})$, and we have also employed the wide-band approximation as well.

We stress that only the central component ($n=0$) of $\mathbf{g}^{<}$ is required to evaluate the time average of $\langle \hat{n} \rangle(t)$ over one period.

Within the M-like Floquet NEGF framework, the time average of the number operator  over one period reads as

\begin{eqnarray}
\label{eq:26}
\begin{aligned}
\overline{ n}
\sum_{n} \int_0^{\omega}
\frac{d \mathcal{E}}{2\pi}
\operatorname{Tr}(g^<_{n n}(\mathcal{E}))
\\&
\int_0^{ \omega}
\frac{d \mathcal{E}}{2\pi}
\operatorname{Tr}
\Big(
\Big).
\end{aligned}
\end{eqnarray}

Here, $[\bar{\Sigma}^{<F}](\mathcal{E})$ denotes a block diagonal matrix made of $\bar{\Sigma}^{<}(\mathcal{E}+m\omega)$.

Next, we simplify current expressions. In the V-like framework, the time average of the current over a period, $\overline{{J}}$, simplifies as

\begin{eqnarray}
\label{eq:27}
\begin{aligned}
\overline{{J}}_{l}=
\int_{-\infty}^{\infty}
\frac{d \mathcal{E}}{2\pi}
\Big(
\Sigma_{l}^{<} (\mathcal{E})+
\Sigma_{l}^a (\mathcal{E})
\Big).
\end{aligned}
\end{eqnarray}

In the M-like framework, based on Eqs. (<ref>) and (<ref>), the time average of the current over one period read as

\begin{eqnarray}
\label{eq:28}
\begin{aligned}
\overline{{J}}_{l}
\sum_{n} \int_0^{\omega}
\frac{d \mathcal{E}}{2\pi}
\operatorname{Tr}
\\&
\int_0^{ \omega}
\frac{d \mathcal{E}}{2\pi}
\operatorname{Tr}
\Big(
\Big).
\end{aligned}
\end{eqnarray}

In the absence of electron-electron interaction, same as in none-Floquet scenarios, the above relation can be further simplified to a Landauer-like expression for a two-terminal device as

\begin{eqnarray}
\label{eq:29}
\begin{aligned}
\overline{{J}}_{l}
\int_0^{ \omega}
\frac{d \mathcal{E}}{2\pi}
\big(
\big),
\\&
\Big(
[G^{rF}]^\dagger (\mathcal{E})
\Big),
\end{aligned}
\end{eqnarray}

where $T_{ll^\prime}(\mathcal{E})$ is the two terminal transmission function.  Here, $[\bar{\Gamma}_{l/l^\prime}^{rF}]$ is the block diagonal matrix made of $\bar{\Gamma}_{l/l^\prime}^{r}$.

§.§ Floquet Green's function for interacting systems

For interaction systems, we introduce an ansatz for the M-like Floquet Green's function related to the equation-of-motion method <cit.>.

We should consider the spin degree of freedom in the total Hamiltonian, Eqs. (<ref>)-(<ref>), by modifying annihilation (creation) operators $\hat{d}_{i}^{(\dagger)}/\hat{c}_{i}^{(\dagger)}
\rightarrow \hat{d}_{i,\sigma}^{(\dagger)}/\hat{c}_{i,\sigma}^{(\dagger)}$ in which $\sigma=\uparrow,\downarrow$. The (time-independent) electron-electron interaction, $\hat{H}_{ee}$, is given by a Hubbard-like two-body interaction as

\begin{eqnarray}
\label{eq:30}
\begin{aligned}
\hat{H}_{ee} &=
\sum_{i}
u_i \,
\hat{d}_{i,\sigma}^\dagger \hat{d}_{i,\sigma}
\hat{d}_{i,\tilde{\sigma}}^\dagger \hat{d}_{i,\tilde{\sigma}}. \\
\end{aligned}
\end{eqnarray}

Here, $\tilde{\sigma}$ refers to opposite spin of $\sigma$. In case of time-independent $h_{ij}$ (None-Floquet case) and with truncating to double-particle Green's function, $G^{(2)r}(\mathcal{E})$, one can obtain a close form for the retarded Green's function as

\begin{eqnarray}
\label{eq:31}
\mathcal{E}[I]
-[h]-[\bar{\Sigma}^{r}_{\sigma}] \Big)
\\
\label{eq:32}
\mathcal{E}[I]
-[h]-[U]-[\bar{\Sigma}^{r}_{\sigma}] \Big)
[G^{(2)r}_{\sigma}](\mathcal{E})=[\langle n_{\tilde{\sigma}}\rangle],~~~~~
\end{eqnarray}

where, $[U]$ refers to a diagonal matrix consist of the elements $\{u_i\}$. The matrix $[\langle n_{\tilde{\sigma}}\rangle]$ is a diagonal matrix consist of occupation numbers $\langle n_{i\tilde{\sigma}}\rangle$. Although, Eqs. (<ref>) and (<ref>) suggest a self-consistent procedure but $\langle n_{i\tilde{\sigma}}\rangle=1/2$ gives a good quantitative prediction on the onsets of current plateaus (charging effects) related to the coulomb blockade phenomenon.

In previous sections where we discussed the Floquet transformation, we saw that the none time-dependent components of Green's function equation turn into fixed diagonal blocks. Here, we extend this observation to the interacting case, and hence, we expect the following M-like Floquet version for the Eqs. (<ref>)-(<ref>) as

\begin{eqnarray}
\label{eq:33}
\mathcal{E}[I^F]
-[h^F]-[\bar{\Sigma}^{rF}_{\sigma}] \Big)
\\
\label{eq:34}
\mathcal{E}[I^F]
-[h^F]-[U^F]-[\bar{\Sigma}^{rF}_{\sigma}] \Big)
[G^{(2)rF}_{\sigma}]=[\langle n_{\tilde{\sigma}}\rangle ^F].~~~~~~~~~~
\end{eqnarray}

In addition, based on the EOM technique and the above argument, we would have also the following expressions for the M-like Floquet lesser Green's functions

\begin{eqnarray}
\label{eq:35}
\begin{aligned}
\mathcal{E}[I^F]
\Big)
\quad \quad \quad \quad
\\&
\end{aligned}
\\
\label{eq:36}
\begin{aligned}
\mathcal{E}[I^F]
\Big)
\\&
\end{aligned}
\end{eqnarray}

Note that, the lesser Green's function for the interacting case is different than the none-interacting one. Consequently, the expectation value of the number operator and the average of the number operator over one period should follow Eq. (<ref>). However, the expression for the average current per cycle is identical to the non-interacting case. It is worth mentioning that the Landauer formula is not correct in the case of interacting systems and hence there would be no Landauer-like expression for the average current per cycle.

§.§ Floquet Quantum Master Equation

In the following we explore two types of Floquet Quantum Master Equation (FQME) whereby we named 1) Hilbert-Space FQME (HS-FQME), and 2) Floquet-Space FQME (FS-FQME).

Our starting point for deriving HS-FQME is the well-known Redfield QME in the Schrödinger picture whereas we will start from the Floquet Liouville von-Neumann (FLvN) equation for deriving HS-FQME.

§.§.§ Hilbert-Space FQME

When the Hamiltonian of the System is time-dependent (not restricted to the time-periodic), the dynamics of the System's density operator can be notified by the following Redfield equation

\begin{eqnarray}
\label{eq:37}
\begin{aligned}
\frac{\partial \hat{\rho}_S(t)}{\partial t}
\hat{\rho}_S(t)\right]-
\\&
\int_0^{\infty} \operatorname{Tr}_B
\left[
\hat{H}_{SB},
\left[\tilde{\tilde{H}}_{SB}(t, \tau), \hat{\rho}_S(t) \otimes \hat{\rho}_B
\right]\right] d \tau,
\end{aligned}
\\
\label{eq:38}
\begin{aligned}
\tilde{\tilde{H}}_{SB}(t, \tau)=
\sum_{i}
\tilde{\tilde{C}}_i^{\dagger}(\tau) \tilde{\tilde{d}}_i(t, \tau)
\tilde{\tilde{d}}_i^{\dagger}(t, \tau)
\tilde{\tilde{C}}_i(\tau)~~~~~~~~~~~~
\end{aligned}
\end{eqnarray}

where $\hat{\rho}_S(t)$ and $\hat{\rho}_B$ are the System and bath many-body density operators.
$\operatorname{Tr}_B$ indicates tracing over the bath degrees of freedom.
In the,
Eq. (<ref>), $\tilde{\tilde{C}}_i(\tau)=\sum_{l k} V_{l k,i}^{\ast} \hat{c}_{l k} e^{i \, \epsilon_{lk}  \tau}$ and it determines how the bath part of the $\hat{H}_{SB}$ evolve in the bath's time scale, $\tau$.

On the other hand, $\tilde{\tilde{d}}_i(t, \tau) =\hat{U}_S(t, t-\tau) \hat{d}_i \hat{U}_S^{\dagger}(t, t-\tau)$ describes the evolution of the System part of the $\hat{H}_{SB}$ in which the time evolution operator is $\hat{U}_S(t, t-\tau)=\mathcal{T} exp\big({-i \int_{t-\tau}^t \hat{H}_S(s) d s}\big)$  ($\mathcal{T}$ refers to the time-ordering operator). Same definition holds for $\tilde{\tilde{d}}^{\,{\dagger}}_i(t, \tau)$. It was Shirley who first formalized the time evolution operator for the time-periodic Hamiltonian. The time evolution operator can be expressed as

\begin{eqnarray}
\label{eq:39}
\begin{aligned}
\hat{U}_S(t, t_0)=
\sum_n\langle n|e^{-i \hat{H}_S^F\left(t-t_0\right)}| 0\rangle e^{i n \omega t},
\end{aligned}
\end{eqnarray}

where $| n \rangle$ ($| 0 \rangle$) refers to the $n$th ($0$th) state in the Fourier basis set. Here, the Floquet many-body Hamiltonian defines identical to the one-body counterpart as $\hat{H}_S^F=\sum_{n,m}\hat{H}_{S\, m n}-\hat{I} (m \omega) \delta_{mn}$ in which $\hat{H}_{S\,m n}$ defines similar to Eq. (<ref>).

In general, disassembling the double commutator in Eq. (<ref>) results in eight major integrands. In addition, having two parts in Eq. (<ref>) results in a total of sixteen integrands. From these, only eight integrands should be considered (terms in which one of the operators multiples to the conjugate transpose of the other). For example, the first and second non-vanishing terms (Dissipators) can be given as

\begin{eqnarray}
\label{eq:40}
\begin{aligned}
%\int_0^{\infty} \operatorname{Tr}_B
%\Big(\hat{H}_{SB}\,\tilde{\tilde{H}}_{SB}(t, \tau)
%\,\hat{\rho}_S(t) \otimes\hat{\rho}_B\Big)d \tau=
\sum_{i,j}
\int_0^{\infty}
\operatorname{Tr}_B
\Big(
\big(
\hat{C}_{i}^\dagger \hat{d}_i+
\hat{d}_i^\dagger \hat{C}_{i}
\big)
\big(
\tilde{\tilde{C}}_j^{\dagger}(\tau) \tilde{\tilde{d}}_j(t, \tau)
\tilde{\tilde{d}}_j^{\dagger}(t, \tau)
\\&
\tilde{\tilde{C}}_j(\tau)
\big)
\,
\hat{\rho}_S(t) \otimes \hat{\rho}_B
\Big)
d \tau=
\sum_{i,j}
\int_0^{\infty}
\operatorname{Tr}_B
\Big(
\hat{C}_{i}^\dagger \hat{d}_i
\tilde{\tilde{d}}_j^{\dagger}(t, \tau)~~~~~~
\\&
\tilde{\tilde{C}}_j(\tau)
\hat{\rho}_S(t) \otimes \hat{\rho}_B
\hat{d}_i^\dagger \hat{C}_{i}
\tilde{\tilde{C}}_j^{\dagger}(\tau)
\tilde{\tilde{d}}_j(t, \tau)\hat{\rho}_S(t) \otimes \hat{\rho}_B
\Big)d \tau.
\end{aligned}
\end{eqnarray}

For the sake of brevity, we only focus on the first Dissipator which can be further simplified to

\begin{eqnarray}
\label{eq:41}
\begin{aligned}
\sum_{i,j}
\int_0^{\infty}
\operatorname{Tr}_B
\Big(
\hat{C}_{i}^\dagger \tilde{\tilde{C}}_j(\tau)
\hat{\rho}_B
\otimes
\hat{d}_i
\tilde{\tilde{d}}_j^{\dagger}(t, \tau)
\hat{\rho}_S(t)
\Big)d \tau= \\
\sum_{i,j,kl}
\int_0^{\infty}
\big(
V_{l k,i} V_{l k,j}^{\ast}
e^{i \,\epsilon_{lk} \tau}
f(\epsilon_{lk}, \mu_l)
\hat{d}_i
\tilde{\tilde{d}}_j^{\dagger}(t, \tau)
\hat{\rho}_S(t)
\big)
d \tau.
\end{aligned}
\end{eqnarray}

Note that
\Big(\hat{c}_{lk}^{\dagger}
\hat{c}_{l^\prime k^\prime} \hat{\rho}_B(\mu)
\Big)=f(\epsilon_{lk}, \mu_l)\delta_{k,k^\prime}\delta_{l,l^\prime}$. Until now, the derivation procedure is not far from the non-Floquet QME. However, complexity in the time dependency of the $\tilde{\tilde{d}}_j^{\,{\dagger}}(t, \tau)$ requires Shirley's Floquet Formalism, Eq. (<ref>), as

\begin{eqnarray}
\label{eq:42}
\begin{aligned}
\tilde{\tilde{d}}_j^{\,\dagger}(t, \tau)=
\sum_{n, m}
\langle n|
e^{-i \hat{H}^F_S\tau}
\hat{d}_j^{\,\dagger}
\langle 0|
e^{i \hat{H}^F_S\tau}
| m\rangle
e^{i(n-m) \omega t}.~~~~
\end{aligned}
\end{eqnarray}

To proceed, we shall use the eigenbasis of the Floquet System Hamiltonian $\hat{Y}^{\dagger} \hat{H}^F_S \hat{Y}=\hat{\Lambda}^F$ in which $\hat{Y}$ is the rotating operator. Then, the above relation can be simplified to

\begin{eqnarray}
\label{eq:43}
\begin{aligned}
\tilde{\tilde{d}}_j^{\,\dagger}(t, \tau)=
\sum_{n, m}
\langle n|
\hat{Y}
e^{-i \hat{\Lambda}^F\tau}
\hat{\mathcal{D}}_{j}^{{o}\,\dagger}
e^{i \hat{\Lambda}^F\tau}
\hat{Y}^{\dagger}
| m\rangle
e^{i(n-m) \omega t},~~~~
\end{aligned}
\end{eqnarray}

where $\hat{\mathcal{D}}_{j}^{{o}\,\dagger}=
\hat{Y}^{\dagger}
\big(|0\rangle
\hat{d}_j^{\,\dagger}
\langle 0|\big)
\hat{Y}$ and we used the property $\hat{Y}^{\dagger}\hat{Y}=\hat{I}$ multiple times.

Next, we should consider a basis set $\{|a\rangle \}$ for the Hilbert space which allows us to present elements of density operator as $\rho_{ab}(t)= \langle a |\hat{\rho}(t)|b\rangle$ and turn operators in Eq. (<ref>) to matrices. In particular, it allows us to work with the (hybrid) Floquet basis, $\{| \gamma \rangle \}$ such that we can simplify the $e^{-i\hat{\Lambda}^F\tau}
\hat{\mathcal{D}}_{j}^{{o}\,\dagger}
e^{i \hat{\Lambda}^F\tau}$. A matrix element of this term is given by

\begin{eqnarray}
\label{eq:44}
\begin{aligned}
\Big(
e^{-i \hat{\Lambda}^F\tau}
\hat{\mathcal{D}}_{j}^{{o}\,\dagger}
e^{i \hat{\Lambda}^F\tau}
\Big)_{\gamma,\nu}=
e^{-i \Omega_{\gamma \nu} \tau}
\big(
\hat{\mathcal{D}}_{j}^{{o}\,\dagger}
\big)_{\gamma,\nu},
\end{aligned}
\end{eqnarray}

where $\Omega_{\gamma \nu}=(\mathcal{E}_{\gamma} -\mathcal{E}_{\nu})$. Clear separation of the two time scales ($t$ and $\tau$) allows us to perform the time integration over $\tau$ as

\begin{eqnarray}
\label{eq:45}
\begin{aligned}
e^{i \,(\epsilon_{lk}-\Omega_{\gamma \nu}) \tau}
d \tau=
\pi \delta(\epsilon_{lk}+\Omega_{\gamma \nu})
\big(
\frac {1}{\epsilon_{lk}+\Omega_{\gamma \nu} }
\big).~~
\end{aligned}
\end{eqnarray}

Neglecting the Cauchy's principle, $P$, and employing the energy-independent wide-band approximation, $\Gamma^l_{ij}=2\pi \sum_kV_{l k,i} V_{l k,j}^{\ast}  $ results in the following matrix expression for the first Dissipators

\begin{eqnarray}
\label{eq:46}
\begin{aligned}
\sum_{i,j,l}
\frac {\Gamma^l_{ij} }{2}
\sum_{n, m}
\langle n|
\circ
| m\rangle
e^{i(n-m) \omega t}
\end{aligned}
\end{eqnarray}

where $\circ$ refers to the Hadamard product. We should repeat this procedure for the other seven non-vanishing Dissipators. A more compact matrix form of the Floquet QME in the Hilbert space can then be given as

\begin{eqnarray}
\label{eq:47}
\begin{aligned}
\begin{gathered}
\frac{\partial \rho_S(t)}{\partial t}=
\left[
{H}_S(t), {\rho}_S(t)
\right]
\frac{\Gamma^l_{ij}}{2}
\Big(d_i\widetilde{d}_{jl}^{\, \dagger}(t)\rho_S(t)+~~~
\\
{d}_i^{\, \dagger} \overline{d}_{jl}(t) \rho_S(t)-d_i \rho_S(t) \widetilde{d}_{jl}^{\, \dagger}(t)
-{d}_i^{\dagger} \rho_S(t) \overline{d}_{jl}(t)
\Big)+\text { h.c. },
\end{gathered}
\end{aligned}
\end{eqnarray}

where $\widetilde{d}_{jl}(t)$ and $\overline{d}_{jl}(t)$ defined as

\begin{eqnarray}
\label{eq:48}
\widetilde{d}_{jl}^{\, \dagger}(t)
\sum_{n, m}
\langle n|
\,
\circ
\,
| m\rangle
e^{i(n-m) \omega t},~~~~
\\
\label{eq:49}
\overline{d}_{jl}(t)
&=&\sum_{n, m}
\langle n|
\,
\circ
\,
| m\rangle
e^{i(n-m) \omega t}.
\end{eqnarray}

We stress that within the Hilbert-Space FQME, $\widetilde{d}_{jl}(t)$ and $\overline{d}_{jl}(t)$ (and their hermitian conjugates) are time-dependent matrices. Finally, we can simplify the notation of Eq. (<ref>)  as ${\partial \rho_S(t)}/{\partial t}=-i
\left[
{H}_S(t), {\rho}_S(t)
\right]+\sum_{l} \mathcal{L}_l^{HS}(t)\rho_S(t)$.

§.§.§ Floquet-Space FQME

LvN equation is the starting point for varieties of QMEs (e.g. the Redfield QME). For a periodic Hamiltonian, the LvN equation can be transformed into the Floquet LvN equation which is very general and has a long history in the analysis of nuclear magnetic resonance signals (closed systems) <cit.>. The Floquet LvN equation is given by

\begin{eqnarray}
\label{eq:50}
\frac{\partial \hat{\rho}^F(t)}{\partial t}
\left[
\hat{H}^F, \hat{\rho}^F(t)
\right]
\\
\label{eq:51}
\hat{H}^F
\sum_{n}
\widehat{L}_n \otimes \hat{H}^{(n)}
+\widehat{N} \otimes \hat{\mathbf{1}},
\\
\label{eq:52}
\hat{H}^{(n)}
\frac{1}{T} \int_0^T dt\,  e^{i n \omega t} \hat{H}(t) \\
\label{eq:53}
\hat{\rho}^F(t)
\sum_{n}
\widehat{L}_n \otimes \hat{\rho}^{(n)}(t).
\end{eqnarray}

Here, $n$ is an integer in the space $[-N, N]$, $\widehat{L}_{n}$ is the Ladder operator of the order $n$, $\widehat{N}$ is the Number operator in the Fourier space, and $\hat{\mathbf{1}}$ is the identity operator in the Hilbert space. In the matrix form, $\widehat{L}_{n}$ represents a $N\times N$ off-diagonal matrix with ones in its off-diagonal and $\widehat{N}$ represents a diagonal matrix with the integer vector $\{-N,..., N\}$ placed on it's diagonal.

The major benefit of the Floquet LvN equation is that $\hat{H}(t)$ turns into a time-independent Hamiltonian $\hat{H}^F$.

Note that the operator coefficient $\hat{\rho}^{(n)}(t)$ is time dependent whereas $\hat{H}^{(n)}$ is not. The definition of $\hat{H}^{(n)}$ is consistent with Eq. (<ref>) and hence the $\hat{H}^F$ given in Eq. (<ref>) is effectively similar to those discussed before.

Assuming that $\hat{H}(t)$ is the total time-periodic Hamiltonian, one can partition the $\hat{H}^F$ as

\begin{eqnarray}
\label{eq:54}
\hat{H}^F
\hat{H}^F_S \otimes \hat{I}_B+
\hat{I}^F_S \otimes \hat{H}_B+
\hat{H}^F_{SB}, \\
\label{eq:55}
\hat{H}_{SB}^F
\sum_{i}
\hat{C}_{i}^\dagger \hat{\mathcal{D}}^{}_i +\mathrm{H.c.}
\end{eqnarray}

where $\hat{I}^F_S=\widehat{L}_0 \otimes \hat{I}_S$ and $\hat{\mathcal{D}}^{}_i= \widehat{L}_0 \otimes \hat{d}_i$. Then,  we can define the following rotation protocols

\begin{eqnarray}
\label{eq:56}
\tilde{H}^F_{BS}(t)
e^{i \hat{H}_B t}
e^{i \hat{H}^F_S t}
\hat{H}^F_{BS}
e^{-i \hat{H}_B t}
e^{-i \hat{H}^F_S t}, \\
\label{eq:57}
\tilde{\rho}(t)
e^{i \hat{H}_B t}
e^{i  \hat{H}^F_S t}
\hat{\rho}(t)
e^{-i \hat{H}_B t}
e^{-i \hat{H}^F_S t}.
\end{eqnarray}

The above rotations follow the transformation to the interaction picture in the context of conventional QME but, here within this version of Floquet QME, we have used $\hat{H}^F_{S/SB}$. Next, we can take the (initial factorization) assumption $\tilde{\rho}^{F}(t)=\tilde{\rho}^F_S(t) \otimes \tilde{\rho}_B$ and proceed just as the conventional QME. This eventually leads to the following Floquet Markovian relation in the Schrödinger frame

\begin{eqnarray}
\label{eq:58}
\begin{aligned}
\frac{\partial \hat{\rho}_S^F(t)}{\partial t}
&=-i\left[\hat{H}_S^F, \hat{\rho}_S^F(t)\right]-\int_0^{\infty} \operatorname{Tr}_B
\\
\quad
\left[\hat{H}^F_{SB},
\left[\tilde{\tilde{H}}^F_{S B}(\tau), \hat{\rho}_S^F(t) \otimes \hat{\rho}_B\right]\right] d \tau,
\end{aligned}
\end{eqnarray}

where $\tilde{\tilde{H}}^F_{S B}(\tau)=e^{-i \hat{H}_B \tau} e^{-i \hat{H}_S^F \tau} \hat{H}^F_{S B} e^{i \hat{H}_B \tau} e^{i \hat{H}_S^F \tau}$.

Because $\hat{\rho}_S^F(t)$ is in the Floquet space, we can employ Floquet basis, $\{| \gamma \rangle \}$, to present elements of the reduced Floquet density operator (Floquet Hamiltonian) as $\hat{\rho}^F_{S\, \gamma \nu}= \langle \gamma |\hat{\rho}_S^F(t)| \nu \rangle$

($H^F_{S\, \gamma \nu}(t)= \langle \gamma | \hat{H}_S^F | \nu \rangle$).
After taking the integrations over $\tau$ and employing Wide-band approximation, we arrive at the following matrix expression for the dynamic of the Floquet density matrix

\begin{eqnarray}
\label{eq:59}
\begin{aligned}
\begin{gathered}
\frac{\partial \rho^F_S(t)}{\partial t}=
\left[
{H}^F_S, {\rho}^F_S(t)
\right]
\sum_{i,j,l}
\frac{\Gamma_{ij}^{l}}{2}
\Big(\mathcal{D}_i\widetilde{\mathcal{D}}_{jl}^{\, \dagger}\rho^F_S(t)+~~~
\\
\mathcal{D}_i^{\, \dagger} \overline{\mathcal{D}}_{jl}\rho^F_S(t)
-\mathcal{D}_i \rho^F_S(t) \widetilde{\mathcal{D}}_{jl}^{\, \dagger}
-\mathcal{D}_i^{\dagger} \rho^F_S(t) \overline{\mathcal{D}}_{jl}
\Big)+\text { h.c. },
\end{gathered}
\end{aligned}
\end{eqnarray}

where $
\widetilde{\mathcal{D}}_{jl}^{\, \dagger} = {Y} \,f_l(-\Omega) \circ{\mathcal{D}}_{j}^{{}\,\dagger}\,{Y}^{\dagger}$
and  $\overline{\mathcal{D}}_{jl}={Y} \,f_l(\Omega) \circ{\mathcal{D}}_{j}^{} \,{Y}^{\dagger}$.
The Floquet-Space FQME is thus easier to be implemented as compared to the Hilbert-Space FQME. The disadvantage of Floquet-Space FQME is that it requires more memory since all its components are in the Floquet space. Finally, we can simplify the notation of Eq. (<ref>)  as ${\partial \rho^F_S(t)}/{\partial t}=-i
\left[
{H}^F_S, {\rho}^F_S(t)
\right]+ \sum_ {l}\mathcal{L}_l^{FS}\rho^F_S(t)$. Note that $\mathcal{L}_l^{FS}$ is time-independent.

§.§ Observable in Floquet QME

Within the Hilbert-Space FQME, it is straightforward to obtain expectation values of the particle number, $\langle \hat{n} \rangle$ and current passing through the contact $l$, $\langle {\hat{J}_l}\rangle$, <cit.> as

\begin{eqnarray}
\label{eq:60}
\langle \hat{n} \rangle(t)
Tr\big(\hat{n} \hat{\rho}(t)\big), \\
\label{eq:61}
\langle \hat{J}_l \rangle(t)
Tr\big(\hat{n} \mathcal{L}_l^{HS}(t) \hat{\rho}(t)\big).
\end{eqnarray}

To evaluate the expectation values of an observable in the Floquet-Space FQME, we need to project the operators from the Floquet space to the Hilbert space <cit.> as

\begin{eqnarray}
\label{eq:62}
\langle \hat{n} \rangle(t)
\sum_m
\langle m|
\hat{n}^F \hat{\rho}^F(t)
|0\rangle e^{i m \omega t}
\big),
\\
\label{eq:63}
\langle \hat{J}_l \rangle(t)
\sum_m
\langle m|
\hat{n}^F
\mathcal{L}_l^{FS}
\hat{\rho}^F(t)
|0\rangle e^{i m \omega t}
\big).
\end{eqnarray}

To evaluate $\overline{n}$ and $\overline{J}_l$ in both Floquet QME methods, we must evaluate the time-average of observables, $\hat{O}(t)$, numerically by using $\overline{O}=1/T \int_0^T \langle \hat{O} \rangle(t) dt$.

§ NUMERICAL RESULTS

To demonstrate the implications of light-matter interactions on the transport characteristics of an open quantum system, and to examine consistency between the discussed methods (V-like Floquet NEGF, M-like Floquet NEGF, Hilbert-Space FQME, and Floquet-Space FQME), we have first taken a spinless two-level quantum dot placed between two metallic contacts. This two-level dot then interacts with a monochromatic light. The one-body Hamiltonian of a such quantum dot can be given by

\begin{eqnarray}
\label{eq:64}
\begin{aligned}
[h] (t)= \left[
\begin{array}{cc}
\varepsilon_{1} & A \cos (\omega t) \\
A \cos (\omega t) & \varepsilon_{2}
\end{array}
\right],
\end{aligned}
\end{eqnarray}

where $\omega$ is the frequency of light and $A$ represents the light intensity. Off-diagonal elements on the above Hamiltonian represent the dipole approximation in the context of light-matter interaction (Rabi model in quantum optics).
This Hamiltonian is connected to the left (right) contacts with small $\Gamma_{L(R)}$ (weak coupling regime). The schematic of such a configuration is given in Fig. <ref>.

(a) Schematics of a quantum transport system comprised of a two-level System under light illumination. The two-level System with monochromatic light can roughly be thought of as multiple quasi-levels.

§.§ Results for Floquet NEGF methods

As the first example, we calculate the $\overline{n}$ and $\overline{J}_L$ (current passing through the left contact) using Eqs. (<ref>-<ref>) for a range of external bias, $V_{SD}=\mu_R-\mu_L$ and multiple light frequencies.

Hereafter, we will fix the right contact (Drain) at $\mu_R=-0.4$ and sweep over $\mu_L$.

In the absence of the light, and when off-diagonals of $[h]$ are zero, there should be only two plateaus on the $I-\mu_L$ curve such that the onset of these two plateaus are at $\varepsilon_{1}$ and $\varepsilon_{2}$.

When the $A$ is large enough and the frequency of light is tuned to the difference between two levels, $\hbar\omega=\varepsilon_{2}-\varepsilon_{1}$, then maximum alteration of the quantized conductance will happen such that the number of plateaus is doubled. This has been shown in Figs. <ref>(a) and <ref>(b). The extra quantized plateaus can be the hallmark of an electron-photon hybrid state (Floquet states).

Floquet NEGF: (a) One-period time-average of the number operator, $\overline{n}$, for multiple driving frequencies. Parameters are: $\varepsilon_{1}=-0.1~eV$, $\varepsilon_{2}=0.1~eV$, $A=0.1~eV$,  $N=3$, $T=4.2~K$, $\Gamma_{L}=\Gamma_{R}=0.0025~eV$, $\mu_R=-0.4~V$.

In extreme off-tuning situations (associated with $\hbar\omega=0.05$ and $\hbar\omega=0.4$), the transport characteristics should be almost identical to the non-Floquet transport characteristics. We can see it from the first and last curves in Figs. <ref>(a) and <ref>(b).

At low frequencies, the M-like Floquet NEGF does not provide accurate results with only $N=3$. Hence, V-like Floquet NEGF delivers better numerical performance in low frequencies.

Next, we fixed the driving frequency at $\hbar \omega=0.2~(eV)$ and calculated the two observables for a range of external bias, $\mu_L$, and multiple driving amplitudes, $A$. Results are plotted in Figs. <ref>(a) and <ref>(b). Here, we can see, that the light intensity, $A$, shifts the onset of plateaus symmetrically around $\mu_L=0$.

Floquet NEGF: (a) One-period time-average of the number operator, $\overline{n}$, and (b) one-period time-average of the left contact's current operator, $\overline{J}_L$, for multiple driving amplitudes at $\hbar \omega=0.2~(eV)$. All other parameters are the same as Fig. <ref>.

§.§ Results for Floquet QME methods

We now use the spinless many-body operators $\hat{d}/\hat{d}^\dagger$ to make the time-dependent many-body Hamiltonian and its Floquet counterpart.

Then, by identifying the many-body basis set, an initial state (initial many-body density matrix), and the Fourier space (hence the many-body Floquet space), we will be able to perform dynamical solutions to Eqs. (<ref>) and Eq. (<ref>) and afterward evaluate time-averages of the observables, Eqs. (<ref>)-(<ref>). Note that, we have used the Runge-Kutta 4th-order method to solve Floquet QMEs. The time step has to be adjusted carefully. The rest of the parameters are the same as what is given in Fig. <ref>.

The numerical cost of these solutions is much higher than the Floquet NEGF methods and hence, we only plotted the time averages of the number/current operator versus voltage for the perfect tuning case, (i.e., $\hbar\omega=0.2~eV$) with the amplitude $A=0.1~(eV)$. Results are plotted in Fig. <ref>. To perform these calculations, we have considered the two-level system to be empty at $t=0$ such that there will be no Rabi oscillation as time goes on.

Floquet QME: (a) time-dependent expectation value of the number operator evaluated by the Hilbert space Floquet QME. (b) Same as (a) evaluated by the Floquet space Floquet QME. (c) and (d) One-period time-average of the number operator, and the left contact's current operators, $\overline{n}$, $\overline{J}_L$, at $\hbar \omega=0.2~(eV)$ for a range of external bias. All other parameters are the same as Fig. <ref>.

Firstly, we can see no pronounced oscillation in the $\langle \hat{n} \rangle(t)$ of the Hilbert-Space FQME in the short-time limit  [Fig. <ref> (a)] whereas $\langle \hat{n} \rangle(t)$ of the Floquet-Space FQME, shows oscillation feature [Fig. <ref> (b)] in the early times. Such oscillations became steady as the time went by.
Secondly, the two Floquet QME methods are identical in the long-time limit. Hence, calculated $\overline{n}$, and $\overline{J}_L$ by the two Floquet QME is almost identical, see the solid and dashed lines in Figs. <ref> (c) and (d). Thirdly, $\overline{n}$, and $\overline{J}_L$ is also identical to the Floquet NEGF methods, compare Figs. <ref> (c) and (d) with the central curves in Figs. <ref> (a) and (b). Note that, we have found both Floquet QMEs require much smaller time steps as compared to the none Floquet QME. However, in comparison, the Hilbert-Space FQME is less sensitive to time step than the Floquet-Space FQME.

§.§ Results for Floquet driven interacting system

Finally, we present a few numerical calculations performed for the spin case with and without interactions using the Floquet QME.

Floquet QME for the interacting case: (a) $\overline{n}$ vs the external bias when light is absent ($A=0$). $[U]\neq 0$ gives rise to the extra quantized plateaus. We take $u_1=u_2=0.1$, for the case where $[U]\neq 0$. All other parameters are the same as Fig. <ref>. (b) $\overline{n}$ vs the external bias in the presence of light at $\hbar \omega=0.2~(eV)$. (c) and (d) show $\overline{n}$ vs the external bias when only dipole approximation is used in the time-periodic Hamiltonian.

When $[U]$ and $A$ (none-Floquet case) are both zero, having two spins, obviously, should result in the degeneracy factor $2$ on the observed plateaus. In the presence of electron-electron interaction, $[U] \neq 0$, the number of plateaus increases such that the onset of the spin-up plateaus shifts with respect to the spin-down plateaus by the values $\{u_i\}$ or vice versa. This observation is verified with our Floquet QME when $A=0$, see the dashed red line in Fig. <ref> (a). The correct height ratio for the interacting case should be 2:1 <cit.>. Here in this section, we only show $\overline{n}$ vs the external bias because $\overline{J}$ follows the same pattern.

Note that, the equation of motion technique in the NEGF, can reproduce the correct onsets for the plateaus however it does not often deliver the correct height for the plateaus. For a spin system and within the NEGF, we are required to include two identical time-dependent Hamiltonian, one for spin $\downarrow$ and another for $\uparrow$, whereas, within the QME, we should only employ the spin version of many-body operators $\hat{d}/\hat{d}^\dagger$.

In Fig. <ref> (b) we examined what will happen to the well-quantized plateaus of a two-level spin system subjected to a Floquet deriving with ($[U]\neq 0$) and without two-body interaction ($[U]= 0$). Here, we include the dipole approximation for both spins [as indicated by $H(t)_{\uparrow \downarrow}$]. We can see there is an interplay between the interaction and the Floquet driving such that even in the case of $[U]= 0$, the quantized Floquet plateaus of the spinless case do not simply multiply to the degeneracy factor $2$. We can perceive that the electron-electron interaction in the presence of light, overall, lowers the population probability.

To further speculate this interplay, we have evaluated $\overline{n}$ vs the external bias when the dipole approximation is taken into account only for one of the spins with and without interactions. Results of $[U]= 0$ and $[U] \neq 0$ are plotted in Figs. <ref> (c) and (d) respectively. Firstly, we can see the non-interacting spin system shows different quantization when only one spin is allowed to be excited. Secondly, electron-electron interaction makes the finer quantized population of the spin system smoother [compare Figs. <ref> (c) with (d)].

§ CONCLUSION

In the present work, we have derived four full Floquet-based formulations for the quantum transport through a time-periodic dot. In particular, we have shown the two Floquet Green's function formalisms rely on the expansion of Green's function, the first  Floquet quantum master formalism, HS-FQME, relies on the Floquet time evolution operator, and the second Floquet quantum master formalism, FS-FQME relies on expanding the Hilbert space to the Floquet space.

In addition, we have connected the equation of motion for each formulation to the time-averaged observables.

Through a quantum dot model Hamiltonian, we quantitatively demonstrated that all Floquet quantum transport methods provide consistent results for non-interacting Hamiltonian in the weak coupling regime. The frequency and amplitude of time-periodic driving can significantly manipulate the time-average of observables.

We have found an intense light can shift the one-set and the number of quantized plateaus in the voltage-current curve provided choosing an appropriate light frequency.

As a drawback (even with parallel computation), we found Floquet quantum master equation frameworks to be computationally expensive.

We also found that Floquet QME methods need much smaller time steps compared to the conventional QME. As usual with QME,
when the lead-dot coupling or temperature decrease, the number of time step has to be increased. While the QME is the most practical method for interacting systems in the weak regime, it is computationally challenging to employ the full Floquet quantum master equation formalisms for large many-body systems. We believe the power of Floquet NEGF and Floquet  QME methods goes well beyond the considered example and one needs to combine Floquet QMEs with newly developed machine learning methods to fully benefit from them.

This material is based upon work supported by the National Natural Science Foundation of China and (NSFC No. 22273075). V.M. acknowledges funding from the Summer Academy Program for International Young Scientists (Grant No. GZWZ[2022]019). W.D. acknowledges the startup funding from Westlake University.

§ CONFLICT OF INTEREST

There are no conflicts to declare.

§ DATA AVAILABILITY STATEMENT

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