Asymptotic Floquet states of open quantum systems: The role of interaction

We investigate the asymptotic state of a periodically driven many-body quantum system which is weakly coupled to an environment. The combined action of the modulations and the environment steers the system towards a state being characterized by a time-periodic density operator. To resolve this asymptotic non-equilibrium state at stroboscopic instants of time, we introduce the dissipative Floquet map, evaluate the stroboscopic density operator as its eigen-element and elucidate how particle interactions affect properties of the density operator. We illustrate the idea with a periodically modulated Bose-Hubbard dimer and discuss the relations between the interaction-induced bifurcations in a mean-field dynamics and changes in the characteristics of the genuine quantum many-body state. We argue that Floquet maps provide insight into the system relaxation towards its asymptotic state and may help to understand whether it is possible (or not) to construct a stroboscopic time-independent generator mimicking the action of the original time-dependent one.


Introduction
Many-body effects in combination with a coupling to an environment and in the presence of time variations of parameters give rise to a variety of phenomena useful for quantum technologies. Interactions sculpt the spectrum of different collective states and moderate transitions between them [1]. Effects of the system-environment coupling, however weak they are, play a decisive role in out-shaping the system's asymptotic state. Indeed, such effects may not necessarily present a nuisance but can be of beneficial use. Particularly, they can be exploited to steer the system towards desired states, including pure and high-entangled ones [2,3]. This recent idea of engineering by dissipation [4,5,6,7,8] has promoted a dissipative time evolution to the same level of importance as that of a unitary evolution. Finally, time periodic modulations can strongly modify the state of a quantum system. In the coherent limit, modulations implicate an explicit timeperiodicity of the Hamiltonian, H(t + T ) = H(t + 2π/ω) = H(t), where T is the driving period and ω is the frequency of modulations. The system dynamics is determined by the Floquet states [9,10,11], the eigenstates of the unitary Floquet propagator U T = T exp − i T 0 H(τ )dτ , where T is the time-ordering operator. The particular structure of the Floquet propagator, and thus the properties of the Floquet states, depend on modulation parameters. This allows to grasp effects [11,12,13,14,15,16,17,18,19,20] which are out of reach of experimentally available time-independent Hamiltonians in atom optics, optomechanics, and solid state physics .
The combined effect of all three factors, (i) many-body interactions, (ii) coupling to an environment, and (iii) periodic external driving, provides the motivation for the investigation we undertake with this study. We start by introducing the notion of Floquet maps which constitute an extension to the dissipative case of the Floquet propagator [9,10,11], and demonstrate how those can be used to obtain non-equilibrium asymptotic states. We next study a many-body model, a driven Bose-Hubbard dimer, and use both, a full quantum mechanical treatment and a mean-field description, to gain insight into the properties of the time-periodic asymptotic states.

