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

Many-body effects in combination with a coupling to an environment give rise to a variety of phenomena which are of beneficial use 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 as well of practical 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–8] has promoted a dissipative time evolution of the system dynamics to the same level of importance as that obtained with a unitary evolution.

Time periodic modulations can also strongly modify the state of a quantum system. In the coherent limit, time-periodic modulations implicate an explicit time-periodicity of the Hamiltonian, i.e., $H(t+T)=H(t+2\pi /\omega )=H(t)$, where T is the driving period and ω is the frequency of modulations. The system dynamics is governed by the Floquet states [9–11]; i.e., the eigenstates of the unitary Floquet propagator ${U}_{T}={ \mathcal T }\exp \left[-\tfrac{{\rm{i}}}{{\hslash }}{\int }_{0}^{T}H(\tau ){\rm{d}}\tau \right]$, where ${ \mathcal 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–20] which are out of reach of experimentally available time-independent Hamiltonians in atom optics, optomechanics, and solid state physics.

Our key objective here is to investigate the combined effect of all three factors, namely (i) many-body interactions, (ii) coupling to an environment, and (iii) periodically varying external driving. We start by introducing the notion of Floquet maps which constitute an extension to the dissipative case [11] of the unitary Floquet propagator [9, 10], 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.

We consider the dynamics of a general M-dimensional system modeled with a quantum master equation whose generator ${ \mathcal L }$ is of Lindblad form [21–24]

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\varrho ={{ \mathcal L }}_{t}(\varrho )=-{\rm{i}}[H(t),\varrho ]+{{ \mathcal D }}_{t}(\varrho ).\end{eqnarray} \tag{ 1 }$$

The first term on the rhs describes the unitary evolution of the system's density operator ϱ, governed by the time-periodic Hamiltonian H(t). The dissipator

$$\begin{eqnarray}&&{{ \mathcal D }}_{t}(\varrho )=\displaystyle \sum _{l,k=1}^{{M}^{2}-1}{\gamma }_{{kl}}(t)\left[{V}_{k}\varrho {V}_{l}^{\dagger }-\displaystyle \frac{1}{2}\{{V}_{l}^{\dagger }{V}_{k},\varrho \}\right],\end{eqnarray} \tag{ 2 }$$

is built from the set of operators $\{{V}_{k}\}$, which, together with the normalized identity, ${V}_{0}={\mathbb{1}}/M$, span the Hilbert–Schmidt space ${{ \mathcal B }}_{M}$ of the operators acting in the M-dimensional Hilbert space [25]. Note that all parameters of the system are scaled with respect to the Planck constant ${\hslash }$.

