## 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

and

${\{\mid {E}_{i}\u3009\}}_{i=1}^{N}$, respectively. For simplicity, we assume that the eigenenergies are not degenerate and *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

(1)We assume that the driving term is periodic with period *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,

. Thus, the Fourier series of *H*_{ext}(*t*) can be written as

(2)

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)

(3)

$$\mathrm{D}[\mathrm{\rho}(t)]\equiv \sum _{i,j}{\mathrm{\Gamma}}_{\mathit{ij}}[{L}_{\mathit{ij}}\mathrm{\rho}(t){L}_{\mathit{ij}}^{\mathrm{\u2020}}-\frac{1}{2}\{{L}_{\mathit{ij}}^{\mathrm{\u2020}}{L}_{\mathit{ij}},\mathrm{\rho}(t)\left\}\right]$$(4)Here, *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}(5)where β is the inverse temperature of the bath coupled to the system. We also assume that the matrix Γ* _{ij}* is a nonnegative irreducible matrix (

*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):

. The time-independent part ℒ_{eff} is represented by the effective Hamiltonian

(6)with *H*_{eff} = *H*_{0} + ω^{−1}∑_{n > 0}[*H*_{−n}, *H _{n}*]/

*n*+

*O*(ω

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

is the so-called micromotion operator periodic in time

$\mathrm{G}(t+T)=\mathrm{G}(t)$and is given by

$\mathrm{G}(t)(\mathrm{\rho})={\mathrm{\omega}}^{-1}{\sum}_{m\ne 0}[{H}_{m},\mathrm{\rho}]{e}^{-\mathit{\text{im}}\mathrm{\omega}t}+O({\mathrm{\omega}}^{-2})$. Without loss of generality, we suppose the initial time to be *t*′ = 0, having

(7)with ρ(0) being our initial state.

To obtain the asymptotic behavior of ρ(*t*), we focus on the first part

. Note that this is the solution of the time-independent Lindblad equation *d*ρ′(*t*)/*dt* = ℒ_{eff} ρ′(*t*) from the initial state

. Under our assumptions on D, ρ′(*t*) approaches, irrespective of the initial state, the unique state

characterized by

${\mathcal{L}}_{\text{eff}}{\mathrm{\rho}}_{\infty}^{\prime}=0$. Thus, we come to the first main result, obtaining the asymptotic behavior

$$\mathrm{\rho}(t)\to {\mathrm{\rho}}_{\text{ness}}(t)={e}^{\mathrm{G}(t)}{\mathrm{\rho}}_{\infty}^{\prime}\text{as}t\to \infty $$(8)Because

$\mathrm{G}(t)=\mathrm{G}(t+T)$, this NESS is also periodic in time. Focusing on the leading-order contribution, we have a simple explicit formula for ρ_{ness}(*t*)

(9)in which both σ_{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

(10)We have named it after the following two properties of σ_{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

(11)and 〈*E _{k}*∣σ

_{FE}∣

*E*〉 = 0 for all

_{k}*k*, where Δ

*H*

_{eff}≡

*H*

_{eff}−

*H*

_{0}=

*O*(ω

^{−1}),

is the Boltzmann weight, and

${\mathrm{\gamma}}_{\mathit{kl}}\equiv {\displaystyle \sum _{i}}({\mathrm{\Gamma}}_{\mathit{ik}}+{\mathrm{\Gamma}}_{\mathit{il}})/2$ represents the symmetric transition-rate matrix (see the Supplementary Materials for the generalization to degenerate *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 γ* _{ij}* → 0, where ρ

_{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*)

(12)where *B _{s}* is the static Zeeman field,

and

${\mathrm{N}}_{\mathrm{xy}}$are the coupling constants of magnetic anisotropic terms, and

${H}_{\text{ext}}^{\text{circ}}(t)\equiv -{B}_{d}({S}_{x}\text{cos}\mathrm{\omega}\mathit{t}+{S}_{y}\text{sin}\mathrm{\omega}\mathit{t})$represents the coupling to the circularly polarized ac magnetic field. We note that the energy eigenstates of the time-independent part of Eq. 12 are analytically obtained, and our formula can be computed almost analytically in this model. Because

${\mathrm{N}}_{z}\gg {\mathrm{N}}_{\mathrm{xy}}$ in the NV centers (*35*), we set

and

${\mathrm{N}}_{\mathrm{xy}}=0.05$in our analysis. The ac-field irradiation leads to spin dynamics, as schematically illustrated in Fig. 1A.

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}/(

${\mathit{e}}^{-\mathrm{\beta}{\mathit{E}}_{\mathit{i}}}$+

${\mathit{e}}^{-\mathrm{\beta}{\mathit{E}}_{\mathit{j}}}$) for *i* ≠ *j*, with rate constant γ = 0.2 and inverse temperature β = 3, and Γ* _{ii}* = 0 for all

*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

$\overline{A}(\mathrm{\omega})={T}^{-1}{\int}_{t}^{t+T}\mathit{ds}\text{tr}[\mathrm{\rho}(s)A]$ for *t* ≫ γ^{−1}. For this test, we also define

, where ρ_{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

and

${\overline{A}}^{\text{HF}}(\mathrm{\omega})$ [recall that the micromotion part σ_{MM}(*t*) does not contribute to the one-cycle averages]. At high-frequency ω ≳ 10, the difference

decreases quite well. In Fig. 2 (A and B), we plot Δ*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)

${\mathrm{\Sigma}}_{A}(\mathrm{\omega})={[{T}^{-1}{\int}_{t}^{t+T}\mathit{ds}{\{\text{tr}[\mathrm{\rho}(s)A]-\overline{A}(\mathrm{\omega})\}}^{2}]}^{1/2}$, which quantifies the micromotion amplitude. Similarly to

${\overline{A}}^{\text{HF}}(\mathrm{\omega})$, we also introduce the one-cycle SD based on our high-frequency formulas,

${\mathrm{\Sigma}}_{A}^{\text{HF}}(\mathrm{\omega})={[{T}^{-1}{\int}_{t}^{t+T}\mathit{ds}{\{\text{tr}[{\mathrm{\rho}}^{\text{HF}}(s)A]-{\overline{A}}^{\text{HF}}(\mathrm{\omega})\}}^{2}]}^{1/2}$. These quantities are suitable for testing Eq. 10 because they include the contribution only from the micromotion part σ_{MM}(*t*). Because Σ* _{A}*(ω) is an

*O*(ω

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

is *O*(ω^{−2}). This criterion is satisfied, as shown in Fig. 2 (C and D) for *A* = *S _{z}* and {

*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} has the independent solutions ∣ψ* _{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}(13)where

${Z}_{\text{FG}}={\displaystyle \sum _{i}}{e}^{-\mathrm{\beta}{\mathrm{\u03f5}}_{i}}$. To obtain the Floquet states and quasienergies in practice, the common method, which we use here, is to calculate the one-cycle unitary evolution

$U(T)=\mathbb{T}\text{exp}[-i{\int}_{0}^{T}\mathit{ds}H(s)]$, where

$\mathbb{T}\text{exp}$ denotes the time-ordered exponential, by numerical integrations of the time-dependent Schrödinger equation. The eigenvectors and eigenvalues of *U*(*T*) correspond to ∣*u _{i}*(0)〉 and

, which give us ϵ* _{i}* and ∣

*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

${\overline{A}}^{\text{FG}}(\mathrm{\omega})={T}^{-1}{\int}_{t}^{t+T}\mathit{ds}\text{tr}[{\mathrm{\rho}}_{\text{FG}}(s)A]$and the SD

${\mathrm{\Sigma}}_{A}^{\text{FG}}(\mathrm{\omega})={[{T}^{-1}{\int}_{t}^{t+T}\mathit{ds}{\{\text{tr}[{\mathrm{\rho}}_{\text{FG}}(s)A]-{\overline{A}}^{\text{FG}}(\mathrm{\omega})\}}^{2}]}^{1/2}$. In Fig. 2 (A and D), we plot the difference

$\mathrm{\Delta}{A}^{\text{FG}}=\mid {\overline{A}}^{\text{FG}}(\mathrm{\omega})-\overline{A}(\mathrm{\omega})\mid $ for the two representative observables *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

$\overline{\{{S}_{x},{S}_{y}\}}(\mathrm{\omega})=0$ for any ω, while it is not true in the actual dynamics. This is due to an antiunitary dynamical symmetry constraining the Floquet states and, hence, the FGS. We take an antiunitary operator *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

is also a Floquet state with quasienergy ϵ* _{i}*. Assuming that quasienergies are not degenerate as in our examples, we have that

and ∣*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

and ∣*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

$\overline{\{{S}_{x},{S}_{y}\}}(\mathrm{\omega})$. We can show that this average vanishes by taking the limit γ* _{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

(14)Note that the CFSS is similar to the FGS (Eq. 13) and obtained by replacing the quasienergy ϵ* _{i}* with the real energy

*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

characterized by

$-i[{H}_{0},\tilde{\mathrm{\rho}}]+\mathrm{D}(\tilde{\mathrm{\rho}})=0$. Correspondingly, our results for the NESS (Eqs. 9 to 11) hold true with the following replacements:

${\mathrm{\rho}}_{\text{can}}\to \tilde{\mathrm{\rho}}$and

${p}_{\text{can}}^{(k)}\to {\tilde{p}}^{(k)}=\u3008{E}_{k}\mid \tilde{\mathrm{\rho}}\mid {E}_{k}\u3009$. Thus, our formulas apply to any periodically driven dissipative systems as long as the dissipation is of Lindblad type and irreducible. Therefore, our formulas are useful for a broad class of systems in exploring generic features of the NESS and in estimating Floquet-engineered physical quantities.

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

(15)We introduce the Fourier series for the Liouvillian as ℒ(*t*) = Σ* _{m}*ℒ

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

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

_{ij}(16)

$${\mathcal{L}}_{m}\mathrm{\rho}=-i[{H}_{m},\mathrm{\rho}](\text{for}m\ne 0)$$(17)The formal solution of the Floquet-Lindblad equation is obtained as

$\mathrm{\rho}(t)=\mathrm{V}(t,t\prime )\mathrm{\rho}(t\prime )$, with the propagator

$\mathrm{V}(t,t\prime )=\mathbb{T}\text{exp}[{\int}_{t\prime}^{t}\mathcal{L}(s)\mathit{ds}]$, where

$\mathbb{T}\text{exp}$denotes the time-ordered exponential. The determining equations for V are

${\mathrm{\partial}}_{t}\mathrm{V}(t,t\prime )=\mathcal{L}(t)\mathrm{V}(t,t\prime )$and

$\mathrm{V}(t\prime ,t\prime )=1$.

The high-frequency expansion in terms of the van Vleck approach makes the following ansatz

$$\mathrm{V}(t,t\prime )={e}^{\mathrm{G}(t)}{e}^{(t-t\prime ){\mathcal{L}}_{\text{eff}}}{e}^{-\mathrm{G}(t\prime )}$$(18)where

$\mathrm{G}(t)$ is periodic in time and ℒ_{eff} is time independent. This ansatz satisfies one determining equation

for any choices of

$\mathrm{G}(t)$ and ℒ_{eff}, and what determines these two is

. As we will see below, this equation only determines the derivative of

$\mathrm{G}(t)$, and thus, we further impose

${\int}_{0}^{T}\mathrm{G}(t)\mathit{dt}=0$to fix the constant of integration.

To obtain the determining equations for

$\mathrm{G}(t)$ and ℒ_{eff}, we substitute Eq. 18 into

, having

${\mathrm{\partial}}_{t}({e}^{\mathrm{G}(t)})+{e}^{\mathrm{G}(t)}{\mathcal{L}}_{\text{eff}}=\mathcal{L}(t){e}^{\mathrm{G}(t)}$. With some algebra (see the Supplementary Materials), we rewrite this equation, obtaining

$${\mathrm{\partial}}_{t}\mathrm{G}(t)=\sum _{k=0}^{\infty}\frac{{B}_{k}}{k!}{({\text{ad}}_{\mathrm{G}})}^{k}[\mathcal{L}(t)+{(-1)}^{k+1}{\mathcal{L}}_{\text{eff}}]$$(19)where

${\text{ad}}_{\mathrm{G}}$is defined by

${\text{ad}}_{\mathrm{G}}\mathrm{\rho}=[\mathrm{G}(t),\mathrm{\rho}]$, and *B _{k}* denotes the Bernoulli number (

*B*

_{0}= 1,

*B*

_{1}= − 1/2,

*B*

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

We determine

$\mathrm{G}(t)$ and ℒ_{eff} from Eq. 19 by the series expansions

(20)We substitute these expansions into Eq. 19 and find the order-by-order solutions, where we assign an order 1 for ℒ(*t*) and *k* for

and

${\mathcal{L}}_{\text{eff}}^{(k)}$ [see (*6*) for the case of unitary dynamics]. With some straightforward calculations (see the Supplementary Materials), we obtain the following low-order solutions

(21)

$${\mathrm{G}}^{\left(1\right)}(t)=\frac{i}{\mathrm{\omega}}{\displaystyle \sum _{m\ne 0}}\frac{{e}^{-\mathit{\text{im}}\mathrm{\omega}t}}{m}{\mathcal{L}}_{m}$$(22)

$${\mathcal{L}}_{\text{eff}}^{(2)}=-\frac{i}{\mathrm{\omega}}{\displaystyle \sum _{m>0}}\frac{[{\mathcal{L}}_{m},{\mathcal{L}}_{-m}]}{m}$$(23)Note that higher-order components are *O*(ω^{−2}), and

is equivalent to that used in the main text. We also note that

${\mathcal{L}}_{\text{eff}}={\mathcal{L}}_{\text{eff}}^{(1)}+{\mathcal{L}}_{\text{eff}}^{(2)}$up to this order can be represented by the effective Hamiltonian

$${H}_{\text{eff}}={H}_{0}+\frac{1}{\mathrm{\omega}}{\displaystyle \sum _{m>0}}\frac{[{H}_{-m},{H}_{m}]}{m}+O({\mathrm{\omega}}^{-2})$$(24)By invoking the Jacobi identity [*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

${\mathcal{L}}_{\text{eff}}{\mathrm{\rho}}_{\infty}^{\prime}=0$for

${\mathrm{\rho}}_{\infty}^{\prime}$ at the leading order of ω^{−1}. It is convenient to work in the energy eigenbasis and separate the diagonal and off-diagonal parts

(25)

$${\mathrm{\rho}\prime}_{\infty}^{(\mathrm{d})}=\sum _{k}{\mathrm{\rho}\prime}_{\infty}^{\mathit{kk}}\mid {E}_{k}\mathrm{\u3009}\mathrm{\u3008}{E}_{k}\mid $$(26)

$${\mathrm{\rho}\prime}_{\infty}^{(\text{od})}=\sum _{\begin{array}{c}k,l\\ (k\ne l)\end{array}}{\mathrm{\rho}\prime}_{\infty}^{\mathit{kl}}\mid {E}_{k}\mathrm{\u3009}\mathrm{\u3008}{E}_{l}\mid $$(27)

First, we consider the off-diagonal elements of both sides of

${\mathcal{L}}_{\text{eff}}{\mathrm{\rho}}_{\infty}^{\prime}=0$, having, for *k* ≠ *l*

(28)where

${\mathrm{\gamma}}_{\mathit{kl}}\equiv \sum _{i}({\mathrm{\Gamma}}_{\mathit{ik}}+{\mathrm{\Gamma}}_{\mathit{il}})/2$ and Δ*H*_{eff} ≡ *H*_{eff} − *H*_{0} = *O*(ω^{−1}). Equation 28 is transformed as

(29)Note that the denominators (*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} are *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

, and Eq. 29 can be solved recursively. This yields the ω^{−1} expansion for the off-diagonal elements

, whose leading order contribution is given by

$${\mathrm{\rho}\prime}_{\infty}^{\mathit{kl}}=\frac{\mathrm{\u3008}{E}_{k}\mid \mathrm{\Delta}{H}_{\text{eff}}\mid {E}_{l}\mathrm{\u3009}}{({E}_{k}-{E}_{l})-i{\mathrm{\gamma}}_{\mathit{kl}}}({\mathrm{\rho}\prime}_{\infty}^{\mathit{kk}}-{\mathrm{\rho}\prime}_{\infty}^{\mathit{ll}})+O({\mathrm{\omega}}^{-2})$$(30)

Next, we consider the diagonal elements of both sides of

${\mathcal{L}}_{\text{eff}}{\mathrm{\rho}}_{\infty}^{\prime}=0$, having

$$\u3008{E}_{k}\mid {\mathcal{L}}_{\text{eff}}{\mathrm{\rho}}_{\infty}^{\prime}\mid {E}_{k}\u3009=-i\u3008{E}_{k}\mid [\mathrm{\Delta}{H}_{\text{eff}},{{\mathrm{\rho}}^{\prime}}_{\infty}^{(\text{od})}]\mid {E}_{l}\u3009+{\mathrm{\Sigma}}_{l}({\mathrm{\Gamma}}_{\mathit{kl}}{{\mathrm{\rho}}^{\prime}}_{\infty}^{\mathit{ll}}-{\mathrm{\Gamma}}_{\mathit{lk}}{{\mathrm{\rho}}^{\prime}}_{\infty}^{\mathit{kk}})=0$$(31)We note

$\u3008{E}_{k}\mid [\mathrm{\Delta}{H}_{\text{eff}},{{\mathrm{\rho}}^{\prime}}_{\infty}^{(\text{od})}]\mid {E}_{l}\u3009=O({\mathrm{\omega}}^{-2})$ because both Δ*H*_{eff} and

are *O*(ω^{−1}). Thus, the diagonal elements

are determined up to *O*(ω^{−1}) by the equation

(32)According to the irreducibility and the detailed balance condition of Γ* _{kl}*, we have the unique solution as

(33)where the error is *O*(ω^{−2}). This result means

(34)Because these diagonal elements

${{\mathrm{\rho}}^{\prime}}_{\infty}^{\mathit{kk}}$ are *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

(35)which implies

$${{\mathrm{\rho}}^{\prime}}_{\infty}^{(\text{od})}={\mathrm{\sigma}}_{\text{FE}}+O({\mathrm{\omega}}^{-2})$$(36)We remark that tr(σ_{FE}) = 0 because each of the diagonal elements of σ_{FE} vanishes.

Now that we have obtained the leading-order expression for

${\mathrm{\rho}}_{\infty}^{\prime}={{\mathrm{\rho}}^{\prime}}_{\infty}^{(\mathrm{d})}+{{\mathrm{\rho}}^{\prime}}_{\infty}^{(\text{od})}$, let us calculate the time-dependent density matrix by

$\mathrm{\rho}(t)={e}^{\mathrm{G}(t)}{\mathrm{\rho}}_{\infty}^{\prime}$. By noticing that

${{\mathrm{\rho}}^{\prime}}_{\infty}^{(\mathrm{d})}$ is *O*(ω^{0}) and

is *O*(ω^{−1}) and using the Taylor expansion

, we obtain

$$\mathrm{\rho}(t)={\mathrm{\rho}}_{\text{can}}+{\mathrm{\sigma}}_{\text{MM}}(t)+{\mathrm{\sigma}}_{\text{FE}}+O({\mathrm{\omega}}^{-2})$$(37)where we have defined

$${\mathrm{\sigma}}_{\text{MM}}(t)=\mathrm{G}(t)[{\mathrm{\rho}\prime}_{\infty}^{(\mathrm{d})}]=\frac{1}{\mathrm{\omega}}\sum _{m\ne 0}\frac{{e}^{-\mathit{\text{im}}\mathrm{\omega}t}}{m}[{H}_{m},{\mathrm{\rho}}_{\text{can}}]$$(38)Thus, we have derived the Eqs. 9 to 11 here. We remark tr[σ_{MM}(*t*)] = 0, which follows from the cyclic property of trace, and hence tr[ρ(*t*)] = 1, at least up to this order.

**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.