Floquet stroboscopic divisibility in non-Markovian dynamics

We provide a general discussion of the Liouvillian spectrum for a system coupled to a non-Markovian bath using Floquet theory. This approach is suitable when the system is described by a time-convolutionless master equation with time-periodic rates. Surprisingly, the periodic nature of rates allow us to have a stroboscopic divisible dynamical map at discrete times, which we refer to as Floquet stroboscopic divisibility. We illustrate the general theory for a Schr\"odinger cat which is roaming inside a non-Markovian bath, and demonstrate the appearance of stroboscopic revival of the cat at later time after its death. Our theory may have profound implications in entropy production in non-equilibrium systems.

We provide a general discussion of the Liouvillian spectrum for a system coupled to a non-Markovian bath using Floquet theory. This approach is suitable when the system is described by a time-convolutionless master equation with time-periodic rates. Surprisingly, the periodic nature of rates allow us to have a stroboscopic divisible dynamical map at discrete times, which we refer to as Floquet stroboscopic divisibility. We illustrate the general theory for a Schrödinger cat which is roaming inside a non-Markovian bath, and demonstrate the appearance of stroboscopic revival of the cat at later time after its death. Our theory may have profound implications in entropy production in non-equilibrium systems. Introduction:-No single system is completely isolated from environment. Hence, an appropriate understanding and modeling of system-environment interaction is necessary to effectively predict the reduced system dynamics. In this way, the memory effect of the bath in open quantum systems is a fundamental and important concept. When a quantum system is in contact with a Markovian (memory-less) bath, the information flows from the system to the bath. However, if the bath has memory, the situation changes dramatically and the dynamics is known as non-Markovian. Recently, the dynamics of systems coupled to non-Markovian reservoirs has been the focus of both theoretical as well as experimental investigations [1][2][3][4][5][6][7][8]. The state of the art quantum technologies allow one to manipulate quantum systems with unprecedented precision and control. In the context of reservoirs engineering [9], for instance, one can manipulate the spectral density of the bath and consider structured reservoirs [8,10,11] that exhibit non-Markovian dynamics [12][13][14]. In the context of quantum simulation [15], there are recent experimental realizations of dynamical maps both in the Markovian regime by using trapped ions [16], as well as in the non-Markovian regime [17] via an all-optical setup. In optomechanical systems, non-Markovian behavior has also been observed experimentally [18].
In the case of Markovian dynamics, the master equation for the reduced density matrix of the system can be written formally in terms of the Liouvillian opera- * Electronic address: cqtvmbv@nus.edu.sg † Electronic address: dimitris.angelakis@qubit.org tor (LO), which is an operator acting on density matrices [19][20][21]. In most of the cases, the LO is timeindependent or has an adiabatic dependence on time, which enables one to find its spectrum. This is of utmost importance because the imaginary parts of the Liouvillian spectrum are related to coherent evolution and the real parts to decay and decoherence. In addition, the kernel of the LO determines the steady state of the system. In the context of non-equilibrium quantum phase transitions, the Liouvillian spectrum plays a fundamental role because by defining the gap between the real parts of the eigenvalues of the LO, one can extend the definition of quantum critical points under the effect of dissipation [22]. Spectral properties of the Liouvillian also determine the onset of bistability in driven-dissipative resonators [23]. For a long time, the theoretical understanding of time-dependent Liouvillians is an open problem [24,25]. These kind of Liouvillians appear naturally in the context of pure-dephasing channels [26,27] and lead to time-convolutionless master equations, which can be non-Markovian [28]. A recent work has shown that in the case of a Markovian master equation with timeperiodic Liouvillian, one can use Floquet theory [29][30][31] to obtain the asymptotic steady state [32].
In this letter, we use Floquet theory to define the Liouvillian spectrum for non-Markovian evolution. The latter is generated through a time-convolutionless master equation, i.e., a Lindblad-type master equation with time-periodic rates that can be negative in certain intervals. In addition, we prove that the dynamical map is divisible at stroboscopic times, which we term as "Floquet stroboscopic divisibility". We illustrate the general theory by considering a quantum harmonic oscillator coupled to a structured bath with an engineered spectral density that defines time-periodic rates. Due to the non-Markovian nature of the bath, we observe the preservation of the Schrödinger cat coherences at stroboscopic times.
Floquet stroboscopic divisibility:-To make a direct connection between the dynamics of an open quantum system and Floquet theory, we consider a Lindblad master equation dρS(t) dt =L(t)ρ S (t) [19,20,32] with time periodic Liouvillian operatorL(t + T ) =L(t). Here, ρ S (t) denotes the reduced density matrix of the system. The master equation is defined via the Liouvillian operatorL(t)(·) = −i[Ĥ S (t), (·)] + l γ l (t)D(Ô l )(·), whereĤ S (t) is the time-periodic Hamiltonian of the system, γ l (t) are time-periodic rates, andD(Ô l )(·) = O l (·)Ô † l − 1 2 {Ô † lÔ l , (·)} are the dissipators. The Floquet theorem (see Supplemental Material [33]) ensures that there exists a dynamical mapΦ(t)-or fundamental matrix-with the formΦ(t) =P (t) exp(L F t), wherê P (t + T ) =P (t) [29,30]. In addition, the dynamical map is the propagator for the density matrix, i.e., ρ S (t) =Φ(t)ρ S (0). The eigenvalues λ α = e LαT of the matrixΦ(T ) are the characteristic multipliers and the eigenvalues {L α } ofL F are the Floquet exponents. As the evolution of the system is non-unitary, the Floquet exponents {L α } are in general complex numbers. The real part of the Floquet exponents are called Lyapunov exponents and they are responsible for decoherence and relaxation. On the other hand, their imaginary components govern the coherent dynamics. Now we have all the elements we need to define stroboscopic Floquet stroboscopic divisibility, which is a fundamental result of our work. This concept is a consequence of Floquet theorem because at stroboscopic times, we haveΦ(mT ) = [Φ(T )] m , which means the dynamical map is divisible. So far, we have not touched upon the Markovian or non-Markovian nature of the dynamical map yet. Furthermore, the Floquet theorem establishes that stable solutions are possible when the Lyapunov exponents are smaller or equal than zero [29,30]. That implies the constraint | detΦ(T )| ≤ 1, which can be derived from the general formula detΦ(T ) = exp{´T 0 Tr[L(τ )]dτ } (see Supplemental Material [33] and Refs. [29,30]). We are interested in the case where the time average of all the rates in one period is positive or zero, in order to satisfy the stability constraint discussed above. We also note that the previous statement does not restrict the rates to be positive values at all times. Based on the geometrical characterization of non-Markovianity [34], the absolute value of detΦ(T ) can be reinterpreted as the volume of the accessible states at time T . Here, the dynamical evolution we consider is discrete due to the stroboscopic nature of the Floquet theorem. A dynamical map is divisible if the rate of change of the volume is smaller than zero [35,36]. In our case, the dynamical map is stroboscopically divisible if the finite difference Depicts a device to demonstrate pendulum waves. In this mechanical device, the system comes back to its initial configuration after one time period T . (b) Quantum evolution of the mean photon number N k (t) = a † k a k (t) of N = 60 modes of the non-Markovian bath we consider in the manuscript (the density plot depicts log N k ). Similarly to the pendulum waves, at a time T , the whole system comes back to its initial configuration. For the coupling g k = he −zk to the modes of the bath we used z = 0.1 and hT = 1.0. For convenience, we consider a zero temperature bath at the initial time with N k (0) = 0 and we prepare the resonator in a cat state Interestingly enough, ∆ m goes to zero when m goes to infinity, and for finite m, when | detΦ(T )| = 1. In the former case, this implies that the system reaches a periodic asymptotic state and the latter case means that the system is purified stroboscopically. In contrary to the results presented in Refs. [34][35][36], ∆ m is a measure of how the volume of accessible states is contracted stroboscopically. The stroboscopic divisibility can also be explained by observing properties of the instantaneous spectra of dynamical maps when the rates are time-dependent [35,36]. In the following, we will apply the general theory presented so far to a simple example: a resonator couples to a non-Markovian bath. Non-Markovian dynamics of a resonator and the dephasing environment:-We begin by considering a timeconvolutionless master equation for a resonator. This master equation is obtained from an exact microscopic derivation for pure dephasing type coupling (see Supplemental Material [33]) wheren = b † b and the coherent evolution is governed byĤ S (t) = ωn − g(t)n 2 . This solution was inspired by previous works on phase damping [37] and dynamics of cavities coupled to moving mirrors [38]. The dephasing rate is given by where β is the inverse temperature. Without loss of generality, we consider a zero-temperature bath throughout the paper. We also consider a structured bath whose spectral density has peaks at frequencies ω k = k s Ω, where s is a positive integer, and Ω = 2π/T . Besides its dissipative effect, the bath also has an effect on the coherent dynamics of the system via the nonlinear driving For the purposes of this work, the bath frequencies are chosen as ω k = kΩ (s = 1) with coupling strengths g k = he −zk/2 , where z > 0 is a real number determining how many modes couple to the resonator. We note that our results are valid for any value of s, and we take s = 1 case for simplicity. In the Supplemental Material [33], we present the exact form of the functions g(t) and γ(t) in the case of a bath with infinite number of modes and propose an implementation of the system in cirquit QED. As we have complete knowledge over the density matrix of the total system, we can also get information of the dynamics of the modes of the bath. In particular, one can derive very simple expressions for the mean photon number of any tal Material [33]). Figure 1(a) depicts a mechanical analogue of the bath we are considering in the manuscript, which is referred to as pendulum waves [39]. One can prepare the ensemble of oscillators in a given configuration and after some time T it will be back to the initial configuration. Similarly, Figure 1(b) shows the dynamics of the bath with N = 60 modes. In this case, the period T = 2π/Ω is determined by the fundamental frequency Ω and one can see that the dynamics of mean photon number of the modes N k (t) is reversed at time t = T /2, exactly as in the mechanical pendulum waves.
Our choice for the frequencies of the bath has dramatic consequences for the time evolution of the total system. In particular, from the expressions for the dephasing rate γ(t) and the bath-induced nonlinearity g(t), we find that these functions turn out to be periodic with period T = 2π/Ω. Besides this, the integral of the dephasing rate over one period is zero. An immediate consequence of this is that when the rates are positive, there is dephasing of the resonator. As the average of the rates in one period is zero, the rates are also negative and there are intervals of time where the coherences are build up again. In short, due to the choice of the frequencies, the information flow between the system and the bath is periodic and the dynamics is non-Markovian.
Properties of the dynamical map:-One advantage of the exact nature of the master equation (2) is that the resulting dynamical map is valid for any strength of the system-bath coupling and for arbitrary spectral densities of the bath. For our choice of the bath frequencies, the Liouvillian is periodic and´T 0 Tr[L(τ )]dτ = 0, which implies a maximally non-Markovian dynamics of the bath. In addition, the dynamical mapΦ(T ) is diagonal and has eigenvalues λ m,n = e −i(Em−En)T e iG(T )(m 2 −n 2 ) , where E n = nω, G(t) is the integral of the function g(t), and  [33]). From this expression one can extract the corresponding Floquet exponents, which are imaginary. As the consequence, we have detΦ(T ) = 1, from which we can see that the system is stroboscopically purified at times t l = lT with l a positive integer. The physical interpretation is that the information trade-off between the system and the bath is balanced. In other words, the information that goes away from the system is recovered after one period T .
The property detΦ(T ) = 1 of the dynamical map has interesting features, because the volume of accessible states [35,36] is conserved stroboscopically. As the Floquet-Liouvillian spectrum is imaginary, it implies a symmetry Im[L (m,n) ] = −Im[L (n,m) ] of imaginary part of the Floquet exponents In the strong coupling regime (hT = 1.0), we plot the imaginary part of the Floquet-Liouville spectrum in Fig. 2(a) and the characteristic multipliers in Fig. 2(b). The bath induces a nonlinearity proportional to G(t) in the coherent evolution of the system. This has an influence in the spectrum because it is not harmonic and its nonlinear behavior is proportional to G(T ), which ex-  Fig. 2(c). This is related to the clustering of levels appearing in the Fig. 2(d).
Return probability and linear entropy:-As we have previously discussed, the dynamics of the bath is periodic and this is related to the time-periodicity of the dephasing rates in the master equation (2). Furthermore, we have also discussed that the bath introduces a time dependent nonlinearity proportional to g(t)n 2 which influences the coherent evolution of the resonator. The interplay between these two effects has nontrivial physical consequences. For example, one can ask the question: Does the system comes back to the initial configuration as the bath does? To measure this, one can calculate the density-density correlation R(t) = tr S [ρ S (t)ρ S (0)] which is also known as the return probability. We have calculated this and the results are shown in Fig. 3(a). One can observe that the system does not comeback to the initial configuration, but it is close to it at a time t = 3T . This is due to the nonlinearity induced from the bath.
As a measurement of bipartite entanglement between the system and the environment, we consider the linear entropy L(t) = 1 − tr S [ρ 2 S (t)]. Figure 3(b) shows that the system is purified in a stroboscopic manner at times t l = lT with integer l > 0. This reinforces our previous discussion because the dynamical map creates mixing of the state of the system when the rates are positive and has the contrary effect when the rates are negative. As the time-average of the rates is zero in a period T , the information lost is retrieved in the system. To visualize the nonlinearity due to the coupling to the bath, we also calculate the Wigner function W (Q, P ) = 1 π Tr ΠD † (α)ρ S (t)D(α) of the resonator, whereD(α) andΠ are displacement and parity operators, respectively [40]. By using the canonical coordinates Q and P one can define α = 1 √ 2 (Q + iP ). The stroboscopic dynamics of the Wigner function is depicted in Figs. 3(c),(d),(e), and (f). After one period of evolution, the Wigner function reveals signatures of the nonlinearity, as the system is not anymore in a cat state. Strikingly, after three periods of the evolution, the systems refocuses in a configuration which is close to the initial state. This behavior is intimately related to the entanglement revival in non-Markovian evolution [6].
Conclusions and outlook:-We have investigated the Liouvillian spectrum of a non-Markovian master equation which is local in time. Based on Floquet theory, we have shown that even though the dynamics is fully non-Markovian, the dynamical map is divisible at stroboscopic times. In addition, we have proven that spectral properties of the Liouvillian determine the contraction of the volume of accessible states at stroboscopic times. To substantiate the general theory, we present a timeconvolutionless master equation derived microscopically for a maximally non-Markovian bath. We show that in this example, the volume of the accessible states [34] is stroboscopically conserved, because detΦ(T ) = 1, which implies that the system is purified at stroboscopic times. We also show the intriguing relation between the mechanical pendulum waves and the dynamics of the coloured bath of our choice. In this particular example, to our surprise, we find the revival of a Schrödinger cat at later time after its death, prompting us to question whether it is possible to conserve the system entropy in periodic manner or even to decrease the system entropy stroboscopically. Possible directions in the future include the theoretical investigation of environments that exhibit phase transitions [41].
Note Added: While this work was near completion, we became aware of an interesting work [42]. The authors define the Liouvillian spectrum for a master equation with a memory kernel which is periodic in time, and study for the asymptotic system behaviour, which is different from what we have presented here.
Acknowledgements:-V.M.B. and T.H.K. acknowledge fruitful discussions with A. Chia and E. Munro, S. Restrepo, and S. Vinjanampathy. The authors are grateful for the financial support through the National Research Foundation and Ministry of Education Singapore (partly through the Tier 3 Grant "Random numbers from quantum processes"); and travel support by the EU IP-SIQS. Equations in the main paper are denoted by Eq. [*].

