Coarse-graining master equation for periodically driven systems
CCoarse-graining master equation for periodically driven systems
Ronja Hotz and Gernot Schaller , ∗ Institut f¨ur Theoretische Physik, Technische Universit¨at Berlin,Hardenbergstr. 36, 10623 Berlin, Germany and Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstraße 400, 01328 Dresden, Germany (Dated: February 8, 2021)We analyze Lindblad-Gorini-Kossakowski-Sudarshan-type generators for selected periodicallydriven open quantum systems. All these generators can be obtained by temporal coarse-grainingprocedures, and we compare different coarse-graining schemes. Similar to for undriven systems, wefind that a dynamically adapted coarse-graining time, effectively yielding non-Markovian dynamicsby interpolating through a series of different but invididually Markovian solutions, yields the bestresults among the different coarse-graining schemes, albeit at highest computational cost.
I. INTRODUCTION
While the propagation of undriven closed quantumsystems via Schr¨odingers equation may be challengingfor many degrees of freedom, it can be calculated withstandard methods. This changes when the Hamiltonianis subject to an external driving, when time-orderingbecomes relevant. Even for the simple case of periodicdriving, where Floquet theory [1, 2] applies, the exactcalculation of the unitary propagator may be notori-ously difficult [3–5]. In addition, it is then often hiddenthat the time-dependent variation of external parame-ters requires that the system is coupled in some wayto an outside world, which to be complete would actu-ally require to model the system as open. While openquantum systems in absence of driving are well-studiedand interesting-subject on their own [6–8], the joint dis-cussion of the effects of system-reservoir coupling andperiodic driving is challenging.However, with in particular periodic types of driv-ing being experimentally quite feasible with currentstandards [9–13], the derivation of proper dissipatorsfor periodically driven open systems is of great con-cern. In principle, it is formally possible to per-form the same approximations that lead to Lindblad-Gorini-Kossakowski-Sudarshan (LGKS) master equa-tions [14, 15] for an undriven system also for its drivenversion [3, 16–19]. In the interaction picture, the ap-plication of Born-, Markov-, and secular approxima-tions generically leads to LGKS master equations thatcome handy due to their stability and – despite the lackof thermalization [20, 21] – their transparent thermo-dynamic interpretation [22]. However, it has also fre-quently been noted that in particular the involved secu-lar approximation may lead to strong artifacts [23, 24].An intuitive explanation for this is that the level spacingin the extended (Floquet) space may become very small,which conflicts with the secular approximation [25].Therefore, traditional Floquet-LGKS approaches havebeen questioned in many works [23, 26].Already for undriven systems, it can be argued that ∗ [email protected] the (non-LGKS) Redfield equation applied in the properregime only leads to small violations of density matrixproperties in trade for a closer description of physicalreality [27]. As a practical benefit, improper densitymatrices can then be used as an indicator for leaving theregion of validity of such perturbative schemes, whereassuch a witness is missing for LGKS-based approaches. tt =0 t +T t +Tt →∞ ∞ TimeSignal
FIG. 1. Sketch illustrating the different coarse-grainingtimescales to derive a master equation in the weak cou-pling limit. For CG1 (orange) and CG2 (green) the coarse-graining extends over one period of the driving field τ = T ,where for CG1 t = 0 and for CG2 the initial coarse-grainingtime is taken to infinity t → ∞ . For dynamical coarse-graining (DCG) (purple) the coarse-graining time is chosendynamically as the physical time. For τ → ∞ all methodsrecover the Born-Markov-secular (BMS) limit (red). Another viewpoint therefore is to preserve the formalGKSL form of the generator while improving in someways the involved approximations [28]. There are var-ious routes to achieve this, e.g. formally restoring theLGKS form by discarding its subspace with negativeeigenvalues [29, 30], by performing only a partial sec-ular approximation [31, 32], or by taking the effects ofthe reservoir on the direction of the system pointer basisinto account [33, 34]. While these approaches apply toundriven systems, we would like to focus on an approachthat can be easily combined with periodic driving. Suchan approach that can also be well-motivated using a mi-croscopic derivation is the temporal coarse-graining ap-proach [27, 35–46]. Intuitively, it can be understood as a r X i v : . [ qu a n t - ph ] F e b an approximate representation of the reduced exact so-lution over an appropriately chosen coarse-graining timeinterval. The freedom in choosing the coarse-graininginterval is exemplified in Fig. 1 and shall be the subjectof this paper with a focus on periodically driven systems.We begin with a brief introduction into the method inSec. II and then compare various dissipators for threeexamples of driven open systems in Secs. III, IV, and V. II. COARSE-GRAINING METHODS FORPERIODICALLY DRIVEN OPEN SYSTEMS
The starting point of our considerations is a repre-sentation of the Hamiltonian of the full universe in theform H ( t ) = H S ( t ) + (cid:88) α A α ⊗ B α + H B (1)with periodically driven system Hamiltonian H S ( t + T ) = H S ( t ), (not necessarily individually Hermitian)system and bath coupling operators A α and B α , re-spectively, and reservoir Hamiltonian H B . Since systemand bath operators act on different Hilbert spaces, weassume from the beginning that [ H S , H B ] = 0 and also[ A α , B β ] = 0. In an interaction picture with respect to H = H S ( t )+ H B (bold symbols), we can formally solvethe von-Neumann equation ˙ ρ = − i (cid:34)(cid:88) α A α ( t ) ⊗ B α ( t ) , ρ (cid:35) (2)by employing the time-evolution operator from time t to time t + τ ρ ( t + τ ) = U ( t + τ, t ) ρ ( t ) U † ( t + τ, t ) . (3)The coarse-graining approach originates from the at-tempt to match the reduced evolution of an open quan-tum system with a time-local master equation e L τ τ ρ S ! = Tr B (cid:8) U ( t + τ, t ) ρ S ⊗ ρ B U † ( t + τ, t ) (cid:9) . (4)Here, we have already put an index L τ to the dissi-pation superoperator to highlight the fact that with aconstant L superoperator it will in general not be possi-ble to capture the exact reduced dynamics for all times,which would require a Kraus representation [47]. Forweak system-reservoir interaction (or small propagationtimes τ ) we can use that U ( t + τ, t ) is close to theidentity and L τ is small. One can then expand the ex-pressions on both sides and assuming furthermore thatTr (cid:8) B α ρ B (cid:9) = 0, one has a defining equation for the coarse-graining generator L τ ρ S = − i 12i τ t + τ (cid:90) t dt dt (cid:88) α ¯ α C α ¯ α ( t , t ) × sgn( t − t ) [ A α ( t ) A ¯ α ( t ) , ρ S ]+ 1 τ t + τ (cid:90) t dt dt (cid:88) α ¯ α C α ¯ α ( t , t ) × (cid:104) A ¯ α ( t ) ρ S A α ( t ) − { A α ( t ) A ¯ α ( t ) , ρ S } (cid:105) , (5)with bath correlation functions C α ¯ α ( t , t ) =Tr B { B α ( t ) B ¯ α ( t ) ρ B } = Tr (cid:8) B α ( t − t ) B β ρ B (cid:9) ,where the last equality holds when [ H B , ρ B ] = 0.Here, the bold symbols denote the interaction pic-ture, and for system operators we have A ( t ) = U † S ( t, AU S ( t,
0) with U S ( t,
0) = T (cid:110) e − i (cid:82) t H S ( t (cid:48) ) dt (cid:48) (cid:111) and T denoting the time ordering operator. Particu-larly, for a periodically driven system it can be expressedas a product of a unitary kick operator U kick ( t,
0) withthe period of the driving and the evolution under aneffective system Floquet Hamiltonian ¯ H : U S ( t,
0) = U kick ( t, e − i ¯ Ht U kick ( t,
0) = U kick ( t + T, , U kick ( nT,
0) = , (6)with n ∈ Z , such that the stroboscopic evolution be-tween full periods is given by the Floquet Hamiltonian U ( nT,
0) = e − i ¯ HnT . However, it should be stressedhere that the above decomposition is not unique. In-troducing the eigenbasis of the Floquet Hamiltonianvia ¯ H | ¯ a (cid:105) = ¯ E a | ¯ a (cid:105) , we can see that the time evolu-tion operator remains invariant under the simultaneoustransformations ¯ H → ¯ H + m Ω | ¯ a (cid:105) (cid:104) ¯ a | and U kick ( t, → U kick ( t, e +i m Ω t | ¯ a (cid:105)(cid:104) ¯ a | , which maintains the periodicityof the kick operator. This gauge degree of freedom al-lows one to deliberately choose the convention that alleigenenergies of the Floquet Hamiltonian should lie inthe first Brillouin zone ¯ E a ∈ [ − Ω / , +Ω / A α ( t ) can be expressed as: A α ( t ) = U † S ( t ) A α U S ( t ) = e i ¯ Ht U † kick ( t ) A α U kick ( t ) e − i ¯ Ht = + ∞ (cid:88) n = −∞ e i ¯ Ht ˆ A nα e i n Ω t e − i ¯ Ht , (7)where due to the periodicity of the kick operator wehave expanded U † kick ( t ) A α U kick ( t ) in a Fourier seriesˆ A nα = Ω2 π T/ (cid:90) − T/ U † kick ( t ) A α U kick ( t ) e − i n Ω t dt, T = 2 π Ω . (8)Now, in the eigenbasis of the Floquet Hamiltonian wecan also write this as A α ( t ) = (cid:88) ab (cid:88) n (cid:104) ¯ a | ˆ A nα | ¯ b (cid:105) e i( ¯ E a − ¯ E b + n Ω) t | ¯ a (cid:105) (cid:104) ¯ b |≡ (cid:88) ab (cid:88) n A + nα,ab e +i( ¯ E a − ¯ E b + n Ω) t L ab , A ¯ α ( t ) = (cid:88) cd (cid:88) n (cid:48) (cid:104) ¯ c | ˆ A n (cid:48) ¯ α | ¯ d (cid:105) e i( ¯ E c − ¯ E d + n (cid:48) Ω) t | ¯ c (cid:105) (cid:104) ¯ d |≡ (cid:88) cd (cid:88) n (cid:48) A − n (cid:48) ¯ α,dc e − i( ¯ E c − ¯ E d + n (cid:48) Ω) t L dc , (9)where L ab ≡ | ¯ a (cid:105) (cid:104) ¯ b | .Additionally, the reservoir correlation functions forreservoir at equilibrium depend only on the differenceof time arguments, and we can introduce their even andodd Fourier transforms C α ¯ α ( t − t ) = 12 π (cid:90) dωγ α ¯ α ( ω ) e − i ω ( t − t ) ,C α ¯ α ( t − t )sgn(t − t ) = 12 π (cid:90) dωσ α ¯ α ( ω ) e − i ω ( t − t ) . (10)The temporal integrals can be analytically performedleading to f τt ( α, β, ω ) ≡ πτ t + τ (cid:90) t t + τ (cid:90) t e − i ω ( t − t ) e i αt e − i βt dt dt = τ π e i( α − β ) t e i( α − β ) τ/ sinc (cid:104) ( α − ω ) τ (cid:105) × sinc (cid:104) ( β − ω ) τ (cid:105) . (11)In particular, for discrete α and β one finds thatlim τ →∞ f τt ( α, β, ω ) = δ αβ δ ( α − ω ) . (12)Inserting it all, we can generally write the dissipator as ˙ ρ S = − i (cid:2) H τ LS , ρ S (cid:3) + D τ ρ S (13)with the single-integral expressions H τ LS = (cid:90) dω (cid:88) α ¯ α σ α ¯ α ( ω ) 12i (cid:88) nn (cid:48) (cid:88) abcd × f τt ( ¯ E a − ¯ E b + n Ω , ¯ E c − ¯ E d + n (cid:48) Ω , ω ) × A nα,ab A − n (cid:48) ¯ α,dc L ab L dc , D τ ρ S = (cid:90) dω (cid:88) α ¯ α γ α ¯ α ( ω ) (cid:88) nn (cid:48) (cid:88) abcd × f τt ( ¯ E a − ¯ E b + n Ω , ¯ E c − ¯ E d + n (cid:48) Ω , ω ) × A + nα,ab A − n (cid:48) ¯ α,dc × (cid:20) L dc ρ S L ab − { L ab L dc , ρ S } (cid:21) . (14) Here, the first term defines a (Hermitian) correction tothe system Hamiltonian (which for large τ convergesto the LGKS Lambshift term), and the second termdenotes the dissipative influence of the reservoir.Since by convention we chose the quasienergies inthe first Brillouin zone ¯ E a ∈ [ − Ω / , +Ω / E a − ¯ E b ∈ ( − Ω , +Ω). The important observation so faris that for any choice of t and τ , the coarse-grainingscheme always yields LGKS-type dissipators [48]. Nev-ertheless, the different coarse-graining schemes sketchedin Fig. 1 may lead to different dissipators. For laterreference, we coin the coarse-graining schemes dynami-cal coarse-graining [37] (DCG) with t = 0 and τ = t , initial period-coarse-graining (CG1) with t = 0 and τ = T , long-term period-coarse-graining with t → ∞ and τ = T (CG2).Finally, we would like to address the limit of verylarge coarse-graining times. Technically, this limit al-lows – based on Eq. (12) – to eliminate the remainingintegration H ∞ LS = 12i (cid:88) α ¯ α (cid:88) nn (cid:48) (cid:88) abcd σ α ¯ α ( ¯ E a − ¯ E b + n Ω) × δ ¯ E a − ¯ E b + n Ω , ¯ E c − ¯ E d + n (cid:48) Ω A nα,ab A − n (cid:48) ¯ α,dc L ab L dc , D ∞ ρ S = (cid:88) α ¯ α (cid:88) nn (cid:48) (cid:88) abcd γ α ¯ α ( ¯ E a − ¯ E b + n Ω) × δ ¯ E a − ¯ E b + n Ω , ¯ E c − ¯ E d + n (cid:48) Ω A + nα,ab A − n (cid:48) ¯ α,dc × (cid:20) L dc ρ S L ab − { L ab L dc , ρ S } (cid:21) . (15)The evaluation of the resonance condition resulting fromthe Kronecker- δ ¯ E a − ¯ E b + n Ω = ¯ E c − ¯ E d + n (cid:48) Ω (16)however requires some care. The typical argument isthat for fast driving, the above resonance can only bemet if separately n (cid:48) = n and ¯ E a − ¯ E b = ¯ E c − ¯ E d .However, the applicability of this argument criticallydepends on the Floquet spectrum. If for a two-levelsystem we by chance have Floquet energies well withinthe first Brillouin zone ¯ E a ∈ {− Ω / , +Ω / } , this gen-erates energy differences ¯ E a − ¯ E b ∈ {− Ω / , , +Ω / } .In this case, the above resonance condition can also bemet with ¯ E a − ¯ E b = +Ω /
2, ¯ E c − ¯ E d = − Ω / n (cid:48) = n + 1, demonstrating that in general, the long-term limit of periodically driven coarse-graining masterequation need not coincide with the significantly sim-pler Floquet-BMS master equation that results fromdemanding the resonance separately (i.e., by setting n (cid:48) = n and ¯ E a − ¯ E b = ¯ E c − ¯ E d ) H BMSLS = 12i (cid:88) α ¯ α (cid:88) n (cid:88) abcd σ α ¯ α ( ¯ E a − ¯ E b + n Ω) × δ ¯ E a − ¯ E b , ¯ E c − ¯ E d A nα,ab A − n ¯ α,dc L ab L dc , D BMS ρ S = (cid:88) α ¯ α (cid:88) n (cid:88) abcd γ α ¯ α ( ¯ E a − ¯ E b + n Ω) × δ ¯ E a − ¯ E b , ¯ E c − ¯ E d A + nα,ab A − n ¯ α,dc × (cid:20) L dc ρ S L ab − { L ab L dc , ρ S } (cid:21) . (17)Even then, we note that the remaining Kronecker- δ leaves many terms that are often neglected. For exam-ple, the above resonance can always be trivially fulfilledwith ¯ E a = ¯ E b and ¯ E c = ¯ E d , even for the undriven case.For many models, one has that γ α ¯ α (0) →
0, such thatsuch terms would not contribute anyways for undrivensystems, where Ω →
0. But in general they will haveto be kept. For comparison we therefore also state theultra-secular approximation (BMU), where only termswith a = c and b = d are kept H BMULS = 12i (cid:88) α ¯ α (cid:88) n (cid:88) ab σ α ¯ α ( ¯ E a − ¯ E b + n Ω) × A nα,ab A − n ¯ α,ba L ab L ba , D BMU ρ S = (cid:88) α ¯ α (cid:88) n (cid:88) ab γ α ¯ α ( ¯ E a − ¯ E b + n Ω) A + nα,ab A − n ¯ α,ba × (cid:20) L ba ρ S L ab − { L ab L ba , ρ S } (cid:21) . (18)In what follows, we will compare the solutions to vari-ous periodically driven problems that are based on thesedifferent coarse-graining approaches, i.e., for the DCGapproach we evaluate Eq. (14) with τ = t and t = 0,for the CG1 approach we use Eq. (14) with τ = T and t = 0, for the CG2 approach we use Eq. (14) with τ = T and t → ∞ , whereas for the BMS and BMU re-sults we use Eqns. (17) and (18), respectively. Further-more, in our calculations we will for simplicity neglectthe Hermitian correction term to the Hamiltonian H ... LS throughout. Whereas for the particular pure-dephasingmodel considered below this term has no effect anyways,in general it can only be neglected in the weak-couplingregime where | H ... LS | (cid:28) | H S | . III. PURE-DEPHASING MODELS
By the term pure-dephasing models we summarizemodels where the interaction commutes with the sys-tem Hamiltonian at all times, i.e., system and reservoircannot exchange energy. The Hamiltonian of a generaldriven-system pure-dephasing model with a reservoir of bosonic oscillators is then given by: H ( t ) = H S ( t ) + A ⊗ (cid:88) k (cid:16) h k b k + h ∗ k b † k (cid:17) + (cid:88) k ω k b † k b k , [ H S ( t ) , A ] = 0 . (19)The assumption of commuting coupling operator A andsystem Hamiltonian H S ( t ) allows to solve the dynamicsexactly. A. Exact solution
Specifically for a driven two level system with H S = σ z (cid:2) ∆2 + λ cos(Ω t ) (cid:3) and coupling A = σ z we ob-tain that the populations in the σ z eigenbasis remainconstant ρ S, ( t ) = ρ S, , ρ S, ( t ) = ρ S, , (20)whereas the coherences evolve according to ρ S, ( t ) = e +2i ( ∆2 t + λ Ω sin(Ω t ) ) (21) × exp − π ∞ (cid:90) −∞ dωγ ( ω ) sin (cid:0) ωt (cid:1) ω ρ S, . compare App. A. Here, γ ( ω ) = Γ( ω )[1 + n B ( ω )] withspectral density Γ( ω ) ≡ π (cid:80) k | h k | δ ( ω − ω k ) (analyti-cally continued via Γ( − ω ) = − Γ( ω ) to the complete realaxis) and Bose distribution n B ( ω ) = [ e βω − − . In ab-sence of driving ( λ → B. Master equation solutions
We can now compare this exact solution with the var-ious approximate approaches discussed before by iden-tifying the single coupling operator A = σ z . Partic-ularly, with an odd spectral density, the Fourier trans-form of the only reservoir correlation function becomes γ ( ω ) = Γ( ω )[1 + n B ( ω )] and the matrix elements ofthe coupling operators become rather trivial A nα,ab = δ α, δ n, δ ab (cid:104) ¯ a | σ z | ¯ a (cid:105) . To compute them, it suffices to re-alize that for this simple problem the system FloquetHamiltonian is just ¯ H = ∆2 σ z , such that its eigenstatesare trivial. C. Comparison
We analytically find that the DCG approach (adap-tively choosing τ = t and t = 0 in Eq. (14)) reproducesthe exact solution. Furthermore, we also analyticallyfind that CG1 and CG2 approaches ( τ = T in E. 14))both agree for this model class. We also find analytically dimensionless time ω t -0,4-0,200,20,4 c oh e r e n ce R e [ ρ S , ] = < σ x > / exact solution / DCGCG 1 / CG 2 for τ = TBMS/BMU FIG. 2. Plot of (cid:104) σ x (cid:105) / τ = T = π Ω (yellow) and the BMS/BMU solution (red)to the exact analytical solution (blue). Here, CG1 and CG2as well as BMS and BMU solutions are identical, whereasDCG ( τ = t ) reproduces the exact solution. By construction,CG1/CG2 (yellow) and DCG solutions (blue) intersect at t = T (vertical dashed line). Parameters: Ω = 10 ∆, λ = ∆, Γ( ω ) = Γ ωe −| ω | /ω c with Γ = 0 . ω c = 20 ∆, β ∆ = 1, ρ = 1 / + σ x ). from the trivial structure of the model that BMS (17)and BMU (18) solutions must concide. In Fig. 2 onecan see that an additional driving modifies the decayof coherences by super-imposing additional oscillations.Apart from that, we see that the exact reduced dynam-ics cannot be reproduced by a single Markovian gener-ator (like CG1, CG2, BMS or BMU). However, sincepure-dephasing-type models are a very specific modelclass, we consider models that allow the exchange ofenergy between system and reservoir below. IV. CIRCULAR DRIVING
As a second model, we consider a periodically drivenopen two level system of the form H ( t ) = ∆2 σ z + P σ + e − iΩ t + P ∗ σ − e +iΩ t + (cid:88) k ( σ + h k b k + σ − h ∗ k b † k ) + (cid:88) k ω k b † k b k , (22)where P denotes a driving amplitude. In absence ofdriving and for a vacuum reservoir, the model is ex-actly solvable [7], but we are here interested in the caseof finite driving amplitudes but weak system-reservoircoupling. By employing the unitary transform V ( t ) = exp (cid:34) − i (cid:32) Ω t σ z + Ω t (cid:88) k b † k b k (cid:33)(cid:35) (23) one can see from the observations V † ( t ) σ z V ( t ) = σ z , V † ( t ) σ ± V ( t ) = σ ± e ± iΩ t ,V † ( t ) b k V ( t ) = e − iΩ t b k , V † ( t ) b † k V ( t ) = e +iΩ t b † k , (24)that the total Hamiltonian becomes time-independentunder the gauge transformation V † ( t ) H ( t ) V ( t ) = ∆2 σ z + P σ + + P ∗ σ − (25)+ (cid:88) k ( σ + h k b k + σ − h ∗ k b † k ) + (cid:88) k ω k b † k b k . Accordingly, by transforming into the correspondingpicture, the expectation value of an observable O canbe written as (cid:104) O (cid:105) = Tr (cid:8) OU ( t ) ρ U † ( t ) (cid:9) ≡ Tr (cid:110) ˜ O ( t ) ˜ U ( t ) ρ ˜ U † ( t ) (cid:111) (26)with the effective time evolution operator˜ U ( t ) = V † ( t ) U ( t ) and transformed observable˜ O ( t ) = V † ( t ) OV ( t ). The effective time evolutionoperator evolves according to a time-independentpicture ddt ˜ U ( t ) = − i ˜ H ˜ U ( t ) , (27)˜ H = V † ( t ) H ( t ) V ( t ) + i ˙ V V = ∆ − Ω2 σ z + P σ + + P ∗ σ − + (cid:88) k ( σ + h k b k + σ − h ∗ k b † k ) + (cid:88) k ( ω k − Ω) b † k b k , such that we can write ˜ U ( t ) = e − i ˜ Ht . A. Benchmark solution
Therefore, under the gauge transform (25), the timeevolution is just that of a time-independent open two-level system. In this picture, we can set up the Born-Markov secular approximations and derive the BMSmaster equation for an undriven system. Formally,this looks like Eq. (17) with keeping only the n = 0term and instead of the Floquet Hamiltonian eigenstatesand eigenvalues using the eigenstates and eigenvalues of˜ H S = ∆ − Ω2 σ z + P σ + + P ∗ σ − . Due to a shifted KMS re-lation, this master equation does not thermalize even inthe time-independent picture but reaches some nonequi-librium steady state, which – after transforming back tothe original picture – exhibits a time-dependence withthe period of the driving.In the long-term limit, we can apply an even sim-pler analysis: For non-degenerate eigenvalues of theeffective system Hamiltonian ˜ H S , the populations de-couple from the coherences and the evolution is givenby a simple rate equation: ˙˜ ρ S,aa = (cid:80) b γ ab,ab ˜ ρ S,bb − (cid:80) b γ ba,ba ˜ ρ S,aa with transition rates from energy level b to a given by γ ab,ab . The steady state solution is thengiven by ¯˜ ρ S = P − |−(cid:105) (cid:104)−| + (1 − P − ) | + (cid:105) (cid:104) + | ,P − = γ − + , − + γ − + , − + + γ + − , + − , (28)where |−(cid:105) and | + (cid:105) denote the ground and excited stateof ˜ H S , respectively. For our example the two transitionrates read γ − + , − + = ˜ γ ( E + − E − ) | (cid:104)−| σ − | + (cid:105) | + ˜ γ ( E + − E − ) | (cid:104)−| σ + | + (cid:105) | ,γ + − , + − = ˜ γ ( E − − E + ) | (cid:104) + | σ − |−(cid:105) | + ˜ γ ( E − − E + ) | (cid:104) + | σ + |−(cid:105) | , (29)with ˜ γ ( ω ) = Θ( ω + Ω)Γ( ω + Ω)[1 + n B ( ω + Ω)] and˜ γ ( ω ) = Θ( − ω + Ω)Γ( − ω + Ω) n B ( − ω + Ω) exemplify-ing the broken Kubo-Martin-Schwinger [49–51] relation ˜ γ ( − ω )˜ γ (+ ω ) = e − β ( ω +Ω) . In this picture, system observ-ables can be computed via (cid:104) O (cid:105) → Tr S (cid:110) e +iΩ t/ σ z Oe − iΩ t/ σ z ˜ ρ S ( t ) (cid:111) , (30)where for large times we can insert ˜ ρ ( t ) → ¯˜ ρ S as givenby Eq. (28).The blue curves in Figs. 3 and 4 provide the full time-dependent BMS solution in the time-independent pic-ture, which in the long-term limit coincides with thelimit discussed above. B. Master equation solutions
We now go back to the original Schr¨odinger picturerepresentation of Eq. (22), where we can identify sys-tem coupling operators A = σ + and A = σ − . Ad-ditionally, the non-vanishing correlation functions be-come γ ( ω ) = Θ( ω )Γ( ω )[1 + n B ( ω )] and γ ( ω ) =Θ( − ω )Γ( − ω ) n B ( − ω ). Furthermore, one finds thatkick-operator and Floquet Hamiltonian are given by¯ H = ∆ − Ω2 σ z + P σ + + P ∗ σ − − Ω2 ,U kick ( t,
0) = exp (cid:18) − i Ω t + σ z ) (cid:19) . (31)This allows us to compute the coarse-graining dissipa-tors CG1, CG2, and DCG as well as BMS and BMUdissipators. C. Comparison
We compare the resulting dynamics with the previ-ously discussed asymptotic solution in Figs. 3 and 4. dimensionless time ω t -0,4-0,200,20,40,6 e xp ec t a ti on v a l u e < σ x > benchmark solutionDCG : τ = tCG 1 : τ = TCG 2 : τ = TBMSBMU
82 85 88 91 94-0,08-0,0400,040,080 3,1 6,3 9,4-0,200,20,4
FIG. 3. Plot of the time-dependent expectation value ofPauli matrix σ x according to the asymptotic solution in thetime independent frame from Eq. (30) (blue), CG1 (orange),CG2 (green), BMS/BMU (red/red dashed) and DCG (pur-ple crosses). The gridlines mark multiples of the driving pe-riod T . At long times, all DCG, BMS, and BMU solutionsall show a perfect numerical agreement with the long-termbenchmark solution (right inset). For short times, the DCGsolution agrees better with the benchmark than the BMSand BMU solutions, which also differ from each other then(left inset). Parameters: Ω = 2 ∆, P = ∆ (with Floquetenergies ( ¯ E ± = ± −√ ∆) in the first Brillouin zone), β ∆ =0 .
1, Γ = 0 . ω c = 15 ∆, ρ S = 1 / + 0 . σ x + 0 . σ z ]. dimensionless time ω t -0,4-0,200,2 e xp ec t a ti on v a l u e < σ z > benchmark solutionDCGCG 1 : τ = TCG 2 : τ = TBMSBMU
82 85 88 91 94-0,04-0,0200 3,1 6,3 9,4-0,2-0,100,10,20,30,4
FIG. 4. Plot of the time-dependent expectation value ofPauli matrix σ z . By construction, for short times the DCGsolution is superior to the other perturbative methods, butagain DCG, BMS, and BMU solution converge to the steady-state long-term solution for large times. Color coding andparameters identical to Fig. 3. Again, we see a convincing agreement of the DCG ap-proach with the steady-state benchmark solution (blue)in the long-time limit, but in this limit, the computa-tionally much simpler BMS and BMU approaches per-form equally well. Here one finds strong differences onlytransiently, where for short times the DCG approach byconstruction will approximate the exact solution andmust therefore also be close to the benchmark solutionfor the chosen parameters.
V. FAST-DRIVING SOLUTION
As a last example, we consider systems where thedriving breaks the pure-dephasing character of the sys-tem H ( t ) = H S + λ cos(Ω t ) C + A ⊗ (cid:88) k ( h k b k + h ∗ k b † k )+ (cid:88) k ω k b † k b k , (32)with arbitrary system Hamiltonian H S and system cou-pling operator A = A † fulfilling [ H S , A ] = 0 and systemdriving operator C = C † with in general [ H S , C ] (cid:54) = 0.Under a naive RWA approximation, the contribution ofthe driving vanishes, and the exact solution from theprevious section in absence of driving ( λ = 0) appliesfor all values of the system-reservoir coupling strength.To go beyond the naive RWA, we transform into an in-teraction picture with respect to the driving U ( t ) = e − i λ/ Ω sin(Ω t ) C , (33)and in this picture (marked by a tilde e.g. via ˜ A ( t ) = U † ( t ) AU ( t )) the transformed Hamiltonian assumes atime-dependent pure-dephasing form˜ H ( t ) = ˜ H S ( t ) + ˜ A ( t ) ⊗ (cid:88) k ( h k b k + h ∗ k b † k ) + (cid:88) k ω k b † k b k , (34)where we still have (cid:104) ˜ H S ( t ) , ˜ A ( t ) (cid:105) = U † ( t ) (cid:2) H S , A (cid:3) U ( t ) = 0. Nevertheless, we cannotapply the naive polaron treatment of App. A, sincenow also the coupling operator has picked up a periodictime-dependence. However, depending on the problemone may have the situation that also the time averagesof the operators¯ H S ≡ T (cid:90) T ˜ H S ( t ) dt , ¯ A ≡ T (cid:90) T ˜ A ( t ) dt (35)commute with each other [ ¯ H S , ¯ A ] = 0, and in this casethe Hamiltonian after the RWA in the interaction pic-ture¯ H tot = ¯ H S + ¯ A ⊗ (cid:88) k ( h k b k + h ∗ k b † k ) + (cid:88) k ω k b † k b k (36)is of pure-dephasing type and can be solved with stan-dard methods, see App. B for details. The transfor-mation into this frame is of course not equivalent tothe interaction picture representation. However, underthe RWA in this transformed frame, the time evolutionoperator of the system becomes U S ( t ) ≈ U ( t ) e − i ¯ H S t ,such that we can identify approximations to kick oper-ator U kick ( t ) ≈ U ( t ) and system Floquet Hamiltonian¯ H ≈ ¯ H S , respectively. A. Fast driving benchmark
Specifically, for a two level system with H S = ∆2 σ z , A = σ z , C = σ x , (37)we obtain with the Fourier decomposition˜ A ( t ) = U † ( t ) σ z U ( t )= σ z ∞ (cid:88) n = −∞ J n (cid:18) λ Ω (cid:19) e i2 n Ω t − i σ y ∞ (cid:88) n = −∞ J n +1 (cid:18) λ Ω (cid:19) e i(2 n +1)Ω t (38)that this method – under the RWA in the transformedframe – yields the following dynamics (cid:104) σ x (cid:105) = (cid:2) cos ( µ ∆ t ) (cid:104) σ x (cid:105) − sin ( µ ∆ t ) (cid:104) σ y (cid:105) (cid:3) Σ( t ) , (cid:104) σ y (cid:105) = − sin (cid:0) λ Ω sin(Ω t ) (cid:1) (cid:104) σ z (cid:105) + cos (cid:0) λ Ω sin(Ω t ) (cid:1) × (cid:2) sin ( µ ∆ t ) (cid:104) σ x (cid:105) + cos ( µ ∆ t ) (cid:104) σ y (cid:105) (cid:3) Σ( t ) , (cid:104) σ z (cid:105) = + cos (cid:0) λ Ω sin(Ω t ) (cid:1) (cid:104) σ z (cid:105) + sin (cid:0) λ Ω sin(Ω t ) (cid:1) × (cid:2) sin ( µ ∆ t ) (cid:104) σ x (cid:105) + cos ( µ ∆ t ) (cid:104) σ y (cid:105) (cid:3) Σ( t ) , (39)with Σ( t ) = exp (cid:40) − π ( µ ) ∞ (cid:82) −∞ γ ( ω ) sin ( ωt ) ω dω (cid:41) and µ = J (cid:0) λ Ω (cid:1) being given by the Bessel function of thefirst kind [52]. Details on this can be found in App. Bby inserting the expressions for a two level system inEq. (B4). In absence of of driving λ → µ → ρ ( t ) = 12 ( (cid:104) σ z (cid:105) + 1) → ρ , (40) ρ ( t ) = (cid:10) σ + (cid:11) → e +i∆ t e − π ∞ (cid:82) −∞ dωγ ( ω ) sin2 ( ωt ) ω ρ , as J (cid:0) λ Ω (cid:1) → λ →
0. Furthermore, at vanishingsystem-reservoir coupling we have Σ( t ) →
1, where onecan immediately verify that the purity of an initial stateis preserved. Finally, we can see that for reasonablespectral densities we find that lim t →∞ Σ( t ) = 0, suchthat the long-term asymptotics of (39) is given by (cid:104) σ y (cid:105) → − sin (cid:18) λ Ω sin(Ω t ) (cid:19) (cid:104) σ z (cid:105) , (cid:104) σ z (cid:105) → + cos (cid:18) λ Ω sin(Ω t ) (cid:19) (cid:104) σ z (cid:105) . (41)Importantly, we stress that the fast-driving solution isvalid also for stronger system-reservoir couplings [54]. dimensionless time ω t -0,2-0,100,10,2 e xp ec t a ti on v a l u e < σ x > DCGFD solutionCG 1 : τ = TCG 2 : τ = TBMSBMU
30 33 35 38-0,04-0,0200,020,04
FIG. 5. Plot of the time-dependent expectation value ofPauli matrix σ x comparing the analytical fast-driving solu-tion (blue), CG1 (yellow), CG2 (green), BMS (red), BMU(red dashed) and DCG (purple cross). The gridlines markmultiples of the driving period T . A higher order approx-imation for the driving yields similar results (not shown).Parameters: Ω = 25 ∆, λ = ∆, β ∆ = 1, Γ = 0 . ω c = 15 ∆, ρ = 1 / + 0 . σ x + 0 . σ z ). B. Master equation solutions
Since we are interested in the fast drivingregime here, we consider kick operator U kick ( t ) ≈ exp {− i λ/ Ω sin(Ω t ) σ x } and system Floquet Hamiltonian¯ H ≈ ∆2 J (cid:0) λ Ω (cid:1) σ z also by applying the RWA in thetransformed frame. Considering the scaling of theBessel functions at small arguments, we therefore con-sistently only keep the lowest order terms with n ∈{− , , +1 } in Eq. (38). Then, the coupling operatorin the interaction picture can be approximately writtenas A ( t ) ≈ J (cid:18) λ Ω (cid:19) σ z + 2 sin(Ω t ) J (cid:18) λ Ω (cid:19) e +i ¯ Ht σ y e − i ¯ Ht + O (cid:26) λ Ω (cid:27) , (42)which suffices to set up all master equations. We havenumerically checked convergence in the shown parame-ter regime by considering also the case n ∈ {− , . . . , +2 } (not shown). C. Comparison
In the limit of both fast driving and weak coupling wecan compare the benchmark solution with the pertur-bative CG1/CG2/DCG/BMS/BMU approaches, whichleads to Fig. 5 and Fig. 6.There, we can observe that the DCG approach per-forms best. It has however the disadvantage that thegenerator has to be re-computed for any desired time.Still, all perturbative approaches fail to correctly cap- dimensionless time ω t e xp ec t a ti on v a l u e < σ z > DCGFD solutionCG 1 : τ = TCG 2 : τ = TBMSBMU FIG. 6. Plot of the time-dependent expectation value ofPauli matrix σ z . Color coding and parameters analogous toFig. 5. ture the long-term dynamics of the fast-driving solu-tion (39), which is not surprising as the latter is non-perturbative in the system-reservoir coupling strength. VI. CONCLUSIONS
Our intention in this study was to find an improvedbut simple (Markovian) master equation applicable toperiodically driven systems. For this, we analyzeddriven qubit systems coupled to thermal reservoirs withcoarse-graining approaches. All methods used had theformal advantage of being in LGKS form.The negative result is that the analyzed coarse-graining approaches over one period of the driving CG1and CG2 did not match the expectations, i.e., they wereinferior to the DCG approach, even in the special caseof pure dephasing. On the positive side, we found thatthe BMS approach can yield quite reliable results in thelong term limit provided the secular approximation isperformed in a proper way. The DCG approach per-formed at least as good as the BMS variant in the long-term limit and approaches the exact solution by con-struction in the short-term limit, such that we considerit as the most accurate second-order perturbative so-lution among the schemes tested. Its drawback is thecomputational cost, since for each time a suitable dissi-pator needs to be calculated numerically. To avoid this,we note that by using spectral densities with a simplerpolynomial structure, the involved integrals in (14) canbe performed analytically. The fact that it is not possi-ble to capture the full dynamics faithfully with a singlecoarse-graining time could be taken as a hint that therestriction to LGKS dynamics is too severe such thatone should rather aim at deriving Kraus representationsfrom first principles.
ACKNOWLEDGMENTS
R.H. acknowledges financial support by the JapaneseStudent Services Organization JASSO. G.S. grate-fully acknowledges discussions with J. Ablaßmayer, T.Becker, A. Eckardt, and A. Schnell.
Appendix A: Exact pure-dephasing solution
For a driven Hamiltonian of pure dephasing form (19),the exact dynamics can be calculated by using a polaron(or Lang-Firsov [55]) transform: U p = e A ⊗ (cid:18) hkωk b k − h ∗ kωk b † k (cid:19) . (A1)Applying this transformation to system and bath oper-ators yields: U † p H S ( t ) U p = H S ( t ) , U † p AU p = A , (A2) U † p b k U p = b k − h ∗ k ω k A , U † p b † k U p = b † k − h k ω k A , which leads to decoupled system and bath Hamiltonians U † p H ( t ) U p = H S ( t ) − (cid:32)(cid:88) k | h k | ω k (cid:33) A + (cid:88) k ω k b † k b k ≡ ¯ H S ( t ) + (cid:88) k ω k b † k b k , (A3)with ¯ H S ( t ) denoting an effective system Hamiltonian.Hence, the time evolution operator in the polaron pic-ture ˜ U ( t ) is given by the system and bath evolutionseparately˜ U ( t ) = T (cid:110) e − i (cid:82) t ¯ H S ( t (cid:48) ) dt (cid:48) (cid:111) e − i (cid:80) k ω k b † k b k t , (A4)where T denotes the time-ordering operator of the sys-tem. For a general unitary transformation V ( t ), thetime evolution operator in the original frame U ( t ) canbe expressed in terms of the time evolution operator inthe new frame ˜ U ( t ) as follows: U ( t ) = V ( t ) ˜ U ( t ) V † (0) . (A5)Thus, for the polaron transform this relation yields: U ( t ) = U p ˜ U ( t ) U † p . (A6)From this, the expectation value of any system observ-able O S , (cid:104) O S (cid:105) = Tr { ( O S ⊗ B ) ρ ( t ) } = Tr (cid:8) U † ( t )( O S ⊗ B ) U ( t ) ρ (cid:9) , (A7)can be calculated as: (cid:104) O S (cid:105) = Tr (cid:110) U p ˜ U † ( t ) U † p ( O S ⊗ B ) U p ˜ U ( t ) U † p ρ (cid:111) . (A8) For the example of a two level system, H S ( t ) = ∆2 σ z , A = σ z , (A9)the populations ρ S, ( t ) = (cid:104) + σ z (cid:105) and ρ S, ( t ) = (cid:104) − σ z (cid:105) are constant as and σ z commute with U p and ˜ U ( t ): ρ S, ( t ) = ρ S, , ρ S, ( t ) = ρ S, . (A10)For the coherences ρ = (cid:104) σ − (cid:105) and ρ = (cid:104) σ + (cid:105) , thecalculations are more involved. Using that U † p σ ± U p = e ± (cid:80) k (cid:18) h ∗ kωk b † k − hkωk b k (cid:19) σ ± , (A11)and organising bath and system parts together, one gets (cid:10) σ ± (cid:11) = Tr (cid:110) U p U † B ( t ) e ± (cid:80) k (cid:18) h ∗ kωk b † k − hkωk b k (cid:19) U B ( t ) U † S ( t ) σ ± × U S ( t ) U † p ρ (cid:111) , (A12)with U B ( t ) = e − i (cid:80) k ω k b † k b k t and U S ( t ) = e − i ( ∆2 t + λ Ω sin(Ω t ) ) σ z . Now, one can use that e i (cid:80) k ω k b † k b k t e ± (cid:80) k (cid:18) h ∗ kωk b † k − hkωk b k (cid:19) e − i (cid:80) k ω k b † k b k t = e ± (cid:80) k (cid:18) h ∗ kωk b † k e i ωkt − hkωk b k e − i ωkt (cid:19) ,e i xσ z σ ± e − i xσ z = e ± x σ ± , (A13)which leads by inserting the identity U † p U p = to: (cid:10) σ ± (cid:11) = e ± ( ∆2 t + λ Ω sin(Ω t ) ) (A14) × Tr (cid:40) U p e ± (cid:80) k (cid:18) h ∗ kωk b † k e i ωkt − h . c . (cid:19) U † p U p σ ± U † p ρ (cid:41) . Next, the polaron transformation is applied to the sys-tem and bath parts separately: U p e ± (cid:80) k h ∗ kωk b † k e i ωkt − h . c . U † p = e ± (cid:80) k h ∗ kωk (cid:16) b † k + hkωk σ z (cid:17) e i ωkt − h . c . ,U p σ ± U † p = e ∓ (cid:80) k (cid:18) h ∗ kωk b † k − hkωk b k (cid:19) σ ± . (A15)Separating system and bath parts and inserting the ini-tial condition ρ = ρ S ⊗ ¯ ρ B yields (cid:10) σ ± (cid:11) = e ± ( ∆2 t + λ Ω sin(Ω t ) ) × Tr (cid:40) e ± (cid:80) k | hk | ω k sin( ω k t ) σ z σ ± ρ S (cid:41) B ± ( t ) ,B ± ( t ) = Tr (cid:110) e ± (cid:80) k (cid:18) h ∗ kωk b † k e i ωkt − hkωk b k e − i ωkt (cid:19) × e ∓ (cid:80) k (cid:18) h ∗ kωk b † k − hkωk b k (cid:19) ¯ ρ B (cid:111) . (A16)0The trace over the bath parts B ± ( t ) gives: B ± ( t ) = exp (cid:40) ∓ (cid:88) k | h k | ω k sin( ω k t ) (cid:41) (A17) × exp (cid:40) − (cid:88) k | h k | ω k [1 − cos( ω k t )][1 + 2 n B ( ω k )] (cid:41) . Here, n B ( ω k ) = e βωk − denotes the Bose distribution.Plugging the expression for B ± ( t ) into Eq. (A16) andusingTr (cid:40) e ± (cid:80) k | hk | ω k sin( ω k t ) σ z σ ± ρ S (cid:41) = e ± i (cid:80) k | hk | ω k sin( ω k t ) × Tr (cid:8) σ ± ρ S (cid:9) , (A18)eventually yields for the expectation value of σ ± : (cid:10) σ ± (cid:11) ( t ) = e ± ( ∆2 t + λ Ω sin(Ω t ) ) (A19) × e − π ∞ (cid:82) dω Γ( ω ) sin2 ( ωt ) ω [1+2 n B ( ω )] (cid:10) σ ± (cid:11) , where the spectral coupling density Γ( ω ) =2 π (cid:80) k | h k | δ ( ω − ω k ) has been inserted and it hasbeen used that 1 − cos( ωt ) = 2 sin (cid:0) ωt (cid:1) . Appendix B: Analytic fast-driving approximation
We consider the simplified form of Eq. (32) with H S = ∆2 A . First, we move into an interaction picture withrespect to the driving by performing the transformation U ( t ) = e − i λ Ω sin(Ω t ) C , which yields after applying therotating wave approximation:˜ H RWA = ˜ A (cid:32) ∆2 + (cid:88) k ( h k b k + h ∗ k b † k ) (cid:33) + (cid:88) k ω k b † k b k , ˜ A = Ω2 π (cid:90) π Ω U † ( t ) AU ( t ) dt . (B1) This is the Hamiltonian of a simple pure dephasing sys-tem, see Eq. (19), with time independent system Hamil-tonian.Applying the polaron transform (A1), system andbath part of the Hamiltonian can be decoupled and (A3)becomes H p = U † p ˜ H RWA U p = ∆2 ˜ A − (cid:88) k | h k | ω k ˜ A + (cid:88) k h k b † k b k , (B2)which yields for the time evolution operator in this pic-ture:˜ U ( t ) ≈ e − i (cid:32) ∆2 ˜ A − (cid:80) k | h k | ω k ˜ A (cid:33) t e − i (cid:80) k h k b † k b k t . (B3)In an anologous way as for the pure dephasing system,we can write the expectation value of an observable as: (cid:104) O S (cid:105) = Tr (cid:110) U p ˜ U † ( t ) U † p U † ( O S ⊗ B ) U U p ˜ U ( t ) U † p ρ (cid:111) . (B4)For a two level system Eq. (37), we get˜ A = J (cid:0) λ Ω (cid:1) σ z , (B5)with J n ( x ) denoting the Bessel function of firstkind [56], which is defined as the solution of the differ-ential equation x J (cid:48)(cid:48) n ( x ) + x J (cid:48) n ( x ) − ( z + n ) J n ( x ) = 0. [1] G. Floquet, Annales scientifiques de l’´Ecole NormaleSup´erieure , 47 (1883).[2] J. H. Shirley, Phys. Rev. , B979 (1965).[3] M. Grifoni and P. H¨anggi, Physics Reports , 229(1998).[4] A. Eckardt and E. Anisimovas, New Journal of Physics , 093039 (2015).[5] M. Bukov, L. D’Alessio, and A. Polkovnikov, Advancesin Physics , 139 (2015).[6] U. Weiss, Quantum Dissipative Systems , Series of Mod-ern Condensed Matter Physics, Vol. 2 (World Scientific,Singapore, 1993). [7] H.-P. Breuer and F. Petruccione,
The Theory of OpenQuantum Systems (Oxford University Press, 2002).[8] M. Schlosshauer,
Decoherence and the Quantum-To-Classical Transition (Springer-Verlag Berlin Heidelberg,2007).[9] S. Chu, Science , 861 (1991).[10] S. Rice and M. Zhao,