Dissipative Floquet maps
We consider the dynamics of a general M -dimensional system modeled with a quantum master equation whose generator L is of Lindblad form [21,22,23,24] The first term on the r.h.s. describes the unitary evolution of the system's density operator , governed by the time-periodic Hamiltonian H(t). The dissipator is built from the set of operators {V k }, which, together with the normalized identity, V 0 = 1/M , spans the Hilbert-Schmidt space B M of the operators acting in Mdimensional Hilbert space [25]. Note that all parameters of the system are scaled with respect to the Planck constant . Strictly speaking, modulations affect both the unitary and the dissipative part of the generator L t , and in general both the Hamiltonian H(t) and the dissipative matrix Γ(t) = γ kl (t), have to become time-dependent [12,23,26]. The complete theoretical foundation of a time-dependent master equation of the type (1) still remains an open problem [27,28]. However, if the matrix Γ(t) is positive semidefinite at any instant of time, the propagator P s,t = T exp( t s L τ dτ ) is completely positive and trace-preserving. In this case a master equation in the form of (1) is meaningful [23]. This is a set-up which is now frequently used to model quantum systems far from equilibrium; see, e.g., [8,29,30,31].
When the generator L t is time dependent, the propagator P s,t depends on both parameters, the starting time s and the final time t. The closure of the set of propagators for different times is lost and they no longer form a semi-group. It is stated by Lendi [28] that "the best chance to find a solution to a master equation [with a time-dependent generator] is only offered by a possible existence of transformations which eliminate the time dependence". Consistently, most studies until now have focused on removing the time-dependence when dealing with time-periodic generators, either by (i) finding a proper gauge, which makes the original time-periodic Hamiltonian time-independent [23,32] and then assuming that the dissipator remains time-independent in the new frame (this is often a good approximation in quantum optics, where frequencies of modulations are much higher than the decay rates), or (ii) by changing to the Floquet basis of the Hamiltonian H(t) and then performing an additional secular approximation [8], or (iii) by trying to construct a Magnus expansion-like approximations [33]. All these strategies result in deriving an effective time-independent generator L eff of the Lindblad form. Once the time-dependence is removed, one has to calculate the kernel of L eff to find the asymptotic state ∞ of the system, i.e., L eff ∞ = 0. Under fairly general conditions [23,34], the time-homogeneous propagator P t = exp(L eff t) relaxes towards a unique attractor ∞ of the dissipative quantum evolution. However, the above discussed approximations cannot always be justified. Below we propose an approach which does not demand the reduction to a time-independent form and avoids corresponding approximation schemes. Note also that it is not necessary to switch to the Floquet basis of the system Hamiltonian, a step often performed; see Refs. [8,12,31].
Because the master equation (1) is linear, in the case of a time-periodic generator L t we can resort to Floquet theory [11,35]. We concentrate on the propagator P F ≡ P 0,T = T exp( T 0 L τ dτ ) which we refer to as the Floquet map. The Floquet map has at least one eigenvalue 1 and all other eigenvalues lie inside the unit circle. Assuming that the Floquet map is irreducible [2], the attractor 1 of the map is its fixed point, i.e., the eigen-operator corresponding to the eigenvalue 1, P F 1 = 1 . The uniqueness of the attractor solution is related to the symmetries of the generator L t [36]. After a sufficiently large time any initial density operator 0 (t 0 ) will converge to the time-periodic asymptotic state a (t), i.e., 0 (t 0 ) → a (mT +t 0 ) for an integer m 1.
The operator 1 = a (0) = a (T ) is the asymptotic density operator of the system at the stroboscopic instants of time. Because a (t) is periodic in time, a (mT + s) = a (s), the asymptotic density matrix for any instance of time t, can be calculated via propagating a (0) up to time s.