A. Floquet theory for open quantum systems
In the systems of linear ordinary differential equations with periodic coefficients, there is a very powerful tool, which is referred to as the Floquet theorem [29,30]. In this part of the supplemental material, we discuss briefly the Floquet theorem and its application in the theory of open quantum systems [31,32].

A quick review of the Floquet theorem
Let us consider a linear system of M coupled ordinary differential equations dX(t) is a vector and A(t + T ) = A(t) is a M times M matrix, which is periodic in time [29,30]. In general, the solutions of this system of differential equations are not periodic, but the Floquet theorem enables us to obtain important information about their behaviors. One can define a square matrix M(t) whose columns are M solutions X 1 (t), . . . , X M (t) of the system of differential equations. This matrix is called the fundamental matrix. In particular, it can be shown that M(t + which is derived from the Liouville-Jacobi formula [30]. The complex numbers e C l T are the eigenvalues of B = e CT , which are also referred to as characteristic multipliers. The arguments C l of the exponential are called Floquet exponents. Given that M(0) = I, the fundamental matrix evaluated at one period M(T ) = B = e CT is known as the monodromy matrix [30]. With these tools we can understand the main features of the Floquet theorem. This theorem establishes that the fundamental matrix of the system of differential equations with periodic coefficients can be written as M(t) = P(t)e Ct , where P(t + T ) = P(t). Due to the form of the monodromy matrix, the Floquet exponents are not uniquely defined, i.e., if C l is a Floquet exponent, C l + 2πni/T with integer n, is also a Floquet exponent. The fundamental matrix M(t) also plays the role of a propagator in physical applications. This means that given an initial condition X(0), the solution at a later time can be written as such that X l (t + T ) = X l (t). The vectors X l (T ) are eigenvectors of the monodromy matrix B = e CT with eigenvalues e C l t . The coefficients a l = X T l (T ) · X(0) are the overlap between the eigenmodes X l and the initial condition X(0).

Floquet theorem in open quantum systems
Given a master equation of the form dρS(t) dt =L(t)ρ S (t) [19][20][21], whereL(t + T ) =L(t) is the Liouvillian, we can apply the Floquet theorem for ordinary differential equations with periodic coefficients [29,30]. Based on the previous discussion, given a basis for the Hilbert space, one can always construct a matrix representation of the Liouvillian operator. In this basis, the master equation turns out to be just a system of coupled ordinary differential equations with periodic coefficients and one can translate one to one all the elements presented above. For example, if one represents the density matrixρ S (t) as a vector, it will play the role of X(t) and the matrix representation of the Liouvillian will be time-periodic matrix A(t).
The Floquet theorem ensures that the propagator-or dynamical map-has the formΦ(t) =P (t) exp(L F t), wherê P (t+T ) =P (t). To make the connection with the theory presented above, the matrix representation of the dynamical mapΦ(t) is the fundamental matrix M(t). Correspondingly, the matrix representation of the operator exp(L F T ) plays the role of the monodromy matrix B = e CT . The Liouville-Jacobi formula implies that the determinant of the dynamical map satisfies The discussion presented in the main text was based on spectral properties of the dynamical map, but the Floquet theorem gives us more information. For example, the solution of the master equation can be written asρ S (t) = α c α e Lαtρ α (t), whereρ α (t + T ) =ρ α (t) andΦ(T )ρ α (T ) = e LαTρ α (T ). One should also take into account that the Floquet exponents are not uniquely defined because one can always add to them a complex phase 2πni/T with integer n such that one gets the same characteristic multiplier, i.e, e (Lα+2πni/T )T = e LαT .

B. Microscopic derivation of the master equation
We consider the following Hamiltonian in the example presented in the main text (onwardsh = 1.)- whereĤ B = N k=1 ω kâ † kâ k the bath Hamiltonian, and the system-bath couplingĤ SB =n N k=1 g k â † k +â k . The effect of the bath on the system is to create pure dephasing [26,27]. A related problem was discussed in the context on phase damping [37] and dynamics of cavities coupled to moving mirrors [38].
The system-bath evolution at time t is obtained bŷ One and only assumption we make in the entire derivation is that we assumeρ(0) =ρ S (0) ⊗ρ B (0), whereρ B = e −βĤB /Z with Z = tr(e −βĤB ) is a thermal state of the bath. A direct consequence from the above assumption is that we can write the reduced system density matrix aŝ Now, the entire problem reduces to finding the influence functional F nm (t) analytically [27]. Its expression is whereĤ Here, we note that n is the eigenvalue ofn in the Fock state basis. The propagator of this Hamiltonian in the interaction picture satisfies the Schrödinger equation whereĤ (n) To make progress, we make an ansatz for U (n) where f n,k and Λ n,k are obtained from solution of the Schrödinger equation, Eq. 12. By inserting the ansatz equation into Eq. 12, and equating the LHS and RHS, we find Since we find the exact expression of the propagator in the interaction picture, the propagator in the Schrödinger picture is given byÛ (n) =Û BÛ , which can then be expressed aŝ whereD is the displacement operator. At last, we invoke a beautiful mathematical equality where we assume the initial bath density matrix to be a thermal state and β is the inverse temperature [27]. With the above equality, we then obtain When we take the time derivation of the reduced system density matrix, we have Here,Ĥ 1. Explicit expressions for the observables of the bath and the functions g(t) and γ(t) As we have an exact solution for the density matrix of the total system, we can explicitly calculate observables of the bath. For example, the mean photon number N k = a † k a k for the kth mode reads Similarly, the expectation values of the quadraturesX k = 1 The physical intuition behind this solution is that the bath is out of equilibrium due to the coupling between the system and the bath and its time evolution is affected by the number of photons in the system. This is a total opposite to a Markovian dynamics, where the bath is not influenced by the system. Now we provide the explicit expressions for the bath-induced non-linearity and the dissipation rates γ(t) The previous analytical expressions were obtained for g k = he −zk/2 , where z > 0 and ω k = kΩ.

C. Circuit QED implementation
Circuit quantum electrodynamics (QED) has emerged as a promising platform to engineer strongly-correlated states of quantum matter, where "particles" arise from excitations of low-temperature electrical circuits [43]. In this section we provide a circuit design that implements the system-bath Hamiltonian discussed in the main text, i.e., The modesâ 0 andâ † 0 will play the role ofb andb † in the main text, respectively. In addition, the frequency of the zeroth mode is the frequency of the resonator ω 0 = ω in the main text. To avoid nonlinear coupling between too many pairs of resonators, we apply the decomposition in eigenmodesâ k = N m=1 V k,mâm , where V k,m are the elements of a M times M square matrix. In terms of the new bosonic modes we obtain where the couplingg 1 = N k=1 V m,k is related to the eigenmode decomposition discussed above. To substantiate the structure of the mode decomposition, we assume periodic boundary conditionsâ 1 =â † N +1 for the bosonic modes. This can be achieved by introducing a capacitive coupling between the nodes 1 and N in Fig.4. In this case, the coefficients appearing in the mode decomposition read V m,k = 1 √ N e ikm . The latter equation is much less demanding as nonlinear coupling between only one pair of resonators is required. The circuit diagram for implementing the latter is shown in Fig.4. As will be shown below, each LC circuit forms a resonator with the frequencyω m ≈ 1/ √ L m C m . These LC circuits may as well be replaced with transmission lines which can be fabricated with higher precision [44]. However, the calculation for the latter is more troublesome, so we restrict ourself to the LC circuits instead for simplicity without compensating the physics. Nonlinear coupling comes from the use of the Josephson junctions with an external magnetic driving field.
Following the standard circuit quantization procedure [45], we first write down the circuit's Lagrangian as where C m , C m,m+1 are capacitance, L m are inductance, Φ 0 =h/2e is the flux quantum, E J is the Josephson energy, Φ b (t) = π + φ b (t) is a flux bias and φ m = −´V m dt is a flux variable, with V m being a voltage at the corresponding position. We choose the flux bias field φ b (t) to be an oscillating field with the frequency ω p , which can be implemented using an external AC magnetic field [46]. The drive frequency will be chosen to eliminate undesired terms in the cosine expansion using the rotating wave approximation (RWA). The Hamiltonian can be obtained using the Legendre transformation [47], where q m = C m ∂L/∂φ m is a conjugate momentum of φ m ,C m = C m,m−1 + C m,m+1 + C m is an effective capacitance.
Here we have assumed that C m /C m 1. We then quantized φ m and q m by defining ladder operatorsâ m ,â † m according toφ m = (L m /4C m ) 1/4 (â m +â † m ) andq m = i(C m /4L m ) 1/4 (−â m +â † m ). It follows that whereω m = 1/ L mCm and J m = − √ω mωm+1 C m,m+1 /2C m . Here we have assumed that J m ω m and hence the rotating term,â † mâ † m+1 + H.c., can be ignored with RWA. When working in the low excitation regime justified by a weak driving φ b (t) ω m , the expansion of the cosine function term can be kept up to the fourth order (η = 0, 1, 2). The quadratic term (η = 1) will simply renormalize the frequency of the resonators. This lefts us with only the fourth-order terms, We then choose the coherent drive with the frequency ω p =ω 1 , i.e. φ b (t) = Ω(be −iω1t + H.c.) whereb is promoted to a c-number. RWA can be applied for a weak driving Ω ω m . The only non-rotating term in Eq.33 that survives after the RWA is thenâ † 0â 0 (â † 1 +â 1 ) as desired. As we discussed before, by engineering the energiesω m and couplings J m , one can obtain a linear dispersion for the frequencies ω k and the desired couplings g k .