General description for nonequilibrium steady states in periodically driven dissipative quantum systems.

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 (–). 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 () and the ensuing theoretical developments (–). 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 () in solids (), ultracold atomic gases (), and in photonic wave guides (), Floquet time crystals () in nitrogen-vacancy (NV) centers () and trapped ions (), and control of quantum magnets () and their interactions ().

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 (–). 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 (, ), periodic thermodynamics (), dynamical localization (), and generalized Bose-Einstein condensation (). 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 (–). 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) (, , ), i.e., the Gibbs state with the effective Hamiltonian, but the conditions for the FGS being realized have shown quite restrictive (). 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 (, , –), to the Lindblad equation with periodic Hamiltonians. As exemplified in an effective model for the NV center in diamonds (), 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 H0 denote the time-independent Hamiltonian, which describes our system in the absence of driving. The eigenenergies and eigenstates of H0 are denoted by and , respectively. For simplicity, we assume that the eigenenergies are not degenerate and E1 < E2 < ⋯ < E (the generalization to degenerate H0 is formulated in the Supplementary Materials). The effect of the driving is represented by a time-dependent part Hext(t) of the total HamiltonianWe assume that the driving term is periodic with period T: Hext(t + T) = Hext(t) and hence H(t + T) = H(t). Without loss of generality, the decomposition (Eq. 1) is defined so that the time average of Hext(t) vanishes, . Thus, the Fourier series of Hext(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 (–) (we set ℏ = 1 throughout this paper)Here, L ≡ ∣E〉〈E∣ is the time-independent Lindblad operator describing the transition from the j-th to the i-th eigenstates of the undriven Hamiltonian H0. Each Lindblad operator L represents a decay (excitation) process for i < j (i > j). The real number Γ ( ≥ 0) denotes the rate for the corresponding process, and we set Γ = 0 for each 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 d tr[ρ(t)]/dt = 0, and thus, we use the normalization tr[ρ(t)] = 1. We assume that the transition rates Γ satisfy the detailed balance conditionwhere β is the inverse temperature of the bath coupled to the system. We also assume that the matrix Γ is a nonnegative irreducible matrix (). 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−β/Z of H0 with Z = tr(e−βH). We note that the Lindblad operators L may depend on the driving in general if we consider more microscopic theories of dissipation (). 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 (). Among several formulations, we adopt the van Vleck perturbation theory (, ), which leads to the following propagation for ρ(t) (see Materials and Methods): . The time-independent part ℒeff is represented by the effective Hamiltonianwith Heff = H0 + ω−1∑[H−, H]/n + O(ω−2). The time-dependent part is the so-called micromotion operator periodic in time and is given by . Without loss of generality, we suppose the initial time to be t′ = 0, havingwith ρ(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 , ρ′(t) approaches, irrespective of the initial state, the unique state characterized by . Thus, we come to the first main result, obtaining the asymptotic behaviorBecause , this NESS is also periodic in time. Focusing on the leading-order contribution, we have a simple explicit formula for ρness(t)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 byWe 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 byand 〈E∣σFE∣E〉 = 0 for all k, where ΔHeff ≡ Heff − H0 = O(ω−1), is the Boltzmann weight, and represents the symmetric transition-rate matrix (see the Supplementary Materials for the generalization to degenerate H0). 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 ρ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(ρcanA) = 0, but nonzero value for the NESS, tr(σFEA) ≠ 0. This situation means that one can implement an appropriate periodic driving Hext(t) and hence ΔHeff, 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〉} of the time-independent Hamiltonian H0. 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 N2-dimensional vector and the superoperator ℒ as an N2 × N2 matrix. Thus, the computational complexity for each time step is O(N4) in general. On the other hand, the complexity for our formula is one order smaller and given by O(N3), which derives from the exact diagonalization of H0. 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 O(N0) nonzero elements, the computational complexity for one time step is O(N2). 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 ()where B is the static Zeeman field, and are the coupling constants of magnetic anisotropic terms, and 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 in the NV centers (), we set and in our analysis. The ac-field irradiation leads to spin dynamics, as schematically illustrated in Fig. 1A.

Fig. 1

NESSs of an NV center in diamond.

(A) Schematic illustration of the NV-center spin driven by laser magnetic field. (B and C) Time evolution of (B) S and (C) {S, S} = S + S of the spin described by the Floquet-Lindblad equation (Eq. 3) with the Hamiltonian (Eq. 12) (see text for the parameters). The dashed line shows the one-cycle average at a sufficiently long time. (D and E) One-cycle average for (D) A = S and (E) A = {S, S} calculated from the dynamics simulation (solid line), our formula (Eq. 9) (circle), and the FGS (Eq. 13) (square). All the other parameters except ω are the same as in the top panels.

Typical time evolutions A(t) ≡ tr[ρ(t)A] are shown for two representative observables A = S and {S, S} ( = S + S) in Fig. 1 (B and C). In this calculation, we take the thermal state for the time-independent part of HNV at t = 0 and let the state evolve according to the Floquet-Lindblad equation (Eq. 3) with H(t) being HNV(t). The static Zeeman field is B = 0.3, and the driving parameters are B = 0.1 and ω = 1. As for the Lindblad operators, we take Γ according to the heatbath method as Γ = γ/( + ) for i ≠ j, with rate constant γ = 0.2 and inverse temperature β = 3, and Γ = 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, S} is initially zero for a symmetry reason (e.g., S → − S) but becomes nonzero on time average. Namely, this observable is engineered by the periodic drive .