Strictly speaking, modulations affect both the unitary and the dissipative part of the generator ${{ \mathcal L }}_{t}$, and in general both the Hamiltonian H(t) and the dissipative matrix ${\rm{\Gamma }}(t)={\gamma }_{{kl}}(t)$, become time-dependent [11, 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 ${\rm{\Gamma }}(t)$ is positive semi-definite at any instant of time, the propagator ${{ \mathcal P }}_{s,t}={ \mathcal T }\exp ({\int }_{s}^{t}{{ \mathcal L }}_{\tau }{\rm{d}}\tau )$ is completely positive and trace-preserving. In this case a master equation in the form of (1) is meaningful [23]. This in turn provides a set-up which is frequently employed to model quantum systems operating far from equilibrium; see, e.g., [8, 29–31].

When the generator ${{ \mathcal L }}_{t}$ is time dependent, the propagator ${{ \mathcal P }}_{s,t}$ depends on both 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 constructing Magnus expansion-like approximations [33]. All these strategies result in deriving an effective time-independent generator ${{ \mathcal L }}_{\mathrm{eff}}$ of the Lindblad form. Once the time-dependence is removed, one has to calculate the kernel of ${{ \mathcal L }}_{\mathrm{eff}}$ to find the asymptotic state ${\varrho }_{\infty }$ of the system, i.e., ${{ \mathcal L }}_{\mathrm{eff}}{\varrho }_{\infty }={\mathbb{0}}$. Under fairly general conditions [23, 34], the time-homogeneous propagator ${{ \mathcal P }}_{t}=\exp ({{ \mathcal L }}_{\mathrm{eff}}t)$ relaxes towards a unique attractor ${\varrho }_{\infty }$ of the dissipative quantum evolution. However, the above discussed approximations cannot always be justified away from the case of high frequency driving. Below we propose an approach which does not demand the reduction to such a time-independent form and thus avoids those corresponding approximation schemes. Note also that it is not necessary to switch to the Floquet basis of the driven system Hamiltonian.

Because the master equation (1) is manifestly linear, we can in the case of a time-periodic generator ${{ \mathcal L }}_{t}$ readily resort to the Floquet theorem [11, 35, 36]. We concentrate next on the one period propagator ${{ \mathcal P }}_{{\rm{F}}}\equiv {{ \mathcal P }}_{0,T}={ \mathcal T }\exp ({\int }_{0}^{T}{{ \mathcal L }}_{\tau }{\rm{d}}\tau )$ which we refer to as the Floquet map. The Floquet map possesses at least one (possibly degenerate) eigenvalue 1 and all other eigenvalues lie inside the unit circle. Assuming that the Floquet map is irreducible [2], the attractor ϱ₁ of the map is given by its fixed point, i.e., the eigen-operator corresponding to the eigenvalue 1, ${{ \mathcal P }}_{{\rm{F}}}{\varrho }_{1}={\varrho }_{1}$. More generally, the number of different attractor solutions is directly related to the symmetries of the generator ${{ \mathcal L }}_{t}$ [37]. Particularly, in absence of such additional symmetries the resulting asymptotic attractor assumes the unique fixed point solution. Therefore, after a sufficiently large time span any initial density operator ${\varrho }_{0}({t}_{0})$ will converge to the time-periodic asymptotic state ${\varrho }_{{\rm{a}}}(t)$, i.e., ${\varrho }_{0}({t}_{0})\to {\varrho }_{{\rm{a}}}({mT}+{t}_{0})$ for an integer $m\gg 1$. The operator ${\varrho }_{1}={\varrho }_{{\rm{a}}}(0)={\varrho }_{{\rm{a}}}(T)$ is the asymptotic density operator of the system at the stroboscopic instants of time. Because ${\varrho }_{{\rm{a}}}(t)$ is periodic in time, ${\varrho }_{{\rm{a}}}({mT}+s)={\varrho }_{{\rm{a}}}(s)$, the asymptotic density matrix for any instance of time t can be calculated via propagating ${\varrho }_{{\rm{a}}}(0)$ up to time s.

To exemplify the ample physics expected to emerge from the interplay of many-body interactions, dissipation and periodic driving, we consider a system composed of N interacting bosonic atoms hopping over a dimer which is subjected to periodic driving. We consider the system Hamiltonian

$$\begin{eqnarray}&&H(t)=-J({b}_{1}^{\dagger }{b}_{2}+{b}_{2}^{\dagger }{b}_{1})+\displaystyle \frac{U}{2}\displaystyle \sum _{j=1,2}{n}_{j}({n}_{j}-1)+\varepsilon (t)({n}_{2}-{n}_{1}),\end{eqnarray} \tag{ 3 }$$

where J denotes the tunneling amplitude, U is the interaction strength, and $\varepsilon (t)$ represents the modulation of the local potential. In particular we choose $\varepsilon (t)=\varepsilon (t+T)={\mu }_{0}+{\mu }_{1}\sin (\omega t)$, where ${\mu }_{0}$ presents a static and ${\mu }_{1}$ models a dynamic energy offset between the two sites. Here, b_j and ${b}_{j}^{\dagger }$ are the annihilation and creation operators of an atom at site j, and ${n}_{j}={b}_{j}^{\dagger }{b}_{j}$. This Hamiltonian has been previously studied theoretically in [38–41] and has been implemented in several recent experimental studies [42, 43]. However, to the best of our knowledge, the joint action of all three ingredients—interaction, dissipation and temporal driving—has not been addressed before.

With the coupling constant γ taken to be time-independent, we use ${{ \mathcal D }}_{t}={ \mathcal D }$ with the single jump operator [3, 44]

$$\begin{eqnarray}&&V=({b}_{1}^{\dagger }+{b}_{2}^{\dagger })({b}_{1}-{b}_{2}).\end{eqnarray} \tag{ 4 }$$

This dissipator tends to 'synchronize' the dynamics on the dimer sites by constantly recycling anti-symmetric out-phase modes into the symmetric in-phase ones. Note that our particular setup serves as an illustration only. The Floquet map approach applies equally well to other cases, e.g., when both parts of the generator, i.e., the unitary and dissipative parts both are time-periodic or when there are several jump operators acting on the system. Because the jump operator (4) is non-Hermitian, the propagators ${{ \mathcal P }}_{s,t}$ are not unital and the attractor solution is not the maximally mixed state, ${\varrho }_{{\rm{a}}}\ne {\mathbb{1}}/M$.

To gain additional insight into the physics of the model, we derive a set of mean-field equations and compare its attractor solutions with those of the quantum Floquet map ${{ \mathcal P }}_{{\rm{F}}}$. For the dimer problem, it is convenient to recast the master equation (1) in terms of the spin operators

$$\begin{eqnarray}&&{{ \mathcal S }}_{x}=\displaystyle \frac{1}{2N}({b}_{1}^{\dagger }{b}_{2}+{b}_{2}^{\dagger }{b}_{1}),\quad {{ \mathcal S }}_{y}=-\displaystyle \frac{{\rm{i}}}{2N}({b}_{1}^{\dagger }{b}_{2}-{b}_{2}^{\dagger }{b}_{1}),\quad {{ \mathcal S }}_{z}=\displaystyle \frac{1}{2N}({n}_{1}-{n}_{2}),\end{eqnarray} \tag{ 5 }$$

and then study their evolution in the Heisenberg picture [24]. For a large number of atoms $N\gg 1$, the commutator $[{{ \mathcal S }}_{x},{{ \mathcal S }}_{y}]={\rm{i}}{{ \mathcal S }}_{z}/N$ becomes negligibly small and similarly for other cyclic permutations. Replacing operators with their expectation values, $\langle {{ \mathcal S }}_{k}\rangle =\mathrm{tr}[\varrho {{ \mathcal S }}_{k}]$, and denoting $\langle {{ \mathcal S }}_{k}\rangle$ by S_k, we end up with

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{S}_{x}}{{\rm{d}}t}=2\varepsilon (t){S}_{y}-2{{UNS}}_{z}{S}_{y}+8\gamma N({S}_{y}^{2}+{S}_{z}^{2}),\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{S}_{y}}{{\rm{d}}t}=-2\varepsilon (t){S}_{x}+2{{UNS}}_{x}{S}_{z}+2{{JS}}_{z}-8\gamma {{NS}}_{x}{S}_{y},\end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{S}_{z}}{{\rm{d}}t}=-2{{JS}}_{y}-8\gamma {{NS}}_{x}{S}_{z},\end{eqnarray} \tag{ 6c }$$

