SciPost Phys. Core 4, 033 (2021)

Nonequilibrium steady states in the Floquet-Lindblad systems:
van Vleck’s high-frequency expansion approach
Tatsuhiko N. Ikeda1? , Koki Chinzei1 and Masahiro Sato2
1 Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan
2 Department of Physics, Ibaraki University, Mito, Ibaraki 310-8512, Japan? tikeda@issp.u-tokyo.ac.jp

Abstract
Nonequilibrium steady states (NESSs) in periodically driven dissipative quantum systems are vital in Floquet engineering. We develop a general theory for high-frequency
drives with Lindblad-type dissipation to characterize and analyze NESSs. This theory isbased on the high-frequency (HF) expansion with linear algebraic numerics and without
numerically solving the time evolution. Using this theory, we show that NESSs can deviate from the Floquet-Gibbs state depending on the dissipation type. We also show thevalidity and usefulness of the HF-expansion approach in concrete models for a diamondnitrogen-vacancy (NV) center, a kicked open XY spin chain with topological phase transition under boundary dissipation, and the Heisenberg spin chain in a circularly-polarizedmagnetic field under bulk dissipation. In particular, for the isotropic Heisenberg chain,
we propose the dissipation-assisted terahertz (THz) inverse Faraday effect in quantum
magnets. Our theoretical framework applies to various time-periodic Lindblad equations
that are currently under active research.
Copyright T. N. Ikeda et al.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 26-07-2021
Accepted 09-11-2021
Check forupdatesPublished 24-12-2021
doi:10.21468/SciPostPhysCore.4.4.033

Contents
1 Introduction

2

2 Floquet–Lindblad Equation (FLE)

43 van Vleck high-frequency expansion of Liouvillian

5

4 General aspects of NESS solution

8

5 Phenomenological time-independent dissipators
5.1 High-frequency expansion
5.2 Example 1: Open XY Chain with Boundary Dissipation
5.3 Example 2: Three-level system

10
11
11
15

6 Microscopically-derived dissipators by weak thermal contact
6.1 Floquet-Lindblad equation obtained by RWA

18
20

1

SciPost Phys. Core 4, 033 (2021)

6.2
6.3

6.1.1 Floquet-Lindblad equation and its characteristic properties
6.1.2 Conditions for Floquet-Gibbs state (FGS)
6.1.3 Implementation of high-frequency expansion
Example 3: Three-level system revisited and deviation from FGS
Example 4: Inverse Faraday Effect in Heisenberg chain and dissipation-assisted
Floquet engineering

20
23
24
24
29

7 Conclusion

33

A van Vleck high-frequency expansion in isolated systems

34

B Commutator calculations

36

C Majorana representation of the XY chain

36

D Derivation of time-dependent FLE for weak thermal contact
D.1 Interaction picture
D.2 Schrödinger picture
D.3 Time evolution of density matrix
D.3.1 Nondegenerate quasienergies
D.3.2 Degenerate quasienergies
D.4 Conditions for FGS when quasienergies are degenerate

38
38
41
42
42
44
44

References45

1

Introduction

Periodically driven quantum systems have seen a resurgence of interest motivated by laser
technology advancement and theoretical developments [1–5]. Such systems are theoreticallydescribed by the Floquet theorem [6, 7], which enables us to study nonequilibrium states ofmatter systematically [8, 9]. An important application is the Floquet engineering, i.e., artificially creating useful functionalities of physical systems by designing an appropriate driving
protocol. Floquet engineerings for various systems have been proposed theoretically and realized experimentally: the dynamical localization [10–14], Floquet topological states [15–19],
Floquet time crystals [20–22], the inverse Faraday effect [23–25], etc.
Despite those extensive studies, most of their theoretical analyses neglect the dissipation
effect and focus mainly on well-controlled artificial systems and ultraclean materials. Considering dissipation is important in two ways. First, most physical systems such as genericmaterials contact their environment, and dissipation is not negligible. Second, without dissipation, the Floquet-engineered states are eventually broken by injected heat, and the system
becomes the featureless infinite-temperature state [26–28]. Thus, in isolated systems, the Floquet engineering is usually considered in the Floquet prethermalization regime of finite time
window [29–33]. In dissipative systems, nonequilibrium steady states (NESSs), where the energy injection by the periodic driving is balanced with the energy dissipation, are establishedand give us the opportunity to the long-lived Floquet engineering.
There have been several approaches to analyze NESSs in periodically driven dissipative
systems. The nonequilibrium Green function method is a successful approach and particularly
useful in electron systems coupled to several leads [34–37]. Several systems have been ana2

SciPost Phys. Core 4, 033 (2021)

lyzed with this approach, such as the Floquet topological insulators [38], strongly-correlated
electrons [39], and so on. Quantum master equation [40, 41] is another useful formulation to
analyze various periodically-driven dissipative systems and has been applied to periodic thermodynamics [42–44], Bose-Einstein condensation [45], Floquet-band occupation [46], and soon. Fermi’s golden rule is a similar technique to define the NESS [47].
Among quantum master equations, the Floquet-Lindblad equation (FLE) is a powerful andflexible description for periodically-driven systems subject to Markovian dissipation [48–51].
The FLE is an extension of the Lindblad (or Gorini-Kossakowski-Sudarshan-Lindblad [52, 53])
for time-periodic systems and has nice properties such as the complete positivity and trancepreserving condition. However, it is generally hard to solve each FLE and find the NESS except
for some special cases [54–56]. Recently, based on the high-frequency expansion method, apart of the present authors [57] have analytically solved NESS in a phenomenological FLE.
However, there are two major remaining issues to be solved. First, that work assumed thatthe dissipation is time-independent and satisfies the detailed-balance condition, but dissipation without these properties has also attracted attention in the literature. Second, the NESS
was solved at the leading-order approximation for 1/ω (ω is the driving frequency), and thesystematic extension to higher orders has not been studied.
In this paper, we extensively study the NESS in the FLE for high-frequency drives. For thispurpose, we formulate how to characterize and obtain the NESS using the high-frequency (HF)
expansion for the Liouvillian (Lindbladian) from the van Vleck point of view. Compared to the
Floquet–Magnus viewpoint [50], the van Vleck approach involves fewer terms and enables us
to study NESSs, including the micromotion. Following the general formulation, we apply thismethod to example models: an effective three-level model for a nitrogen-vacancy (NV) centerin diamond, a periodically driven topological XY spin chain, and the Heisenberg spin chain
under a circularly polarized magnetic field. The HF-expansion approach enables us to analyze
the NESS systematically by numerically obtaining an eigenstate of the effective Liouvillian
without resorting to the direct numerical integration of the FLE. From our formalism, we also
reveal the condition that the NESS approaches or deviates from the Floquet-Gibbs state.
The rest of this paper is organized as follows. In Sec. 2, we introduce the FLE and formulate
our problem of finding the NESSs that we address in this paper. In Sec. 3, we develop the vanVleck HF expansion for Liouvillians. Importantly, we generally prove that a zero-mode exists
for the effective Liouvillian at each order of the HF expansion. The zero mode is shown tocorrespond to the NESS in Sec. 4, where general aspects of the NESS are discussed within the
FLE. In Sec. 5, we apply the HF expansion methods to time-independent dissipators, including two models. We first discuss an open XY spin chain under a periodic drive and boundary
dissipation and analyze its topological phase diagram together with its stability against dissipation. We second analyze an effective model for the NV center in diamonds, demonstrating
how the HF-expansion method gives the NESS accurately. In Sec. 6, we apply the HF expansion to a class of time-dependent dissipators derived from the system-bath coupling and therotating-wave approximation (RWA). We derive such FLEs with a slight generalization of previous studies in that one allows the Floquet quasienergies to be degenerate, and we reveal
the conditions for the NESS being approximated by the Floquet-Gibbs state. We then apply
these general discussions to two concrete example models: the effective model for the NV
center and the Heisenberg spin chain under a circularly polarized magnetic field. For the NV
center model, unlike with the time-independent dissipator, we show that the relationship between the driving frequency and the bath spectral cutoff brings about a nontrivial effect on
the NESS. Namely, when the photon energy ħhω is below the cutoff, some photon-exchange
processes are active, and the NESS cannot be described by the Floquet-Gibbs state. Even inthis case, the HF expansion is valid and enables us to obtain the NESS with the systematic
improvement of accuracy. For the Heisenberg chain, we first provide a microscopic theory for3

SciPost Phys. Core 4, 033 (2021)

the dissipation-assisted inverse Faraday effect in magnetic insulators, i.e., the emergence oftime-averaged magnetization due to the circularly polarized ac magnetic field. Unlike related
studies, the magnetization is activated not by magnetic anisotropy of the Hamiltonian (theHeisenberg spin chain is isotropic) but by dissipation, i.e., the system-bath coupling. Finally,in Sec. 7, we summarize our work and present some outlooks.

2

Floquet–Lindblad Equation (FLE)

In this section, we introduce the most general form of the Floquet–Lindblad equation (FLE)that we study in this work and formulate our problem.
Let us consider the quantum master equationdρ(t) = L t [ρ(t)] .
dt(1)Here, ρ(t) is the density operator describing the quantum state of the system of interest, andL t is the Liouvillian superoperator generating dynamics. We use the term “Liouvillian” forthe right-hand side of Eq. (1) by analogy with the Liouville equation in classical mechanics.
Throughout this paper, we use calligraphic symbols for superoperators. We suppose that the
Liouvillian is periodic in time: L t+T = L t , where T is the period. In the following, we shalluse the corresponding angular frequency ω ≡ 2π/T .
We suppose that the Liouvillian is of Lindblad (or GKSL [52, 53]) form. Namely, the Liouvillian consists of the Hamiltonian and dissipator parts:

L t (ρ) = H t (ρ) + D t (ρ) ,(2)

H t (ρ) ≡ −i[H(t), ρ] .

(3)

with

Here, H(t) = H(t + T ) is the time-dependent Hamiltonian (involving the Lamb shift contribution [40]) in general). The dissipator is described by the jump operators Lα (t) as

X
1 †
D t (ρ) =Lα (t)ρ Lα† (t) −Lα (t)Lα (t), ρ .
(4)2αWe assume that D t is also periodic: D t+T = D t . Under this assumption, we further assume
that each L k (t) is also periodic.1 In the following, we call L t as the Lindbladian when weemphasize that it is a Liouvillian of Lindblad form. We call the quantum master equation (1)generated by a time-periodic Lindbladian the FLE.
The periodicity in time of the FLE enables systematic analysis by Fourier expansions. We
Fourier-expand the Hamiltonian and dissipator superoperators as
XX
Ht =Hm e−imωt ; D t =Dm e−imωt .
(5)
m

m

Here, the Fourier components Hk are simply given by

Hm (ρ) = −i[H m , ρ] ,(6)

This assumption can be slightly relaxed as follows: For each α, there exists θα (t) ∈ R such thatLα (t) = e iθα (t) LαP (t) with periodic part LαP (t + T ) = LαP (t). In such cases, we can replace Lα (t) by LαP (t) without changing D t since the nonperiodic phase factors e iθα (t) cancel between Lα (t) and Lα† (t). Thus, under thisassumption, we can assume Lα (t) are periodic without loss of generality.
1

4

SciPost Phys. Core 4, 033 (2021)

P
where H m are defined in H(t) = m H m e−imωt . In contrast, Dm are a little more complicatedbecause it is not linear in Lα (t) and Lα† (t). Fourier expanding them asXX†Lα (t) =Lα,m e−imωt ; Lα† (t) =Lα,me+imωt ,(7)mm

we haveDm (ρ) =X
α,n†Lα,m−n ρ Lα,n
−

©
1¦ †
Lα,n Lα,m−n , ρ .
2

(8)

†
We note that Lα (t) is not necessarily Hermitian, and Lα,−m 6= Lα,min general.
The formal solution of the FLE (1) is represented by the propagator V (t, t 0 ) as

ρ(t) = V (t, t 0 )ρ(t 0 ) ,

(9)

with
0

V (t, t ) = exp+

Z t



Ls ds ,

(10)t0