To test our formula (Eq. 9) quantitatively, we first focus on the one-cycle average 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 [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 ΔAHF(ω) against ω for A = S and {S, S}, respectively. We stress that the difference ΔAHF(ω) decreases more rapidly than O(ω−1) [O(ω−2) for S and O(ω−3) for {S, S}]. This means that Eqs. 9 and 11 perfectly describe the actual NESS at the level of O(ω−1). As shown in the Supplementary Materials, ΔAHF(ω) = 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.

Fig. 2

Asymptotic errors of our formula and the Floquet-Gibbs state.

(Top) Difference ΔA(ω) between the one-cycle average calculated from the dynamics simulation and that from our formula (Eq. 9) (circle) or the FGS (Eq. 13) (square) plotted against ω for (A) A = S and (B) A = {S, S}. (Bottom) Difference ΔΣ(ω) between the one-cycle SD Σ(ω) calculated from the dynamics simulation and that from our formula (Eq. 9) (circle) or the FGS (Eq. 13) (square) plotted against ω for (C) A = S and (D) A = {S, S}. In all the panels, the solid and dashed lines are guides for the eye, showing ∝ω−2 and ∝ω−1, respectively.

For the complementary test of our formula (Eq. 9), we consider the one-cycle standard deviation (SD) , which quantifies the micromotion amplitude. Similarly to , we also introduce the one-cycle SD based on our high-frequency formulas, . These quantities are suitable for testing Eq. 10 because they include the contribution only from the micromotion part σMM(t). Because Σ(ω) 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 and {S, S}, respectively. We remark that our formula leads to Σ(ω) = 0 at O(ω−1) for these observables, which can be analytically shown by noticing H±1 ∝ S± = S ± iS. Thus, the plotted data correspond to Σ(ω) itself for the actual dynamics, and ΔΣ(ω) could be reduced by dealing with the higher-order terms in Eq. 9. In any case, the fact that ΔΣ(ω) is 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 (, , ). To define the FGS, we introduce the Floquet state ∣u(t)〉 and its quasienergy ϵ. According to the Floquet theorem, the time-dependent Schrödinger equation has the independent solutions ∣ψ(t)〉 = ∣u(t)〉 (i = 1,2, …, N), with periodicity ∣u(t + T)〉 = ∣u(t)〉 and 〈u(t)∣u(t)〉 = δ. In terms of the Floquet states, the FGS is defined bywhere . To obtain the Floquet states and quasienergies in practice, the common method, which we use here, is to calculate the one-cycle unitary evolution , where 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(0)〉 and , which give us ϵ and ∣u(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 and the SD . In Fig. 2 (A and D), we plot the difference for the two representative observables A = S and {S, S}. Remarkably, the difference is 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 weak point of the FGS is highlighted in Fig. 1E, in which the FGS gives 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† = − S and VS† = S (α = x and z). Then, we notice the dynamical symmetry VHNV(T − t)V† = HNV(t), which implies that is also a Floquet state with quasienergy ϵ. Assuming that quasienergies are not degenerate as in our examples, we have that and ∣u(t)〉 are equivalent up to an overall phase shift. Because of VAV† = − A with A = {S, S}, the one-cycle averages of A calculated for and ∣u(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 . We can show that this average vanishes by taking the limit γ → 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 γ 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 byNote that the CFSS is similar to the FGS (Eq. 13) and obtained by replacing the quasienergy ϵ 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 γ → 0 (see the Supplementary Materials for details).

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 Γ (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).

Fig. 3

Validity of NESS formulas in parameter space.

Our formula is valid when the driving frequency ω is larger than energy scales of the system. When the dissipation strength is smaller than the (nonzero) minimum energy gap, our formula reduces to the CFSS. At lower frequency, the NESS description remains an open question.

One should note that the Floquet-Lindblad equation (Eq. 3) becomes invalid when Γ is too large. The Lindblad-type dissipation is derived from several approximations such as the Born-Markov approximation (). 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 Γ 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.

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 Γ is irreducible. In this case, the solution of the Lindblad equation without driving is not the canonical ensemble ρcan but another state characterized by . Correspondingly, our results for the NESS (Eqs. 9 to 11) hold true with the following replacements: and . 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), where the time-dependent Liouvillian ℒ(t) is defined byWe introduce the Fourier series for the Liouvillian as ℒ(t) = Σℒ−. Because the Lindblad operators L are time independent in this work, each Fourier component is given as followsThe formal solution of the Floquet-Lindblad equation is obtained as , with the propagator , where denotes the time-ordered exponential. The determining equations for are and .

The high-frequency expansion in terms of the van Vleck approach makes the following ansatzwhere is periodic in time and ℒeff is time independent. This ansatz satisfies one determining equation for any choices of and ℒeff, and what determines these two is . As we will see below, this equation only determines the derivative of , and thus, we further impose to fix the constant of integration.

To obtain the determining equations for and ℒeff, we substitute Eq. 18 into , having . With some algebra (see the Supplementary Materials), we rewrite this equation, obtainingwhere is defined by , and B denotes the Bernoulli number (B0 = 1, B1 = − 1/2, B2 = 1/6, …).

We determine and ℒeff from Eq. 19 by the series expansionsWe 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 [see () for the case of unitary dynamics]. With some straightforward calculations (see the Supplementary Materials), we obtain the following low-order solutionsNote that higher-order components are O(ω−2), and is equivalent to that used in the main text. We also note that up to this order can be represented by the effective HamiltonianBy 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 for at the leading order of ω−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 , having, for k ≠ lwhere and ΔHeff ≡ Heff − H0 = O(ω−1). Equation 28 is transformed asNote that the denominators (E − E) − iγ do not vanish because γ > 0 is ensured by the nonnegativity and irreducibility of Γ. Now, as a working hypothesis, we suppose that the diagonal elements are O(ω0) as verified later. Then, the first term on the right-hand side of Eq. 29 is O(ω−1) because ΔHeff = 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

Next, we consider the diagonal elements of both sides of , havingWe note because both ΔHeff and are O(ω−1). Thus, the diagonal elements are determined up to O(ω−1) by the equationAccording to the irreducibility and the detailed balance condition of Γ, we have the unique solution aswhere the error is O(ω−2). This result meansBecause these diagonal elements 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 elementswhich impliesWe remark that tr(σFE) = 0 because each of the diagonal elements of σFE vanishes.

Now that we have obtained the leading-order expression for , let us calculate the time-dependent density matrix by . By noticing that is O(ω0) and is O(ω−1) and using the Taylor expansion , we obtainwhere we have definedThus, 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.