where we have neglected terms proportional to γ of lower order in N. The replacement of operators by their expectation values is justified provided that $\langle {AB}{\rangle }_{t}\approx \langle A{\rangle }_{t}\langle B{\rangle }_{t}$. This is not guaranteed a priori, and, for a dissipative system, the commutator behaves differently compared to the unitary setup. A necessary favorable comparison with the results of the exact quantum analysis then justifies the validity of this mean-field approximation.

The structure of the mean-field equations in (6a)–(6c) implies that $\tfrac{{\rm{d}}}{{\rm{d}}t}{S}^{2}=0$. Therefore the quantity ${S}^{2}={S}_{x}^{2}+{S}_{y}^{2}+{S}_{z}^{2}$ is a constant of motion. This is consistent with the preservation of the total number of bosons N; see the definitions given by (5). We therefore can reduce the mean-field evolutions to the surface of a Bloch sphere; i.e.,

$$\begin{eqnarray}&&({S}_{x},{S}_{y},{S}_{z})=\displaystyle \frac{1}{2}[\cos (\varphi )\sin (\vartheta ),\sin (\varphi )\sin (\vartheta ),\cos (\vartheta )],\end{eqnarray} \tag{ 7 }$$

yielding the equations of motion

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\vartheta & = & 2J\sin (\varphi )+4\gamma N\cos (\varphi )\cos (\vartheta ),\\ \displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\varphi & = & 2J\displaystyle \frac{\cos (\vartheta )}{\sin (\vartheta )}\cos (\varphi )-2\varepsilon (t)+{UN}\cos (\vartheta )-4\gamma N\displaystyle \frac{\sin (\varphi )}{\sin (\vartheta )}.\end{eqnarray} \tag{ 8 }$$

