## Abstract

Laser technology has developed and accelerated photo-induced nonequilibrium physics, from both the scientific and engineering viewpoints. Floquet engineering, i.e., controlling material properties and functionalities by time-periodic drives, is at the forefront of quantum physics of light-matter interaction. However, it is limited to ideal dissipationless systems. Extending Floquet engineering to various materials requires understanding of the quantum states emerging in a balance of the periodic drive and energy dissipation. Here, we derive a general description for nonequilibrium steady states (NESSs) in periodically driven dissipative systems by focusing on systems under high-frequency drive and time-independent Lindblad-type dissipation. Our formula correctly describes the time average, fluctuation, and symmetry properties of the NESS, and can be computed efficiently in numerical calculations. This approach will play fundamental roles in Floquet engineering in a broad class of dissipative quantum systems from atoms and molecules to mesoscopic systems, and condensed matter.

## INTRODUCTION

State-of-the-art laser technology has opened new research fields in physics: the Floquet science and engineering (*1*–*3*). The main focus of these fields is the nonequilibrium states driven periodically by external fields, e.g., intense laser fields. Physical properties of the nonequilibrium states are mainly understood by the so-called effective Hamiltonian, which reflects the periodic driving, according to the Floquet theorem (*4*) and the ensuing theoretical developments (*5*–*8*). Conversely, designing a suitable driving protocol, one can engineer the effective Hamiltonian, which enables us to have desirable properties and functionalities of physical systems. Various exotic states and useful manipulation of matter have been theoretically proposed, and some of them have been experimentally realized: Floquet topological states (*9*) in solids (*10*), ultracold atomic gases (*11*), and in photonic wave guides (*12*), Floquet time crystals (*13*) in nitrogen-vacancy (NV) centers (*14*) and trapped ions (*15*), and control of quantum magnets (*16*) and their interactions (*17*).

However, these Floquet-theoretical predictions based on the effective Hamiltonian are quantitatively reliable only in ultraclean materials or well-designed artificial systems, where dissipation is negligible. For the Floquet science and engineering in a variety of materials, it is indispensable to understand the nonequilibrium steady state (NESS), which emerges in a balance of the energy injection by the periodic driving and the energy dissipation (*18*–*21*). For individual systems, by considering specific sources of dissipation, i.e., system-bath couplings, one can calculate physical quantities in the NESS and predict interesting phenomena such as Floquet topological insulators (*22*, *23*), periodic thermodynamics (*24*), dynamical localization (*25*), and generalized Bose-Einstein condensation (*26*). In this research direction, the Floquet-Green function approach has developed and enables us to calculate various physical effects dependent on the type of the system-bath coupling (*27*–*29*). Another research direction, which we address here, is to seek for a universal characterization for the NESS. We could imagine that there exists a simple and general expression for the NESS when the dissipation is weak and featureless. An attempt is to conjecture that the NESS is generally described by the Floquet-Gibbs state (FGS) (*21*, *30*, *31*), i.e., the Gibbs state with the effective Hamiltonian, but the conditions for the FGS being realized have shown quite restrictive (*30*). Hence, despite its importance, the general formula for the NESS has been still an elusive problem.

Here, in exchange for restricting ourselves to the high-frequency drivings, we deal with generic systems and driving protocols, obtaining simple and general formulas for the NESS (Eqs. 9 to 11 below). We derive these formulas by applying the high-frequency expansion technique, which has been recently developed (*5*, *6*, *32*–*34*), to the Lindblad equation with periodic Hamiltonians. As exemplified in an effective model for the NV center in diamonds (*35*), our formulas correctly describe both the time average and fluctuation of the NESS at the leading order of ω^{−1} (ω denotes the driving frequency). These formulas also capture nontrivial behaviors of physical quantities due to the dynamical-symmetry breaking that cannot be described by the effective Hamiltonian or the FGS and will thereby play critical roles in the Floquet science and engineering in dissipative quantum systems.

## RESULTS

### Floquet-Lindblad equation with detailed balance condition

