Response functions as quantifiers of non-Markovianity

Quantum non-Markovianity is crucially related to the study of dynamical maps, which are usually derived for initially factorized system-bath states. We here demonstrate that linear response theory also provides a way to derive dynamical maps, but for initially correlated (and in general entangled) states. Importantly, these maps are always time-translational invariant and allow for a much simpler quantification of non-Markovianity compared to previous approaches. We apply our theory to the Caldeira-Leggett model, for which our quantifier is valid beyond linear response and can be expressed analytically. We find that a classical Brownian particle coupled to an Ohmic bath can already exhibit non-Markovian behaviour, a phenomenon related to the initial state preparation procedure. Furthermore, for a peaked spectral density we demonstrate that there is no monotonic relation between our quantifier and the system-bath coupling strength, the sharpness of the peak or the resonance frequency in the bath.

Introduction.-A central problem of non-equilibrium statistical mechanics is to obtain a closed dynamical description for some "relevant" degrees of freedom without the need to explicitly model the remaining "irrelevant" degrees of freedom. Within the theory of open quantum systems, the complete system state ρ S (t) is usually regarded as relevant while the bath is traced out [1,2]. Using the Nakajima-Zwanzig projection operator formalism, this can be done in a formally exact way, but unfortunately, initial system-bath correlations prevent the reduced dynamics from being closed due to the appearence of an inhomogeneous term.
We here show that within linear response theory it is possible (under certain conditions stated below) to obtain a reduced dynamical description for a set of system observables even in presence of an initially entangled system-bath state. Our findings will allow us to define a rigorous, yet very simple quantifier of non-Markovianity, which we can even express analytically for the Caldeira-Leggett model -a result which is very demanding to derive based on previous approaches [3,4].
Linear response theory.-We consider the standard system-bath setup and assume a global equilibrum state ρ SB (t 0 ) ∼ e −β(H S +H I +H B ) (where H S/I/B denotes the system/interaction/bath part of the Hamiltonian) prior to the "initial" time t 0 . We then suddenly perturb the system part of the Hamiltonian such that where the A i are system observables and the a i ∈ R (assumed to be sufficiently small) describe the respective strengths of the delta-kick δ(t − t 0 ). The purpose of the delta-kick is to generate a local unitary transformation U 0 = exp( i i a i A i ) ⊗ 1 B , which prepares the system in a nonequilibrium state at t 0 . The initial state has then, to linear order, expectation values Here, . . . β denotes an expectation value with respect to the global equilibrium state. Furthermore, we introduced the skew-symmetric matrix χ + with entries denotes the commutator. We remark that the bath state does not change during this preparation procedure and the system-bath correlations (as measured by the mutual information) also remain the same. In the following we consider only centered observables such that A i β = 0 without loosing generality.
The expectation value of A i at a later time t ≥ t 0 is then connected to the response function χ ij (t) ≡ i Θ(t) [A i (t), A j ] β via the Kubo formula (for an overview of linear response theory see Ref. [5]). In matrix notation we have where we introduced the mean value propagator G(t) ≡ χ(t)χ −1 + , which will be the central object of interest in what follows. Note that lim t 0 χ(t) = χ + . Eq. (3) amounts to our fundamental assumption in this paper as it is not guaranteed that the inverse of χ + exists (this is more thoroughly discussed in the supplementary material). If it exists, Eq. (3) describes a closed evolution equation for the mean values of the set of observables A i for all times t ≥ t 0 . Two properties of G(t) will be very important in the following. First, G(t) is independent of the initial state as it does not depend on any of the a i . Second, the propagator G(t) depends only on the elapsed time, which follows from the fact that the response function is expressed in terms of time-translationally invariant equilibrium correlation functions (CFs).
Therefore, if the system behaves Markovian, the mean value propagator must obey a condition which is also called divisibility. Equivalently, this implies for the response functions Finally, for later use we note that the response function also determines the temporal behaviour of the equilibrium CFs due to the fluctuation dissipation theorem (FDT). Out of the many possible forms of the FDT, we will only need where, in general, which also depend only on the time difference. Furthermore, we introduced the Fourier transformf (ω) ≡ ∞ −∞ dte iωτ f (τ ). Note that the FDT also fixes the real part of the response function via the Kramers-Kronig relation.
Comparison with previous approaches.-Before proceeding, let us contrast our approach with the conventional one. Arguably the common considered scenario in the theory of open quantum systems starts with an initial product state ρ S (t 0 ) ⊗ ρ B (t 0 ) [1][2][3][4] (for an exception see Ref. [6]). This bears the advantage that the inhomogeneous term in the Nakajima-Zwanzig equation disappears and the reduced dynamics of the system is described by a completely positive and trace preserving (CPTP) dynamical map where U is a unitary evolution operator acting on the joint system-bath state. Although it has been recently studied in greater generality whether it is possible to relax the initial product state assumption [7][8][9][10], the family of initially entangled states considered above will in general not give rise to a CPTP map. Therefore, there is no direct connection between our approach and previous results, although we can draw analogies. Indeed, while Φ(t, t 0 ) as G(t) is independent of the initial system state, the former does not propagate mean values but the complete system state ρ S (t 0 ) to arbitrary later times t ≥ t 0 . It is interesting to ask whether G(t) can be extended to a dynamical map for the entire system density matrix by looking at a complete set of system observables {A i }, whose expectation values are isomorphic to ρ S (t). In the supplementary material we demonstrate that this is not possible because χ + in Eq. (3) becomes non-invertible.
Finally, to characterize non-Markovianity within the standard approach based on Eq. (7), the concept of CP divisibility is important. A CP divisible quantum stochastic process is characterized by a family {Φ(t 2 , t 1 )|t 2 ≥ t 1 ≥ t 0 } of CPTP maps, which obeys analogous to the classical Chapman-Kolmogorov equation. Consequently, if a process is CP divisible, then the evolution of the density operator is Markovian (although there seems to be less agreement about the converse statement [3,4]). Based on this concept or a related notion, various quantifiers of non-Markovianity have been recently put forward [11][12][13][14][15][16] and direct experimental evidence is also accumulating [17,18].
Unfortunately, evaluating non-Markovianity for time evolutions generated by Eq. (7) is demanding as it requires, e.g., optimization procedures, the inversion of dynamical maps or the integration over complicated disconnected domains. In part, this problem is caused by the fact that the CPTP map Φ(t, t 0 ) has a complicated time dependence: even if the unitary U in Eq. (7) is generated by a time-independent Hamiltonian, the dynamical map is not time-translational invariant, i.e., Φ(t, t 0 ) = Φ(t − t 0 ) in general. This is in strong contrast to our result in the linear response regime, where G(t) always depends only on the elapsed time and which allows us to check the simpler condition (4) instead of Eq. (8).
Distance quantifier.-To introduce new quantifiers of non-Markovianity within our approach, we need to quantify the distance between two functions f (t) and g(t). We use the standard L 2 scalar product f, g = ∞ −∞ dtf (t)g * (t) and the induced norm f = f, f , where it is tacitly assumed that the integrals are converging. We then define the distance By Cauchy-Schwarz' inequality 0 ≤ D(f, g) ≤ 1 and D(λf, λg) = D(f, g) for any λ ∈ C, i.e., the difference has the favourable properties that it is positive, bounded and independent of any global scaling. By analogy with the Euclidean scalar product, D(f, g) = | sin(φ)| can be seen as quantifying the "angle" φ between the two vectors f (t) and g(t). Most importantly for our applications, by Parseval's theorem we can deduce that D(f, g) = D(f ,g), where the right hand side is computed by using the L 2 scalar product in Fourier space, . New quantifiers of non-Markovianity.-It will be advantageous to work in Fourier space in the following. In the supplementary material we show that integrating Eq. (5) over s from zero to t implies in Fourier space Then, to measure violations of Eq. (10) as a consequence of the (assumed) divisibility property, we propose the As a second quantifier of non-Markovianity, we also check the validity of the regression theorem (RT) [19,20], which allows us to relate the evolution of CFs to the evolution of mean values. Within our setting the Markovian assumption enters here by using that Eq. (3) holds for all initial states and that there exists a dynamical map Φ(t, t 0 ), which is independent of ρ S (t 0 ). It is worth emphasizing that the validity of the RT does not a priori rely on an initial product state assumption or on the property of CP divisibility. It merely signifies that it is possible to find for any initial system state a map G(t) to propagate the mean values and -in addition to what is required to evaluate N Here, we have added the superscript "RT" to emphasize that this is the predicted CF assuming the validity of the RT. Note that denotes in general an out-of-equilibrium CF, but we will be only interested in equilibrium CFs which we denote with a calligraphic C. For them we can deduce in Fourier space that (see supplementary material) We add that the behaviour of CFs (often in relation with the validity of the RT) has played an important role historically to define a quantum Markov process [19][20][21][22][23] and was also investigated in the recent debate about non-Markovianity in Refs. [24][25][26]. However, its use in the linear response regime has not been noted before, although it is well-known that all quantum systems violate the RT in that regime [27]. Then, based on the comparison of the exact equilibrium CFs [obtained from the FDT (6)] and their Markovian prediction [obtained from the RT (13)], we propose We here assume that the equilibrium covariance matrix C(0) is exactly known such that the prediction (13) uses the correct initial value.
To conclude, the magnitude of both, N ij and N (2) ij , measures by how much we fail by naively assuming that the process is Markovian. They can be computed without the need to a priori derive any quantum master equation -only the knowledge of the linear response functions or the equilibrium CFs is required.
We will now treat an important class of open system models with Gaussian dynamics exactly, i.e., without any approximation about the temperature of the bath, the system-bath coupling strength or the spectral features of the bath. We also remark that for this class our results are valid beyond linear response. Related studies about non-Markovianity of Gaussian dynamics based on different approaches and various approximations can be found in Refs. [14,18,[28][29][30][31].
Quantum Brownian motion.-We consider the standard Caldeira-Leggett model with Hamiltonian (in suitable mass-weighted coordinates) Here, q and p refer to the position and momentum of the system with frequency ω 0 > 0, whereas the bath oscillators with frequencies ω k > 0 are specified with an additional index k. Furthermore, c k denotes the coupling strength between the system and the k'th oscillator. Of central importance is the spectral density (SD) It characterizes the coupling between system and bath and it is assumed to be a continuous function of ω in the limit of a large bath fulfilling J(0) = 0 = J(ω → ∞).
A great benefit of the Brownian motion model is that almost all quantities of interest are computable in closed form [32,33], e.g., the matrix of response functions reads (see supplementary material for a derivation) whereγ(ω) is the Fourier transform of the memory kernel γ(t) = Θ(t) 2 π ∞ 0 dω J(ω) ω cos(ωt). In view of the general theory outlined above, our set of system observables will be the position and momentum of the system, {A 1 , A 2 } = {q, p}, and one easily verifies that χ + is symplectic with (χ + ) pq = 1 = −(χ + ) qp . The delta-kick now creates the unitary which shifts the position and momentum operators, thereby shifting the mean values but leaving the covariances unchanged. Furthermore, the equilibrium covariance matrix is diagonal with entries We now have all quantities at hand to compute our quantifiers. For the rest of the paper we will set t 0 = 0. Classical Ohmic limit.-We consider the simplest case of a classical particle ( = 0) coupled to an Ohmic bath, which corresponds to a memory kernel of the form γ(t) = Dδ(t). This follows from a linear SD J(ω) = Dω in the limit of an infinitely high cutoff frequency. The resulting Langevin equation for the system reads (see supplementary material for a detailed derivation) Here, the noise obeys ξ(t) = 0 and ξ(t)ξ(s) = γ(t − s)/β with the crucial requirement that . . . refers to an average over an initial conditional equilibrium state of the bath [33][34][35][36] Here, the position q(0 − ) of the Brownian particle prior to the delta-kick is a random variable distributed according to a Gaussian P [q(0 − )] ∼ e −βω 2 0 q(0 − ) 2 /2 such that, shortly before the unitary kick, the global system-bath state is in equilibrium.
If we would not disturb the state, a q = 0 and a p = 0, and Eq. (21) reduces to the standard Langevin equation. However, the presence of the unitary kick results is an initial system state described by a shifted Gaussian P [q(0)] ∼ e −βω 2 0 [q(0)−aq] 2 /2 while the bath still resides in the state (22). The fact that the bath has no time to adapt to a new conditional equilibrium state causes non-Markovian behaviour as we can rigorously show with our quantifier. For instance, in view of Eq. (10) we find that which is clearly non-zero and only becomes negligible in the weak coupling regime, see Fig. 1. The subtle importance of the initial state preparation procedure for the validity of the Langevin equation was already noted in Ref. [34][35][36], but it had not been rigorously quantified.
We remark that it is a special property of the Caldeira-Leggett model that the first moments do not depend on . This changes for CFs, which depend on and have a non-trivial dependence on the inverse bath temperature, see Fig. 1 again.
Peaked SD.-We now turn to a non-trivial case described by the SD This corresponds to the SD felt by a system, which is coupled with strength D to another harmonic oscillator of frequency Ω, which is in turn coupled to an Ohmic bath with SD Γω [37]. Note that the parameter Γ controls the structure of the SD: a small Γ corresponds to a sharp peak around the frequency Ω whereas a larger Γ smears out the peak resulting in an increasingly flat SD. Furthermore, the real part of the Fourier transformed memory kernel is [γ(ω)] = J(ω)/ω and the imaginary part becomes (see supplementary material for more details) In practical considerations, non-Markovian behaviour is often associated with a strong system-bath coupling and a structured SD [2]. Thus, one would intuitively expect that the degree of non-Markovianity increases for larger D and smaller Γ and that it reaches a maximum, if the system is on resonance with the oscillator in the bath, i.e., if ω 0 ≈ Ω. As Fig. 2 demonstrates, this intuition is not always correct. We observe that there is no simple (i.e., monotonic) relation between our quantifier of non-Markovianity and the parameters D, Γ and |ω 0 − Ω|. In fact, one could ask whether this results from the particular definition (9) and (11) which we have used and which always entails a certain level of arbitrariness. Therefore, we have also plotted the time-evolution of the observable p(t) in Fig. 2 (bottom row), whose deviation from an exponentially damped oscillation seems to be roughly in agreement with our quantification scheme.
Summary.-This work shows that it is possible to quantify non-Markovianity in the linear response regime in a rigorous and straightforward manner. Since we can only treat initially correlated states, our approach is rather "orthogonal" to previous ones, but for many scenarios of experimental interest this might be indeed a more realisitic assumption. Furthermore, for the Caldeira-Leggett model our quantifier is valid beyond linear response and can be expressed analytically in terms of an integral over known functions. We have then shown that even a classical particle coupled to an Ohmic bath can behave non-Markovian depending on the initial state preparation procedure and that one should not expect a simple relation between our quantifier of non-Markovianity and parameters in the SD or the temperature of the bath.

SUPPLEMENTARY MATERIAL
This appendix contains in the following order: • A detailed discussion about the existence of the mean value propagator G(t).
• A detailed discussion and derivation of the RT (12).
• A derivation ot the matrix of linear response functions (17).
• A derivation of the generalized Langevin equation (21).
• A derivation of a general relation between the Fourier transformed memory kernel and the spectral density as well as a derivation of Eq. (25).
For notational convenience, we will often denote Existence of G(t) The existence of G(t) is guaranteed if the matrix χ + can be inverted. To approach the problem, we first consider the simplest case where the perturbation is caused only by a single observable A 1 . Then, χ + is the scalar Clearly, in this case χ + can never be inverted.
To study the Caldeira-Leggett model in the main text, a perturbation caused by two observables A 1 and A 2 turns out to be sufficient. In fact, one can in general hope that this is always the case because Hence, if we choose A 1 and A 2 such that [A 1 , A 2 ] β = 0, the existence of the inverse of χ + is guaranteed. As an additional example we consider an open two-level system (TLS) such as the spin-boson model, but the particular form of the environment is unimportant for the present reasoning. Let us choose a basis such that the reduced equilibrium state of the TLS is aligned along the Paulimatrix σ z such that σ z β = 0. Then, if we decide to perturb the system using the observables A 1 = σ x and A 2 = σ y , we obtain which is invertible. Finally, it is interesting to ask what happens if we choose a basis set of N 2 − 1 observables in case of an N -dimensional quantum system. As the expectation values of these observables is isomorphic to the complete system density operator, this would provide us with a full dynamical map for the system and questions related to its complete positivity would then be relevant. Let us start with the TLS from before and let us add A 3 = σ z to the set of observables {A 1 , A 2 }. Then, It is easy to confirm that the determinant of this matrix is zero, and hence χ + is not invertible. Indeed, this will always be the case. To demonstrate this we show that there are in case of a complete set of N 2 − 1 observables multiple a's, which give rise to the same state of the system and hence, χ + maps many to one and is not invertible. For this purpose we introduce the effective Hamiltonian H eff which describes the reduced Gibbs state of the system prior to the perturbation, where Z is the global partition function. As the set A i is assumed to be a complete set of observables, we can actually expand the effective Hamiltonian, too, for some real-valued parameters b 0 , b 1 , . . . b N . Let us now choose a = (b 1 , . . . , b N ) where is some constant choosen small enough such that a can be safely treated within linear response. Then, the net effect will be that U 0 and H eff commute such that we obtain the same initial state ρ S (t 0 ) = tr B {U 0 e −β(H S +H I +H B ) U † 0 }/Z = e −βH eff independent of the choice of . Since there are multiple possible, this proves that χ + can not be invertible.

Derivation of Eq. (10)
The assumed divisibility property (5) depends explicitly on the choice of s ∈ [0, t]. To get a simple and average quantifier, we decide to integrate Eq. (5) over s from zero to t. This yields Now, we recognize that due to the Heaviside step functions involved in the definition of χ(t), we can write the right hand side as Fourier transformation then immediately yields as claimed in the main text.

Derivation of the regression theorem
We here state, prove and discuss the quantum regression theorem following a rather general approach. The discussion of the RT in the current literature is in fact often motiviated by a particular physical situation.
Theorem. If there exists a set of system observables {A i } such that their dynamics is closed, i.e., if for all times t ≥ t 0 and for some propagator G(t, t 0 ) independent of the initial state ρ S (t 0 ), and if there exists a dynamical map Φ(t, t 0 ) independent of the initial state ρ S (t 0 ), then for all times t ≥ t 0 .
Proof. According to quantum mechanics, the expectation value of a system observable is given by 1 We add that the existence of some CPTP map Φ(t, t 0 ), which maps the initial state ρ S (t 0 ) to the final state ρ S (t), is always guaranteed [38], but in general its construction is highly non-unique and depends on ρ S (t 0 ). We now use our two assumptions. Namely, by comparing Eq. (38) with Eq. (36) and by using that Φ(t, t 0 ) is the correct dynamical map for any initial system state, we conclude that must hold for any ρ S (t 0 ). By using that where Φ * (t, t 0 ) is the adjoint dynamical map in the Heisenberg picture, we obtain with the help of Eq. (39) the chain of equalities which is the RT (37). QED.
Remarks. From the proof above it becomes evident that we neither had to use the quantum Chapman-Kolmogorov equation (8) nor any product state assumption for the system-bath state. Physically, of course, we know that deriving a dynamical map Φ(t, t 0 ) is often admissible only for factorized initial state. However, also purely classically correlated initial states yield to a CPTP map Φ(t, t 0 ) [7] and even in presence of quantum correlations it is sometimes possible to derive such maps [8][9][10]. Therefore, it is important to distinguish the question When can we derive Φ(t, t 0 ) in a physical setting? from What does an assumed existence of such a Φ(t, t 0 ) imply?
Within our context the validity of the RT and the value of the corresponding quantifier N (2) ij in Eq. (14) give us key information about the question whether knowledge of the initial system state and its evolution suffices to infer the evolution of the (equilibrium) CFs, but it does not reveal any more insights.
Moreover, as Figs. 1 and 2 explicitly demonstrate, the RT can even fail in the classical limit β → 0. The reason for that can be traced back to the fact that the propagator G(t) in Eq. (3) is only well-defined for a i = 0 and cannot be used at equilibrium. Thus, our first assumption that G(t) is the correct propagator for the mean values for any initial system state ρ S (t 0 ) was clearly wrong.

Equilibrium correlation functions from the regression theorem
At equilibrium we have the symmetry A i (−t)A j β = A i A j (t) β and using the RT we obtain C RT (−t) = C(0) T G(t) T for t > 0 where T denotes the transpose and C(0) the initial equilibrium covariance matrix. In terms of the linear response function, we can also write C RT (−t) = −C(0) T χ −1 + χ(t) T where we used the skew symmetry χ + = −χ T + . Therefore, Fourier transformation yields For Hermitian observables it now holds true that χ(−ω) =χ * (ω) such that we obtaiñ as claimed in the main text. More explicitly, using Eq. (17), we get Note the symmetry relationC RT qp (ω) =C RT pq (ω) * . This has to be compared with the true CF, which follows from Eq. (20) and the fluctuation dissipation theorem (6): (47)

Response function of the Brownian motion model
We here derive Eq. (17) using the well-known result for the position-position response function [32,33] χ qq (ω) = 1 To derive the other response functions, one uses (i) that where H denotes the unperturbed Hamiltonian (15) of the Caldeira-Leggett model, (ii) partial integration and where we explicitly wrote out the "dependence" of the operators on the initial time t 0 = 0 in the Heisenberg picture. Then, in detail the other response functions can be derived as follows: where we used lim t→∞ [q(t), q(0)] β = 0 and [q(0), q(0)] = 0 at the end. Analogously, one also derivesχ pq (ω) = −iωχ qq (ω). To derive the corresponding result forχ pp (ω), two partial integrations are necessary: Generalized Langevin equation We here go through the derivation of the generalized Langevin equation for the Caldeira-Leggett model (15) with an arbitrary linear system force added to it, Our treatment is only slightly more general than the one in Ref. [33], where a momentum-dependent force is not considered. Below we will assume that right before t 0 = 0 the system and the bath is in a global equilibrium state ω(β) and we explicitly require that the perturbation is switched on at times t ≥ 0, i.e., f q (t) = 0 = f p (t) for t < 0. The equations of motion for the positions and momenta of the system and the bath become (note that there form is identical classically and quantum mechanically) p k (t) = −ω 2 k q k (t) + c k q(t).
The last two equations are formally solved by q k (t) = q k (0) cos(ω k t) + p k (0) ω k sin(ω k t) It turns out to be convenient to rewrite the integral using partial integration: Note that in presence of f p (t) it is in general not valid to replaceq by p under the integral. With the help of this formal solution we obtaiṅ p(t) = − ω 2 0 q(t) + f q (t) + k c k q k (0) cos(ω k t) + k which obeys in the classical case ξ(t) = 0 and ξ(t)ξ(s) = γ(t − s)/β, where the double bra-ket notation refers to an average over the conditionally equilibrated bath state in Eq. (22). Using also the definition of the SD and the memory kernel, we can write which allows us to arrive at the more compact expressioṅ This equation generalized Langevin equation reduces to Eq. (21) in the limit of an Ohmic SD.