We next analyze the quantum dynamics by using both the Floquet map computed via (1)–(4) and contrast the results with the mean-field equations (8).

To construct the Floquet map, we use the standard scheme of a vectorization of the density matrix, which allows to transform (1) into a system of linear differential equations with time-periodic coefficients. The Floquet map is obtained by propagating the δ-Kronecker basis [45] over the full period T. Finally, the asymptotic density matrix is given as the eigen-element of the map corresponding to the unique eigenvalue one.

To extract 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. The so obtained bifurcation diagram is presented in figure 1(a). As the interaction UN varies, we detect regions containing limit cycles of different periods, chaotic attractors, and transitions between them [46]. We anticipate that different dynamical regimes of the mean-field description are characterized by significantly different properties of the system in the quantum limit for $N\gg 1$. To check this hypothesis, we calculate the time-averaged purity

$$\begin{eqnarray}&&\bar{P}=\displaystyle \frac{1}{T}{\int }_{0}^{T}\ \mathrm{tr}\ [{\varrho }_{a}{(t)}^{2}]{\rm{d}}t\end{eqnarray} \tag{ 9 }$$

and also the time-averaged negativity,

$$\begin{eqnarray}&&\bar{{ \mathcal N }}=\displaystyle \frac{1}{T}{\int }_{0}^{T}{ \mathcal N }[{\varrho }_{a}(t)]{\rm{d}}t.\end{eqnarray} \tag{ 10 }$$

Here, ${ \mathcal N }$ represents the negativity [47],

$$\begin{eqnarray}&&{ \mathcal N }[\varrho ]=\displaystyle \frac{1}{2}\displaystyle \sum _{k\ne l}| {\varrho }_{k,l}| ,\end{eqnarray} \tag{ 11 }$$

which characterizes the degree of entanglement in a two-mode system of N indistinguishable bosons. Figures 1(c) and (d) show the dependence of the two quantities as functions of the interaction strength. It is interesting that, as the number of bosons increases, changes of the time-averaged purity and negativity become more pronounced in the vicinity of bifurcations of the mean-field equations.

Zoom In Zoom Out Reset image size

The inverse tangent bifurcation [46] near ${UN}/J=2$ (the transition from chaos to a period-one limit cycle) is striking: both the negativity and the purity of the asymptotic state move to higher values at this point. This relates to the concept of 'dissipative engineering' used to shape a stationary many-body system into a pure highly-entangled equilibrium state with the help of specially designed dissipative operators [4]. Here, we observe a trend towards a pure highly-entangled non-equilibrium state upon increasing the particle number N. An intriguing question arises as to which values both characteristics saturate in the thermodynamic limit $N\to \infty$. Will, for example, the purity value approach unity? If 'yes' then we would have a first example of dissipative engineering of a time-periodic quantum state. Unfortunately, it was not possible for us to go beyond $N\approx 300$ by using the numerical spectral methods.