Model study: Bose-Hubbard dimer
Evidently, the proposed approach, being general and mathematically formal, does not differentiate between single-particle and of many-particle systems. To exemplify the ample physics expected to emerge from the interplay of many-body interaction, dissipation and periodic driving, we consider a system composed of N interacting bosonic atoms hopping over a dimer which is subjected to periodic driving. The considered system Hamiltonian reads where J denotes the tunneling amplitude, U is the interaction strength, and ε(t) presents the modulation of the local potential. In particular we choose ε(t) = ε(t + T ) = µ 0 +µ 1 sin(ωt), where µ 0 presents a static and µ 1 models a dynamic energy offset between the two sites. Here, b j and b † j are the annihilation and creation operators of an atom at site j, and n j = b † j b j . This Hamiltonian has been previously studied theoretically [37,38,39,40] and additionally has been implemented in several recent experimental studies [41,42]. However, to the best of our knowledge, the joint action of all three ingredients -interaction, dissipation and temporal driving -has not been addressed before.
For the dissipator D we use the single jump operator [3,43], which tends to "synchronize" the dynamics on the dimer sites, by constantly recycling anti-symmetric out-phase modes into the symmetric in-phase ones. The coupling constant γ is taken to be time-independent. This particular set-up serves as an illustration only. Note that the Floquet map approach applies equally well to other cases, e.g., when both parts of the generator, unitary and dissipative ones, are timeperiodic or when there are several jump operators acting on the system. Because the jump operator (4) is non-Hermitian, the propagators P s,t are not unital and the attractor solution is not the maximally mixed state, i.e. a = 1/M . To gain additional insight into the physics of the model, we derive a set of meanfield equations and compare its attractor solutions with those of the quantum Floquet map P F . For the dimer problem, it is convenient to recast the master equation (1) in terms of the spin operators and then study their evolution in the Heisenberg picture [24]. For a large number of atoms N 1, the commutator [S x , S y ] = iS z /N becomes negligible small and similarly for other cyclic permutations. Replacing operators with their expectation values, S k = tr[ S k ], and denoting S k by S k , we end up with where we have neglected terms proportional to γ of lower order in N . The replacement of operators by their expectation values is justified provided that AB t ≈ A t B t . This is not guaranteed a priori, and, for a dissipative system, the commutator behaves differently compared to the unitary setup. A necessary comparison with the results of the exact quantum analysis then does justify the validity of the approximation. As the quantity S 2 = S 2 x + S 2 y + S 2 z is a constant of motion, we can reduce the mean-field evolutions to the surface of a Bloch sphere; i.e., (S x , S y , S z ) = 1 2 [cos(ϕ) sin(ϑ), sin(ϕ) sin(ϑ), cos(ϑ)] , yielding the equations of motioṅ ϑ = 2J sin(ϕ) + 4γN cos(ϕ) cos(ϑ), We next analyze the quantum dynamics by using both the Floquet map computed via equations (1)(2)(3)(4) and contrast the results with the mean-field equations (8). To find the classical attractor solution of the mean-field system, we evolve (8) from randomly chosen initial conditions, and, after a transient time 10 4 T , record the value of S z at the next 250 stroboscopic instants of time, cf. figure 1(a). As the interaction U N varies, we detect regions containing limit cycles of different periods, chaotic attractors, and transitions between them. Figures 1(c,d) show that the different dynamical regimes of the mean-field description are characterized by significantly different properties of the system in the quantum limit. To quantify this observation, we calculate the timeaveraged puritȳ and also the time-averaged negativity [44], which characterizes the degree of entanglement in a two-mode system of N indistinguishable bosons. It is noteworthy that, as the number of bosons increases, changes of the time-averaged purity and negativity become more pronounced near the vicinities of the mean-field bifurcation transitions. Spectral properties of a Floquet map may also provide insight into the relaxation towards the corresponding quantum attractor. A typical spectrum of a map is shown in the inset of figure 1(b). It has a shape inherent to the spectra of completely positive trace-preserving maps [45]. Namely, it has the spectral radius 1, includes the single eigenvalue λ 1 = 1, and is invariant under complex conjugation. The spectral gap, ∆ = 1 − |λ 2 |, where λ 2 is the second largest eigenvalue (by absolute value), can be used to estimate the inverse relaxation time from a randomly chosen initial state [46,47,48]. We find that the spectral gap also exhibits a strong dependence on the interaction strength; see figure 1 Different mean-field regimes can be visualized by plotting stroboscopic Poincaré sections on the plane {ϑ, ϕ}. In Figure 2, classical Poincaré sections are compared with the Poincaré-Husimi distributions p(ϑ, ϕ) of the quantum asymptotic state obtained by projecting the density operator a (0) on the set of the generalized SU(2) coherent states [49]. For U N/J = 0.2(0.8), the mean-field model predicts two(six) points on the Poincaré section, corresponding to period-two (period-two plus period-four) attractor(s); see symbols in Fig. 2(a,b). The Poincaré-Husimi distributions, Fig. 2(a-b), reveal a concentration of p(ϑ, ϕ) near these points (we attribute the minor mismatch to still present finite-size effects). For U N/J = 1, the mean-field system (8) exhibits a chaotic attractor, Fig. 2(d), and the Poincaré-Husimi distribution, Fig. 2(c), fits the structure of this classical attractor for N = 250. Fig. 3 shows three-dimensional plots of the quantum attractors super-imposed on the classical Poincaré sections both for the case in which the mean-field equations predict two points (from a period-two limit cycle) or a chaotic attractor [corresponding respectively to Fig. 2(a) and Fig. 2(c,d)]. It is noteworthy that the inverse particle number 1/N can be thought of as an effective Planck constant, thus allowing for the comparison with the results obtained for single-particle models [50,51,52,53]; note in addition those cited in the mini-review [52].
The Floquet map is a completely-positive and trace-preserving map which belongs, following the nomenclature introduced in Ref. [54], to the class of time-dependent Markovian channels. It is an interesting question whether it is possible to find an effective time-independent generator L eff , which is of Lindblad form, Eqs. (1-2), and can mimic the action of the original generator at stroboscopic instants of time, such that P F = exp(L eff T ). The formal recipe to get an answer -"yes" or "no" -for a particular P F is given in Ref. [55]. However, it is hardly realizable in practice when M > 2 because it involves repeated solution of an O(M 4 ) optimization problem within the mixed-integer semidefinite programming framework [56].
There is a necessary condition a map P F should fulfill in order to have a generator L eff of Lindblad form, which is much easier to check. Namely, the spectrum {λ j } of the map should not contain real negative eigenvalues (strictly speaking, of odd algebraic multiplicity); otherwise it is impossible fulfill the condition of the symmetry under complex conjugation for the spectrum of the generator [45]. Figure 4 shows the number of negative real eigenvalues [57] as a function of the driving frequency ω. The dependence reveals that the condition is not fulfilled in the most interesting case of nonadiabatic driving, when the asymptotic state of the dimer is sculpted by the modulations. Apparently, an effective stroboscopic time-independent Lindblad generator does not exist in this regime.

