Towards Laser Control of Open Quantum Systems: Memory Effects
Laser control of Open Quantum Systems (OQS) is a challenging issue as compared to its counterpart in isolated small size molecules, basically due to very large numbers of degrees of freedom to be accounted for. Such a control aims at appropriately optimizing decoherence processes of a central two-level system (a given vibrational mode, for instance) towards its environmental bath (including, for instance, all other normal modes). A variety of applications could potentially be envisioned, either to preserve the central system from decaying (long duration molecular alignment or orientation, qubit decoherence protection) or, to speed up the information flow towards the bath (efficient charge or proton transfers in long chain organic compounds). Achieving such controls require some quantitative measures of decoherence in relation with memory effects in the bath response, actually given by the degree of non-Markovianity. Characteristic decoherence rates of a Spin-Boson model are calculated using a Nakajima-Zwanzig type master equation with converged HEOM expansion for the memory kernel. It is shown that, by adequately tuning the two-level transition frequency through a controlled Stark shift produced by an external laser field, non-Markovianity can be enhanced in a continuous way leading to a first attempt towards the control of OQS.
33.80.-b, 03.65.Yz, 42.50.Hz
pen Quantum Systems, Measure of Non-Markovianity, Coherence, Spin-Boson Model, Laser Control.
The theory of Open Quantum Systems (OQS) deals with non-unitary dynamics among quantum states of a sub-system coupled with some unobservable degrees of freedom in their surroundings [1, 2, 3]. This is actually a very common situation since no physical system can truly be considered as isolated. The additional degrees of freedom that make up the environment can take many physical forms such as phonons in solids, photons in cavities, charge fluctuations, or molecular vibrations. They are responsible for a wide variety of basic processes and, in particular, thermalization, energy and charge transfers, decoherence. As a consequence, the scope of OQS encompasses several disciplines of both fundamental and technological importance, from struggling against decoherence in superconducting qubits of quantum information, to electronic and proton transfers in flexible proteins. This is why OQS theory has attracted considerable attention both from physics and chemistry communities.
Understanding non-equilibrium dynamics in OQS where dissipation and coherence evolve simultaneously is by itself a great challenge, but even more important is the depiction of feasible quantum control strategies to optimize physical observables such as decoherence rates or efficient and fast charge transfers over large molecules. As an illustrative example, non-linear ultra-fast optics experiments have shown the possibility to protect from environmental fluctuations, on picosecond time scales, photosynthetic organisms in molecular light harvesting materials [4, 5]. This has even led to the emergence of a new interdisciplinary field: Quantum Biology . In this work we are concerned by a strong laser control scheme addressing OQS described in a standard topology of a central system of few degrees of freedom coupled to an external and much larger set of environmental modes, namely, the residual bath which damps the dynamics. This problem could naively be approached by a time-dependent wave function evolution incorporating both the central system and the bath. But its feasibility would be very soon affected with the increasing number of degrees of freedom. It is important to have in mind that even referring to high performance variational codes like MCTDH with multi-layer, time-adapted basis set expansions, but only a few hundreds of such vibrational modes can reasonably be taken into account . Stochastic methods would also be limited to comparable numbers of degrees of freedom . A strong field Floquet type of approach, as has been applied to H photodissociation by André Bandrauk and coworkers, would even be worse, as the resulting number of close coupled equations to describe the multiphoton excitation process would lead to computationally non tractable issues . One of the goals of OQS theory is precisely to avoid a full integration of all the degrees of freedom by describing the dynamics in a reduced Hilbert space, through a reduced density matrix. More specifically, such a reduction is obtained by tracing out the bath degrees of freedom from the full density matrix. In some cases the damping could be an almost pure, memory-less dissipation.The dynamics is then approximately described referring to Markovian master equations in the so-called Lindblad form [10, 11] extensively used in atomic physics, quantum optics, semiconductors. This approximation assumes an instantaneous recovery of the bath from the interaction with a continuous flow of information from the central system to its environment. But, when aiming at control OQS at different times and length scales, or different energy and temperature ranges, dissipation and decoherence dynamics could be such that the bath response to the central system can no longer be neglected, leading thus to non-Markovian behaviours characterized by a back-flow of information.
Two classes of control strategies can be devised: either acting on the central system alone, or on both the central system and its environmental bath. The present work is concerned by the simplest scenario where a strong static field produces a Stark shift affecting the eigen-energies of a two-level central system leading to different coupling schemes with the bath. A more complete scenario would be the one which takes advantage of the non-Markovian response of the bath to control the full system dynamics . Such an ultimate goal could presumably be reached by referring to some collective modes which guide the flow of information from the central system to the bath in a reversible manner . But we still expect that the Stark shift affecting the central system can be used as a first step for enhancing non-Markovian behaviours. Obviously a requisite of any control strategy exploiting non-Markovianity is a quantitative measure of such a behaviour. Several measures characterizing non-Markovianity have been proposed in the literature [14, 15, 16, 17], among them the volume of accessible states  and more recently the appearance of negative decoherence rates from several possible decoherence channels . The relative merits of such measures based on the quantum dynamical map (with a matrix representation ) of the reduced density matrix elements deserve interest. More specifically, the volume of accessible states within the Bloch sphere, given by the determinant of , evolves in time with a total rate . The time-dependence of , as opposite to a constant value, is a signature of memory effects. But an even more sensitive signature would be reached through the partial decoherence rates towards the different decoherence channels of the dynamics. Their calculation is based on, roughly speaking, a logarithmic derivative of , diagonalizing . These partial rates sum up to but could individually and temporarily be negative, as signatures of information back-flow from the bath to the central system. We are postponing their detailed analysis to a forthcoming paper.
The article is organized as follows. Section 2 is devoted to the theory describing the dynamics of a two-level system coupled to a bosonic bath, the so-called Spin-Boson model. A canonical (system Hamiltonian-independent) form for a time-local master equation is introduced, for a generic mapping matrix. This is done by closely following the derivation of Ref.. The system-specific reduced density matrix evolution and the resulting physical is calculated by numerically solving a Nakajima-Zwanzig type of evolution using hierarchical equations of motion (HEOM) up to convergence, for the memory kernel [18, 19, 20, 21, 22]. The results are discussed in Section 3 in terms of the time-evolution of an observable taken as the volume of states decaying with . The most important observation is the enhancement of the non-Markovian behaviour when the two-level system internal transition frequency is off-resonance with respect to the maximum of the spectral density measuring the frequency distribution of the system-bath couplings. The control knob is then the external field strength inducing the Stark shift that monitors the system transition frequency (from on- to off-resonance). At this stage, the control targets the optimization of decoherence rates in order to slow down the overall decay, which by itself opens important applications in protecting quantum information lost or obtaining long-lived alignment-orientation in molecular dynamics (although not described by a Spin-Boson Hamiltonian), for instance. Future perspectives on other control goals aiming at adequately designing collective modes in the bath could open the way for efficient charge or proton transfer processes in large molecular systems. These are mentioned in the conclusion, Section 4.
2 Theory and Method
2.1 Canonical Master Equation
In many situations, in particular in liquid and solid state physics, memory effects due to the non-Markovian character of the dynamics have to be taken into account [23, 24]. At this point, a natural question about a given experimental system is to which extent the relaxation effects lead to non-Markovian dynamics. In the past few years, the problem of how to measure non-Markovianity has sparked remarkable interest, starting with Refs. [25, 14, 15, 25]. Number of quantitative measures have been proposed since these initial attempts [16, 17, 23, 26, 27, 28, 28, 29]. In addition to the theoretical studies, recent experimental works [30, 31] have shown that it is by now possible to engineer OQS and to drive it from the Markovian to the non-Markovian regimes. However, the measures recently proposed in the literature are well-suited to theoretical or abstract situations, where for example the dynamics can be solved analytically. Very few papers have explored more complicated model systems where the time evolution of the OQS can only be achieved numerically [32, 33, 34]. A benchmark example in this category is the spin boson model which will be analyzed in this paper. We will consider for this model system a specific measure of non-Markovianity, namely the volume of accessible state space [17, 35]. The goal of this paragraph is to briefly summarize the computation of this measure. The reader is referred to the original papers for technical details.
The dynamics of a Markovian open quantum system is typically ruled by a quantum dynamical semi-group . The most general representation of such semigroup is given by the Lindblad-Kossakowski equation [10, 11] which can be expressed as:
where is the density matrix of the central system under study and a generator in Lindblad form. This latter can be written as:
being the Hamiltonian of the system, the coefficients the relaxation rates and the the time-independent Lindblad operators. Hereafter, we assume that is a density matrix of a finite-dimensional Hilbert space, typically of dimension two for a spin system. In this case, the dynamics is characterized by three relaxation rates, .
This formalism can be generalized to a time-local master equation where the generator depends on time while preserving its general form:
We stress here that both the relaxation rates and the Lindblad operators may now depend on time. The process remains Markovian if the rates are positive for any time , and presumably becomes non-Markovian otherwise .
Recent studies have investigated the issue of the description of the dynamics of non-Markovian systems which can be given either by a non-local master equation with memory kernel obtained, e.g., from the Nakajima-Zwanzig technique , or by a local in time equation [17, 36]. We give here a simple argument showing that any time evolution of the density matrix of the system can be ruled by a time-local equation. We introduce the quantum dynamical map of the system which satisfies:
Differentiating Eq. (4) with respect to time leads to:
If we assume that the dynamical map is not singular at time , i.e. that the determinant of the corresponding matrix is different from 0, then the inverse of can be defined and Eq. (5) can be rewritten as follows:
which is a master equation local in time. The non-Markovian character of the dynamics is then associated with the values of the relaxation rates of the generator .
is also referred to when mapping the density matrix on its corresponding Bloch ball using Pauli matrices together with the identity as an orthogonal basis set. The time evolution of the volume of this Bloch ball , namely the volume of accessible states, is given by the determinant of :
which constitutes another signature of non-Markovian dynamics . More precisely, the dynamics is said to be non-Markovian if:
for a given time . In this work we will merely focus on this signature, leaving the calculation of partial relaxation rates and their possible negativity for a forthcoming paper.
2.2 Hierarchical equations of motion
We consider a spin-boson Hamiltonian currently used to model a lot of processes such as electron or proton transfer, excitation energy transfer and qubit in a surrounding
are Pauli matrices ( and ) written in a zeroth order, so called diabatic basis labelled . is the central system interstate coupling in a diabatic representation and the vibrational modes are written in mass weighted coordinates. Adopting an adiabatic representation for the central system, by diagonalizing the system Hamiltonian results into the central system transition frequency , which is hereafter taken as the single parameter defining the central system. The system-bath coupling takes the form
where is a collective bath coordinate which couples to the system Hamiltonian by inducing fluctuations of the energy gap and is the generic central system coordinates. is the renormalization energy shifting the system energies. The key tool in dissipative dynamics is the reduced density matrix which is the partial trace of the full density matrix over the bath degrees of freedom.
Quantum information is exchanged between the system and the environment causing energy relaxation and decoherence.
Well known projection techniques leading to non-Markovian master equations either time non local in the Nakajima-Zwanzig approach [37, 38] or time local in the Hashitsumae-Shibata-Takahashi formalism [39, 40] provide effective equations for the relevant central system part
At the initial time (), the following factorization is assumed:
where is the density matrix of the bath at equilibrium.
The exact reduced density matrix can in principle be obtained from the Feynman-Vernon influence functional  and has been implemented in some simple cases [42, 43]. However the evaluation of the memory kernel , in practical applications, remain cumbersome for complex systems. For more numerical efficiency different strategies belonging to two main classes, projection techniques or hierarchy equations of motions, have addressed differential equations of motion. HEOM originally proposed by Tanimura and Kubo [18, 44, 19] is a powerful method providing a non-perturbative calculation of in Markovian or non-Markovian baths. Details about the derivation of the HEOM equations can be found in the original papers and in different reviews [45, 45]. It is born from the path integral method in the case of a harmonic bath allowing for analytical treatment of the influence functional, in particular the derivation of a hierarchy of time local coupled equations when the correlation function of the collective bath mode can be expanded as a sum of exponential functions. The correlation function of the collective coordinate is a Boltzmann average over the equilibrium bath ensemble defined by
where is the Heisenberg representation of the bath coordinate and the bath density matrix. The fluctuation-dissipation theorem relates the correlation function to a spectral density 
where and the Boltzman factor.
Different schemes have been proposed to decompose the spectral density  which is nothing but a frequency distribution of the system-bath coupling. We discuss the operational equations in the case where the spectral density is expanded as a sum of Lorentzian functions with Ohmic behaviour at low frequencies [47, 48]
being the coupling amplitude of the central system to the th labelled collective-bath mode. Referring to Cauchy residues theorem, the correlation function takes the form of an exponential series involving both the poles of and those of the Bose function through Matsubara frequencies  with complex and
The complex conjugate correlation function can be written by using the same but different expansion coefficients 
where being the total number of dissipative modes The reduced system density matrix is then the first element of a chain of auxiliary density matrices . The evolution is driven by time local coupled equations:
The hierarchy is built with levels. Each level corresponds to a given sum of the occupation numbers , of the dissipation modes. Level corresponds to the system matrix with the vector , i.e. all ’s are zero. Each density matrix of level is coupled to matrices of level . This very efficient algorithm only uses a single correlation function to describe the system bath interaction. Moreover, the structure is well suited for implementation on parallel computers . In the following Eq. (14) has been approximately solved at a given level of hierarchy corresponding to perturbation order, by an appropriate truncation of Eq. (22)
2.3 The Model.
As been discussed in the previous paragraph the central part of the Spin-Boson system is basically modeled by a single parameter two-level system, namely its internal transition frequency defining the energy gap between states and . Note that in some cases, this could result from the diagonalization of a primary so-called diabatic two levels and directly interacting through a potential coupling . The bosonic bath is described in terms of harmonic oscillators of frequency (). As for the coupling of the central system to the bath, it is taken into account referring to a spectral density, . It is worthwhile noting that, in the absence of an external field, the individual levels making up the central system are only indirectly coupled through their environmental bath. By analogy with a standard Fano model of two discrete levels facing and interacting with a continuum, one can refer to as a frequency representation of an energy-dependent discrete-continuum coupling scheme appropriately averaged on the density of levels of the discretized continuum. The dynamics using a full Fourier transform of the spectral density also implies an anti-symmetrical form for , that is:
Among the different functional forms that have been devised in the literature, we are hereafter referring to the so-called ohmic function, given by Eq. (19), retaining but only two Lorentzians (), enough to provide a well-structured form for this spectral density within the parameter range suggested in Ref.
The corresponding couplings and frequencies are gathered in Table 2.3, whereas the resulting spectral density and the corresponding correlation function, at room temperature are illustrated in Fig. 1(a) and 1(b), respectively.
Two important arguments advocate for this specific choice of :
i) A highly structured is in favor of enhanced memory effects in the bath response. In analogy with the Fano picture, such structures in could be attributed to some Feshbach resonances locally modifying the density of states. Actually they are supporting collective modes of the bath, and due to their possibly long enough lifetimes proportional to could temporarily trap the dynamics, leading to memory effects signatures.
ii) The central system transition frequency can be tuned in such an amount that it matches either one of the two maxima of the spectral density (case referred to as on-resonance) or, as an extreme situation, the minimum between the two peaks (case referred to as off-resonance). The off-resonance case is expected to produce the most important memory effects, as has already been pointed out in the literature  and as could be rationalized in terms of the back flow of information from the bath to the central system has to organized through the two neighbouring resonances. The specific two well separated peaks structure of the spectral density, offers a control flexibility by tuning the central transition frequency from on- to off-resonant cases. More specifically this could be achieved through an external laser field acting on the energy splitting of the central system via its transition dipole producing a controlled Stark shift.
2.4 HEOM convergence.
With the parameters that have finally been retained together with a bath temperature of , the overall dynamics is completely determined through the correlation function displayed in Fig. 1(b). It is interesting to note that shows about three damped oscillators of period , with almost negligible values for times exceeding , showing that memory effects can develop on ultra fast time scales. As has been explained in the previous theory section, this correlation function enters the memory kernel of the Nakajima-Zwanzig equation, which is further expanded in terms of successive HEOM levels (. Our purpose is now to check the numerical convergence of the dynamical calculations as a function of increasing level of this hierarchy. For doing this we retain two observables on the reduced density matrix , namely, its diagonal and off-diagonal terms , and the volume of the accessible measuring the decoherence of the central system. The calculations are carried out fixing which corresponds to an off-resonance case close to the minimum between the two peaks of the spectral density displayed in Fig. 1(a).
Fig. 2 gathers the results as a function of time, starting from a diabatic initial state where the whole population is in the central system ground state in the diabatic representation, all other matrix elements of being zero. The corresponding adiabatic picture, resulting from the diagonalization of , consists in taking all matrix elements . At HEOM level, often referred in the literature to be as converged for these rather low coupling regime (corresponding to the values of listed in Table 1), we obtain a non-monotonic decay of the volume with, at least, two clearly identified bumps occurring at and . Such bumps would be considered as convincing signatures of memory effects leading to non-Markovian behaviour. However, increasing HEOM level up to has as a consequence to partially wash out these bumps, fortunately still maintaining a non-monotonic decay especially observable for times longer than . Finally, convergence is assumed to be satisfactory, as HEOM level gives almost the same results (see Fig. 2) with reminiscences of the level bumps slightly shifted at times and . Concerning the diagonal and off-diagonal matrix elements of displayed in panels (b) and (c) of Fig. 2, the converged calculations differ from the level one by erasing the shoulder at . It is worthwhile noting that since , we have accordingly checked that , at any time . Off-diagonal elements asymptotically go to zero for the converged calculations. As a conclusion, for the forthcoming part of this work, we assume that is the HEOM level for which convergence is reached both for on- and off-resonance cases.
3 Results and Discussion
The volume of the accessible states is decreasing as a function of time according to a time-dependent total decay rate , following an exponential law:
as has previously been discussed . As opposite to a situation where the total decay rate is a constant, the above-mentioned law (Eq. (23)) induces a non-monotonic behaviour that could characterize non-Markovian dynamics. As has already been mentioned in the literature, a characteristic non-Markovian behaviour with a back-flow of information from the bath to the central system can be observed for temporarily negative values . Referring to converged calculations carried at HEOM level, we now examine two system-bath coupling situations, namely on- and off-resonant cases.
3.1 On-Resonant Case
The central system internal transition frequency is tuned such as to match the first peak maximum of the spectral density, that is as indicated in Fig.3(a). The results concerning the dynamical description of the reduced density matrix are depicted as the trace Tr in panel (c) and the real and imaginary parts of the off-diagonal elements in panel (b). The state-space volume is displayed in panel (d). We checked that the two curves resulting from Eq. (7) as the determinant of , or from (Eq. (23)) as the decaying law with , are perfectly superimposed. The initial values are 1 for Tr and the volume , and 0.5 for the off-diagonal elements . Basically two decoherence time scales are observed. The first, of about , concerns the decay of Tr and both showing rather monotonic behaviour, with a shoulder at about and a small amplitude oscillation for the imaginary part of at about . Tr finally reaches its asymptotic value of 0.5. The second time scale, much shorter, characterizes the dynamics of the volume of the accessible states which is decaying almost monotonously within about (panel d). This is not only a signature of a fast dynamical evolution (typically a vibrational period), but also of a presumably Markovian, or close to Markovian behaviour.
For this resonant case our conclusion is that the volume of accessible states space is not by itself a probe sensitive enough to detect a clear non-Markovian behaviour.
3.2 Off-Resonant Case
The central system internal transition frequency in now tuned off-resonance around the minimum of the spectral density (see Fig. 4(a)). As resulting from the discussion of the model, this situation is expected to support better marked non-Markovianity. As in the on-resonant case, two decoherence time scales are also observed here. The long one, characterizing Tr and , is still about . The difference being that both dynamics are not monotonic and the asymptotic value of Tr is now much larger (more than 0.8) showing a better preservation of the central system, and even more interestingly, an increasing behaviour starting from about (Fig. 4(c) and (b) and supporting a back flow of information from the bath. As for the observable taken as the signature of non-Markovianity, the volume of accessible states is displayed in Panel (d), with the second characteristic shorter time-scale. Contrary to the on-resonant case, a clear non-monotonic behaviour is obtained with a small bump at , the overall complete decay occurring at . In other words, the characteristic decoherence time has been slowed down from to .
For this off-resonant case, our conclusion is that, as expected the non-Markovianity is better marked on the volume total decay rate.
Finally, the control strategy, as suggested previously, is to adequately tune the central system internal transition frequency from on- to off-resonant cases, in order to increase the non-Markovian behaviour. This can be achieved using a strong external laser field inducing a Stark shift between the levels of the central system.
4 Conclusion and Outlook
Different laser control scenarios could be looked for when dealing with open quantum systems with the potentiality to be applied to a broad range of processes from biology (proton transfer in large proteins), to chemistry (charge transfers in donor-acceptor systems) and physics (molecular alignment-orientation, quantum information). Such scenarios can roughly be classified in two groups according to the situation where the control field is acting on the central two-level system alone (by modifying its internal characteristic transition frequency) or on the whole system including the bath (by modifying appropriate collective modes responsible for the above mentioned physical processes). In the present work we are only referring to the simplest control scenario of the first class with a control field just producing a Stark shift on the dressed levels of the central system with a consequence on its coupling scheme with the bath. To the best of our knowledge, this is one of the very first attempts toward control in such complex systems.
The model in consideration is a standard Spin-Boson Hamiltonian with a double-peaked spectral density function. The calculations are carried by solving a Nakajima-Zwanzig type master equation with a converged HEOM level for the memory kernel. More precisely, for the chosen parameters of the mode, the central system typically supporting a molecular vibrational mode coupled to a bath at an equilibrium room temperature, fully converged calculations need not less than level 4 for the HEOM hierarchy. The external field intensity is continuously varied to tune the transition frequency of the central system from on- to off-resonant situations, that are reached when this frequency matches respectively the maximum or minimum of the spectral density. The observable measuring the overall decoherence of the two-level system towards the bath taken into account is the time-dependent decay of the accessible states mapping the central system reduced density matrix. A non-monotonic behaviour is considered as a possible signature of non-Markovianity. We numerically show and rationalize that the on-resonant case corresponds to basically memory-less, fast dynamics, whereas the off-resonant configuration may bring into play better marked memory effects and longer decoherence times.
As a word of conclusion, a simple control can be exerted merely by an appropriate tuning of the transition frequency to enhance non-Markovian bath response. The resulting longer decoherence times may advantageously be exploited to enhance the efficiency and long period robustness of molecular alignment-orientation or qubits information preservation. Moreover, very challenging perspectives are now opening with more sophisticated control strategies aiming at optimizing and coherently interfering bath collective modes through appropriate combinations of the eigen-channels of the decoherence dynamics [51, 52, 53, 54, 55, 56]. Ultimately, fast and efficient charge or proton transfers in long protein chains could be considered. We are actively pursuing our investigations along these lines.
We acknowledge Prof. Christiane Koch, Prof. Christoph Meier, Prof. Arne Keller and Prof. Eric Charron for fruitful discussions. O.A. thanks the organizers of the Molecules and Laser Fields Symposium in honour of Andre Bandrauk, in Orford (QC), Canada, May 4-7, 2016, where part of the control scheme has been discussed. We acknowledge support from the ANR-DFG, under Grant No. ANR-15-CE30-0023-01. This work has been performed with the support of the Technische Universität München Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under Grant Agreement No. 291763.
-  H.P. Breuer and F. Petruccione The Theory of Open Quantum Systems Oxford University Press (2002)
-  U. Weiss Quantum Dissipative Systems World Scientific (1999)
-  I. de Vega and D. Alonso In Press: Rev. of. Mod. Phys (2016) arXiv:1511.06994.
-  A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga and M. B. Plenio Nat. Phys., 9 113-118 (2013)
-  S.Gélinas, A.Rao, A. Kumar, S. L. Smith, A. W. Chin, J. Clark, T. S. van der Poll, G. C. Bazan and R. H. Friend Science, 343 512-516 (2014)
-  G. D. Scholes, G. R. Fleming, A. Olaya-Castro and R. van Grondelle Nat. Chem., 3 763-774 (2011)
- Meyer et al.  H.-D. Meyer, F. Gatti and G. A. Worth Multidimensional Quantum Dynamics: MCTDH Theory and Applications Wiley (2009)
- Gaspard and Nagaoka  P. Gaspard and M. Nagaoka J. Chem. Phys., 111 5676-5690 (1999)
-  S. Miret-Artés, O. Atabek and A. D. Bandrauk Phys. Rev. A, 45 8056-8063 (1992)
-  V. Gorini, A. Kossakowski and E. C. G. Sudarshan J. Math. Phys., 17 821 (1976)
-  G. Lindblad Com. Math. Phys., 48 119 (1976)
-  C. P. Koch J. Phys. Cond. Mat., 28 213001 (2016)
-  I. Burghardt, R. Martinazzo and K. H. Hughes J. Chem. Phys., 137 (2012)
-  H.-P. Breuer, E.-M. Laine and J. Piilo Phys. Rev. Lett., 103 210401 (2009)
-  E.-M. Laine, J. Piilo and H.-P. Breuer Phys. Rev. A, 81 062115 (2010)
-  S. Lorenzo, F. Plastina and M. Paternostro Phys. Rev. A, 88 020102 (2013)
-  M. J. W. Hall, J. D. Cresser, L. Li and E. Andersson Phys. Rev. A, 89 042120 (2014)
-  Y. Tanimura and R. Kubo J. Phys. Soc. Japan, 58 101-114 (1989)
-  A. Ishizaki and Y. Tanimura J. Phys. Soc. Japan, 74 3131-3134 (2005)
-  M. Schröder, M. Schreiber and U. Kleinekathöfer J. Chem. Phys., 126 114102 (2007)
-  A. Ishizaki and G. R. Fleming J. Chem. Phys., 130 234111 (2009)
-  A. Shabani, M. Mohseni, H. Rabitz and S. Lloyd Phys. Rev. E, 86 011915 (2012)
-  H.-P. Breuer, E.-M. Laine, J. Piilo and B. Vacchini Rev. Mod. Phys., 88 021002 (2016)
-  Á. Rivas, S. F. Huelga and M. B. Plenio Rep. Prog. Phys., 77 094001 (2014)
-  Á. Rivas, S. F. Huelga and M. B. Plenio Phys. Rev. Lett., 105 050403 (2010)
-  R. Vasile, S. Maniscalco, M. G. A. Paris, H.-P. Breuer and J. Piilo Phys. Rev. A, 84 052118 (2011)
-  J. Liu, X.-M. Lu and X. Wang Phys. Rev. A, 87 042103 (2013)
-  X.-M. Lu, X. Wang and C. P. Sun Phys. Rev. A, 82 042103 (2010)
-  D. Chruściński and S. Maniscalco Phys. Rev. Lett., 112 120404 (2014)
-  J. T. Barreiro, M. Muller, P. Schindler, D. Nigg, T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller and R. Blatt Nature, 470 486-491 (2011)
-  B.-H. Liu, L. Li, Y.-F. Huang, C.-F. Li, G.-C. Guo, E.-M. Laine, H.-P. Breuer and J. Piilo Nat. Phys., 7 931-934 (2011)
-  G. Clos and H.-P. Breuer Phys. Rev. A., 86 012115 (2012)
-  U. Lorenz and P. Saalfrank Eur. Phys. J. D., 69 46 (2015)
-  P. Rebentrost and A. Aspuru-Guzik J. Chem. Phys., 134 101103 (2011)
-  E. Andersson, J. D. Cresser and M. J. W. Hall J. Mod. Opt., 54 1695-1716 (2007)
-  D. Chruściński and A. Kossakowski Phys. Rev. Lett., 104 070406 (2010)
-  S. Nakajima Prog.Theo. Phys., 20 948-959 (1958)
-  R. Zwanzig Physica, 30 1109-1123 (1964)
-  F. Shibata, Y. Takahashi and N. Hashitsume J. Stat. Phys., 17 171-187 (1977)
-  N. Hashitsumae, F. Shibata and M. Shingu J. Stat. Phys., 17 155-169 (1977)
-  R.P Feynman and F.L Vernon Annals of Physics, 24 118-173 (1963)
-  N. Makri J. Math. Phys., 36 2430-2457 (1995)
-  N. Makri J. Chem. Phys., 111 6164-6167 (1999)
-  Y. Tanimura Phys. Rev. A, 41 6676-6687 (1990)
-  R.-X. Xu, P. Cui, X.-Q. Li, Y. Mo and Y. Yan J. Chem. Phys., 122 041103 (2005)
-  Hao Liu, Lili Zhu, Shuming Bai and Qiang Shi J. Chem. Phys., 140 134106 (2014)
-  C. Meier and D. J. Tannor J. Chem. Phys., 111 3365-3376 (1999)
-  A. Pomyalov, C. Meier and D.J. Tannor Chemical Physics, 370 98-108 (2010)
-  A. Nieto Comp. Phys. Comm., 92 54-64 (1995)
-  J. Strümpfer and K. Schulten J. Chem. Theo. Comp., 8 2808-2816 (2012)
-  J.-S. Tai, K.-T. Lin and H.-S. Goan Phys. Rev. A, 89 062310 (2014)
-  D. M. Reich, N. Katz and C. P. Koch Sci. Rep., 5 12430 (2015)
-  S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, D. Sugny and F. K. Wilhelm Eur. Phys. J. D., 69 279 (2015)
-  A. Chenel, G. Dive, C. Meier and M. Desouter-Lecomte J. Phys. Chem. A., 116 11273-11282 (2012)
-  A. Chenel, C. Meier, G. Dive and M. Desouter-Lecomte J. Chem. Phys., 142 024307 (2015)
-  F. F. Floether, P. de Fouquieres and S. G. Schirmer New J. Phys., 14 073023 (2012)