Furthermore, 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 ${\lambda }_{1}=1$, and is invariant under complex conjugation. The spectral gap

$$\begin{eqnarray}&&{\rm{\Delta }}=1-| {\lambda }_{2}| ,\end{eqnarray} \tag{ 12 }$$

where ${\lambda }_{2}$ is the second largest eigenvalue by absolute value, can be used to estimate the inverse relaxation time from a randomly chosen initial state [48–50]. The spectral gap also exhibits a strong dependence on the interaction strength; see figure 1(b).

Different mean-field regimes can be visualized by plotting stroboscopic Poincaré sections on the plane $\{\vartheta ,\varphi \}$. In figure 2, classical Poincaré sections are compared with the Poincaré–Husimi distributions $p(\vartheta ,\varphi )$ of the quantum asymptotic state obtained by projecting the density operator ${\varrho }_{{\rm{a}}}(0)$ on the set of the generalized SU(2) coherent states [51]. For ${UN}/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 figures 2(a, b). The Poincaré–Husimi distributions, figures 2(a) and (b), reveal a concentration of $p(\vartheta ,\varphi )$ near these points. We attribute the minor mismatch to finite-size effects. For ${UN}/J=1$, the mean-field system (8) exhibits a chaotic attractor, figure 2(d), and the Poincaré–Husimi distribution, figure 2(c), fits the structure of this classical attractor for N = 250. Figure 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 figures 2(a) and (c), (d).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 [52–55]; note in addition those cited in the mini-review [54].

The Floquet map ${{ \mathcal P }}_{{\rm{F}}}$ is a completely-positive and trace-preserving map which belongs, following the nomenclature introduced in [56], to the class of time-dependent Markovian channels. It is an interesting question whether it is possible to find an effective time-independent generator ${{ \mathcal L }}_{\mathrm{eff}}$ of Lindblad form (1) and (2) that can mimic the action of the original generator at stroboscopic instants of time, such that ${{ \mathcal P }}_{{\rm{F}}}=\exp ({{ \mathcal L }}_{\mathrm{eff}}T)$. There are three necessary (and altogether sufficient) conditions which any Lindblad generator has to fulfill: (i) trace preservation, $\mathrm{tr}[{{ \mathcal L }}_{\mathrm{eff}}A]=\mathrm{tr}[{\rm{A}}]$, (ii) Hermiticity preservation, ${({{ \mathcal L }}_{\mathrm{eff}}A)}^{\dagger }={{ \mathcal L }}_{\mathrm{eff}}A$, if ${A}^{\dagger }=A$, and (iii) so-called 'conditionally complete positivity' [57].

This is in distinct contrast with the case of a unitary evolution [9, 14]. The effective time-independent Hermitian operator H_eff can always be obtained as the logarithm of the unitary Floquet propagator ${U}_{T}={ \mathcal T }\exp \left[-\tfrac{{\rm{i}}}{{\hslash }}{\int }_{0}^{T}H(\tau ){\rm{d}}\tau \right]$. Moreover, not only its prinicpal branch but any branch of the logarithm yields a valid ${H}_{\mathrm{eff}}$. The non-unitary case is much more restricted: only the branch of the logarithm of a dissipative Floquet map which produces an operator possessing properties (i–iii) yields a legitimate Lindblad generator ${{ \mathcal L }}_{\mathrm{eff}}$ [57].

Condition (i) holds by default if one starts with a trace preserving map (which is our case). Condition (iii) can formally be checked with the algorithm given in [57]. However, it is hardly realizable in practice when $M\gt 2$ because it involves repeated solution of an ${ \mathcal O }({M}^{4})$ optimization problem within the mixed-integer semidefinite programming framework [58]. Condition (ii) is much easier to check.