We begin by considering a quantum system defined on an *N*-dimensional Hilbert space. This system can be single body or many body as long as it satisfies the requirements that will be described below. We let *H*_{0} denote the time-independent Hamiltonian, which describes our system in the absence of driving. The eigenenergies and eigenstates of *H*_{0} are denoted by *E*_{1} < *E*_{2} < ⋯ < *E _{N}* (the generalization to degenerate

*H*

_{0}is formulated in the Supplementary Materials). The effect of the driving is represented by a time-dependent part

*H*

_{ext}(

*t*) of the total Hamiltonian

*T*:

*H*

_{ext}(

*t*+

*T*) =

*H*

_{ext}(

*t*) and hence

*H*(

*t*+

*T*) =

*H*(

*t*). Without loss of generality, the decomposition (Eq. 1) is defined so that the time average of

*H*

_{ext}(

*t*) vanishes,

*H*

_{ext}(

*t*) can be written as

To study driven dissipative systems, we consider the density operator ρ(*t*), whose dynamics is described by the Floquet-Lindblad equation (*36*–*39*) (we set ℏ = 1 throughout this paper)*L _{ij}* ≡ ∣

*E*〉〈

_{i}*E*∣ is the time-independent Lindblad operator describing the transition from the

_{j}*j*-th to the

*i*-th eigenstates of the undriven Hamiltonian

*H*

_{0}. Each Lindblad operator

*L*represents a decay (excitation) process for

_{ij}*i*<

*j*(

*i*>

*j*). The real number Γ

*( ≥ 0) denotes the rate for the corresponding process, and we set Γ*

_{ij}*= 0 for each*

_{ii}*i*. The transition rates Γ

*must be small enough for the Floquet-Lindblad equation to be valid (see Discussion). Note that Eq. 3 is trace-preserving*

_{ij}*d*tr[ρ(

*t*)]/

*dt*= 0, and thus, we use the normalization tr[ρ(

*t*)] = 1. We assume that the transition rates Γ

*satisfy the detailed balance condition*

_{ij}*is a nonnegative irreducible matrix (*

_{ij}*40*). These assumptions ensure that, without driving, the system goes, irrespective of the initial state, to the thermal equilibrium state, or the canonical ensemble ρ

_{can}=

*e*

^{−βH0}/

*Z*of

*H*

_{0}with

*Z*= tr(

*e*

^{−βH0}). We note that the Lindblad operators

*L*may depend on the driving in general if we consider more microscopic theories of dissipation (

_{ij}*21*). However, we neglect this dependence in this work for simplicity.

### General formulas for NESS

The key idea to obtain the NESS is the high-frequency expansion for the Floquet-Lindblad equation (*34*). Among several formulations, we adopt the van Vleck perturbation theory (*5*, *6*), which leads to the following propagation for ρ(*t*) (see Materials and Methods): _{eff} is represented by the effective Hamiltonian*H*_{eff} = *H*_{0} + ω^{−1}∑_{n > 0}[*H*_{−n}, *H _{n}*]/

*n*+

*O*(ω

^{−2}). The time-dependent part

*t*′ = 0, having