Numerical methods
The key technical step of our approach is the construction of the dissipative Floquet map. This task can be performed by unfolding the propagator P F into a matrix form by using any Hilbert-Schmidt basis in B M , {F l } l=0,M 2 −1 , which yields P F,kl = tr[F † k , P F F l ]. Notably, this basis may differ from the basis {V l } used in (2). The elements P F F l need to be calculated numerically by propagating the corresponding basis element over one temporal period T of the periodically varying modulation. The Floquet map P F has the eigenvalue 1 and the corresponding eigenoperator gives after normalization the asymptotic stroboscopic state a (0).
Towards this goal, one may use the Kronecker product to write the Lindblad master equation (1) in the vectorized form [58] with the M 2 × M 2 Liouvillian given bỹ Here, A * denotes complex conjugation of every element of the matrix A, and the linear operator vec(A) denotes a vectorization which transforms an m × m matrix A = (a kl ) into a column vector vec(A) = (a 11 , . . . , a m1 , a 12 , . . . , a m2 , a 13 , . . . , a mm ) T .
Then, the Floquet propagator P F ≡ P 0,T can be evaluated readily by propagating every unit vector e j = (δ 0j , δ 1j , . . . , δ M j ) T over a full period T [59]. Often, the Hamiltonian H and the jump operators V k are given as sparse matrices. For example, in the Fock-basis the Hamiltonian (3) and the jump operator (4) are tridiagonal implying that the number of non-zero matrix elements is of order O(N ) where N is the number of bosons. The Kronecker product of sparse matrices also yields a sparse matrix, and it can be shown that the LiouvillianL possesses O(N 2 ) non-zero matrix elements. Moreover, the action ofL on an arbitrary vector v, i.e., the multiplicatioñ Lv, can be implemented with a O(N 2 ) efforts.
For an efficient numerical calculation of the Floquet propagator it is important to use a method for solving the linear system of differential equations (11) which does exploit the sparseness ofL. One option is to use the standard Runge-Kutta scheme. As the multiplication ofL with a vector needs O(N 2 ) operations, calculating all elements of the Floquet propagator P F costs approximately O(N 4 ). The step size of the Runge-Kutta scheme depends on the parameters µ 0 , µ 1 , J, U, ω and γ, and is also weakly dependent on the number of particles N . Therefore, this algorithm scales slightly worse than O(N 4 ), but is still much faster than a simplistic approach with a O(N 6 ) scaling.
While it may be possible to calculate all elements of the Floquet propagator P F , this is not necessary if one is solely interested in obtaining the asymptotic state and a few eigen-elements of P F . The asymptotic state can be calculated using the so termed power method which corresponds to propagating a random density operator until it converges to the asymptotic state. If one is also interested in a few additional eigenvalues and eigen-elements of P F , one may use the implicitly restarted Arnoldi method [60]. This method does not use the matrix elements of P F directly, however; it merely requires the action of the propagator P F on an arbitrary state v, P F v, to calculate eigenvalues and eigenvectors. As the operation P F v can be implemented in O(N 2 ) operations, using the Runge-Kutta method and exploiting the sparseness ofL; this latter scheme scales approximately as O(N 2 ) only.
Finally, in order to calculate the asymptotic state we employed ARPACK's implementation of the implicitly restarted Arnoldi method [61], using a wrapper from SciPy [62].

Conclusions
We demonstrated that the concept of dissipative Floquet maps provides an operational way to identify 'quantum attractors', i.e., asymptotic time-periodic states of modulated open quantum systems. To illustrate this idea, we have applied the concept to a dissipative and periodically driven many-body model. We study the model both exactly and, in the limit of a large particle number, within a mean-field picture. The latter predicts transitions from regular to chaotic attractor as the interaction strength is varied. The analysis shows a strong dependence on the interaction strength, especially in proximity of the dynamical transitions predicted by mean-field theory.
We conclude by pointing out possible research directions which may benefit from the use of Floquet maps. It has been proposed to use time-periodic driving to create, for the situation with coherent Hamiltonian systems, effective topologically protected states [15,63]. The important problems of the stability of these states against dissipation or their creation with a synthetic dissiaption [64] could be investigated by making use of our concept. Another interesting question is whether the idea of 'engineering by dissipation' [2,6,7,25] can be extended somehow to periodically modulated systems. Finally, recent progress in the field of many-body localization (MBL) inaugurates yet another potential application; the effect of temporal driving on the localization has been addressed in [65,66,67] and, very recently, the dynamics of open MBL systems was considered in [68,69,70]. We expect that the idea to combine the two latter ingredients may soon invigorate the MBL community in pursuing future research in this spirit (see a very recent Ref. [71]).