Hermiticity preservation demands that the spectrum of ${{ \mathcal L }}_{\mathrm{eff}}$ is invariant under complex conjugation; in other words, it should consist of real eigenvalues or/and of complex conjugated pairs of eigenvalues. If there are negative real-valued eigenvalues (strictly speaking, of odd algebraic multiplicity) in the spectrum $\{{\lambda }_{j}\}$ of the map ${{ \mathcal P }}_{{\rm{F}}}$, it is impossible to fulfill the condition of the invariance of the spectrum of ${{ \mathcal L }}_{\mathrm{eff}}$ (which is a logarithmic branch of ${{ \mathcal P }}_{{\rm{F}}}$) under complex conjugation. This is because any branch of the logarithm of a negative real-valued number can neither produce a real number nor a complex conjugated pair. Figure 4 depicts the number of eigenvalues $\{{\lambda }_{j}\}$ with $\mathrm{Re}({\lambda }_{j})\lt -\varepsilon$ and $| \mathrm{Im}({\lambda }_{j})| \lt \varepsilon ,\varepsilon ={10}^{-7}$, as a function of the driving frequency ω. The dependence reveals that the condition is not fulfilled in the most interesting case of non-adiabatic and non-diabatic 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 parameter region.

Zoom In Zoom Out Reset image size

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, and estimate the relaxation time towards them. To illustrate this idea, we have applied the concept to a dissipative and periodically driven many-body model. We have studied the model both numerically exactly and, in the limit of a large particle number, within a mean-field description. The latter predicts bifurcations from regular to chaotic attractors as the interaction strength is varied. The analysis shows a strong dependence of quantum characteristics of the asymptotic non-equilibrium many-body state, such as the purity and the negativity, on the interaction strength, especially in proximity of bifurcations predicted by the mean-field theory.

It is interesting to contrast the idea of Floquet maps produced by time-periodic Lindblad generators, and an approximate Bloch-Redfield master, a well-known alternative to the Lindblad formalism [59]. Typically one starts from a bilinear coupling of the system to a heat bath of harmonic oscillators. The bath is characterized by its spectral properties. When the model Hamiltonian is time-periodic, it is possible, by assuming an Ohmic heat bath and following the Born–Markov ideology, to derive the so-called Floquet-Markov equation [11, 60]. This linear equation, similar to (1), governs the evolution of the system density operator; it is also local in time with a time-periodic generator. Formally it is thus possible to construct a corresponding Floquet map in this case as well. However, this so obtained map does not guarantee completely positivity for the evolution of the reduced density operator [61]; even more problematic is that it may even not necessarily assure the positivity of the reduced density operator on the way to its asymptotic limit, see [62–64] for detailed comparisons.

We conclude by pointing out possible research directions which may benefit from the use of Floquet maps within the Lindblad framework. It has been proposed to use time-periodic driving to create, for the situation with coherent Hamiltonian systems, effective topologically protected states [15, 65]. The important problems of the stability of these states against dissipation or their creation with a synthetic dissipation [66] 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 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 [67–69] and, very recently, the dynamics of open MBL systems was considered in [70–72]. 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 also a very recent [73].

The numerical simulations were supported by the Russian Science Foundation Grant No. 15-12-20029 (MI, SD and PH) and were performed on the Lobachevsky cluster of the University of Nizhny Novgorod. DP and PH acknowledge the support by Singapore MOE Academic Research Fund Tier-2 project (Project No. MOE2014-T2-2-119, with WBS No. R-144-000-350-112). SD and PH acknowledge support by the Deutsche Forschungsgemeinschaft (DFG) via grants (DFG) HA1517/35-1 (PH), DE1889/1-1 (SD). PH also acknowledges support by the Singapore Ministry of Education and the National Research Foundation of Singapore.