where exp+ denotes the time-ordered exponential. One merit of assuming that L t is a Lindbladian at each time t is that the propagator V (t, t 0 ) is guaranteed to be a completely positive
and trace preserving (CPTP) map [52, 53]. Thus, ρ(t) is qualified as a density operator at
any t during time evolution, which is not guaranteed in other master equations such as the
Redfield equation [58, 59].
Finally, we discuss how the NESS is characterized in the FLE, following Ref. [60]. To thisend, we take an initial time t = 0 and consider a long-time evolution to t. To utilize theperiodicity L t+T = L t and hence V (t + T, t 0 + T ) = V (t, t 0 ), it is useful to denote t = t 0 + `T ,where 0 ≤ t 0 < T and ` ∈ Z. Together with the property V (t, t 0 ) = V (t, t 00 )V (t 00 , t 0 ) for any
t 00 , we haveρ(t) = V (t 0 , 0)V F` ρ(0) ,(11)where V F ≡ V (T, 0) is the one-cycle propagator. Being a CPTP map, V F is known to have an
eigenvalue 1, and we let the corresponding eigenstate be η. Supposing, for simplicity, that all
the other eigenvalues have absolute values less than 1 (see Ref. [61] for a sufficient condition
for this), we obtain the long-time behavior ρ(t) → V (t 0 , 0)η as t → ∞. This is the NESS
solution periodic in time, and V (t 0 , 0) is called the micromotion within a period. Although
the NESS is thus obtained from the one-cycle and micromotion parts of the propagator, it isdifficult to obtain them analytically and a hard task to do them numerically.

3

van Vleck high-frequency expansion of Liouvillian

The high-frequency (HF) expansions offer systematic ways for looking into the propagatoranalytically when ω is large enough. As developed for isolated systems, there are, at least,
three versions: the Floquet–Magnus [62], van Vleck [8], and Brillouin–Wigner [9]. These
three look different but are related to each other by appropriate transformations [8, 9]. The
Floquet–Magnus approach has been generalized to the FLE in the literature [50, 63]. In thissection, we generalize the van Vleck approach to the FLE, which will be useful for analyzing
the NESS in the following.
Before discussing the FLE, we briefly review the van Vleck HF expansion in isolated sysd
tems. In isolated systems, one tries to solve the Schrödinger equation dt|ψ(t)〉 = −iH(t) |ψ(t)〉5

SciPost Phys. Core 4, 033 (2021)

for a periodic Hamiltonian H(t). The formal solution is given by |ψ(t)〉 = U(t, t 0 ) |ψ(t 0 )〉 withU(t, t 0 ) being the unitary propagator from t 0 to t. The HF expansion is derived from thefollowing decomposed form of U(t, t 0 ) [8, 9].
00U(t, t 0 ) = e−iK(t) e−iHeff (t−t ) eiK(t ) ,(12)where Heff is a time-independent effective Hamiltonian, K(t) = K(t + T ) is a periodic Hermitian operator, and eiK(t) is called the micromotion operator. With the series expansion,P∞P∞ (k)(k)
K(t) = k=1 K (k) (t) and Heff = k=0 Heff , we obtain each of K (k) and Heff order by order,(k)where K (k) and Heff are O(ω−k ). First few terms are given in Appendix A.
Now we generalize the van Vleck HF expansion to the FLE in open systems. To this end, we
invoke the formal analogy between the FLE (1) and the Schrödinger equation. Although there
is a difference between the operator −iH(t) and the superoperator L t , these are both linear
operators. Thus, the derivation of the HF expansion goes in parallel by the formal substitution
H(t) → iL t . More concretely, we put the following ansatz:
0

V (t, t 0 ) = eG t eLeff (t−t ) e−G t 0 ,(13)

where Leff and G t (= G t+T ) are superoperators corresponding to the effective Liouvillian andthe micromotion. Introducing the series expansion

Leff =∞
X(k)
Leff ;

Gt =

∞
X

k=0
(k)

(k)

with Leff = O(ω−k ) and G t

(k)

Gt ,

(14)

k=1

= O(ω−k ), we obtain(0)

iLeff = iL0 ,X [iL , iL ](1)−mm
iLeff =,2mω
m6=0
X [[iL , iL ] , iL ] X X [[iL , iL
(2)
−m
0
m
−m
m−n ] , iLn ]iLeff =
+
,
2
2
2m ω
3mnω2
m6=0
m6=0 n6=0,m

(15)
(16)
(17)

and
(1)

−G t

=−

X iL

m −imωte,mω
X X [iL , iL

(18)

m6=0

(2)

−G t

=

m−n ] −imωt
e
+
2mnω2
n

m6=0 n6=0,m

X [iL , iL ]m

m6=0m2 ω 20e−imωt ,(19)

up to the second order of 1/ω. Here, Lm denotes the Fourier components of the Liouvillian:X
Lt =Lm e−imωt ,(20)
mand the results have been obtained by the replacements
H m → iLm ,(21)

Heff → iLeff ,(22)

K(t) → iG t ,(23)

in the HF expansion for isolated systems (see Appendix A). One can obtain higher-order terms
by continuing the procedure systematically.
6

SciPost Phys. Core 4, 033 (2021)

We remark the relation to the Floquet–Magnus (FM) approach in the literature [50, 63].
In this approach, we take a reference time t = t 0 and define the one-cycle evolution superoperator V F (t 0 ) ≡ V (t 0 + T, t 0 ). The Floquet Liouvillian L F (t 0 ) is defined, in the FM approach,by

V F (t 0 ) = e T LF (t 0 ) .

(24)

On the other hand, in the van Vleck HF expansion [Eq. (13)], this one-cycle propagator isrepresented as V F (t 0 ) = eG t0 e T Leff e−G t0 . By equating these, we havee T LF (t 0 ) = eG t0 e T Leff e−G t0 ,(25)

where we have used the periodicity G t 0 +T = G t 0 . Thus, the effective Liouvillians L F (t 0 ) and
Leff in the different approaches are related to each other by the similarity transformation eG t0 .
The FM expansion is the series expansion for L F (t 0 ) with 1/ω, and we have L F (t 0 ) = Leff = L0
at the zeroth order of 1/ω. At higher orders, the L F (t 0 ) involves more terms than Leff and thedependence on a reference time t 0 similarly to the case of isolated systems [9]. We note that
L F (t 0 ) contains the micromotion information.
Despite the formal analogy in the derivations of the van Vleck HF expansion in open andisolated systems, we need to be careful about the properties of G t and Leff . First, the micromotion superoperator eG t is not unitary in general unlike eiK(t) . Second, Leff may not be of
Lindblad form. As shown in Refs. [64, 65], even if the one-cycle evolution e T LF (t 0 ) is a CPTP
map, its logarithm T L F (t 0 ) may not be a Lindbladian. One representative situation is when
the one-cycle evolution e T LF (t 0 ) has a negative real eigenvalue [64, 65]. Noting that Eq. (25)
implies that e T LF (t 0 ) and e T Leff have the same eigenvalues,2 it also happens that Leff is not a
Lindbladian. In addition, Mizuta et al. have recently shown that the FM expansion for L F (t 0 )may not be a Lindbladian [67]. In the following, we will show that, in the van Vleck approach,
Leff is of Lindblad form at the first order of 1/ω for physically relevant models. Another recent
work [68] has also pointed out that Leff can be of Lindblad form while L F (t 0 ) is not.
Nevertheless, Leff has a good property for analyzing the dynamics and NESSs: Leff has atleast one eigenvalue equal to zero at every order of the HF expansion. To show this, we firstprove the following lemma.
Lemma 1 At each order of the HF expansion, we havetr[Leff (A)] = 0 ,for any operator A.

(26)

Proof. We begin by proving
tr[Lm (A)] = tr[Hm (A) + Dm (A)] = 0 ,∀A and ∀m .

(27)

To show Eq. (27), we first notice Eq. (6) implies tr[Hm (A)] = −i tr([H m , A]) = 0, where wehave used the cyclic property of the trace tr(H m A) = tr(AH m ). We second notice Eq. (8) leads††††
to tr[D(A)] = tr[Lα,m−n ALα,n− 21 {Lα,nLα,m−n , A}] = tr[Lα,nLα,m−n A − Lα,nLα,m−n A] = 0, where
we have again used the cyclic property. These two steps prove Eq. (27).
Now we prove Eq. (26) using Eq. (27) at each order N of the HF expansion. As
PN
PN
(k)
(k)
Leff = k=0 Leff and hence tr[Leff (A)] = k=0 tr[Leff (A)] at an N -th order, it is sufficient
to prove(k)

tr[Leff (A)] = 0 ,∀A and ∀k .

(28)

Similarity transformations do not change eigenvalues [66]. Equation (25) means that e T LF (t 0 ) and e T Leff areconnected by a similarity transformation eG t0 .
2

7

SciPost Phys. Core 4, 033 (2021)(0)

For k = 0, Eq. (28) is obtained from tr[Leff (A)] = tr[L0 (A)] = 0, where the lastequality follows from Eq. (27) for m = 0. Let us hence prove Eq. (28) for k ≥ 1. As(m)represented in Eqs. (16) and (17), Leff consists of nested commutators between Lm ’s fordifferent m’s. Unraveling all the commutators, we have the following formal expression
P(k)
Leff = m1 ,m2 ,...,mk Cm1 ,m2 ,...,mk Lm1 Lm2 . . . Lmk , where Cm1 ,m2 ,...,mk are complex numbers. This
expression givesX(k)tr[Leff (A)] =Cm1 ,m2 ,...,mk tr[(Lm1 Lm2 . . . Lmk )(A)](29)
m1 ,m2 ,...,mk=X

Cm1 ,m2 ,...,mk tr[Lm1 (Am2 ,m3 ,...,mk )] ,(30)

m1 ,m2 ,...,mk

where we have defined operators Am2 ,m3 ,...,mk ≡ (Lm2 . . . Lmk )(A). Here, substituting m = m1
and A = Am2 ,m3 ,...,mk into Eq. (27), we have tr[Lm1 (Am2 ,m3 ,...,mk )] = 0 holds for every set(k)(m1 , m2 , . . . , mk ), meaning that tr[Leff (A)] = 0. Thus, we have proved Eq. (28) and hence
Eq. (26).
From this lemma, the following theorem holds true.
Theorem 1 Leff has at least one eigenvalue equal to zero at every order of the HF expansion .
Proof. Recall that 〈X , Y 〉 ≡ tr(X † Y ) serves as an inner product for two operators X and Y . We
translate Eq. (26) as 〈I, Leff (A)〉 = 0 = 〈L∗eff (I), A〉, where I is the identity operator and L∗eff
is the adjoint superoperator3 for Leff . Remembering that A can be any, we have L∗eff (I) = 0,
which means that L∗eff has a zero eigenvalue (and I is the left eigenvector of Leff ). Therefore,
Leff also has a zero eigenvalue
As we will see below, the zero eigenvalue corresponds to the NESS.
In this work, we assume that Leff ’s eigenvalues all have nonpositive real parts. Since wehave shown the existence of the zero eigenvalue, it means that the maximum of the eigenvalues
real parts is zero. This property is important to obtain sensible time evolution since if Leff had
an eigenvalues with a positive real part (and hence V F had one with absolute value greaterthan 1), the density operator would blow up in many cycles of evolution. In physically relevant
setups discussed in Secs. 5 and 6, we will see that Leff is of Lindblad form at O(ω−1 ) andindeed has the good property. Thus, if the zero eigenvalue is not degenerate (as is the casein all example models in this paper), the higher-order corrections do not break the property
for large enough ω. If it is degenerate, there may appear small positive real parts at higher
orders, and thus one must be careful about this possibility in general. We leave the general
proof of the nonpositivity as an open problem and assume the nonpositivity throughout this
work.

4

General aspects of NESS solution

With the van Vleck HF expansion discussed in the previous section, we can obtain the NESS
solution, including the micromotion, without explicitly calculating time evolution. We discuss
how it generally works in this section and will analyze physically relevant examples in thefollowing sections.
Let us suppose that we have an arbitrary initial state ρ0 at time t = 0 and consider its
long-time evolution. Equations (9) and (13) give us the quantum state at time t as
ρ(t) = eG t e t Leff e−G0 ρ0 .
3

For clarity, we use A† for the adjoint operator for an operator A and L∗ for that for a superoperator L.

8

(31)

SciPost Phys. Core 4, 033 (2021)

In analyzing the long-time behavior, it is convenient to absorb the micromotion by the similarity
transformation:
ρ 0 (t) ≡ e−G t ρ(t) ,(32)

ρ 0 (t) = e t Leff ρ 0 (0) .

(33)

which satisfiesdρ 0 (t)In other words, ρ 0 (t) obeys the time-independent master equation dt = Leff ρ 0 (t) with theinitial condition ρ 0 (t = 0) = e−G0 ρ0 .
To investigate the long-time behavior of ρ 0 (t), it is useful to introduce the eigenstates ofLeff . For simplicity, we assume that Leff is diagonalizable and introduce the eigenstates of Leff
as follows:

Leff (ηλ,a ) = ληλ,a .

(34)Here, ηλ,a is the eigenstate of Leff (ηλ,a is a complex matrix acting on the Hilbert space)belonging to the complex eigenvalue λ, where a (= 1, 2, . . . , Nλ ) labels the degenerate eigenstates with Nλ being the degree of degeneracy. As discussed in the previous section, there
exists λ = 0 eigenvalue, and we assume that Reλ ≤ 0 for all λ.
We note that the trace preserving nature of the effective evolution e t Leff imposes the following condition,tr(ηλ,a ) = 0 ,(for λ 6= 0) .

(35)

Here the trace preserving nature means tr e t Leff A = tr(A) for any matrix A, which follows
from tr[Leff (A)] = 0 at each order of the HF expansion. To prove Eq. (35), one considers ρ 0 (t)starting from an arbitrary initial states ρ 0 (0), which can be represented as a linear combinationof ηλ,a ’s (see Eqs. (37) and (38) below). Since tr[ρ 0 (t)] is time-independent, one obtains thattr(ηλ,a ) must vanish except for the zero modes λ = 0. For the zero modes, tr(η0,a ) can benonvanishing, and, if so, we impose the following normalization
tr(η0,a ) = 1 .

(36)

By using the eigenstates ηλ,a , we solve the asymptotic behavior of ρ 0 (t). To this end, we
express the initial state by these eigenstates asX
ρ 0 (t = 0) = e−G0 ρ0 =cλ,a ηλ,a ,
(37)
λ,a

which leads with Eq. (34) toρ 0 (t) =X

cλ,a eλt ηλ,a .

(38)

λ,a

Since Reλ < 0 holds true for all the λ’s except λ = 0, these eigenmodes all vanish in thelong-time limit t → ∞, and we haveX0
ρ 0 (t) → ρ∞=
c0,a η0,a .
(39)aThus, ρ 0 (t) approaches a time-independent asymptotic state, which can be calculated withlinear algebras for Leff instead of explicit time-evolution simulations.
9SciPost Phys. Core 4, 033 (2021)

0Once we have ρ∞, we immediately obtain the NESS solution for ρ(t),X
0
ρ(t) → ρness (t) = eG t ρ∞=
c0,a eG t η0,a .

(40)

a

The periodicity of G t = G t+T ensures that the ρness is also periodic
ρness (t + T ) = ρness (t) .

(41)

Depending on whether Nλ=0 = 1 or Nλ=0 ≥ 2, Eq. (40) gives different physical consequences.
When Nλ=0 = 1, the sum over a is absent and we have c0,1 = 1 owing to the trace condition0tr(ρ∞) = 1. Thus, Eq. (40) is further simplified asρness (t) = eG t η0,1

(when Nλ=0 = 1) .

(42)

We emphasize that there is no dependence on the initial state ρ0 in Eq. (42). When Nλ=0 = 1,
there is the unique steady state for Leff , and ρ 0 (t) converges to this special state no matter
what initial state we take. On the contrary, when Nλ=0 ≥ 2, there remain some initial state
dependence in c0,a ’s. To obtain the NESS, we need to solve the linear equations [Eq. (37)] forthe unknown coefficients cλ,a , plugging cλ=0,a into Eq. (40). The multiple zero modes happen
in physical models typically when the model has symmetries (see, e.g., Ref. [60]).
For completeness, we discuss how to calculate expectation values of observables in the
NESS assuming Eq. (42). Using Eq. (14) together with Eqs. (18) and (19), we have
(1)
P
[G ]2
(1)
ρness (t) = (1 + G t + t2 + G (2) . . . )η0,1 = m ρm e−imωt , where ρm is time-independentand obtained up to a desired P
order of ω−1 . For an observable A, we have its expectation
value as A(t) = tr[ρness (t)A] = m tr(ρm A)e−imωt . Therefore, once we have ρm ’s from the HF
expansion, we can immediately obtain expectation values at arbitrary times t without time
integration.
The above argument in deriving the NESS solution is exact once Leff and G t are given at
an arbitrary order of the HF expansion. In this approach, the analytical expressions for theapproximate Leff and G t give us physical intuitions for the effects of the drive and dissipation.
Also, the HF expansion avoids dynamics simulations and enables efficient NESS calculations
and possibly analytical calculations in some problems. For instance, one can use the Lanczos
algorithm to numerically obtain η0,1 more efficiently than the direct dynamics simulations.
In the following sections, we classify the problems and develop and demonstrate the NESS
calculations based on the HF expansion.

5

Phenomenological time-independent dissipators

In this section, we focus on the case in which the dissipator is time-independent D t = D. This
class of problems is widely studied in, e.g., quantum-optic [69], Rydberg atoms [70], cavityQED [71], electronic [72,73], and spin [74–77] systems. Bloch equations used in magnets [54]and semiconductors [78,79] can also be viewed as equations of motion with time-independentdissipators. First, we reduce the HF expansion for the Liouvillian in Sec. 3 with the Hamiltonianand dissipator. Then, we apply the reduced formulas to the example model of a dissipative
spin chain [80–82] and discuss what can be known by the HF expansion approach. We notethat the van Vleck HF expansion for a time-independent dissipator was studied in Ref. [57] atthe leading order O(ω−1 ). Here we study general dissipators at higher orders.

10

SciPost Phys. Core 4, 033 (2021)

5.1

High-frequency expansion

When D t is time-independent as D t = D, we have Dm = δm0 D. Therefore, the Liouvillian
Fourier components are simplified as

Lm = Hm + δm0 D .

(43)

We substitute this special form of Lm into the Liouvillian HF expansion derived in Sec. 3,
obtaining more concrete formulas.
Since the derivation is rather straightforward, we write down the results as follows:(0)

Leff (ρ) = H0 (ρ) + D(ρ) ,
(1)

(44)

(1)

iLeff (ρ) = [Heff , ρ] ,
(2)

(2)

iLeff (ρ) = [Heff , ρ] +(45)
X [[iHm6=0−m , iD ] , iHm ] (ρ),
2m2 ω2

(46)

and
(1)

−G t (ρ) = [iK (1) (t), ρ] ,(2)

−G t (ρ) = [iK (2) (t), ρ] +(47)
X [iH , iD] (ρ)
m

m6=0

m2 ω2

e−imωt .

(48)

(k)

Here, Heff and K (k) (t) are the effective Hamiltonian and micromotion defined in the HF expansion for isolated systems (see Appendix A for their explicit forms). In the derivation, we
have used a useful formula for nested commutators presented in Appendix B. Higher-orderresults are similarly obtained by straightforward calculations.
We remark that, up to the first order, the effective Liouvillian is in Lindblad form:

Leff (ρ) = −i[Heff , ρ] + D(ρ) + O(ω−2 ) ,(49)

which follows from Eqs. (44) and (45). This is a good property particular to the van Vleck
approach and does not hold for L F (t 0 ) [Eq. (24)] in the FM representation [67]. However, at
higher orders, Leff is not necessarily of Lindblad form due to the second term of the right-handside in Eq. (46). These kinds of terms derive from the interplay between the external drive H mand dissipation D. As for G t , Eq. (47) dictates that the micromotion is essentially the same asin isolated systems up to the first order, whereas the interplay sets in at the second order, as
one can see in Eq. (48).

5.2

Example 1: Open XY Chain with Boundary Dissipation

Let us apply the HF expansion to an open XY spin chain subject to dissipation acting on the twoedges of the chain. This class of models is mapped to quadratic Majorana fermions [80] andstudied extensively. Here we discuss such a chain under periodic drive f (t) = f (t + T ) [83]:H(t) =

N
−1 
X
j=1


N
X
1+ g x x
1− g y yσ j σ j+1 +σ j σ j+1 + h f (t)
σzj
2
2
j=1

(50)

and the dissipatorD(ρ) =

4 
X
α=1

1
Lα ρ Lα† −
2

11


 †
Lα Lα , ρ ,(51)

SciPost Phys. Core 4, 033 (2021)

withL1,2 =

Ç

L
σ1± ,
Γ1,2

L3,4 =

Ç

R
Γ1,2
σ±
N.

(52)

In the literature,P this model was first studied analytically for the periodic kicks
f (t) = δ T (t) ≡ T m∈Z δ(t − mT ) [83]. They showed the following two properties:(i) This model exhibits a rich phase diagram in the (T, h)-plane,(ii) the boundary dissipations {Lα }4α=1 do not qualitatively change the phase diagram .
Here, the rich phase diagram was characterized by the spin-spin correlations [83], and, itscorrespondence to the topological edge mode was elucidated afterwords [81, 82]. Our aim
here is to reexamine this model by applying our HF expansion method and to provide further
insights. Note that the HF expansion approach can be used for various driving protocols f (t)other than δ T (t) in exchange for restricting ourselves to the HF regime.
We begin by considering the property (i) described above by the HF expansion. As noted
above, this property is related to the bulk Hamiltonian (the bulk-boundary correspondence wasconfirmed for the topological edge mode [81]). Thus, we neglect the dissipator and consider
the infinite chain in discussing the property (i).
The spin Hamiltonian (50) is many-body and hard to solve analytically, and it is convenient
to map this to Majorana fermions, or the two Hermitian components of the complex JordanWigner fermions. We introduce 2N Majorana fermions w j ( j = 1, 2, . . . , 2N ) for an N -site chain
as
!
!
n−1
n−1
Y
Y
(53)
w2n−1 w2n = iσzn ,σzj σny ,w2n =σzj σnx ,w2n−1 =
j=1

j=1

which satisfy the Majorana commutation relations {w i , w j } = 2δi j and translate Eq. (50) into
the following quadratic Hamiltonian
H(t) = −i

N
−1 
X
j=1


N
X
1+ g
1− g
w2 j w2 j+1 −w2 j−1 w2( j+1) − ih f (t)w2 j−1 w2 j .
2
2
j=1

(54)

The Majorana fermions are interpreted as spinon excitations behaving as magnetic domain
walls. From the Majorana fermions to the spins, we can use the inverse transformation of
Eq. (53) given by!
!
2(n−1)2(n−1)YY
σnx = (−i)n−1
w j w2n−1 , σny = (−i)n−1
w j w2n , σzn = −iw2n−1 w2n .
j=1

j=1

(55)
It is convenient to introduce the two-component fermions W j = t (w2 j−1 , w2 j ) and make the
P −ik jfk ∝Fourier transform W j → W
W j in the limit of N → ∞. Then we have the Heisenje
fk (see Appendix C for derivation):berg equation for W
fk (t)dW
fk (t) ,= −ih(k, t)W
dt

(56)

h(k, t) = −2{(g sin k)τ x + [cos k − h f (t)]τ y } ,

(57)

where

12

SciPost Phys. Core 4, 033 (2021)

where τα (α = x, y, and z) denote the Pauli matrices. Thus, the original spin Hamiltonianhas been mapped to noninteracting two-component Majorana fermions. Although the mapping (53) is nonlinear and does not simply tell us every spin observable, the two eigenvaluesof Eq. (57) describe the instantaneous energy dispersion relations for the elementary spinon(Majorana fermion) excitations.
As shown in Refs. [74, 83], the nontrivial phase is related to the quasienergy bands of theelementary Majorana fermions. Since h(k, t) is Hermitian and traceless, the two eigenvaluesRTof the one-cycle unitary V (k) = exp+ [−i 0 H(k, t)d t] are given as exp[±iε(k)T ], and ±ε(k)determine the two quasienergy dispersion relation for the Floquet modes in k ∈ (−π, π]. The
appearance of the nontrivial phase is signaled by the appearance of the nontrivial solutions
for dε(k)/d k = 0 [74, 83]. Here, the nontrivial solutions mean k 6= 0 or π,Æ
since, without theperiodic drive, or f (t)=0, the two quasienergy bands are given by ±ε = ±2 cos2 k + g 2 sin2 k
and dε(k)/d k = 0 at k = 0 and π.
Let us now derive the phase diagram by applying the HF (or small T ) expansion to Eq. (57) for the periodic kicks f (t) = δ T (t). In this example, the FourierRT
dt h(k, t)e imωt /T are simply given bycomponents of the Hamiltonian hm (k) =0
hm=0 = h(k, t) = −2{g sin kτ x + (cos k − h)τ y } and hm6=0 = 2hτ y .
(0)
der of the HF expansion, we have heff (k)
Æ
(0)
2
2
2

At the zeroth or-= hm=0 , which gives the quasienergy bands±ε (k) = ±2 (cos k − h) + g sin k. In this order, ε(0) (k) does not depend on T , andthe number of solutions in dε(0) (k)/d k = 0 changes whether |h| ≥ 1 or |h| < 1. Whereas,for |h| ≥ 1, the only solutions in k ∈ [0, π] are k = 0 and π, there is an extra solution k = k∗ ∈ (0, π) for |h| < 1. Correspondingly, the spin-spin correlations [83] andthe number of topological edge modes become different in these parameter regions asshown in Fig. 1(a). At the first order of the HF expansion, heff does not acquire a cor(1)rection: heff = 0. This follows from the time-reversal symmetry of the driving protocol,(2n+1)f (−t) = f (t), and heff

= 0 (n ∈ N) more generally. At the second order, we obtain a(2)nonvanishing correction from heff , and the effective Hamiltonian up to this order is given byheff = −2{g[1 − 2(hT )2 /3] sin kτ x + (cos k − h)τ y }. We note that the periodic kick in the
Zeeman coupling has reduced the anisotropy |g|. In this order, heff depends on T and gives
more complex phase diagrams as shown in Figs. 1(b) and (c). We remark that the phase diagrams are consistent with the exact solution [83] at the small-T (i.e., high-frequency) region.
We leave it an open question to examine how the phase diagram changes for higher-order
calculations, which require considerable effort.
One advantage of the HF expansion approach is that we can analyze other driving protocols
rather than periodic kicks. While the problem is exactly solvable for the kicks, it is not for, e.g.,
the harmonic drive h f (t) = b + a cos ωt. To obtain the phase diagram for this drive, we
decompose h(k, t) into its Fourier components: hm=0 is the same as before, hm=±1 = hσ y , andhm = 0 otherwise. Thus, at the zeroth order, we obtain the same phase diagram as before(see Fig. 1(a)). Again, all the odd-order corrections vanish in the harmonic drive as well, andthe first nontrivial correction comes at the second order. The effective Hamiltonian up to thesecond order is given by heff = −2{g[1 − (bT )2 /π2 ] sin kτ x + (cos k − a)τ y }, in which we finda quantitative difference in the change of anisotropy g.
Now, we discuss the other property (ii). As we have seen in Sec. 5.1, the effective Liouvillian involves the contributions from the interplay between the drive and dissipation (see thesecond terms on the RHS of Eqs. (46) and (48)). These terms, in general, act on sites otherthan the two edges even if the bare dissipator D acts only on the edges. As we go higher orders, the dissipators propagate into the bulk and can cause a nonnegligible effect on the bulk
property in general.
Nevertheless, we can show that the present model is so special that all drive-dissipation13

SciPost Phys. Core 4, 033 (2021)

(a)

(b)

(c)

Figure 1: Number of solutions for dε(k)/d k = 0 in k ∈ (0, π) obtained by the HF
expansion. The blue (yellow) region represents the number 0 (1). The zeroth-order(a) and second-order (b) calculation for g = 0.1. (c) The second-order calculation
for g = 0.9.
interplay terms vanish andLeff ρ = −i[Heff , ρ] + D(ρ) ,(58)

holds at each order of the HF expansion. Namely, the drive does not let the dissipator propagate
into the bulk, and the effective dynamics is described by the effective Hamiltonian plus the bareboundary dissipation.
To prove Eq. (58), P
it is enough to show [Hm , D] ∝ [Sz , D] = 0, where Sz denotestheP
superoperator Sz ρ ≡ [ j σzj , ρ]. Note that Hm ∝ Sz for m 6= 0 because H m ∝ j σzj (seeEq. (50)). One can easily show this by straightforward calculations for the superoperator
commutators, but this approach cannot gain physical insights. Here we provide another proof
based on the vectorization of the operator space [66]. In this technique, the density matrix ρ
acting on the Hilbert space H is regarded as a vector ρ~ in H ⊗ H. Correspondingly, a superoperator ρ → LρR with L and R being matrices acting on H is regarded as a matrix L ⊗ R T (Tdenotes the transpose) acting on H ⊗ H. Graphically, each superoperator becomes an operator
acting on two copies of spin chains, as illustrated in Fig. 2. In our model, the dissipator isrepresented as

D = D1 + D2 ,
4
X
+
R −
−
D1 ≡
Lα ⊗ (Lα† ) T = Γ1L σ1+ ⊗ σ1+ + Γ2L σ1− ⊗ σ1− + Γ1Rσ+
N ⊗ σ N + Γ2 σ N ⊗ σ N ,

(59)
(60)

α=1

4

D2 ≡ −


1 X †1
Lα Lα ⊗ 1 + 1 ⊗ (Lα† Lα ) T = − (K ⊗ 1 + 1 ⊗ K) ,
2 α=1
2

(61)

+
R + −
K ≡ Γ1L σ1− σ1+ + Γ2L σ1+ σ1− + Γ1Rσ−
N σ N + Γ2 σ N σ N=

Γ2L − Γ1L
2

σ1z +Γ2R − Γ1R
2

σzN +

Γ1L + Γ2L + Γ1R + Γ2R
2

,

(62)

where we have used σ+j σ−j = (1 + σzj )/2. Although the explicit form involves many terms,
its physical meaning is obvious if we think of D as a virtual Hamiltonian acting on the two
spin chains. First, D1 consists of pairwise spin raising or lowering on the left or right edges, as
shown in Fig. 2(a). Second, D2 gives local Zeeman energy on each corner of the two chains,
where K ⊗1 (1⊗K) acts on the upper (lower) chain as illustrated in Fig. 2(b). Thus, D does not
14

SciPost Phys. Core 4, 033 (2021)

(a)

j=1

(b)

2

N-1

j=1

N

2

N-1

N

Figure 2: Schematic illustration of matrices (a) D1 and (b) D2 in the vectorization
scheme of ρ → ρ~ ∈ H ⊗ H. The upper (lower) chain corresponds to the left (right)
part of H ⊗ H. The colored sites show nontrivial actions of parts of D1 and D2 (seeEqs. (60)–(62)).

conserve the total spin along the z direction but conserves the difference between the in-chain
total spins along z. Now, we recall that Sz is given in the matrix representation as
X
X
(63)
Sz =
σzj ⊗ 1 − 1 ⊗
σzj ,
j

j

which is the in-chain spin difference. Therefore, Sz is a conserved quantity under the “Hamiltonian” D, and we obtain [Sz , D] = 0 that means [Hm , D] = 0 (∀m) leading to Eq. (58).
We remark that the disentangling of the drive and dissipation is particular to the present
dissipator (52), which is absent in generic models (another special case is the dephasing
Lα ∝ σzj ). For example, we can think of other dissipators such as spin raising and lowering
along x instead. In this case, D1 involves the terms like σ1+ ⊗ σ1− , and the difference between
the in-chain total σz is no longer conserved, implying that the drive-dissipation-interplay could
influence the bulk property more.
Let us summarize this subsection. We have studied the prominent dissipative open spin
chain, addressing the issues (i) and (ii) from the HF expansion viewpoint. As for (i), we haveprovided a way to analyze generic driving protocol f (t) at high frequency. As for (ii), we
have found the drive-dissipation disentangling, which means that the effective Liouvillian isof Lindblad form, and the dissipator is confined in the spin chain edges.

5.3

Example 2: Three-level system

Here we demonstrate, more quantitatively, the NESS calculation by the HF expansion. To
this end, we take an effective Hamiltonian for the NV center in diamonds under a circularlypolarized ac magnetic field [84]:circ
HNV (t) = H0 + Hext(t) ,H0 = −Bs Sz + Nz Sz2 + Nx y (S 2x − S 2y ) ,

(64)
circ
Hext
(t) ≡ −Bd (S x cos ωt + S y sin ωt) ,(65)

where S x, y,z are the spin-1 operators, Bs is the static Zeeman field, Nz and Nx y represent thecircmagnetic anisotropy terms, and Hext(t) denotes the coupling to the circularly-polarized acmagnetic field. In our analysis, we set Nz = 1 and Nx y = 0.05 because of Nz  Nx y in NV
centers [84]. Since we consider the spin S = 1, the spin matrices are 3 × 3. This model was
analyzed in Ref. [57] at O(ω−1 ), and we here extend the analysis to O(ω−2 ). We write theH0 ’s eigenstates as H0 |Ei 〉 = Ei |Ei 〉 (i = 1, 2, and 3).
As in Ref. [57], we here consider time-independent dissipation and will analyze timedependent dissipation to compare the two cases in Sec. 6.2. To make the comparison easy,

15

SciPost Phys. Core 4, 033 (2021)

we use the dissipator that represents the coupling of S x to an ohmic bosonic thermal reservoir [40]:

X
©
1¦ †
D(ρ) =
γ(Ei j ) A Ei j ρA†E −A E A Ei j , ρ ,
(66)
ij
ij
2
i, j
where Ei j ≡ Ei − E j , A Ei j = Π Ei S x Π E j (Π Ei ≡ |Ei 〉 〈Ei |), and2− εεe 2Λ2γ(ε) = γ0,1 − e−βε(67)

is the bath spectral function. Here, γ0 gives the system-bath coupling strength, β the inverse temperature, and Λ the spectral cutoff. In the following numerical calculations, we set
γ0 = 0.2, β = 3, and Λ = 10.
We remark that the dissipator (66) is approximate because it is represented by the energy eigenvalues and eigenstates for the undriven H0 . Later in Sec. 6.2, we will analyze amore precise dissipator that is derived microscopically and involves the driving effect, findingthat the dissipator difference only gives small quantitative corrections except some symmetry. Given that similar approximate time-independent dissipators are widely used in NV-centerstudies [84], the comparison between Secs. 5.3 and 6.2 offers a theoretical validation for theapproximation.
The dissipator (66) satisfies the detailed balance condition. To see this, we rewrite Eq. (66)
as
X 
©
1¦ †
†
D(ρ) =
Γi j L i j ρ L i j −L Li j , ρ ,
(68)
2 ij
i, j
where L i j = |Ei 〉 〈E j | are jump operators between the energy eigenstates, andΓi j = γ(Ei j )| 〈Ei |S x |E j 〉 |2 are the corresponding transition rates. Noting that γ(−ε) = e−βε γ(ε),
we have Γi j e−β E j = Γ ji e−β Ei . This detailed balance condition ensures that, without drive
circ
Hext
(t), the stationary solution of the Lindblad equation corresponds to the canonical ensemble ρcan = e−β H0 /Z with Z = tr(e−β H0 ).
In this setup, we calculate the NESS for the Lindblad equation in two ways. Thedρ
first way is the numerical time-integration of the Lindblad equation d t = L t ρ for
L t ρ = −i[HNV (t), ρ] + D(ρ). By solving for independent initial conditions, we numericallyobtain all the matrix elements for the one-cycle superoperator (9 × 9 matrix) V F . Then we
diagonalize V F , finding η as the eigenvector with eigenvalue 1. Finally, the NESS ρness (t) is
obtained by the time-integration starting from η. For clarity, we write ρness (t) thus obtainedexact
as ρness(t). This way of calculation is exact but demands many time integrations.
The second way of calculation is the HF expansion approach. Following Eqs. (44)–(46),
we write out Leff at an order N (= 0, 1, 2, . . . ). Then we numerically find η as the eigenvector(N )of Leff with zero eigenvalue. Then we obtain ρness(t) = eG t η, where G t is calculated from
Eqs. (47) and (48), depending on N . In this approach, we do not need any time-integration(N =1)
but only use linear algebra. At the first order [57], ρness(t) was exactly solved for any γ0and shown to coincide with the canonical Floquet steady state in the limit of γ0 → 0. We notethat the correction terms to Heff in the present model are given by(1)

Heff = −

2Bd2
ω

Sz ,

(2)

Heff = −2Bd2 
ω2


−Bs Sz + 3Nz Sz2 + Nx y (S 2x − S 2y ) .

16

(69)

SciPost Phys. Core 4, 033 (2021)

Whereas the first term represents an effective Zeeman field, the second modifies the nematic
terms Nz and Nx y . The corrections to the micromotion read
2Bd[sin(ωt)S x − cos(ωt)S y ] ,ω
2Bd Bs[sin(ωt)S x − cos(ωt)S y ]K (2) (t) =ω2
2Bd (Nx y + Nz )2Bd (Nx y − Nz )+cos(ωt){S y , Sz } −cos(ωt){Sz , S x } ,2ωω2
K (1) (t) =

(70)

(71)

which consist of circularly polarized ac magnetic field and similar terms for nematics. The
drive-dissipation-interplay terms are not written in simple form and hence treated just numerically.
The one-cycle evolutions for the NESS calculated with these two ways are plotted in Fig. 3,

×10 4(a) Sx

×10 4
5
0
5

5
0
5
0.00 0.25 0.50 0.75 1.00
t /T
8.96808.96758.96708.9665(b) Sy×10 2(c) Sz0.00 0.25 0.50 0.75 1.00
t /T
×10 1

0.00 0.25 0.50 0.75 1.00
t /T
1.485
1.490
1.495
1.500

(e) Sz2×10 2

0.00 0.25 0.50 0.75 1.00
t /T
×10 5

1.25925
1.25900
1.25875(d) Sx2 Sy2(f) {Sx, Sy}

50
50.00 0.25 0.50 0.75 1.00
t /T
×10 2

0.00 0.25 0.50 0.75 1.00
t /T(g) {Sy, Sz}

×10 2(h) {Sz, Sx}

1
0
1

1
0
1
0.00 0.25 0.50 0.75 1.00
t /T

0.00 0.25 0.50 0.75 1.00
t /T

Figure 3: NESS in one cycle calculated by the exact numerical integration (blue) andby the HF expansion approach at order 0 (orange), 1 (green), and 2 (red). All the
independent 8 observables are plotted in each panel. The parameters are chosen as
Bs = 0.3, Nz = 1, Nx y = 0.05, Bd = 0.1, ω = 10, β = 3, γ0 = 0.2, and Λ = 10.
17

SciPost Phys. Core 4, 033 (2021)

where we have set Bs = 0.3, Bd = 0.1, ω = 10, β = 3, γ0 = 0.2, and Λ = 10. For the HF
approach, we plot the results for the zeroth, first, and second orders. At the zeroth order,
Heff = H0 and G (t) = 0, and thus the NESS corresponds to the thermal equilibrium showing
no time dependence. As we increase the order, the HF result tends to approach the exact
numerical integration, which we verify quantitatively below. We note that there are someperiod-T /2 (or the second harmonic) oscillations in, e.g., S 2x − S 2y in Fig. 3. These oscillations
cannot be taken up to the first order since the micromotion K (1) (t) contains only frequency ω.
(N =2)At the second order, ρness(t) = eG t η = (1 + G t + 21 G t2 )η, where 21 G t2 η involves the frequency
2ω.
The HF expansion approach becomes more accurate for higher frequency ω. Figure 4 showthe ω-dependence of one-cycle averages of observables
Z T
dt
O=
tr[ρness (t)O] ,
(72)
0 T

1
((S + )2 − (S − )2 ) in the NESS (S ± = S x ± iS y ). We
for O = Sz , S 2x − S 2y , Sz2 , and {S x , S y } = 2inote that the other observables have zero one-cycle averages for symmetry reasons, as shown
in Ref. [57]. We observe that the higher-order HF expansion reproduces the exact results. To
quantify the goodness of the HF expansion approach more strictly, we introduce the following
measureZ T
1/2
d t (N )
exact
δρN =
,
(73)
(t)k2
kρness (t) − ρnessT
0

where k · · · k denotes the Hilbert-Schmidt norm for matrices. If δρN = 0, it follows that
(N )
exact
ρness
(t) = ρness(t) for any t. Figure 5 shows δρ for the orders 0, 1, and 2 plotted against
ω. We observe1δρN ∝ N +1 ,(74)
ω
for high frequencies ω ¦ 2. This is a clear indication that the N -th order approximation
(N )
exact
ρness
(t) completely describe ρness(t) up to O(ω−N ).
Let us summarize this subsection. We have used the HF expansion approach to calculate
the NESS quantitatively. This approach does not require numerical time integration but uses
linear algebra. In this approach, the effective Hamiltonian and micromotion are interpretable.
Although the drive-dissipation-interplay terms are difficult to interpret, they are important toreproduce the exact results quantitatively.As discussed at the end of Sec. 4, this approach isnumerically efficient especially for many-body systems. To obtain the NESS, the time integration requires many matrix-vector multiplications (i.e., applications of L t onto ρ(t)) before theNESS is reached. The number of multiplications is estimated by ∼ 1/(γ0 ∆t), where γ0 is thedissipation strength and ∆t is the time stepping. To obtain high-accuracy results, we haveto decrease ∆t and need numerous multiplications. On the other hand, in the HF-expansion
approach, we are to find the eigenstate η0,1 of Leff with the largest real part. For this purpose, we can use the famous Lanczos algorithm, where the required number of matrix-vector
multiplications are greatly suppressed. So we can reach the NESS more efficiently with the HFexpansion approach within the HF approximation. The efficiency also applies to calculating
physical observables as discussed at the end of Sec. 4.

6

Microscopically-derived dissipators by weak thermal contact

In Sec. 5, we have discussed phenomenological Floquet-Lindblad equations, where the dissipators are time-independent. However, those dissipators are only approximate, strictly speaking.
18

SciPost Phys. Core 4, 033 (2021)

9.00 ×10(a) Sz2

1.45 ×10
1.47

8.98

1.48

8.97

1.49
101

100

1.30 ×10

1.50

102

(c) Sz2

1

0.75

1.28

0.50

1.27

0.25

1.26

0.00
101

100

0.25

102

101

100

1.00 ×10

1.29

1.25

(b) Sx2 Sy2

1.46

8.99

8.96

2

100

5

102

(d) {Sx, Sy}

101

102

Figure 4: One-cycle averages of observables in the NESS [Eq. (72)] plotted against
ω. Here ρness (t) is calculated by the exact numerical integration (blue) and by the
HF expansion approach at order 0 (orange), 1 (green), and 2 (red). The parameters
are the same as in Fig. 3 except ω.

100
10 2
order=0
order=1
order=2
slope=-1
slope=-2
slope=-3

10 4
10 6
10 8

101

100

102

exact
(N )
Figure 5: Difference between NESSs ρness(t) and ρness(t) [Eq. (73)] plotted against
ω. The data points are the results for N = 0 (blue), 1 (orange), and 2 (green), andthe dashed lines are guides to the eye for some slopes (see legend).

19

SciPost Phys. Core 4, 033 (2021)

The dissipator originates from the fact that the system of interest is coupled weakly to its environment [40], and thus it can be modified and time-dependent when the system is driven
periodically.
In this section, we discuss time-dependent dissipators that are microscopically derived
for the system-bath coupling setup. The most well-known example is the quantum masterequation obtained by the rotating wave approximation (RWA) [41, 85]. Also, there are othertime-dependent Lindblad equations that have recently been derived without the RWA [86–88].
In these FLEs with or without the RWA, one can apply the HF expansion for the Liouvillian
derived in Sec. 3 to the expansions for the Hamiltonian, Lamb shift, and dissipator. Unlike
time-independent ones, time-dependent dissipators have multiple Fourier components, which
induce more terms in the HF expansion. One subtlety is that the Lamb shift and dissipator cannot be Taylor-expanded for ω−1 in general because the bath spectral function γ(ε) cannot be
Taylor-expanded (see, e.g., Eq. (67)) but depend on ω in it if the dissipator is affected by thedrive. Nevertheless, one can still use the HF expansion if we treat these ω-dependence rigorously while HF-expand other ω−1 dependences (see a related discussion below in Sec. 6.1.3).
Throughout this section, we focus on a widely-accepted special class of FLEs: the quantum master equation obtained by the RWA [41, 85]. As we will see below, these FLEs havethe special property that the time-dependence can be eliminated by an appropriate unitarytransformation. Thus, we will derive a different HF-expansion-based approach to analyze theNESS utilizing this useful special property. This new approach is different from the previous
one derived in Secs. 3 and 4, but one could also use the previous one to the examples in thissection.

6.1

Floquet-Lindblad equation obtained by RWA

6.1.1

Floquet-Lindblad equation and its characteristic properties

Leaving its derivation in Appendix D, we here summarize the FLE obtained microscopically with the RWA. Suppose that we have a system-bath coupled Hamiltonian,
Htot (t) = HS (t) + H B + HSB , where HS (t + T ) = HS (t)P is for the periodically driven system ofinterest, H B for the heat bath (reservoir), and HSB = α Aα ⊗ Bα for the system-bath coupling.
The bath is characterized by the correlator,Z∞1Γαβ (ω) =ds eiωs 〈Bα (s)Bβ (0)〉 = γαβ (ω) + iSαβ (ω) ,(75)
20
where Bα (t) = e iH B t Bα e−iH B t , and γαβ (ε) and Sαβ (ε) are Hermitian matrices. We assume thatthe bath is in thermal equilibrium ρB ∝ e−β H B at inverse temperature β, which implies the
Kubo-Martin-Schwinger (KMS) condition
γαβ (−ω) = e−βω γβα (ω) .

(76)

The system-bath coupling acts on the system through Aα whose natural basis is the Floquet
eigenstates:id|ψm (t)〉 = HS (t) |ψm (t)〉 ;dt|ψm (t)〉 = e−iεm t |um (t)〉 ;|um (t + T )〉 = |um (t)〉 ,(77)

where εm is called the quasienergy.
Under appropriate Born and Markov approximations [41, 85] (see also Appendix D), we
obtain the Floquet-Lindblad equation with the Lindbladian

L t (ρ) = −i[HS (t) + ΛLS(t), ρ] + D t (ρ) ,
20

(78)

SciPost Phys. Core 4, 033 (2021)

with the Lamb shift ΛLS(t) and the dissipator D t given byX
βΛLS(t) =Sαβ (ε)Aα†
ε (t)Aε (t) ,(79)

1Aβε (t)ρAα†ε (t) −

(80)

α,β,ε

D t (ρ) =X

γαβ (ε)



α,β,ε

2


 α†
β
Aε (t)Aε (t), ρ .

Here, the jump operators are given byAαε (t) = e−iεtXAαmn (ε) |ψm (t)〉 〈ψn (t)| ,

(81)

m,n

where the matrix elements Aαmn (ε) are defined by the Fourier expansionX〈ψm (t)|Aα |ψn (t)〉 =Aαmn (ε)e−iεt .

(82)

ε

and periodic in time as Aαε (t) are periodic. Since 〈ψm (t)|Aα |ψn (t)〉 = e i(εm −εn )t 〈um (t)|Aα |un (t)〉and |un (t + T )〉 = |u(t)〉, the sums over ε in the above equations are taken forε = εnm;k ≡ εn − εm + kω(k ∈ Z) .

Thus, one can also rewrite Eq. (82) asX
X
Aαmn;k e−i(εn −εm +kω)t .
〈ψm (t)|Aα |ψn (t)〉 =Aαmn (εn − εm + kω)e−i(εn −εm +kω)t ≡

(83)

(84)

k

k

As shown in Sec. D, the jump operator Aαε (t) (Aα†ε (t)) lowers (raises) quasienergy by ε.
One remarkable property of this FLE is the existence of a reference frame in which the
Lindbladian is time-independent. This frame is the interaction picture as is evident in thederivation of the FLE (see Appendix D). Also, as
Pshown in Appendix D, the Lamb shift can be
written, without loss of generality, as ΛLS(t) = n λn |ψn (t)〉 〈ψn (t)|, where λn are the (real)Pβeigenvalues of the Hermitian matrix α,β,ε,l Sαβ (ε)Aα∗(ε)Aln (ε). Using this expression formlΛLS(t) and writing ρ(t) asX
ρ(t) =|ψm (t)〉 σmn (t) 〈ψn (t)| ,
(85)
m,n

we have the Lindblad equation for the matrix σ(t) (in the interaction picture):


X
dσ(t)
1  α†βα†β= −i[λ, σ] +γαβ (ε) A (ε)σ(t)A (ε) −A (ε)A (ε), σ(t) ,
dt
2
α,β,ε

(86)

where (λ)mn ≡ δmn λm . Thus, written for σ(t), the Lindbladian is time-independent. This is a
special property of the FLE of the RWA and not shared with other FLEs in general.
In the rest of this subsection 6.1.1, we assume that the quasienergies are not degenerate
and analyze the NESS using the FLE in the time-independent frame (86). Then, the diagonal
elements, Pn (t) ≡ σnn (t), form a closed set of differential equations (see Appendix D.3 for thederivation)d Pn (t) X=[Wnm Pm (t) − Wmn Pn (t)] ,
dt
m

21

(87)

SciPost Phys. Core 4, 033 (2021)

bath spectrum

k-photon process(k = 1)E2E1

ϵ2 + ω
ϵ1 + ω

∼Λϵ2ϵ1
ϵ2 − ω
ϵ1 − ω

Figure 6: Schematic illustration of the k-photon process for k = 1 in a periodicallydriven two-level system. The curved arrow indicates a dissipation-induced transition between the Floquet states m = 1 and n = 2 entailing the quasienergy change
ε2 + ω − ε1 . Upon this process, the bath’s energy changes by ε1 − ε2 − ω, which
compensates the system’s quasienergy change.

with
Wnm ≡

X
α,β,ε

=

X
α,β,k

X

γαβ (ε)[Aβ (ε)]nm [Aα† (ε)]mn =βγαβ (ε)Aα∗nm (ε)Anm (ε)

(88)

α,β,ε
β

γαβ (εm − εn − kω)Aα∗
nm;k Anm;k ,(89)

where we have used Eq. (84) to obtain the last line. As discussed below, the off-diagonal
elements, forming another closed set of equations, do not contribute to the NESS in mostcases, and we ignore them at this moment (see Appendix D.3 for detail). Since γ(ε) is apositive Hermitian matrix, Wnm is nonnegative for any m and n and hence can be regarded as
the transition rates from the Floquet state |ψm (t)〉 to another |ψn (t)〉. Equation (87) dictatesthat the Floquet-state population Pn (t) obeys the classical master equation with the transition
rates Wmn .
Here we make a physical interpretation of the transition rate Wnm . Equation (89) dictates that the transition rate Wnm consists of contributions labeled by integer k. For eachβk, in addition to the transition amplitude Anm;k , there appears the bath’s spectral weightγαβ (εm − εn − kω). For k = 0, the argument of the spectral weight is εm − εn corresponding to the quasienergy difference between the initial |ψm (t)〉 and final |ψn (t)〉 Floquet statesin the transition. For k 6= 0, the argument involves an extra energy kω, which corresponds
to the energy of k photons. From these observations, we interpret that k in Eq. (89) stands
for the number of photons involved in the Floquet-state transitions, and the quasienergy difference and the photon energy are exchanged between the system and bath (see Fig. 6 forillustration). Thus, in the following, we will refer to each k contribution in the sum of Eq. (89)as the k-photon process.
To study the NESS, we assume that Wmn is irreducible so that the classical master equation (87) has the unique steady-state solution Pnss , which Pn (t) approaches as t → ∞. This
22

SciPost Phys. Core 4, 033 (2021)

solution is characterized by Eq. (87) asX
Wnm Pmss − Wmn Pnss = 0 ,

(90)

m

which is equivalent toX

Snm Pmss = 0 ,Snm ≡ Wnm − δmn

m

X

Wmn .

(91)

m

Equation (91) means that Pmss is the zero-eigenvalue eigenvector for the matrix Smn . Let us also
assume that the NESS density matrix is unique,4 in which case all the off-diagonal elements
σm6=n are guaranteed to vanish in the NESS (see Appendix D.3 for the proof). Once having
Pnss , we obtain the NESS density matrix asXρness (t) =Pnss |ψn (t)〉 〈ψn (t)| .
(92)n

It is noteworthy that Pnss and hence ρness (t) are independent of the system-bath-couplingstrength. Suppose that we multiply a common real number c onto Aα (or Bα ). This results inthe change Wnm → c 2 Wnm . However, this change does not affect Eq. (91), meaning that Pnss
is independent of c. This independence is a special property of the FLE obtained by the RWA
and does not hold for general FLEs such as time-independent ones discussed in Sec. 5.
6.1.2

Conditions for Floquet-Gibbs state (FGS)

Before analyzing concrete examples, we generally study whether the NESS coincides with theso-called Floquet-Gibbs state (FGS)ρFG (t) =X

PnFG |ψn (t)〉 〈ψn (t)| ;PnFG ≡n

1 X −βεn
;
e
Z n

Z≡

X

e−βεn ,

(93)n

for high-frequency drives. The difference between the exact NESS (92) and FGS (93) isonly the Floquet-state population: While Pnss are determined microscopically by Wmn as in
Eq. (91), PnFG are the simple canonical distribution for the quasienergies. As discussed above,
the quasienergies εm are defined modulo ω, and the FGS is ill-defined in general. However, when ω is much larger than the system’s energy scale, we have a natural choice in that−ω/2 < εn < ω/2 and εn approaches, in the limit of ω → ∞, each eigenenergy of the timeaveraged Hamiltonian. When we discuss the FGS in this paper, we implicitly assume thosesituations.
The FGS well approximates the NESS if the driving frequency ω is not only larger than
the system’s energy scale but also than the high-frequency cutoff Λ of the bath spectral function. For simplicity, we here discuss the case in that the quasienergies εm are not degenerate and make use of the classical master equation (87). However, we can generalize thearguments to the case of degenerate quasienergies (see Appendix D.4 for detail). When
ω is much larger than the bath spectral cutoff Λ, for, e.g., the ohmic bath (67), we have
(ε −ε +kω)222γ(εm − εn + kω) ∝ exp − m 2Λn 2≈ exp − k2Λω2 ≈ 0 (similar arguments hold for any
other spectral functions). Thus, Eq. (89) can be approximated as
X
β
Wnm ≈
γαβ (εm − εn )Aα∗
(94)
nm;0 Anm;0 .
α,β4This assumption corresponds to the one that ensures the uniqueness of the NESS eigenvector η mentioned inSec. 4.

23

SciPost Phys. Core 4, 033 (2021)

Within this approximation, we obtain the detailed balance condition Wnm = eβ(εm −εn ) Wmn
with the help of the KMS condition (76). The detailed balance condition means that the
FGS population PnFG = e−βεn /Z satisfies Eq. (91), meaning that the FGS (93) is the NESS
ρness (t) ≈ ρFG (t).
Physically speaking, all k-photon processes except k = 0 are negligible in this limiting caseof ω  Λ. As discussed in Sec. 6.1.1, for a k-photon process to occur, the accompanied photon
energy kω has to be compensated by the bath. However, if the bath spectral cutoff Λ is muchsmaller than the photon energy ω, this compensation is impossible.
However, in more realistic situations with ω ® Λ, k-photon (k 6= 0) processes can berelevant. As one can check easily, if we do not ignore k 6= 0 terms in Wmn (89), the detailedbalance condition is not satisfied, and the FGS is not necessarily the NESS except for somespecial cases [89].
6.1.3

Implementation of high-frequency expansion

Let us now discuss how we use the HF expansion approach to obtain the NESS (92). As thestandard HF expansion applies to |ψn (t)〉, we consider how to calculate Pnss . Since Pnss are
given as the zero-eigenvalue eigenvector for Smn as in Eq. (91), it is sufficient to have Wmn by
the HF expansion.
First, we consider the case of ω  Λ. In this case, as we have discussed in Sec. 6.1.2,
we can ignore k 6= 0 terms in Eq. (89), having Eq. (94). Once we insert the standard ω−1expansion for the quasienergies εm and Aαnm , we obtain the HF expansion for Wnm and hence
Pnss . As shown in Sec. 6.1.2, the first-order approximation in this method leads to the FGS as
the NESS, and higher-order calculations give more accurate NESSs.
Second, we consider the other case ω ® Λ. In contrast to the first case, in Eq. (89), k’s
with |k|ω ® Λ give nonnegligible contributions to Wmn that cannot be expanded for ω−1 ingeneral. Thus, we have
Wnm ≈

kc
X X

β
γαβ (εm − εn + kω)Aα∗nm (εm − εn + kω)Anm (εm − εn + kω) ,(95)

α,β k=−kcwhere kc > 0 is a cutoff for the photon number k and depends on the accuracy which we
require for the results. Within the approximation (95), we can use the HF expansion for thequasienergies εm and Aαnm , obtaining a good approximation for Wnm and hence Pnss .
In the following sections, we demonstrate that this method quantitatively works in concretemodels and show that the FGS may not be the NESS in some cases.

6.2

Example 3: Three-level system revisited and deviation from FGS

Now we study concrete example models and their NESSs for time-dependent dissipators. Asa first example, we revisit the three-level system studied in Sec. 5.3, making the dissipatortime-dependent within the RWA. For simplicity, we neglect the Lamb shift. In the absence ofnematic terms Nz and Nx y , this model has recently been analyzed for a general spin S, andthe environment-controlled Floquet-state paramagnetism has been proposed [90].
Following Eq. (80), we write down the time-dependent dissipator. For concrete calculations, we use |um (t)〉 instead of |ψm (t)〉. In doing so, we choose εm so that εm satisfy
Em − Ω/2 ≤ εm < Em + Ω/2 for m = 1, 2, and 3. Note that this choice is possible only forhigh-enough frequencies that we consider here. We also define the set of possible values forquasienergy differences asSδε = {E = εn − εm | 1 ≤ m, n ≤ 3} .
24

(96)

SciPost Phys. Core 4, 033 (2021)

Note that this set has 7 elements since εm are not degenerate in our model. Then the dissipatoris given by

D t (ρ) =X∞Xγ(E + kω)E ∈Sδε k=−∞

1
AE ,k (t)ρA†E ,k (t) −
2

¦
©
†
AE ,k (t)AE ,k (t), ρ ,

(97)

where
X

AE ,k (t) = e ikωtAmn;k |um (t)〉 〈un (t)| ,

(98)

m,n
(εn −εm =E )Amn;k =Z T
0

dt ikωte〈um (t)|S x |un (t)〉 .
T(99)As remarked in Sec. 6.1, we can utilize, for this class of problems, a convenient frame in which
the Lindbladian is time-independent. Thus, we make use of it and consider Eq. (88) that leads
to
Wmn =

∞
X

γ(εn − εm + kω)|Amn;k |2 ,

(100)

k=−∞

which gives the Floquet-state population Pnss by Eq. (92).
Our aim here is to show how well the HF-expansion approach formulated in Sec. 6.1.3
gives a systematic approximation for ρness (t). Let us suppose that we try to obtain Pnss up
to O(ω−N ) for an N (≥ 0). For this, we invoke the HF expansion for isolated systems tohave the Floquet states |un (t)〉 and their quasienergies εm up to this order, which give Amn;kup to the desired order O(ω−N ). More concretely, we calculate Heff and K(t) up to O(ω−N )by the van Vleck expansion (12), which give |un (t)〉 = e−iK(t) |Eneff 〉 and εn = Eneff , where
Heff |Eneff 〉 = Eneff |Eneff 〉. Although these calculations provide appropriate approximations foreach term in the summation of Eq. (100), we need to take a large-enough cutoff kc for
Wmn ≈

kc
X

γ(εn − εm + kω)|Amn;k |2 ,

(101)

k=−kc

so that the truncation error can be negligible in our O(ω−N ) calculation. For the present
model, the appropriate choice is kc = b(N + 1)/2c (b· · · c denotes the floor function) becauseof the following reason. As one can check easily in our model, |Amn;k |2 = O(ω−2|k| ), whereas
γ(εn − εm + kω) is at most O(ω1 ) for k 6= 0, and their product is of O(ω−2|k|+1 ). Thus, once
we choose kc = b(N + 1)/2c, the truncation error (i.e., the difference between Eqs. (100) and(101)) is negligible in our O(ω−N ) calculation. Once we have Wmn , we obtain the steady-state
population Pnss from Eq. (91) and ρness (t) by Eq. (92).
The one-cycle evolutions for the NESS calculated in this way and in the exact numericaltime integration are plotted in Fig. 7. As in Sec. 5.3, we use the bath parameters as β = 3, andΛ = 10 (the NESS does not depend on γ0 in the present case). For the HF approach, we plot
the results for the zeroth, first, and second orders, which are obtained by expanding Heff andK(t) up to O(ω−N ) for N = 0, 1, and 2, respectively. It is worth noting the similarity between
Figs. 7 and 3. This means that the approximate bath model in Sec. 5.3 is quantitatively goodfor high-frequency drives.
Similarly to the time-independent dissipator studied in Sec. 5.3, at the zeroth order, theNESS corresponds to the thermal equilibrium showing no time dependence. As we increase
the order, the HF result tends to approach the exact numerical integration, which we verify quantitatively below. We note that there are some period-T /2 (or the second harmonic)
25

SciPost Phys. Core 4, 033 (2021)

×10 4(a) Sx

×10 4
5
0
5

5
0
5
0.00 0.25 0.50 0.75 1.00
t /T
×10 2

0.00 0.25 0.50 0.75 1.00
t /T(c) Sz×10 2(d) Sx2 Sy2

1.490
1.495
1.500

8.98
8.97
0.00 0.25 0.50 0.75 1.00
t /T
1.261

(b) Sy×10 1

0.00 0.25 0.50 0.75 1.00
t /T(e) Sz2

×10 5(f) {Sx, Sy}

5
0
5

1.260
1.259
0.00 0.25 0.50 0.75 1.00
t /T
×10 2

0.00 0.25 0.50 0.75 1.00
t /T(g) {Sy, Sz}

×10 2(h) {Sz, Sx}

1
0
1

1
0
1
0.00 0.25 0.50 0.75 1.00
t /T

0.00 0.25 0.50 0.75 1.00
t /T

Figure 7: NESS in one cycle calculated by the exact numerical integration (blue) andby the HF expansion approach at order N = 0 (orange), 1 (green), and 2 (red). All
the independent 8 observables are plotted in each panel. For each N , we expand
Heff and K(t) up to O(ω−N ) and set kc = b(N + 1)/2c. The parameters are chosen as
Bs = 0.3, Nz = 1, Nx y = 0.05, Bd = 0.1, ω = 10, β = 3, γ0 = 0.2, and Λ = 10.

26

SciPost Phys. Core 4, 033 (2021)

9.4 ×10(a) Sz2

1.46 ×10(b) Sx2 Sy2

2

1.30 ×10

9.3

1.47

1.29

9.2

1.48

1.28

9.1

1.49

9.0

1.50

100

101

102

(c) Sz2

1

1.27
1.26
101

100

1.25 0
10

102

101

102

Figure 8: One-cycle averages of observables in the NESS [Eq. (92)] plotted against
ω. Here ρness (t) is calculated by the exact numerical integration (blue) and by the
HF expansion approach at order N = 0 (orange), 1 (green), and 2 (red). For eachN , we expand Heff and K(t) up to O(ω−N ) and set kc = b(N + 1)/2c. The parameters
are the same as in Fig. 7 except ω.

(a)

10 2

P ss

10 5
10 8
10 11

10 5

order=0, kc = 0
order=1, kc = 1
order=2, kc = 1
order=1, kc = 0
order=2, kc = 0

10 8
101

103

10 11

105

(c)

10 2

10 8

101

103

105

(d)

10 2

P ss

10 5

10 11

(b)

10 2

10 5

order=0, kc = 0
order=1, kc = 1
order=2, kc = 1
order=1, kc = 0
order=2, kc = 0

10 8
101

103

105

10 11

101

103

105

Figure 9: Accuracy of HF expansion approach for the NESS density matrix δρN (73)(left panels) and the Floquet-state population δPNss (105) (right panels) plottedagainst ω for the bath spectral cutoffs Λ = 10 (upper panels) and 106 (lower panels).
Each data set corresponds to a pair of (N , kc ) as shown in legend, for which we use
Heff and K(t) expanded up to O(ω−N ). The other parameters are chosen as Bs = 0.3,
Nz = 1, Nx y = 0.05, Bd = 0.1, γ0 = 0.2, and β = 3.

27

SciPost Phys. Core 4, 033 (2021)

oscillations in, e.g., S 2x − S 2y in Fig. 3. These oscillations cannot be taken up to the firstorder since the micromotion K (1) (t) contains only frequency ω. At the second order,(N =2)
ρness(t) = eG t η = (1 + G t + 12 G t2 )η, where 21 G t2 η involves the frequency 2ω.
The HF expansion approach becomes more accurate for higher frequency ω. Figure 8
shows the ω-dependence of one-cycle averages (72) of observables O = Sz , S 2x − S 2y , and Sz2
in the NESS as the other observables give vanishing averages due to symmetry reasons. Fora more strict measure of the accuracy of the HF approach, we again consider δρN definedin Eq. (73). Figure 9(a) shows δρ for the orders 0, 1, and 2 plotted against ω. We here
observe δρN ∝ ω−(N +1) for high frequencies. This is a clear indication that the N -th order
(N )
exact
approximation ρness(t) completely describe ρness(t) up to O(ω−N ).
We remark the qualitative difference between the time-independent and time-dependentdissipators found in the one-cycle average in a quadrupolar (spin-nematic) operator
1{S x , S y } = 2i((S + )2 − (S − )2 ) [91–94]: While it does not vanish for the time-dependent dissipator, it does for the time-independent one. As shown in Ref. [57], one can understand thisdifference by the following antiunitary symmetry operator V : V S y V † = −S y and V Sα V † = Sα(α = x and z). One can easily check that our Hamiltonian satisfies the following dynamical
symmetry HNV (t) = V HNV (T − t)V † . This means that |ũn (t)〉 ≡ V |un (T − t)〉 is also a Floquet
state having the same quasienergy with |un (t)〉. If there is no degeneracy in quasienergies asis the case of our model, |ũn (t)〉 and |un (t)〉 are the same state: there exists a phase factor e iθn(θn ∈ R) such that |ũn (t)〉 = e iθn |un (t)〉. Noticing V † {S x , S y }V = −{S x , S y }, we have〈un (t)|{S x , S y }|un (t)〉 = 〈ũn (t)|{S x , S y }|ũn (t)〉(102)

= −〈un (T − t)|{S x , S y }|un (T − t)〉(103)= −〈un (t)|{S x , S y }|un (t)〉 ,(104)

which means 〈un (t)|{S x , S y }|un (t)〉=0. Therefore, from the definition of ρness (t) in Eq. (92),P
we have {S x , S y } = tr[ρness (t){S x , S y }] = n Pnss 〈un (t)|{S x , S y }|un (t)〉 = 0. In words, theHamiltonian’s antiunitary symmetry is not broken by the dissipation if it is treated within the
RWA. This is a special property of this class of FLEs that become time-independent in the interaction picture. In general, antiunitary symmetries can be broken as in the time-independentdissipator that we analyze in Sec. 5.3.
Let us discuss the cutoff Λ of the bath spectral function. We have thus far set Λ = 10
and confirmed that the HF approach works well for ω ¦ Λ = 10. For ω ¦ Λ, the subtlety
about the cutoff kc discussed in Sec. 6.1.3 has not been problematic, and kc has not played
important roles. In fact, in Fig. 9(a), we also plot δρN for N = 1 and 2 with kc = 0 that turnout to coincide with those with kc = 1. As another measure of the HF approach’s accuracy, we
introduce
v
u 3
uX
ss
δPN ≡ t (PNss,n − Pnss )2 ,
(105)
n=1

where PNss,n denotes Pnss obtained by our HF approach at the N -th order whereas Pnss does thenumerical exact solution. In Fig. 9(b), we plot δPNss for N = 0, 1 and 2 with kc = 0 and 1
and find that the choice of kc is not relevant. We do not have to set kc to be larger than zero
because each contribution from k 6= 0 in Eq. (100) is negligibly small as γ(εn − εm + kω) ≈ 0
if ω  Λ.
A more nontrivial situation happens when ω  Λ. To verify how our HF approach works
in this case, we perform a similar numerical analysis for Λ = 106 with all the other parameters
being the same. Figures 9(c) and (d) show δρN and δPNss for orders N = 0, 1, and 2 with
28

SciPost Phys. Core 4, 033 (2021)

10 5

FG

10 7
10 9
10 11

= 102
= 104
= 106
slope=-1

10 13
10 15

101

102

103

104

105

Figure 10: Difference between the FGS (93) and the exact NESS (92) for bath spectral cutoff Λ = 102 (circle), 104 (square), and 106 (diamond). The dashed line is
the guide to the eye showing the slope of −1. The other parameters are chosen as
Bs = 0.3, Nz = 1, Nx y = 0.05, Bd = 0.1, γ0 = 0.2, and β = 3.

cutoff kc = 0 and 1. For the correct cutoffs kc = 0 for N = 0 and kc = 1 for N = 1 and 2, we
obtain δρN ∝ δPNss ∝ ω−N −1 , meaning that the N -th order HF expansion correctly describes
the NESS up to O(ω−N ). However, when we set the smaller cutoff kc = 0 for N = 1 and 2,
δρN and δPNss eventually scales like ω−1 as ω increases. For this range of ω, by taking kc
appropriately, we can obtain the NESS at a desired precision with the HF expansion.
Finally, we comment on the accuracy of the FGS, which we quantify byδρFG ≡ kρFG (t) − ρness (t)k .

(106)

qP
FG
ss 2
Equations (92) and (93) lead to δρFG =n (Pn − Pn ) , which is time-independent. We
plot δρFG in Fig. 10 for different bath spectral cutoffs Λ. We note that, to calculate the FGS
exactly, we need numerical integration to obtain the quasienergies εn and time-dependent
Floquet states |ψn (t)〉. As discussed in Sec. 6.1.2, the FGS becomes accurate rapidly when ω
exceeds Λ, where k-photon (k 6= 0) processes are suppressed. Meanwhile, for ω ® Λ, the FGS’serror is as large as O(ω−1 ), which derives from the 1-photon processes of O(ω−1 ) neglectedin the FGS. Thus, for ω ® Λ, the HF expansion approach with appropriately taking k-photon
process gives better description for the NESS without numerical integration. In the following
section, we further investigate the NESS and FGS in another model.

6.3

Example 4: Inverse Faraday Effect in Heisenberg chain and dissipationassisted Floquet engineering

Magneto-optics has been long studied actively [25, 95] and various ultrafast methods of controlling magnetism with light have been discussed. Among them, the inverse Faraday effect [23–25] is a representative phenomenon and it means a magnetization change by applying circularly polarized light to magnetic materials. The sign of varied magnetization canbe controlled by switching the polarization direction (i.e., left- or right-handed polarization).
This effect can be viewed as one of typical Floquet engineerings. The inverse Faraday effecthas been indeed observed in various magnetic materials [96–100] with application of visible
or infrared light whose photon energy is comparable to the electron band energy.

29

SciPost Phys. Core 4, 033 (2021)

(a)

(b)

× 10−5

Circularly polarized laser
Magnetization

M

Dissipationt
Figure 11: (a) Illustration of inverse Faraday effect in the presence of dissipation.
(b) Time evolution of magnetization along z-axis. The blue (orange) line showsa result for a right (left)-handed circularly polarized laser. The parameters are
B = 0.1, ω = ±5, β = 10, Λ = 5, and L = 8.
Its microscopic theory [24] was first developed by Pershan et al., who predicted circularlypolarized light-driven Zeeman coupling appears when the light is applied to conducting electron systems with spin-orbit (SO) coupling. We emphasize that an SO coupling is necessary to
generate the effective Zeeman interaction between electron spins and the ac electric field of
the light.
Recently, in magnetic insulators, the THz-light-driven (or shortly THz) inverse Faraday
effect has been investigated theoretically [101–103]. The THz light is suitable in these systems as its photon energy is comparable to the magnetic excitation energy of, e.g., magnons
and spinons, especially in antiferromagnetic systems. As in conducting electrons, magnetic
anisotropy originated from the SO coupling is necessary for the THz inverse Faraday effect.
Namely, the spin rotation symmetry has to be broken [102] to induce a THz-laser-driven magnetization.
On the other hand, instead of an SO coupling or a magnetic anisotropy, we expect thateven dissipation can also assist the emergence of the THz inverse Faraday effect because the
Lindblad type dissipation can break spin rotation symmetry. In fact, small SO coupling or
spin-rotation-symmetry breaking interactions (such as spin-phonon coupling and dipole interaction) are always present even in ideal magnetic materials. Such weak but finite interactions
may be viewed as a source of dissipation for the systems we consider.
To demonstrate this dissipation-assisted THz inverse Faraday effect, we consider an SU(2)symmetric Heisenberg spin-1/2 chain driven by a circularly polarized laser. The Hamiltonianand dissipator are given by
H(t) = H0 + V (t) = J

LX

S~ j · S~ j+1 + B

L 
X

j=1

D t (ρ) =yS jx cos(ωt) + S j sin(ωt) ,

(107)

j=1



1
γ(ε) Aεj (t)ρAεj† (t) − {Aεj† (t)Aεj (t), ρ} ,
2
j,ε

X

(108)

where S~ j is the spin- 21 operator on jth site and J > 0 is the antiferromagnetic exchange coupling, and B is the coupling constant of the ac Zeeman interaction of the circularly polarizedj
THz laser with frequency ω. In the dissipator D t (ρ), Aε (t) is defined in Eqs. (81) and (82)with A j = S jx . Here, as we did in other sections, we suppose that γ(ε) is the spectral functionof the ohmic boson bath with a Gaussian cutoff Λ and, for convenience, redisplay its form:2− εεe 2Λ2.
γ(ε) = γ01 − e−βε30(109)SciPost Phys. Core 4, 033 (2021)(a) β = 10, B = 0.1(b) β = 10, B = 1(c) β = 1, B = 0.1(d) β = 1, B = 1

M

ω

ω

ω

ω

z
Figure 12: ω-dependence of total magnetization M =j 〈S j 〉 /L for (β, B) =(a)(10, 0.1), (b)(10, 1), (c)(1, 0.1), and (d)(1, 1). The dashed and dotted lines are
the first-order Floquet-Gibbs state (FGS), and the blue and orange solid lines denotethe exact NESS for L = 6 and 8. Note that the lines for L = 6 and 8 are overlappedin (c) and (d).

P

For simplicity, we have ignored the Lamb shift ΛLS(t) here, which does not change the NESS
because both the Lamb shift and NESS are diagonal in the time-independent frame (see, forexample, Eqs. (86) and (87)). In this subsection, we calculate the magnetization dynamics andits value in the NESS by solving the time-dependent Floquet-Lindblad equation of Eqs. (107)and (108)
P L with the forth-order Runge-Kutta method. The initial state is the ground state of
H0 = J j=1 S~ j · S~ j+1 , which is a spin singlet for even L’s. From now on, we fix the parameters
as J = 1 and γ0 = 0.1.
In the absence of the dissipation (i.e., γ0 = 0), the inverse Faraday effect does not occur because of the spin SU(2) symmetry of H0 , which leads to a conservation law
2[H(t), Stot
] = 0 [102]. The symmetry prohibits the magnetization amplitude from growing
up, although its direction changes due to the precession around the circularly polarized laser.
On the other hand, in the presence of the dissipation (i.e., γ0 6= 0), the conservationj2
law breaks down due to the dissipator, [Aε (t), Stot
] 6= 0, and the total magnetization cangrow up. Figure
11
shows
the
inverseFaradayeffectby dissipation. The total magnetiP zzation M =j 〈S j 〉 /L starts from zero and grows up, approaching a nonzero value. We
note that the right- and left-handed circularly polarized lasers give rise to the opposite directions of the magnetization because they are transformed by the π-rotation around the x-axis,yy
S jx → S jx , S j → −S j , and S zj → −S zj .
As discussed in Sec. 6.1, if the driving frequency is larger not only than the system’s energy
scale but also than the bath spectral cutoff Λ, the NESS in this model is approximately the
Floquet-Gibbs state (93) with the effective Hamiltonian
Heff = H0 −B2 X zS + O(1/ω2 ) .
2ω j j(110)

The second term on the right-hand side in Eq. (110) is an effective magnetic field induced by
the circularly polarized laser, giving rise to the inverse Faraday effect. Recall again that, inthe absence of dissipation, the effective magnetic field only cannot induce magnetization, as
illustrated in Fig. 11(b). Figures 12 (a)-(d) show the ω-dependence of the magnetizations in
the exact NESS and FGS for various β and B. For any β and B, the NESS is well described by
the FGS in the high-frequency region, where the magnetization decays as M ∝ ω−1 .
We briefly discuss the system size dependence although the numerically accessible systemsizes are strongly limited by computational complexity in the dissipative quantum many-body
31

SciPost Phys. Core 4, 033 (2021)

(a)

(b)

(c) γ(ϵ)/γ0Λ = 2012108

Λ = 10

6
4

Λ=5

2

M

5

(d)

3.5

10

15

20

25

30

ϵ

10

15

20

25

30

ϵ

β = 0.5

3.0
2.5
2.0
1.5

β=1

1.0
0.5

Λ

β

β = 205

Figure 13: (a)Λ-dependence of total magnetization M . The circles, triangles, andsquares denote the magnetization in the NESS for ω = 10, 20, and 100, respectively.
The dashed lines are the magnetization in the FGS for each ω. The other parameters
are B = 1, β = 10, and L = 8. (b) β-dependence of total magnetization M . The
circles, triangles, and squares denote the magnetization in the NESS for ω = 10, 20,
and 100, and the dash-dotted, dashed, and dotted lines are the magnetization in the
FGS for corresponding ω. The other parameters are B = 1, Λ = 5, and L = 8. (c,d)γ(ε)/γ0 are plotted for (c) Λ = 5, 10, and 20 with β = 10 and (d) β = 0.5, 1, and20 with Λ = 5. In γ(ε), the y-intercept is 1/β and the position of the peak is O(Λ).
Therefore, for example, |γ(0)/γ(ω = 20)| becomes smaller for the larger Λ or thelarger β.
system. Figure 12 shows that the high-frequency expansion predicts the NESS for each of L = 6
and 8. Meanwhile, we observe that the magnetization shows slow system-size convergencefor lower temperatures. The reason why the system-size dependence for β = 10 is larger thanthat for β = 1 is that the low-temperature state is sensitive to the fictitious energy gap dueto the finite-size effect (the Heisenberg chain is known to be gapless in the thermodynamic
limit [104, 105]).
Another subtlety in approaching the thermodynamic limit L → ∞ is the heating and thebreakdown of the HF expansion. However, we do not expect that this subtlety is a serious
problem in our model. This is because our Hamiltonian H(t) becomes time-independent in thezrotating frame |ψ〉 → e−iS ωt |ψ〉, and thus the system does not heat up. We note that, in thegeneric many-body systems, it is known that the high-frequency expansion does not converge
in the thermodynamic limit, at least in the isolated systems [2,31,32]. Understanding the highfrequency expansion in the generic many-body systems with dissipation is an open question,
and we do not go into detail further in this work.
Let us now study the breakdown of the FGS. As discussed in the previous section, theNESS is not necessarily well approximated by the FGS due to the spectral function γ(ε). The
sufficient condition that the FGS is a good approximation is that the contributions of k 6= 0
in Wmn of Eq. (89) are much smaller than that of k = 0 within the desired order of ω. For
example, when the cutoff Λ in γ(ε) of Eq. (109) is much larger than ω, we need to include thecontribution of k = 1 to obtain the accurate first-order HF expansion because γ(kω) ∼ O(ω)β
and |Aα∗A| ∼ O(ω−2|k| ). Then, the approximated NESS is not the FGS. Conversely,
nm;k nm;kwhen the cutoff Λ is much smaller than ω, the contribution of k 6= 0 is negligible because of|γ(0)|  |γ(kω)| (k 6= 0) (see Fig. 13 (c) for the configuration of γ(ε) with varied Λ). Figure 13(a) shows the Λ-dependence of the magnetization in the NESS and the FGS for various ω.
32

SciPost Phys. Core 4, 033 (2021)

We observe that the FGS indeed becomes less accurate for larger Λ and the threshold of the
breakdown is O(ω).
Besides, the FGS becomes less accurate for lower temperature. This is because γ(kω)
becomes relatively larger in comparison with γ(0) = 1/β for the larger β (see Fig. 13 (d)),
meaning that we cannot neglect the contribution of k 6= 0 in Wmn of Eq. (89) in the lowtemperature. Figure 13 (b) shows the β-dependence of the magnetization in the NESS andthe FGS for various ω. Indeed, the FGS becomes less accurate for larger β and the range
where the FGS is valid becomes wider as ω increases.
Finally, we emphasize that the concept of the dissipation-assisted inverse Faraday effectwould be valid even in real materials as well as nanodevices consisting of a few spins although
the analysis can not be quantitatively correct. The FLE with the RWA that we have usedin this section is known not to work well in many-body systems in general since we cannot
use the rotating wave approximation in the derivation. However, even in real materials, thedissipation can break the spin symmetry, giving rise to the inverse Faraday effect, although
the detailed form of the dissipation can be different from the Lindblad-type dissipation. The
detailed analysis for real materials that are not described by the Lindblad equation is future
work. We also note that dissipation-assisted engineered observables can change, dependingon the type of dissipation. As we discussed in Sec. 6.2 and Ref. [57], in the driven single NVcenter model with the time-independent dissipator, we can find a dissipation-assisted Floquet
engineering of a quadrupolar moment {S x , S y }, which is odd for the antiunitary symmetry
operation V .

7Conclusion

In this paper, we have theoretically studied the nonequilibrium steady states (NESSs) in thetime-periodic quantum master equation of Lindblad (or GKSL) form. Considering the highfrequency regime, we have developed a systematic high-frequency (HF) expansion for Liouvillians. One important consequence is that, although the effective Liouvillian is not necessarilyof Lindblad form (Lindbladian), it is still useful to analyze the NESSs. This is mainly because
the effective Liouvillian is trace-preserving (Lemma 1) and the NESS is guaranteed to exist at
each order of the HF expansion (Theorem 1). With the effective Liouvillian and micromotion
superoperators at high-enough expansion order, we can obtain the NESS density matrix bylinear algebra without directly solving the time evolution obeying the Lindblad equation.
We have first applied this theory to concrete models for one-body and many-body quantum
systems subject to phenomenological time-independent dissipators. With the HF-expansion approach, the effects of external drives are manifest in the effective Liouvillian and micromotion
superoperators, which are represented by the commutators of Hamiltonians and dissipators.
This representation, without solving the Lindblad equation, has enabled us to grasp physicalconsequences of the drive and dissipation, such as the effective Hamiltonian in an NV center(Sec. 5.3) the robustness against boundary dissipation in the XY spin chain (Sec. 5.2).
We have then analyzed the microscopically-derived time-dependent dissipators within therotation wave approximation (RWA). Besides the general theoretical framework in Secs. 3 and4, we have developed a slightly different HF-expansion method utilizing a special property ofthis class of dissipators. For these dissipators realized in a weak contact to heat baths, wehave shown that the cutoff Λ of the bath spectral function serves as an important energy scale.
Within these Floquet-Lindblad equations, we have generally shown that the NESS is welldescribed by the Floquet-Gibbs state [Eq. (93)] irrespective of model details as long as thedriving frequency (ħh)ω is higher than both the system’s energy scale and the cutoff Λ. Meanwhile, when (ħh)ω is higher than the system’s energy scale but smaller than the cutoff Λ, richer

33

SciPost Phys. Core 4, 033 (2021)

physics, such as the breakdown of the Floquet-Gibbs state, happen due to dissipation processes
accompanying photon exchange between the system and baths. In these cases, however, theHF-expansion approach, being implemented appropriately, enables a systematic analysis of
the NESS. We have exemplified these in an NV center (Sec. 6.2) and the dissipation-assisted
inverse Faraday effect in the isotropic Heisenberg spin chain.
One important open problem is whether the HF expansion makes sense in generic manybody systems, where the so-called heating problem matters. In isolated many-body systems,
the periodic drive is known to heat the system to the featureless infinite-temperature state, andthus the HF-expansion description is valid at finite time range [2, 31, 32]. This phenomenon isrelated to the fact that the HF expansion for generic many-body systems is not a convergent
series. The example systems that we have analyzed in this paper are so special that the heating
does not occur. However, it is of great interest how the heating is compensated by dissipationand the convergence property of the effective Liouvillian in generic many-body systems. To
address these systems, we may need to consider Floquet-Lindblad equations beyond the RWA,
which is another important open problem.

Acknowledgements
Fruitful discussions with M. Holthaus, K. Mizuta, T. Mori, F. Nathan, A. Polkovnikov, M. Rudner,
T. Satoh, and H. Tsunetsugu are gratefully acknowledged.
Funding information T.N.I was supported by JSPS KAKENHI Grant No. JP21K13852. K.C.
was supported by JSPS KAKENHI Grant No. 21J11245 and Advanced Leading Graduate Course
for Photon Science at the University of Tokyo. M.S. was supported by JSPS KAKENHI (Grants
No. 17K05513 and No. 20H01830) and a Grant-in-Aid for Scientific Research on Innovative
Areas “Quantum Liquid Crystals” (Grant No. JP19H05825).

Avan Vleck high-frequency expansion in isolated systems

We recap the results of Refs. [8, 9]. In an isolated system, the time evolution is given by the
Schrödinger equation.
i∂ t |Ψ(t)〉 = H(t)|Ψ(t)〉 .
(111)
The propagator is decomposed as00
U t, t 0 = e−iK(t) e−iHeff (t−t ) e iK(t ) ,

34

(112)

SciPost Phys. Core 4, 033 (2021)

Heff =∞X(n)

Heff ,

(113)

n=0
(0)

Heff = H0 ,(114)
X [H , H ]
(1)
−m
m
,
(115)
Heff =
2mω
m6=0
X [[H , H ] , H ] X X [[H , H
(2)
−m
0
m
−m
m−n ] , H n ]Heff =
+
,
(116)
2
2
2m ω
3mnω2
m6=0
m6=0 n6=0,m
X [[[H , H ] , H ] , H ] X X [[[H , H ] , H
(3)
−m
0
0
m
−m
0
m−n ] , H n ]Heff =
+
3 ω3
2 nω3
2m
3m
m6=0
m6=0 n6=0,m
X X [[[H , H
X [[[H , H ] , H ] , H ]
]
,
H
]
,
−m
m−n
0 Hn]
−m
m
−n
n
+
−
(117)
2 ω3
2 ω3
4mn
12mn
m6=0 n6=0,mm,n6=0X X [[H , H ] , [H
X X [[[H , H ] , H ] , H ]−m0m−n , H n ]]
−m
m−l
l−n
n
+
+
2 nω3
3
12m
6l
mnω
m6=0 n6=0,mm,n6=0 l6=0,m,n
X X [[[H , H
X X [[H , H ] , [H , H ]]−mm−n−l ] , H l ] , H n ]
−m
m−l
l−n
n
+
+
,
3
3
24l
mnω
24l
mnω
m,n6=0 l6=0,m−nm,n6=0 l6=0,m,nK(t) =∞
X

K (n) (t) ,

(118)

n=0

iK (1) (t) = −X Hme−imωt ,
(119)
mω
m6=0
X X [H , H
X [H , H ]nm−n ] −imωtm
0 −imωtiK (2) (t) =
e
+
e
,
(120)
2
2 ω2
2mnω
m
m6=0 n6=0,mm6=0
X [[H , H ] , H ]X X [H , [H , H ]]
m
−n
n
m
0
0 −imωtiK (3) (t) = −e+e−imωt
3 ω3
2 nω3
m
4m
m6=0
m6=0 n6=0
X X [[H , H ] , H
X X [[H , Hn0m−n ] −imωtnm−n ] , H 0 ] −imωt
−
e
−
e
(121)
2
3
2mn ω
2m2 nω3
m6=0 n6=0,mm6=0 n6=0,m
X X X [[H , H ] , H ]X X X [H , [H , H
n
l−n
m−l
n
l
m−n−l ]] −imωt
−imωt
−
e
−
e
,
3
3
4mnlω
12mnlω
m6=0 n6=0 l6=0,n,mm6=0 n6=0 l6=0,m−n35SciPost Phys. Core 4, 033 (2021)

B

Commutator calculations

Lemma 2 Suppose iLm ρ = [H m , ρ]. We define for N ≥ 2

KmN ,...,m2 ,m1 = [iLmN , [· · · , [iLm2 , iLm1 ]]] ,(122)

H mN ,...,m2 ,m1 = [H mN , [· · · , [H m2 , H m1 ]]] .

(123)

KmN ,...,m2 ,m1 (ρ) = [H mN ,...,m2 ,m1 , ρ] ,(124)

Then, it follows thatfor any ρ .
Proof. We invoke the mathematical induction. First, for N = 2, the statement follows fromthe Jacobi identity [A, [B, C]] + [B, [C, A]] + [C, [A, B]] = 0. In fact,

Km2 ,m1 (ρ) = [H m2 , [H m1 , ρ]] − [H m1 , [H m2 , ρ]](125)= [H m2 , [H m1 , ρ]] + [H m1 , [ρ, H m2 ]](126)= −[ρ, [H m2 , H m1 ]](127)= [H m2 ,m1 , ρ] .

(128)

Next, we suppose that the statement is true for N = M and prove the statement for N = M +1.
This step is achieved as follows:

KmM +1 ,...,m2 ,m1 (ρ) = iLmM +1 KmM ,...,m2 ,m1 (ρ) − KmM ,...,m2 ,m1 iLmM +1 (ρ)

C

(129)

= [H mM +1 , [H mM ,...,m2 ,m1 , ρ]] − [H mM ,...,m2 ,m1 , [H mM +1 , ρ]](130)= [[H mM +1 , H mM ,...,m2 ,m1 ], ρ](131)= [H mM +1 ,...,m2 ,m1 , ρ] .

(132)

Majorana representation of the XY chain

In Sec. 5.2, we have used the Majorana-fermion representation for the infinite XY chain. Here,
we derive the representation and its Fourier transform.
We introduce 2N Majorana fermions w j ( j = 1, 2, . . . , 2N ) for an N -site chain as
w2n−1 =n−1Y!
σzj σnx ,

(133)

j=1

w2n =

n−1
Y

!
σzj σny ,j=1
w2n−1 w2n = iσzn ,

(134)
(135)

which satisfy the Majorana commutation relations {w i , w j } = 2δi j and translate Eq. (50) into
H(t) = −i

N
−1 
X
j=1


N
X
1− g
1+ g
w2 j w2 j+1 −w2 j−1 w2( j+1) − ih f (t)w2 j−1 w2 j .
2
2
j=1
36

(136)

SciPost Phys. Core 4, 033 (2021)

It is convenient to introduce the two-component fermions W j = t (w2 j−1 , w2 j ), which lead to1(w2 j−1 w2 j − w2 j w2 j−1 ) = t W j iτ y W j ,2
1
1
1
w2 j w2 j+1 = (w2 j w2 j+1 − w2 j+1 w2 j ) = t W j τ− W j+1 − t W j+1 τ+ W j ,
2
2
2
1
1t
1
w2 j−1 w2( j+1) = (w2 j−1 w2( j+1) − w2( j+1) w2 j−1 ) = W j τ+ W j+1 − t W j+1 τ− W j ,222
w2 j−1 w2 j =

(137)
(138)
(139)

where τ± ≡ (τ x ± iτ y )/2 and τα (α = x, y, and z) are the Pauli matrices. These equations
give


N
−1
N
X
X

i
t
t
H(t) = − 
W j (gτ x − iτ y )W j+1 − t W j+1 (gτ x + iτ y )W j + 2h f (t)W j iτ y W j 
4 j=1
j=1
(140)
it
WXW ,
(141)
4
where W is the 2N -component vector, and X is the antiHermitian 2N × 2N matrix. Thus, wehave obtained the quadratic Majorana representation of the Hamiltonian, which is translation
invariant except for the boundaries.
To obtain the quasienergy bands, we derive the time-evolution equation in the Majorana
representation. It is convenient to work in the Heisenberg picture to have=−d W (t)= i[H(t), W (t)] = −X W (t) ,(142)
dt
where W (t) is the Heisenberg-picture operator for W and we have used the commutation
relation {w i , w j } = 2δi j . To recover the translation symmetry, we consider the infinite chain
P −ik jfk ∝(N → ∞) and invoke the Fourier transform: W j → W
W j . In the Fourier basis,je
we havedWk (t)= −X (k, t)Wk (t) ,
(143)
dt
where
X (k) ≡ (gτ x − iτ y )e−ik − (gτ x + iτ y )e ik + 2h f (t)iτ y = −2i{g sin kτ x + [cos k − h f (t)]τ y } , (144)is a 2 × 2 matrix for each k. Physically, the 2 × 2 matrix describes the evolution of the quasiparticles.
According to the Floquet theory, the one-cycle evolution
 Z T
V (k) = exp+ −X (k, t)d t

,

(145)

0

determines the Floquet modes and the two eigenvalues define the quasienergy bands. Since
X (k, t) is antiHermitian and traceless, the eigenvalues of V (k) are given as e±iε(k)T , and we
regard ±ε(k) as the quasienergy bands forÆthe Majorana fermions. In the absence of theexternal drive f (t) = 0, we have ε(k) = 2 cos2 k + g 2 sin2 k. Note that dε(k)/d k = 0 at
k = 0 and ±π.
To use the HF expansion for the Hamiltonian, we notice the formal analogy of Eq. (143)to the Schrödinger equation. In fact,h(k, t) ≡ −iX (k, t) = −2{(g sin k)τ x + [cos k − h f (t)]τ y } ,(146)
RTplays the role of Hamiltonian and gives V (k) = exp+ [−i 0 h(k, t)d t], and, hence, the HF
expansion in Appendix A is applicable to h(k, t).
37

SciPost Phys. Core 4, 033 (2021)

D

Derivation of time-dependent FLE for weak thermal contact

Historically, dissipative Floquet systems weakly coupled to thermal reservoirs were studied byquantum master equations [41, 85], which are equivalent to the FLE. Here, for completeness,
we summarize the derivation with slight generalization and emphasis on the Lindblad (GKSL)
form. Our derivation follows a textbook [40] in which the GKSL equation is derived for thetime-independent HS .

D.1

Interaction picture

We begin by considering the system-bath composite system whose Hamiltonian reads
H(t) = HS (t) + H B + HSB ,(147)

where HS (t + T ) = HS (t)Pis for the periodically driven system of interest, H B for the heat bath(reservoir), and HSB = α Aα ⊗ Bα for the system-bath coupling. Here, precisely speaking,HS (t) (H B ) represents HS (t) ⊗ 1 (1 ⊗ H B ), but we omit ⊗1 and 1⊗ for brevity.
Under the Hamiltonian (147), we consider the time evolution of the total system. To thisend, it is useful to work in the interaction picture:dρ I (t)= −i[H I (t), ρ I (t)] ,
dt

(148)

and its integral formII

ρ (t) = ρ (0) − i

Z t[H I (s), ρ I (s)]ds ,

(149)

0

where ρ I (t) is the density matrix for thesystem in the Rtotal interaction picture, andt
H I (t) = U † (t)HSB U(t) with U(t) = exp+ −i 0 (HS (s) + H B )ds . In the following, we drop
the superscript I for density matrices for brevity until we arrive at the final result (170). We assume that, at the initial time t = 0, the system and bath are not entangled: ρ(0) = ρS (0)⊗ρB ,where ρB represents a static state of the bath. We then substitute Eq. (149) into the RHS of
Eq. (148) and take TrB of the both sides, havingdρS (t)=−dt

Z t

dsTrB [H I (t), [H I (s), ρ(s)]] ,

(150)

0

where we have used TrB [ρ(0)Bα ] = 0 for all α without loss of generality. To make this equation
Markovian, we make the replacement ρ(s) → ρS (t) ⊗ ρB , by which the RHS only depends on
ρS (t). We then remove the dependence on the initial time by replacing s in the integral byRtR∞t − s and extending the integration as 0 → 0 . Then, we obtain the following Markovian
equation:dρS (t) = −dt

Z∞dsTrB {[H I (t), [H I (t − s), ρS (t) ⊗ ρB ]]} ,

(151)

0

which is Markovian but not in Lindblad form yet. We remark that these approximations are
based on physical intuitions, but the errors accompanied by them have recently been quantified
in a mathematically rigorous manner [86, 87]. According to the error bounds, the approximations are valid if the system-bath coupling is weak enough and the bath-correlation time isshort enough.
38

SciPost Phys. Core 4, 033 (2021)

We can simplify Eq. (151) by introducing the
P bath spectral function. To this end, we decompose the system-bath interaction as HSB = α Aα ⊗ Bα and correspondinglyX
H I (t) =Aα (t) ⊗ Bα (t) ,
(152)
α


 Rt
where Aα (t) = US† (t)Aα US (t) with US (t) = exp+ −i 0 dsHS (s) , and Bα (t) = eiH B t Bα e−iH B t .
We define the bath spectral function by〈Bα (t)Bβ (t − s)〉 ≡ TrB ρB Bα (t)Bβ (t − s) .
(153)
Although the above arguments apply to any ρB , we here assume ρB = e−β H B /ZB . Then we have〈Bα (t)Bβ (t − s)〉 = 〈Bα (s)Bβ (0)〉 and the so-called Kubo-Martin-Schwinger (KMS) condition:〈Bα (t)Bβ (0)〉 = 〈Bβ (s)Bα (t + iβ)〉. By using the bath spectral function, Eq. (151) is rewritten
as
XZ ∞


d
ρS (t) =ds〈Bα (s)Bβ (0)〉 Aβ (t − s)ρS (t), Aα (t) + H.c.
(154)
dt
0
α,β
Equation (154) can be further simplified once we Fourier expand Aα (t) and Aβ (t − s). For
the conventional setup where HS is time-independent, the Fourier expansion is done simply
by the decomposition into the eigenmodes of HS . However, in our periodically driven systems,
the decomposition into the Floquet states [6, 7] are useful:id|ψm (t)〉 = HS (t) |ψm (t)〉 ;dt|ψm (t)〉 = e−iεm t |um (t)〉 ;|um (t + T )〉 = |um (t)〉 .
(155)

Here, {|ψm (t)〉}Nm=1 (N < ∞ is the Hilbert-space dimension) are independent solutions forthe time-dependent Schrödinger equation without system-bath coupling, and we call εm and|um (t)〉 the quasienergy and Floquet state, respectively. We remark that εm and |um (t)〉 are
not uniquely defined, but we can redefine them by εm → εm + kω and |um (t)〉 → e iεm t |um (t)〉without changing |ψm (t)〉 and the periodicity of |um (t)〉. For concrete calculations, one may
fix a set of εm ’s, but the physical observables must be invariant under those shifts.
These states give
US (t) =NX|ψm (t)〉 〈ψm (0)| ,(156)|ψm (0)〉 〈ψm (t)|Aα |ψn (t)〉 〈ψn (0)| .

(157)

m=1

which leads toAα (t) =Xm,n

We note that 〈ψm (t)|Aα |ψn (t)〉 = ei(εm −εn )t 〈um (t)|Aα |un (t)〉 has a discrete spectrum, in which
nonzero weights lie at ε = −(εm − εn ) + `ω (` ∈ Z), since the Floquet states |um (t)〉 areperiodic. Thus we introduce the spectral decomposition for the matrix elements asXX〈ψm (t)|Aα |ψn (t)〉 =Aαmn (ε)e−iεt =Aαmn (εn − εm + `ω)e−i(εn −εm +`ω)t ,
(158)
ε

`∈Z

and correspondingly that for the operators asXXAα (t) =e−iεt Aαε ;Aαε ≡Aαmn (ε) |ψm (0)〉 〈ψn (0)| .
ε

m,n

39

(159)

SciPost Phys. Core 4, 033 (2021)

For later use, let us show that Aαε and (Aαε )† are operators that lower and raise quasienergyby ε:
U F Aαε U F† = e+iεT Aαε ,†−iεT α†U F Aα†Aε ,ε UF = e

(160)

where
U F ≡ US (T ) =NX|ψm (T )〉 〈ψm (0)| =

m=1

N
X

e−iεm T |ψm (0)〉 〈ψm (0)|

(161)

m=1

denotes the one-cycle unitary evolution. To prove Eq. (160), we first notice thatU F |ψm (0)〉 = |ψm (T )〉 = e−iεm T |ψm (0)〉 ,(162)

meaning that |ψm (0)〉 is an eigenstate of U F with eigenvalue e−iεm T . This equation also implies〈ψn (0)| U F† = e iεn T 〈ψn (0)| .
Using Eqs. (162) and (163) together with Eq. (159), we haveX
Aαmn (ε)e i(εn −εm )T |ψm (0)〉 〈ψn (0)|U F Aαε U F† =

(163)

(164)

m,n

=

X

Aαmn (ε)e iεT |ψm (0)〉 〈ψn (0)|m,n= e iεT Aαε .

(165)
(166)

Here, to obtain the second line, we have used the fact that Aαmn (ε) is nonvanishing only whenεn − εm = ε + `ω for some ` ∈ Z. Equation (166) is the first equation in Eq. (160) that we
wanted to prove, and we obtain the second one by taking the Hermitian conjugation of thefirst one.
We are ready to obtain the Lindblad form with the final approximation, the rotating wave
approximation (RWA). Substituting Eq. (159) into Eq. (154), we have
!
X


d
α†
i(ε−ε0 )t
β
α† β
(167)
ρS (t) =Γαβ (ε)eAε ρS (t)Aε0 − Aε Aε0 ρS (t) + H.c. ,dt
α,β,ε,ε0
where we have defined the bath spectral function
Γαβ (ε) =Z∞

ds eiεs 〈Bα (s)Bβ (0)〉 .

(168)

0

Since Γαβ (ε) is not a Hermitian matrix, it is convenient to decompose it into the Hermitianand antiHermitian parts:
Γαβ (ε) =1
γαβ (ε) + iSαβ (ε) .
2

(169)

The RWA means that we keep ε = ε0 terms only in Eq. (167), by which we arrive at thefollowing equation of Lindblad form:d I
ρ (t) = −i[ΛLS, ρSI (t)] + D I (ρSI (t)) ,
dt S

40

(170)

SciPost Phys. Core 4, 033 (2021)

where the Lamb shift ΛLS and dissipator D I are defined by
X
β
ΛLS =
Sαβ (ε)Aα†ε Aε ,

(171)

α,β,ε

DI =



1  α† β I
γαβ (ε) Aβε ρSI (t)Aα†
−
A
A
,
ρ
(t)
,
ε
2 ε ε S
α,β,ε
X

(172)

Here we have recovered the superscript I on ρS (t) to emphasize that this is in the interaction
picture.
The RWA is valid when quasienergy spacings are large enough but becomes invalid if thequasienergy spacings become infinitely small with the system-bath coupling and the bath correlation time being fixed. In this sense, Eq. (170) cannot apply to large nonintegrable manybody systems, where (quasi)energy spacings decrease exponentially in the system size. Thus,
Eq. (170) should apply, e.g., to two- or few-level quantum systems like an atom, small clustersof spins, or some integrable quantum systems with nonvanishing energy gaps. There have been
new approaches to derive master equations of Lindblad form without using the RWA [86, 87].
We note that our derivation for Eq. (170) is slightly generalized from that ofthe previous studies [41, 85] in the following sense. We recall that ε takes the values in{∆mn (`) ≡ εn − εm + `ω}m,n,k . Here, we emphasize that the correspondence between ε and(m, n, `) is not one-to-one in general, and there are multiple (m, n, `)’s that give the same ε.
This situation occurs, for example, in the following three cases:
1. The quasienergies {εm } are degenerate.
2. The quasienergy differences {εm − εn } are degenerate, such as for equidistant {εm }.
3. There are some pair (m, m0 ) such that εm − εm0 = ω/2 .
Our derivation of Eq. (170) is valid even if the correspondence between ε and (m, n, `) is not
one-to-one in contrast to that in the previous studies [41, 85].
We note that the Lamb shift ΛLS does not change the quasienergy. This follows from theβabove-shown fact that Aα†ε (Aε ) raises (lowers) the quasienergy by ε. As an operator expression, we have
[U F , ΛLS] = 0 ,(173)which follows from Eq. (160). This is a Floquet counterpart of the following property: The
Lamb shift does not change the energy of the system in the time-independent case. As a consequence, when the quasienergies are not degenerate, ΛLS is diagonal in the Floquet eigenbasisXΛLS =λm |ψm (0)〉 〈ψm (0)| ,(174)m

and hence only shifts the quasienergy by λm .

D.2

Schrödinger picture

Equation (170) is written in the interaction picture. To go back to the Schrödinger picture, weremember ρS (t) = US (t)ρSI (t)US† (t). Then we havedρS (t) = −i[HS (t) + ΛLS(t), ρS (t)]dtX
1  α†β+γαβ (ε) Aβε (t)ρS (t)Aα†
(t)
−
A
(t)A
(t),
ρ
(t)
.
S
ε
ε
2 ε
α,β,ε
41

(175)

SciPost Phys. Core 4, 033 (2021)

Here the jump operators are given byAαε (t) = e−iεtX

ΛLS(t) =X

Aαmn (ε) |ψm (t)〉 〈ψn (t)| ,

(176)

m,n

and the Lamb shift is by

α,β,ε

β
Sαβ (ε)Aα†
ε (t)Aε (t) .

(177)

Here, in Aαε (t), we put the oscillating phase factor e−iεt so that Aαε (t) becomes time-periodicβ(this phase factor keeps the Lindbladian invariant as Aα†ε (t) and Aε (t) appear in pairs).
To show the periodicity Aαε (t) = Aαε (t + T ), we write Eq. (176) in terms ofPαi(−ε−εm +εn )t|um (t)〉: Aαε (t) =|um (t)〉 〈un (t)|, in which we need to provem,n Amn (ε)ei(−ε−εm +εn )te
is periodic (|um (t)〉 are periodic by definition). Since Aαmn (ε) is nonvanishing
only when ε = εn − εm + kω for integers k, the oscillating phase factor e i(−ε−εm +εn )t becomese−ikωt , which is periodic.
It is convenient to introduce the dissipator D t , which is the following superoperatorX
1  α†α†ββD t (ρ) ≡γαβ (ε) Aε (t)ρAε (t) −(178)
A (t)Aε (t), ρ .
2 ε
α,β,ε
We also introduce the Lindbladian superoperator by

L t (ρ) = −i[HS (t) + ΛLS(t), ρ] + D t (ρ) .

(179)

One can easily confirm that this FLE reduces to the well-known Lindblad equation for undriven
systems [40], noting that, in undriven systems, each Floquet eigenstate |ψm (t)〉 corresponds
to the time-evolving energy eigenstate e−iEm t |Em 〉 for HS |Em 〉 = Em |Em 〉.
It is noteworthy that the Lindbladian of the FLE is time-independent in the interaction
picture (170) while not in the Schrödinger picture (175). This is an extraordinary propertyof the FLE based on the weak system-bath coupling and the RWA, and we cannot find such anice frame in which the time dependence is eliminated for a general FLE in the Schrödinger
picture. In this sense, the FLE derived in this Appendix is special, making analysis easier.

D.3

Time evolution of density matrix

Here we analyze the time evolution of density matrices, showing that their off-diagonal elements decay whereas the diagonal ones obey the master equation. In Ref. [85], these properties were obtained under the strong assumption: If εn − εm + `ω = εn0 − εm0 + `0 ω, then(m, n, `) = (m0 , n0 , `0 ) holds true. Our argument here relaxes this assumption and also applies
when the quasienergies are nondegenerate as long as the NESS is unique. We also discuss how
the results change for degenerate quasienergies.
D.3.1

Nondegenerate quasienergies

We first discuss the case when the quasienergies are not degenerate. To analyze the time
evolution, it is convenient to work in the interaction picture (170) and represent the density
matrix in the Floquet-state basis:XρSI (t) =σmn (t) |m〉 〈n| ,
(180)
m,n

42

SciPost Phys. Core 4, 033 (2021)

where we have introduced the abbreviation |m〉 = |ψm (0)〉. Substituting Eq. (180) into
Eq. (170), we haveXX
dσmn (t)= −i(λm − λn )σmn (t) +γαβ (ε) 〈m|Aβε |m0 〉 σm0 n0 (t) 〈n0 |Aα†ε |n〉dtm0 ,n0
α,β,ε

1X
1X
α† β
0
0 α† β−〈m|Aε Aε |m 〉 σm0 n (t) −σmn0 (t) 〈n |Aε Aε |n〉 .
(181)
2 m0
2 n0
β

To simplify the sums in Eq. (181), we recall that Aε and (Aαε )† are operators that lower andβraise quasienergy by ε as shown in Eq. (160). Thus, Aα†ε Aε does not raise or lower quasienegies,
β0meaning that 〈m|Aα†ε Aε |m 〉 ∝ δmm0 as the quasienergies are not degenerate. Together with
Eq. (88), we can rewrite Eq. (181) as

dσmn (t)1X= −i(λm − λn )σmn (t) −(Wn0 m + Wn0 n ) σmn (t)
dt
2 n0
X
X
〈m|Aβε |m0 〉 σm0 n0 (t) 〈n0 |Aα†+γαβ (ε)ε |n〉 .
α,β,ε

(182)

m0 ,n0

We now have a set of N 2 differential equations for 1 ≤ m, n ≤ N , part of which are coupled by
the final term on the right-hand side (RHS) of Eq. (182).
Here, we can prove the remarkable property of this set of equations that the diagonal
ones ({σmm (t)}Nm=1 ) and off-diagonal ones ({σmn (t)}m6=n are decoupled with each other. To
show this, we first substitute n by m in Eq. (182) and prove that only m0 = n0 terms appearβ
on the RHS. As Aε (Aα†ε ) lowers (raises) the quasienergy by ε, σm0 n0 (t) appears only whenεm ≡ εm0 − ε and εm ≡ εn0 − ε (mod ω), implying m0 = n0 . Thus, we have obtained the resultthat the diagonal elements {σmm (t)}m form a closed set of equations. Second, we suppose
m 6= n in Eq. (182) and prove that only m0 6= n0 terms appear on the RHS. The proof goessimilarly. One can easily check that the set of equations for the diagonal elements are thesame as Eq. (87).
The off-diagonal elements σm6=n (t) form a complex set of equations except for the special
case studied previously [85], where εn − εm + `ω = εn0 − εm0 + `0 ω means m = m0 , n = n0 , and` = `0 . In this special case, the final term on the RHS of Eq. (182) vanishes for m 6= n, andeach σm6=n (t) obeys a closed differential equation. As Wn0 m + Wn0 n > 0, the equation meansthat each σm6=n (t) exponentially decreases with oscillation to vanish as t → ∞. However,this special case does not necessarily occur even if we assume that the quasienergies are notdegenerate (recall the example cases 2 and 3 discussed in Sec. D.1.
Nevertheless, in the weaker assumption of nondegenerate quasienergies, we can provethat the off-diagonal elements σm6=n vanish in the NESS as long as the NESS is unique. For
this purpose, we invoke a symmetry argument in the Lindbladian systems [106–108]. Let us
symbolically represent the interaction picture (170) as
d I
ρ (t) = L I [ρSI (t)] .
dt S

(183)

Now we notice that the Lindbladian L I has the following weak symmetry [106] associatedwith the quasienergy:[L I , U F ] = 0 ,(184)

where U F is a superoperator defined by U F (ρ) ≡ U F ρU F† . One can easily confirm Eq. (184)by using Eqs. (160) and (173). One important consequence of this weak symmetry is that43

SciPost Phys. Core 4, 033 (2021)

P
there exists a stationary solution of the following form [106]: ρ = µ pµ τµ , where {τµ } arematrices belonging to the zero eigenvalue of U F . For our U F , these zero-eigenvalue states areτµ = |n〉 〈n|. As we have assumed that the NESS is unique, this type of stationary solution is
the NESS, implying that σm6=n (t) → 0 as t → ∞.
D.3.2

Degenerate quasienergies

When the quasienergies εm ’s are degenerate, we cannot, in general, have the decoupling of
the diagonal and off-diagonal elements σmn (t) of the density matrix. This is due to, e.g.,β0〈m|Aα†ε Aε |m 〉 6∝ δmm0 .
Yet, the weak symmetry constrains the form of the NESS density matrix. To see this, we
introduce a new index a that distinguishes the degenerate Floquet states. In this new notation,
Eq. (161) is represented as
UF =

d
X

e

−iεm Tm=1

Nm
X|m, a〉 〈m, a| ,

(185)

a=1

where d is the number of distinct quasienergies, Nm the degree of degeneracy of eachquasienergy εm , and |m, a〉 = |ψm,a (0)〉. Thus, the zero-eigenvalue states for the superoperator
U F (ρ) = U F ρU F† are |m, a〉 〈m, b|, and the weak symmetry tells us that the NESS is represented
as
I
=
ρness

Nm
d X
X

m
cab
|m, a〉 〈m, b| .

(186)

m=1 a,b=1
m
m
The expansion coefficients within each degenerate subspace cab
depend on models. As cab
areHermitian matrices for each m, they are diagonalizable by a unitary basis transformation, soIis ρness. Note, however, that the differential equations for σmn (t) cannot, even in this basis,be separated for the diagonal and off-diagonal elements like in the nondegenerate case.

D.4

Conditions for FGS when quasienergies are degenerate

Here we show that, within the FLE obtained by the RWA, the NESS coincides with the FGS
for an extremely large frequency even if the quasienergies are degenerate. In this subsecItion, we work in the interaction picture (170) and discuss when ρness
sufficiently approaches
I−β H FρFG = e/Z, whereX
HF =εm |ψm (0)〉 〈ψm (0)| .
(187)mAs discussed in Sec. 6.1.1, the quasienergies εm ’s are defined only modulo ω, and the FGS is
ill-defined in general. Here we consider the situation that ω is much larger than the system’senergy scale, i.e., H F becomes close enough to H0 (H0 is the undriven Hamiltonian). In sucha case, we can uniquely determine εm ’s so that εm − Em = O(1/ω) holds for each m (Em are
the eigenenergies of H0 ). Also, we assume that ω is so large that|εm − εn | < ω ,∀m, n.

(188)

I
To judge if ρFG= e−β H F /Z is a steady-state solution we substitute it into the right-handside of the FLE (170) and ask if it vanishes. Assuming the basis on which ΛLS is diagonal (174),
II
we have [ΛLS, ρFG] = 0. Thus, we are to show D I (ρFG
) = 0.

44

SciPost Phys. Core 4, 033 (2021)

I
First, let us argue that D I (ρFG
) = 0 holds if we are allowed to use[H F , Aαε ] = −εAαε ,α†[H F , Aα†ε ] = εAε ,(189)

which are stronger versions of the commutation relations (160). Assuming Eq. (189) as aworking hypothesis, we haveX
1  α† β −β H F
II−1
β −β H F α†
(190)
A A ,eD (ρFG ) = Z
γαβ (ε) Aε e
Aε −
2 ε ε
α,β,ε
X

α† β
= Z −1 e−β H F
e−βε γαβ (ε)Aβε Aα†
−
γ
(ε)A
A
(191)
αβ
ε
ε ε ,
α,β,ε

β

β

β

−β H F
where we have used [H F , Aα†= e−βε e−β H F Aε , which follow from ourε Aε ] = 0 and Aε eworking hypothesis (189) and the Baker–Campbell–Hausdorff formula. Then we rewrite thefirst term in the sum of Eq. (191) by changing the dummy index as ε → −ε, havingXXβe−βε γαβ (ε)Aβε Aα†eβε γαβ (−ε)A−ε Aα†
(192)
−ε
ε =
α,β,ε

α,β,ε

=

X
α,β,ε

=

X
α,β,ε

β

eβε γαβ (−ε)Aε† Aαε

(193)

α
γβα (ε)Aβ†ε Aε ,(194)

where we have used the KMS condition (76) and Aα−ε = (Aαε )† . Substituting Eq. (194) into
Eq. (191) and changing the dummy indices appropriately, we obtainI) = 0.
D I (ρFG(195)

Thus, we have shown that the FGS is a steady-state if our working hypothesis (189) are justified.
Finally, we justify Eq. (189) when the driving frequency is much larger than the bath spectral cutoff. From the definition of Aαε (81), we haveX
[H F , Aαε ] =(εm − εn )Aαmn (ε) |ψm (0)〉 〈ψn (0)| .
(196)
m,n

We recall that Aαmn (ε) is nonvanishing only when ε = εn − εm + kω for integer k’s. However,like in the argument in Sec. 6.1.2, the k 6= 0 contributions are vanishingly small because itis accompanied by the bath spectral function γαβ (ε) = γαβ (εn − εm + kω) ≈ 0 if ω  Λ.
Therefore, we can neglect k 6= 0 contributions and set εm − εn = −ε in Eq. (196), which leads
toX[H F , Aαε ] =(−ε)Aαmn (ε) |ψm (0)〉 〈ψn (0)| = −εAαε .
(197)
m,n

Thus, we have shown that Eqs. (189) are justified when the driving frequency is much larger
I
I
than the bath spectral cutoff and hence ρness≈ ρFG.

References[1] M. Holthaus, Floquet engineering with quasienergy bands of periodically driven optical lattices, J. Phys. B: At. Mol. Opt. Phys. 49, 013001 (2015), doi:10.1088/09534075/49/1/013001.
45

SciPost Phys. Core 4, 033 (2021)[2] M. Bukov, L. D’Alessio and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to Floquet engineering, Adv. Phys. 64,
139 (2015), doi:10.1080/00018732.2015.1055918.
[3] A. Eckardt, Colloquium: Atomic quantum gases in periodically driven optical lattices, Rev.
Mod. Phys. 89, 011004 (2017), doi:10.1103/RevModPhys.89.011004.
[4] T. Oka and S. Kitamura, Floquet engineering of quantum materials, Annu. Rev. Condens.
Matter Phys. 10, 387 (2019), doi:10.1146/annurev-conmatphys-031218-013423.
[5] A. Haldar and A. Das, Statistical mechanics of Floquet quantum matter: The exact
and emergent conservation laws, J. Phys.: Condens. Matter (2021), doi:10.1088/1361648X/ac03d2.
[6] J. H. Shirley, Solution of the Schrödinger equation with a Hamiltonian periodic in time,
Phys. Rev. 138, B979 (1965), doi:10.1103/PhysRev.138.B979.
[7] H. Sambe, Steady states and quasienergies of a quantum-mechanical system in an oscillating field, Phys. Rev. A 7, 2203 (1973), doi:10.1103/PhysRevA.7.2203.
[8] A. Eckardt and E. Anisimovas, High-frequency approximation for periodically driven
quantum systems from a Floquet-space perspective, New J. Phys. 17, 93039 (2015),
doi:10.1088/1367-2630/17/9/093039.
[9] T. Mikami, S. Kitamura, K. Yasuda, N. Tsuji, T. Oka and H. Aoki, Brillouin-Wigner theory
for high-frequency expansion in periodically driven systems: Application to Floquet topological insulators, Phys. Rev. B 93, 144307 (2016), doi:10.1103/PhysRevB.93.144307.
[10] D. H. Dunlap and V. M. Kenkre, Dynamic localization of a charged particle moving under the influence of an electric field, Phys. Rev. B 34, 3625 (1986),
doi:10.1103/PhysRevB.34.3625.
[11] B. G. Klappauf, W. H. Oskay, D. A. Steck and M. G. Raizen, Observation of noiseand dissipation effects on dynamical localization, Phys. Rev. Lett. 81, 1203 (1998),
doi:10.1103/PhysRevLett.81.1203.
[12] Y. Kayanuma and K. Saito, Coherent destruction of tunneling, dynamic localization, and the Landau-Zener formula, Phys. Rev. A 77, 010101 (2008),
doi:10.1103/PhysRevA.77.010101.
[13] N. Tsuji, T. Oka, P. Werner and H. Aoki, Dynamical band flipping in Fermionic latticesystems: An ac-field-driven change of the interaction from repulsive to attractive, Phys.
Rev. Lett. 106, 236401 (2011), doi:10.1103/PhysRevLett.106.236401.
[14] T. Ishikawa et al., Optical freezing of charge motion in an organic conductor, Nat. Commun. 5, 5528 (2014), doi:10.1038/ncomms6528.
[15] T. Oka and H. Aoki, Photovoltaic Hall effect in graphene, Phys. Rev. B 79, 081406 (2009),
doi:10.1103/PhysRevB.79.081406.
[16] T. Kitagawa, E. Berg, M. Rudner and E. Demler, Topological characterization of periodically driven quantum systems, Phys. Rev. B 82, 235114 (2010),
doi:10.1103/PhysRevB.82.235114.
[17] Y. H. Wang, H. Steinberg, P. Jarillo-Herrero and N. Gedik, Observation of FloquetBloch states on the surface of a topological insulator, Science 342, 453 (2013),
doi:10.1126/science.1239834.
46

SciPost Phys. Core 4, 033 (2021)[18] G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif and T. Esslinger,
Experimental realization of the topological Haldane model with ultracold fermions, Nature
515, 237 (2014), doi:10.1038/nature13915.
[19] J. W. McIver, B. Schulte, F.-U. Stein, T. Matsuyama, G. Jotzu, G. Meier and A. Cavalleri, Light-induced anomalous Hall effect in graphene, Nature Physics 16, 38 (2020),
doi:10.1038/s41567-019-0698-y.
[20] D. V. Else, B. Bauer and C. Nayak, Floquet time crystals, Phys. Rev. Lett. 117, 090402(2016), doi:10.1103/PhysRevLett.117.090402.
[21] S. Choi et al., Observation of discrete time-crystalline order in a disordered dipolar manybody system, Nature 543, 221 (2017), doi:10.1038/nature21426.
[22] J. Zhang et al., Observation of a discrete time crystal, Nature 543, 217 (2017),
doi:10.1038/nature21413.
[23] L. P. Pitaevskii, Electric forces in a transparent dispersive medium, J. Exp. Theor. Phys.
12, 1008 (1961).
[24] P. S. Pershan, J. P. van der Ziel and L. D. Malmstrom, Theoretical discussion of the inverse
Faraday effect, Raman scattering, and related phenomena, Phys. Rev. 143, 574 (1966),
doi:10.1103/PhysRev.143.574.
[25] A. Kirilyuk, A. V. Kimel and T. Rasing, Ultrafast optical manipulation of magnetic order,
Rev. Mod. Phys. 82, 2731 (2010), doi:10.1103/RevModPhys.82.2731.
[26] A. Lazarides, A. Das and R. Moessner, Equilibrium states of generic quantum systems subject to periodic driving, Phys. Rev. E 90, 012110 (2014),
doi:10.1103/PhysRevE.90.012110.
[27] L. D’Alessio and M. Rigol, Long-time behavior of isolated periodically driven interacting
lattice systems, Phys. Rev. X 4, 041048 (2014), doi:10.1103/PhysRevX.4.041048.
[28] H. Kim, T. N. Ikeda and D. A. Huse, Testing whether all eigenstates obey
the eigenstate thermalization hypothesis, Phys. Rev. E 90, 052105 (2014),
doi:10.1103/PhysRevE.90.052105.
[29] D. A. Abanin, W. De Roeck and F. Huveneers, Exponentially slow heating inperiodically driven many-body systems, Phys. Rev. Lett. 115, 256803 (2015),
doi:10.1103/PhysRevLett.115.256803.
[30] T. Mori, T. Kuwahara and K. Saito, Rigorous bound on energy absorption and genericrelaxation in periodically driven quantum systems, Phys. Rev. Lett. 116, 120401 (2016),
doi:10.1103/PhysRevLett.116.120401.
[31] T. Kuwahara, T. Mori and K. Saito, Floquet-Magnus theory and generic transient dynamics in periodically driven many-body quantum systems, Ann. Phys. 367, 96 (2016),
doi:10.1016/j.aop.2016.01.012.
[32] D. A. Abanin, W. De Roeck, W. W. Ho and F. Huveneers, Effective Hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems, Phys.
Rev. B 95, 014112 (2017), doi:10.1103/PhysRevB.95.014112.
[33] P. Peng, C. Yin, X. Huang, C. Ramanathan and P. Cappellaro, Floquet prethermalizationin dipolar spin chains, Nat. Phys. 17, 444 (2021), doi:10.1038/s41567-020-01120-z.
47

SciPost Phys. Core 4, 033 (2021)[34] S. Kohler, S. Camalet, M. Strass, J. Lehmann, G.-L. Ingold and P. Hänggi, Charge transport through a molecule driven by a high-frequency field, Chem. Phys. 296, 243 (2004),
doi:10.1016/j.chemphys.2003.09.023.
[35] G. Stefanucci, S. Kurth, A. Rubio and E. K. U. Gross, Time-dependent approachto electron pumping in open quantum systems, Phys. Rev. B 77, 075339 (2008),
doi:10.1103/PhysRevB.77.075339.
[36] N. Tsuji, T. Oka and H. Aoki, Nonequilibrium steady state of photoexcited correlated electrons in the presence of dissipation, Phys. Rev. Lett. 103, 047403 (2009),
doi:10.1103/PhysRevLett.103.047403.
[37] Y. Murakami, N. Tsuji, M. Eckstein and P. Werner, Nonequilibrium steady states and transient dynamics of conventional superconductors under phonon driving, Phys. Rev. B 96,
045125 (2017), doi:10.1103/PhysRevB.96.045125.
[38] H. Dehghani, T. Oka and A. Mitra, Dissipative Floquet topological systems, Phys. Rev. B
90, 195429 (2014), doi:10.1103/PhysRevB.90.195429.
[39] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka and P. Werner, Nonequilibrium
dynamical mean-field theory and its applications, Rev. Mod. Phys. 86, 779 (2014),
doi:10.1103/RevModPhys.86.779.
[40] H.-P. Breuer and F. Petruccione, The theory of open quantum
Oxford University Press,
Oxford,UK, ISBN 9780199213900
doi:10.1093/acprof:oso/9780199213900.001.0001.

systems,(2007),[41] D. W. Hone, R. Ketzmerick and W. Kohn, Time-dependent Floquet theory and absence ofan adiabatic limit, Phys. Rev. A 56, 4045 (1997), doi:10.1103/PhysRevA.56.4045.
[42] W. Kohn, Periodic thermodynamics,
doi:10.1023/A:1010327828445.

J.

Stat.

Phys.

103,

417

(2001),[43] K. Brandner and U. Seifert, Periodic thermodynamics of open quantum systems, Phys.
Rev. E 93, 062134 (2016), doi:10.1103/PhysRevE.93.062134.
[44] H.-J. Schmidt, J. Schnack and M. Holthaus, Periodic thermodynamics of the Rabi modelwith circular polarization for arbitrary spin quantum numbers, Phys. Rev. E 100, 042141(2019), doi:10.1103/PhysRevE.100.042141.
[45] D. Vorberg, W. Wustmann, R. Ketzmerick and A. Eckardt, Generalized Bose-Einstein condensation into multiple states in driven-dissipative systems, Phys. Rev. Lett. 111, 240405(2013), doi:10.1103/PhysRevLett.111.240405.
[46] T. Iadecola, T. Neupert and C. Chamon, Occupation of topological Floquet bands in open
systems, Phys. Rev. B 91, 235133 (2015), doi:10.1103/PhysRevB.91.235133.
[47] K. I. Seetharam, C.-E. Bardyn, N. H. Lindner, M. S. Rudner and G. Refael, Controlled
population of Floquet-Bloch states via coupling to Bose and Fermi baths, Phys. Rev. X 5,041050 (2015), doi:10.1103/PhysRevX.5.041050.
[48] S. Kohler, T. Dittrich and P. Hänggi, Floquet-Markovian description of the parametrically driven, dissipative harmonic quantum oscillator, Phys. Rev. E 55, 300 (1997),
doi:10.1103/PhysRevE.55.300.

48

SciPost Phys. Core 4, 033 (2021)[49] H.-P. Breuer, W. Huber and F. Petruccione, Quasistationary distributions of dissipative
nonlinear quantum oscillators in strong periodic driving fields, Phys. Rev. E 61, 4883
(2000), doi:10.1103/PhysRevE.61.4883.
[50] C. M. Dai, Z. C. Shi and X. X. Yi, Floquet theorem with open systems and its applications,Phys. Rev. A 93, 032121 (2016), doi:10.1103/PhysRevA.93.032121.
[51] M. Hartmann, D. Poletti, M. Ivanchenko, S. Denisov and P. Hänggi, Asymptotic Floquet
states of open quantum systems: the role of interaction, New J. Phys. 19, 083011 (2017),
doi:10.1088/1367-2630/aa7ceb.
[52] G. Lindblad, On the generators of quantum dynamical semigroups, Commun. Math. Phys.
48, 119 (1976), doi:10.1007/BF01608499.
[53] V. Gorini, A. Kossakowski and E. C. G. Sudarshan, Completely positive dynamical semigroups of N-level systems, J. Math. Phys. 17, 821 (1976), doi:10.1063/1.522979.
[54] C. P. Slichter, Principles of magnetic resonance 3rd edn., Springer, ISBN 978-3540501572(1996).
[55] M. Grifoni and P. Hänggi, Driven quantum tunneling, Phys. Rep. 304, 229 (1998),
doi:10.1016/S0370-1573(98)00022-2.
[56] S. Scopa, G. T. Landi, A. Hammoumi and D. Karevski, Exact solution of timedependent Lindblad equations with closed algebras, Phys. Rev. A 99, 022105 (2019),
doi:10.1103/PhysRevA.99.022105.
[57] T. N. Ikeda and M. Sato, General description for nonequilibrium steady states inperiodically driven dissipative quantum systems, Sci. Adv. 6, eabb4019 (2020),
doi:10.1126/sciadv.abb4019.
[58] A. G. Redfield, On the theory of relaxation processes, IBM J. Res. Dev. 1, 19 (1957),
doi:10.1147/rd.11.0019.
[59] A. G. Redfield, The theory of relaxation processes, Adv. Magn. Opt. Reson. 1, 1 (1965),
doi:10.1016/B978-1-4832-3114-3.50007-6.
[60] K. Chinzei and T. N. Ikeda, Time crystals protected by Floquet dynamical symmetry in Hubbard models, Phys. Rev. Lett. 125, 060601 (2020),
doi:10.1103/PhysRevLett.125.060601.
[61] P. Menczel and K. Brandner, Limit cycles in periodically driven open quantum systems, J.
Phys. A: Math. Theor. 52, 43LT01 (2019), doi:10.1088/1751-8121/ab435a.
[62] E. S. Mananga and T. Charpentier, Introduction of the Floquet-Magnus expansion in solidstate nuclear magnetic resonance spectroscopy, J. Chem. Phys. 135, 044109 (2011),
doi:10.1063/1.3610943.
[63] F. Haddadfarshi, J. Cui and F. Mintert, Completely positive approximate solutions of driven open quantum systems, Phys. Rev. Lett. 114, 130402 (2015),
doi:10.1103/PhysRevLett.114.130402.
[64] M. M. Wolf, J. Eisert, T. S. Cubitt and J. I. Cirac, Assessing non-Markovian quantum
dynamics, Phys. Rev. Lett. 101, 150402 (2008), doi:10.1103/PhysRevLett.101.150402.
[65] A. Schnell, A. Eckardt and S. Denisov, Is there a Floquet Lindbladian?, Phys. Rev. B 101,
100301 (2020), doi:10.1103/PhysRevB.101.100301.
49

SciPost Phys. Core 4, 033 (2021)[66] R. A. Horn and C. R. Johnson, Topics in matrix analysis, Cambridge University Press,
Cambridge, UK, ISBN 9780521305877 (2010), doi:10.1017/CBO9780511840371.
[67] K. Mizuta, K. Takasan and N. Kawakami, Breakdown of Markovianity by interactions
in stroboscopic Floquet-Lindblad dynamics under high-frequency drive, Phys. Rev. A 103,
L020202 (2021), doi:10.1103/PhysRevA.103.L020202.
[68] A. Schnell, S. Denisov and A. Eckardt, High-frequency expansions for time-periodic Lindblad generators, Phys. Rev. B 104, 165414 (2021), doi:10.1103/PhysRevB.104.165414.
[69] C. W. Gardiner and P. Zoller, Quantum noise: A handbook of Markovian and nonMarkovian quantum stochastic methods with applications to quantum optics, Springer
Berlin, Germany, ISBN 9783662041055 (2004), doi:10.1007/978-3-662-04103-1[70] Z. Jin, W.-J. Gong, A.-D. Zhu, S. Zhang, Y. Qi and S.-L. Su, Dissipative preparation of
qutrit entanglement via periodically modulated Rydberg double antiblockade, Opt. Express
29, 10117 (2021), doi:10.1364/OE.419568.
[71] H. Walther, B. T. H. Varcoe, B.-G. Englert and T. Becker, Cavity quantum electrodynamics,
Rep. Prog. Phys. 69, 1325 (2006), doi:10.1088/0034-4885/69/5/r02.
[72] S. A. Sato et al., Microscopic theory for the light-induced anomalous Hall effect ingraphene, Phys. Rev. B 99, 214302 (2019), doi:10.1103/PhysRevB.99.214302.
[73] K. Yamamoto, Y. Ashida and N. Kawakami, Rectification in nonequilibrium
steady states of open many-body systems, Phys. Rev. Res. 2, 043343 (2020),
doi:10.1103/PhysRevResearch.2.043343.
[74] T. Prosen and I. Pižorn, Quantum phase transition in a far-from-equilibrium
steady state of an X Y spin chain, Phys. Rev. Lett. 101, 105701 (2008),
doi:10.1103/PhysRevLett.101.105701.
[75] T. N. Ikeda and M. Sato, High-harmonic generation by electric polarization, spin current, and magnetization, Phys. Rev. B 100, 214424 (2019),
doi:10.1103/PhysRevB.100.214424.
[76] M. Sato and Y. Morisaku, Two-photon driven magnon-pair resonance asa signature of spin-nematic order, Phys. Rev. B 102, 060401 (2020),
doi:10.1103/PhysRevB.102.060401.
[77] M. Kanega, T. N. Ikeda and M. Sato, Linear and nonlinear optical responses in Kitaev spin
liquids, Phys. Rev. Res. 3, L032024 (2021), doi:10.1103/PhysRevResearch.3.L032024.
[78] H. Haug and A.-P. Jauho, Quantum kinetics in transport and optics of semiconductors, Springer Berlin Heidelberg, Berlin, Heidelberg, ISBN 9783540735618 (2008),
doi:10.1007/978-3-540-73564-9.
[79] H. Haug and S. W. Koch, Quantum theory of the optical and electronic properties of semiconductors 5th edn., World Scientific (2009), doi:10.1142/7184.
[80] T. Prosen, Third quantization: a general method to solve master equations for
quadratic open Fermi systems, New J. Phys. 10, 043026 (2008), doi:10.1088/13672630/10/4/043026.
[81] M. Thakurathi, A. A. Patel, D. Sen and A. Dutta, Floquet generation of Majorana end modes and topological invariants, Phys. Rev. B 88, 155133 (2013),
doi:10.1103/PhysRevB.88.155133.
50

SciPost Phys. Core 4, 033 (2021)[82] P. Molignini, E. van Nieuwenburg and R. Chitra, Sensing Floquet-Majorana fermions via
heat transfer, Phys. Rev. B 96, 125144 (2017), doi:10.1103/PhysRevB.96.125144.
[83] T. Prosen and E. Ilievski, Nonequilibrium phase transition in a periodically driven X Y spin
chain, Phys. Rev. Lett. 107, 060403 (2011), doi:10.1103/PhysRevLett.107.060403.
[84] L. Rondin, J.-P. Tetienne, T. Hingant, J.-F. Roch, P. Maletinsky and V. Jacques, Magnetometry with nitrogen-vacancy defects in diamond, Rep. Prog. Phys. 77, 056503 (2014),
doi:10.1088/0034-4885/77/5/056503.
[85] R. Blümel, A. Buchleitner, R. Graham, L. Sirko, U. Smilansky and H. Walther, Dynamical
localization in the microwave interaction of Rydberg atoms: The influence of noise, Phys.
Rev. A 44, 4521 (1991), doi:10.1103/PhysRevA.44.4521.
[86] F. Nathan and M. S. Rudner, Universal Lindblad equation for open quantum systems, Phys.
Rev. B 102, 115109 (2020), doi:10.1103/PhysRevB.102.115109.
[87] E. Mozgunov and D. Lidar, Completely positive master equation for arbitrary driving andsmall level spacing, Quantum 4, 227 (2020), doi:10.22331/q-2020-02-06-227.
[88] T. Becker, L.-N. Wu and A. Eckardt, Lindbladian approximation beyond ultraweak coupling, Phys. Rev. E 104, 014110 (2021), doi:10.1103/PhysRevE.104.014110.
[89] T. Shirai, J. Thingna, T. Mori, S. Denisov, P. Hänggi and S. Miyashita, Effective
Floquet-Gibbs states for dissipative quantum systems, New J. Phys. 18, 53008 (2016),
doi:10.1088/1367-2630/18/5/053008.
[90] O. R. Diermann, H.-J. Schmidt, J. Schnack and M. Holthaus, Environmentcontrolled Floquet-state paramagnetism, Phys. Rev. Res. 2, 023293 (2020),
doi:10.1103/PhysRevResearch.2.023293.
[91] N. Shannon, T. Momoi and P. Sindzingre, Nematic order in square lattice frustrated ferromagnets, Phys. Rev. Lett. 96, 027213 (2006), doi:10.1103/PhysRevLett.96.027213.
[92] M. E. Zhitomirsky and H. Tsunetsugu, Magnon pairing in quantum spin nematic, Europhys. Lett. 92, 37001 (2010), doi:10.1209/0295-5075/92/37001.
[93] M. Sato, T. Hikihara and T. Momoi, Spin-nematic and spin-density-wave orders in spatiallyanisotropic frustrated magnets in a magnetic field, Phys. Rev. Lett. 110, 077206 (2013),
doi:10.1103/PhysRevLett.110.077206.
[94] K. Penc and A. M. Läuchli, Spin nematic phases in quantum spin systems, in Introduction to frustrated magnetism, Springer Berlin Heidelberg, Berlin, Heidelberg, ISBN9783642105883, 331 (2011), doi:10.1007/978-3-642-10589-0_13.
[95] M. Sato, Floquet theory and ultrafast control of magnetism, in Chirality, magnetism andmagnetoelectricity, Springer, Cham, Switzerland, ISBN 9783030628444, 265 (2021),
doi:10.1007/978-3-030-62844-4_11.
[96] J. P. van der Ziel, P. S. Pershan and L. D. Malmstrom, Optically-induced magnetization resulting from the inverse Faraday effect, Phys. Rev. Lett. 15, 190 (1965),
doi:10.1103/PhysRevLett.15.190.
[97] T. Makino, F. Liu, T. Yamasaki, Y. Kozuka, K. Ueno, A. Tsukazaki, T. Fukumura, Y. Kongand M. Kawasaki, Ultrafast optical control of magnetization in EuO thin films, Phys. Rev.
B 86, 064403 (2012), doi:10.1103/PhysRevB.86.064403.
51

SciPost Phys. Core 4, 033 (2021)[98] F. Hansteen, A. Kimel, A. Kirilyuk and T. Rasing, Nonthermal ultrafast opticalcontrol of the magnetization in garnet films, Phys. Rev. B 73, 014421 (2006),
doi:10.1103/PhysRevB.73.014421.
[99] T. Satoh et al., Spin oscillations in antiferromagnetic NiO triggered by circularly polarized
light, Phys. Rev. Lett. 105, 077402 (2010), doi:10.1103/PhysRevLett.105.077402.
[100] A. V. Kimel, A. Kirilyuk, P. A. Usachev, R. V. Pisarev, A. M. Balbashov and T. Rasing,
Ultrafast non-thermal control of magnetization by instantaneous photomagnetic pulses,
Nature 435, 655 (2005), doi:10.1038/nature03564.
[101] S. Takayoshi, H. Aoki and T. Oka, Magnetization and phase transition induced bycircularly polarized laser in quantum magnets, Phys. Rev. B 90, 085150 (2014),
doi:10.1103/PhysRevB.90.085150.
[102] S. Takayoshi, M. Sato and T. Oka, Laser-induced magnetization curve, Phys. Rev. B 90,
214413 (2014), doi:10.1103/PhysRevB.90.214413.
[103] M. Sato, S. Takayoshi and T. Oka, Laser-driven multiferroics and ultrafast spin current generation, Phys. Rev. Lett. 117, 147202 (2016),
doi:10.1103/PhysRevLett.117.147202.
[104] M. Takahashi, Thermodynamics of one-dimensional solvable models, Cambridge University Press, Cambridge, UK, ISBN 9780521551434 (1999),
doi:10.1017/CBO9780511524332.
[105] T. Giamarchi, Quantum physics in one dimension, Oxford University Press, Oxford, UK
ISBN 9780198525004 (2003), doi:10.1093/acprof:oso/9780198525004.001.0001.
[106] B. Buča and T. Prosen, A note on symmetry reductions of the Lindblad equation: transport
in constrained open spin chains, New J. Phys. 14, 073007 (2012), doi:10.1088/13672630/14/7/073007.
[107] B. Buča, J. Tindall and D. Jaksch, Non-stationary coherent quantum manybody dynamics through dissipation, Nat. Commun. 10, 1730 (2019),
doi:10.1038/s41467-019-09757-y.
[108] J. Tindall, C. Sánchez Muñoz, B. Buča and D. Jaksch, Quantum synchronisation enabled by dynamical symmetries and dissipation, New J. Phys. 22, 013026 (2020),
doi:10.1088/1367-2630/ab60f5.

52