To obtain the asymptotic behavior of ρ(*t*), we focus on the first part *d*ρ′(*t*)/*dt* = ℒ_{eff} ρ′(*t*) from the initial state *t*) approaches, irrespective of the initial state, the unique state _{ness}(*t*)_{MM}(*t*) and σ_{FE} are *O*(ω^{−1}), and we call σ_{MM}(*t*) and σ_{FE} the micromotion and Floquet-engineering parts, respectively. Equation 9 is our second main result, which we prove in Materials and Methods. Its generalization in the absence of the detailed balance condition is outlined in Discussion. Note that tr[ρ_{ness}(*t*)] = 1 is satisfied, at least, up to this order because both σ_{MM}(*t*) and σ_{FE} are traceless, as will be evident below.

The micromotion part σ_{MM}(*t*) is defined by_{MM}(*t*). First, this part is periodic in time σ_{MM}(*t* + *T*) = σ_{MM}(*t*) and contributes to oscillations of physical observables. Second, σ_{MM}(*t*) does not contribute to the time averages of physical observables for one period of oscillation. For an observable *A*, we have

The Floquet-engineering part σ_{FE} in Eq. 9 is independent of time and given by*E _{k}*∣σ

_{FE}∣

*E*〉 = 0 for all

_{k}*k*, where Δ

*H*

_{eff}≡

*H*

_{eff}−

*H*

_{0}=

*O*(ω

^{−1}),

*H*

_{0}). The reason why we call σ

_{FE}the Floquet-engineering part is that it describes how the effective Hamiltonian changes physical observables from their values in thermal equilibrium. In contrast to the micromotion part, the Floquet-engineering part contributes to the time-averaged quantities. As we will show below, Eq. 11 is regular in the weak dissipation limit γ

*→ 0, where ρ*

_{ij}_{ness}(

*t*) becomes independent of γ

_{ij}and coincides with the canonical Floquet steady state (CFSS) that we define in Eq. 14 below.

Equation 11 serves as the foundation for the Floquet engineering in dissipative quantum systems. Let us imagine, for example, that an observable *A* has zero expectation value at thermal equilibrium, tr(ρ_{can}*A*) = 0, but nonzero value for the NESS, tr(σ_{FE}*A*) ≠ 0. This situation means that one can implement an appropriate periodic driving *H*_{ext}(*t*) and hence Δ*H*_{eff}, thereby activating the observable *A*. Upon this engineering, Eq. 11 tells us how much activation is possible for observables of interest. We will see some examples below.

In addition to their generality, our formulas (Eqs. 9 to 11) are extremely efficient in practical calculations of the NESS. In the straightforward numerical calculation, one integrates the Floquet-Lindblad equation (Eq. 3) with a sufficiently short time step until the system reaches the NESS. In contrast, our formulas enable us to evaluate the NESS at an arbitrary time *t* without numerical integration once we have the energy eigenstates {∣*E _{k}*〉} of the time-independent Hamiltonian

*H*

_{0}. This difference of efficiency becomes more evident when the Hilbert-space dimension

*N*is large. In the straightforward calculation, the density matrix ρ(

*t*) is commonly treated as an

*N*

^{2}-dimensional vector and the superoperator ℒ

*as an*

_{t}*N*

^{2}×

*N*

^{2}matrix. Thus, the computational complexity for each time step is

*O*(

*N*

^{4}) in general. On the other hand, the complexity for our formula is one order smaller and given by

*O*(

*N*

^{3}), which derives from the exact diagonalization of

*H*

_{0}. Thus, our formulas enable us to evaluate the NESS for larger Hilbert-space dimensions that occur, for example, in quantum many-body systems. For special cases where ℒ

*is a sparse matrix and has only*

_{t}*O*(

*N*

^{0}) nonzero elements, the computational complexity for one time step is

*O*(

*N*

^{2}). Nevertheless, even for these cases, we need many time steps typically larger than

*N*in obtaining accurate results, and hence, our formulas require less computational complexity.

### Numerical verification with the NV center in diamonds

By taking a single spin with *S* = 1, we demonstrate how our formula (Eq. 9) works in the quantum dynamics described by the Floquet-Lindblad equation (Eq. 3). We consider an effective Hamiltonian for the NV center in diamonds (*35*)*B _{s}* is the static Zeeman field,

*35*), we set

Typical time evolutions *A*(*t*) ≡ tr[ρ(*t*)*A*] are shown for two representative observables *A* = *S _{z}* and {

*S*,

_{x}*S*} ( =

_{y}*S*+

_{x}S_{y}*S*) in Fig. 1 (B and C). In this calculation, we take the thermal state for the time-independent part of

_{y}S_{x}*H*

_{NV}at

*t*= 0 and let the state evolve according to the Floquet-Lindblad equation (Eq. 3) with

*H*(

*t*) being

*H*

_{NV}(

*t*). The static Zeeman field is

*B*= 0.3, and the driving parameters are

_{s}*B*= 0.1 and ω = 1. As for the Lindblad operators, we take Γ

_{d}*according to the heatbath method as Γ*

_{ij}*= γ*

_{ij}*i*≠

*j*, with rate constant γ = 0.2 and inverse temperature β = 3, and Γ

*= 0 for all*

_{ii}*i*’s. Figure 1 (B and C) shows that after a sufficiently long time

*t*≫ γ

^{−1}, the system reaches the NESS, in which the observables oscillate with period

*T*= 2π/ω. In particular, the observable

*A*= {

*S*,

_{x}*S*} is initially zero for a symmetry reason (e.g.,

_{y}*S*→ −

_{y}*S*) but becomes nonzero on time average. Namely, this observable is engineered by the periodic drive

_{y}To test our formula (Eq. 9) quantitatively, we first focus on the one-cycle average *t* ≫ γ^{−1}. For this test, we also define _{HF}(*s*) is the density matrix calculated from Eqs. 9 to 11 based on the high-frequency expansion. In Fig. 1 (D and E), we compare the one-cycle averages _{MM}(*t*) does not contribute to the one-cycle averages]. At high-frequency ω ≳ 10, the difference *A*^{HF}(ω) against ω for *A* = *S _{z}* and {

*S*,

_{x}*S*}, respectively. We stress that the difference Δ

_{y}*A*

^{HF}(ω) decreases more rapidly than

*O*(ω

^{−1}) [

*O*(ω

^{−2}) for

*S*and

_{z}*O*(ω

^{−3}) for {

*S*,

_{x}*S*}]. This means that Eqs. 9 and 11 perfectly describe the actual NESS at the level of

_{y}*O*(ω

^{−1}). As shown in the Supplementary Materials, Δ

*A*

^{HF}(ω) =

*O*(ω

^{−2}) holds true not only for these two observables but also for all the other observables. Therefore, we have verified Eq. 9 apart from the micromotion part.

For the complementary test of our formula (Eq. 9), we consider the one-cycle standard deviation (SD) _{MM}(*t*). Because Σ* _{A}*(ω) is an

*O*(ω

^{−1}) quantity in general, the accuracy of our formula is verified if the difference

*O*(ω

^{−2}). This criterion is satisfied, as shown in Fig. 2 (C and D) for

*A*=

*S*and {

_{z}*S*,

_{x}*S*}, respectively. We remark that our formula leads to Σ

_{y}*(ω) = 0 at*

_{A}*O*(ω

^{−1}) for these observables, which can be analytically shown by noticing

*H*

_{±1}∝

*S*

_{±}=

*S*±

_{x}*iS*. Thus, the plotted data correspond to Σ

_{y}*(ω) itself for the actual dynamics, and ΔΣ*

_{A}*(ω) could be reduced by dealing with the higher-order terms in Eq. 9. In any case, the fact that ΔΣ*

_{A}*(ω) is*

_{A}*O*(ω

^{−2}) justifies Eqs. 9 and 10.

### Comparison with the FGS

Let us make comparisons with the FGS, which has been a candidate for the NESS description in periodically driven dissipative quantum systems (*21*, *30*, *31*). To define the FGS, we introduce the Floquet state ∣*u _{i}*(

*t*)〉 and its quasienergy ϵ

*. According to the Floquet theorem, the time-dependent Schrödinger equation*

_{i}*(*

_{i}*t*)〉 =

*u*(

_{i}*t*)〉 (

*i*= 1,2, …,

*N*), with periodicity ∣

*u*(

_{i}*t*+

*T*)〉 = ∣

*u*(

_{i}*t*)〉 and 〈

*u*(

_{i}*t*)∣

*u*(

_{j}*t*)〉 = δ

*. In terms of the Floquet states, the FGS is defined by*

_{ij}*U*(

*T*) correspond to ∣

*u*(0)〉 and

_{i}*and ∣*

_{i}*u*(

_{i}*t*)〉. Note that the Floquet states and quasienergies thus obtained are exact and involve all-order contributions in 1/ω.

Quantitative comparisons between the actual dynamics and the FGS are shown in Fig. 2. For the FGS, we define the one-cycle average *A* = *S _{z}* and {

*S*,

_{x}*S*}. Remarkably, the difference is

_{y}*O*(ω

^{−1}), meaning that the FGS cannot reproduce the leading-order contribution of the one-cycle average in contrast to our Eq. 9. As for the micromotion amplitude, or the one-cycle SD Σ

*(ω), Fig. 2 (C and D) shows that the FGS reproduces the actual values better than our Eqs. 9 and 10. This is, in part, because the FGS involve all-order contributions in 1/ω, while our formulas are the leading-order approximation. As stated above, our formula could be improved order by order when we start over from our first main result (Eq. 8).*

_{A}A weak point of the FGS is highlighted in Fig. 1E, in which the FGS gives *V*: *VS _{y}V*

^{†}= −

*S*and

_{y}*VS*

_{α}V^{†}=

*S*(α =

_{α}*x*and

*z*). Then, we notice the dynamical symmetry

*VH*

_{NV}(

*T*−

*t*)

*V*

^{†}=

*H*

_{NV}(

*t*), which implies that

*. Assuming that quasienergies are not degenerate as in our examples, we have that*

_{i}*u*(

_{i}*t*)〉 are equivalent up to an overall phase shift. Because of

*VAV*

^{†}= −

*A*with

*A*= {

*S*,

_{x}*S*}, the one-cycle averages of

_{y}*A*calculated for

*u*(

_{i}*t*)〉 differ by their signs, meaning that the one-cycle average vanishes. Note that similar arguments apply to other observables, satisfying

*VAV*

^{†}= −

*A*.

We remark that the dissipation can break such an antiunitary dynamical symmetry, and this is the origin of the nonzero one-cycle average of * _{ij}* → 0 in Eq. 11. In other words, the NESS in dissipative systems shows richer properties inferred only from the effective Hamiltonian itself. Our Eq. 9 well describes these properties unlike the FGS (Eq. 13), which incorporates no information about the dissipation, or the Lindblad operators.

One might be interested in an approximate description of the NESS independent of the details of γ* _{ij}* for weak dissipation and expect that the FGS serves such a description. This is not true, at least within our formulation of periodically driven dissipative systems described by Eqs. 3 and 4. Instead, the actual NESS coincides with yet another state, which we name the canonical Floquet steady state (CFSS), defined by

*with the real energy*

_{i}*E*. One can analytically show that the CFSS at high frequency coincides with our formulas (Eqs. 9 to 11) in the limit of γ

_{i}*→ 0 (see the Supplementary Materials for details).*

_{ij}## DISCUSSION

We have derived and verified the simple and general formulas (Eqs. 9 to 11) describing the NESS in dissipative quantum systems under high-frequency periodic drive. We have also exemplified the breaking of an antiunitary dynamical symmetry and the possibility of the Floquet engineering in NV centers in diamonds. Being quite general, our formulas should apply to various quantum systems such as atoms and molecules, trapped ions, condensed matters, and so on, and play the fundamental role in understanding and engineering unusual nonequilibrium states in these systems.

The parameter region in which our formulas are valid is depicted in Fig. 3. On the basis of the high-frequency expansion, our formulas are valid when the driving frequency ω (more precisely, the photon energy ℏω) is greater than the energy scales of the system and the system-drive coupling. We note, however, that our formulas hold true for any strength of dissipation Γ* _{ij}* (or γ

*) within the Floquet-Lindblad equation (Eq. 3). As we have shown, our formulas reduce to the CFSS rather than the FGS when the dissipation strength is smaller than the energy gap, i.e., the nonzero minimum difference between eigenenergies (our formulas are generalized for the degenerate Hamiltonian in the Supplementary Materials).*

_{ij}One should note that the Floquet-Lindblad equation (Eq. 3) becomes invalid when Γ* _{ij}* is too large. The Lindblad-type dissipation is derived from several approximations such as the Born-Markov approximation (

*39*). These approximations require the condition that the relaxation time ~1/Γ

*is longer than the time scale of the system’s dynamics and the bath correlation time. Namely, the dissipation rate Γ*

_{ij}*should be smaller than the other relevant energy scales. Note that the high-frequency driving does not break this condition, while lower frequencies may be problematic.*

_{ij}We remark a further generalization of our results Eqs. 9 to 11. Although we have assumed the detailed balance condition (Eq. 5), this condition can be removed as long as the transition-rate matrix Γ* _{ij}* is irreducible. In this case, the solution of the Lindblad equation without driving is not the canonical ensemble ρ

_{can}but another state

It remains an open question to find a simple and general formula for the NESS at lower frequency. The applicability of the CFSS in many-body systems is also a nontrivial issue because the energy gap can be very small in those systems. Addressing these issues will lead us to the complete understanding of the NESS in dissipative Floquet systems.

## MATERIALS AND METHODS

### High-frequency expansion for the Floquet-Lindblad equation

The Floquet-Lindblad equation that we discuss in this work is symbolically represented as ∂* _{t}*ρ(

*t*) = ℒ(

*t*)ρ(

*t*), where the time-dependent Liouvillian ℒ(

*t*) is defined by

*t*) = Σ

*ℒ*

_{m}

_{m}e^{−imtω}. Because the Lindblad operators

*L*are time independent in this work, each Fourier component is given as follows

_{ij}The high-frequency expansion in terms of the van Vleck approach makes the following ansatz_{eff} is time independent. This ansatz satisfies one determining equation _{eff}, and what determines these two is

To obtain the determining equations for _{eff}, we substitute Eq. 18 into *B _{k}* denotes the Bernoulli number (

*B*

_{0}= 1,

*B*

_{1}= − 1/2,

*B*

_{2}= 1/6, …).

We determine _{eff} from Eq. 19 by the series expansions*t*) and *k* for *6*) for the case of unitary dynamics]. With some straightforward calculations (see the Supplementary Materials), we obtain the following low-order solutions*O*(ω^{−2}), and *A*, [*B*, *C*]] + [*B*, [*C*, *A*]] + [*C*, [*A*, *B*]] = 0, we obtain Eq. 6 in the main text.

### Solution for the NESS

Here, we derive Eqs. 9 to 11 from Eq. 8 in the main text. For this purpose, we solve ^{−1}. It is convenient to work in the energy eigenbasis and separate the diagonal and off-diagonal parts

First, we consider the off-diagonal elements of both sides of *k* ≠ *l**H*_{eff} ≡ *H*_{eff} − *H*_{0} = *O*(ω^{−1}). Equation 28 is transformed as*E _{k}* −

*E*) −

_{l}*i*γ

*do not vanish because γ*

_{kl}*> 0 is ensured by the nonnegativity and irreducibility of Γ*

_{kl}*. Now, as a working hypothesis, we suppose that the diagonal elements*

_{ij}*O*(ω

^{0}) as verified later. Then, the first term on the right-hand side of Eq. 29 is

*O*(ω

^{−1}) because Δ

*H*

_{eff}=

*O*(ω

^{−1}). Notice that the second term depends only on the off-diagonal elements

^{−1}expansion for the off-diagonal elements

Next, we consider the diagonal elements of both sides of *H*_{eff} and *O*(ω^{−1}). Thus, the diagonal elements *O*(ω^{−1}) by the equation* _{kl}*, we have the unique solution as

*O*(ω

^{−2}). This result means

*O*(ω

^{0}), the working hypothesis introduced above has been verified. By substituting Eq. 33 into Eq. 30, we have the leading-order expression for the off-diagonal elements

_{FE}) = 0 because each of the diagonal elements of σ

_{FE}vanishes.

Now that we have obtained the leading-order expression for *O*(ω^{0}) and *O*(ω^{−1}) and using the Taylor expansion _{MM}(*t*)] = 0, which follows from the cyclic property of trace, and hence tr[ρ(*t*)] = 1, at least up to this order.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/27/eabb4019/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We are grateful to S. Higashikawa and H. Fujita for collaboration on the early stage of this work, and to K. Chinzei, T. Mori, T. Sagawa, T. Shirai, and H. Tsunetsugu for the fruitful discussions.

**Funding:**This work was supported by JSPS KAKENHI grant nos. JP18K13495 (T.N.I.), JP17K05513, and JP20H01830 (M.S.), and by Grant-in-Aid for Scientific Research on Innovative Area, “Physical Properties of Quantum Liquid Crystals” (grant no. 19H05825) (M.S.).

**Author contributions:**M.S. conceived, initiated, and supervised the project. T.N.I. performed calculations and obtained the main results. T.N.I. wrote the manuscript and drew figures, with feedback from M.S. Both authors discussed the results and approved the manuscript.

**Competing interests:**The authors declare that they have no competing interest.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).