Exploring complete positivity in hierarchy equations of motion

We derive a purely algebraic framework for the identification of hierarchy equations of motion that induce completely positive dynamics and demonstrate the applicability of our approach with several examples. We find bounds on the violation of complete positivity for microscopically derived hierarchy equations of motion and construct well-behaved phenomenological models with strongly non-Markovian revivals of quantum coherence.


Introduction
Every quantum system inevitably interacts with its environment. This interaction contributes to the emergence of classicality, represents the major challenge in the development of quantum information technological hardware [1], but can also be exploited for advanced means to control quantum systems like laser cooling [2] or dissipative state preparation [3].
The central difficulty in describing open quantum systems resides in the mere size of an environment. Most approaches aim at the reduction of such environments to their aspects that are most relevant for the system dynamics. This ranges from effective descriptions in terms of the system's degrees of freedom only [4,5], via approaches that explicitly include the most relevant environmental degrees of freedom [6,7,8] to numerically expensive treatments that aim at a potentially exact description [9,10,11,12,13]. In practice, one needs to take a compromise between accuracy and numerical effort; a comparatively simple effective description can help to develop a more intuitive model, but may also fail to identify more subtle features.
A completely different approach is based on a phenomenological modelling of open quantum dynamics that promises a better understanding and deeper insights into the underlying physics. The hallmark is provided by the Kossakowski-Lindblad equation [4,5] that allows us to construct models that respect the probabilistic interpretation of quantum mechanics: Markovian Lindblad master equations ensure the property of complete positivity such that all quantities that are interpreted as probabilities are always non-negative. Furthermore, it permits an interpretation of the environmental effects in terms of rates and elementary processes. The intuitive comprehension that comes along with this model suffers the disadvantage of a strongly limited applicability. Many systems that are actively investigated, like quantum dots [14], light-harvesting systems [15], or photonic band-gap materials [16], are characterized by a rather strong system-environment coupling giving rise to non-Markovian effects like the back-flow of information that are not captured by this model.
For that reason, it is desirable to find a framework for the description of non-Markovian processes that ensures complete positivity and is in the best case even phenomenologically accessible. Complete positivity is not only crucial for physically sensible dynamics, but understanding the very conditions under which the equation of motion induces completely positive (CP) processes must also be considered a crucial step towards a phenomenological understanding. However, the addition of memory to the system dynamics gives rise to substantial complications [17] in particular with respect to complete positivity. Most approaches to this problem focus on integral equations over a memory kernel and for certain classes like semi-Markov processes [18,19], collisionmodel-based approaches [20] or continuous time quantum random walks [21] conditions for complete positivity have been found. Some of these classes can also be covered by generalized sufficient conditions that have been obtained recently [22]. Eventually, the post-Markovian master equation [23] is capable to describe CP non-Markovian dynamics by interpolating between the generalized measurement interpretation of the Kraus operator sum and the notion of a continuous measurement for Markovian processes but eventually also imposes conditions on the memory kernel.
In this article we propose a different take on the question of CP non-Markovian dynamics. Rather than by integral forms over a memory kernel, our approach is inspired by the framework of hierarchical equations of motion (HEOM) [9] that employ -in addition to the system density matrix 1 (t) -a set of auxiliary operators i (t) (i = 2, . . . , n) that contain information about the environment and system-environment correlation. Equations of motion of the forṁ i (t) = n j=0 L ij j (t) , (i = 1, 2, . . . , n) (1) have been derived based on microscopic models for a wide range of systems including light-harvesting complexes [24,25], molecular transistors [26], as well as electrons in tunneling junctions [27] and organic heterojunctions exposed to a phonon bath [28]. In particular they have proven to be very reliable in the regime where comparable intra-system and system-environment interactions prevent approximations based on separations of time-scales. We will start out assuming the structure (1) with a finite-dimensional system state 1 (t) and target the identification of conditions on L ij such that the dynamical map Λ 1 (t) with t ≥ 0 defined via 1 (t) = Λ 1 (t) 1 (0) is CP. In the case of n = 1, i.e. equations of motion for the system state only, it is well established that this condition is satisfied if and only if L 11 can be written in Lindblad form [4,5]. A generalization of this criterion to hierarchy equations with n > 1 that also permits an incorporation of non-Markovian dynamics is, however, not known.
In this article, we propose a new method by means of which sufficient conditions for CP dynamics can be derived for HEOMs. In some cases, such conditions can be obtained analytically, whereas highly efficient numerical methods are in place whenever an analytical solution is not feasible. In section 2.1, the formal and algebraic framework is set before a very intuitive geometrical interpretation of the procedure is given in section 2.2. In section 2.3, we provide practical tools that ease an analytical examination as it is carried out for a couple of examples in section 2.4. More advanced techniques that aim for a numerical investigation of more complex problems are subject to section 3, where we introduce a transformation that permits a universal applicability of the method (section 3.1) and demonstrate how semi-definite programming can be employed for a highly efficient numerical examination of HEOMs (section 3.2). A truncation of a HEOM is often carried out in order to obtain approximate dynamics but can also give rise to equations that violate complete positivity. That is why in section 4.1 and 4.2, our approach is modified such that the binary classification into CP and non-CP dynamics is replaced by a quantification of the violation of complete positivity. This permits the application of our method to an important example of a truncated HEOM that describes the spin-Boson model in section 4.3, before we conclude our results in section 5.

Deriving conditions for complete positivity
Since complete positivity is a property of the dynamical map rather than of the system state, it is helpful to re-formulate (1) in terms of time-dependent maps Λ i (t) that satisfy i (t) = Λ i (t) 1 (0). The equations of motion for these dynamical maps are obtained from (1) and readΛ This set of equations can be expressed in a very concise form, when we introduce the extended dynamical map Λ = n i=1 Λ i (t) ⊗ |i) with {|i)} being the orthogonal and normalised basis of an n-dimensional Hilbert space. This notion of extended objects that are formed by operators is highly convenient and will in the course of this article be denoted by bold-faced symbols. The scalar product on this extended space reads for any two extended vectors Λ and Γ. In the same spirit, we also introduce the extended state (t) = n i=1 i (t) ⊗ |i) and define the extended generator L = n i,j=1 L ij ⊗ |i)(j|, which permits to express (2) in the compact form By means of this notation we are in the position to concisely describe the algebraic framework for the derivation of conditions for completely positive dynamics. In practice, it is often convenient to work with Bloch-type representations, where the Bloch vector i b i ⊗ |i) of is defined in terms of the Bloch vectors b i of the operators i , i.e. vectors with components [b i ] j = tr(σ j i ) for j = 0, x, y, z, where σ x , σ y and σ z are the Pauli matrices and σ 0 = I, or their generalisation to higher dimensional systems. For extended maps Λ and generators L the corresponding representations are defined analogously.

The algebraic framework
Let us first detail the concept of valid maps. An extended dynamical map Λ(t) is valid if and only if the system map Λ 1 (t) is CP, i.e. it can be written as with a positive semi-definite matrix χ(t) and a set of mutually orthonormal operators µ i . Showing complete positivity of Λ 1 (t) is thus equivalent to proving positive semidefiniteness of χ(t).
The initial map Λ 1 (0) at the time t = 0 is given by the identity map such that χ(0) has rank one. The rank of a matrix provides information about the number of non-vanishing eigenvalues and does therefore also indicate the order of the highest nonvanishing elementary symmetric polynomial. Those elementary symmetric polynomials read where r × r is the dimension of χ, and when χ has rank k, then e j (χ) = 0 for all j > k.
For the initial matrix χ(0) one obtains e 1 (χ(0)) = 1 and e j (χ(0)) = 0 for all j > 1 but as time evolves, the eigenvalues of χ(t) will typically take non-vanishing values such that the rank of χ(t) increases. We will now consider a time t p > 0 at which all non-trivial (i.e. not permanently vanishing) eigenvalues of χ(t) are strictly positive. This is equivalent to χ(t p ) being positive semi-definite and having rank h ≡ max(rank(χ(t)) | t ≥ 0). Any sign-flip of any eigenvalue of χ(t) for t > t p is then necessarily accompanied by a vanishing polynomial e h , for why it is sufficient to ensure that e h (χ(t)) does never become equal to zero for t > t p in order to show complete positivity from time t p on. If no eigenvalue of χ(t) remains equal to zero (i.e. h = r), then e h is given by the determinant but even when one or more eigenvalues of χ(t) vanish permanently, the proposed method will still be applicable with e h being an elementary symmetric polynomial with lower degree. In order to prove e h (χ(t)) > 0 for t > t p , it is assumed (the assumption will be justified in section 3.1) that one can find a non-linear transformation of the coordinates Λ(t) such that e h (χ(t)) is quadratic in the extended map Λ(t) and can be written as with a suitable operator S and satisfy S = S † . As we will show in the following the dynamical map Λ 1 (t) is CP for all times t > t p , if there exists an operator R = R † that (i) fulfils the operator inequality L † R + R L ≤ 0, To prove this statement, we employ (ii) to reformulate (7) to Since e h (χ(t p )) is strictly positive by assumption, this implies the inequality Condition (iii) allows to bound this further to Due to (i) the time derivative of Λ(t), R Λ(t) is always non-positive so that the right-hand-side of (10) grows monotonically, which implies that e h (χ(t)) > 0 for t > t p > 0.
Whereas χ(t p ) is (by assumption) positive semi-definite and of rank h, the ultimate goal is to show complete positivity for all times t ≥ 0 and to also take into account the initial condition e h (χ(0)) = 0. To this end it is necessary to verify that that all non-trivial eigenvalues of χ(t) become positive in leading order in t. If this is the case, and Λ(t), R Λ(t) is not strictly constant (i.e. it decreases in leading order in t following (i)), then (ii) can be replaced by To prove the latter statement, we first find t p > 0 such that χ(t p ) ≥ 0 and e h (χ(t p )) > 0, as well as Λ(t p ), R Λ(t p ) < Λ(0), R Λ(0) , and Λ 1 (t) is CP for all 0 ≤ t ≤ t p . To show complete positivity for t > t p , we rewrite (7) trivially as Following (ii') this implies where the first two terms satisfy This proves strict positivity of e h (χ(t)) as shown above in (9).

Geometric interpretation
The proposed procedure can be intuitively understood in geometrical terms: Just like normalized, Hermitian operators are quantum states (i.e. positive semi-definite operators) of a two-level system if and only if they are represented by points inside the Bloch sphere, also valid extended dynamical maps form a geometrical object which we denote by P. Although the dimension of P is substantially larger than the dimension of the Bloch sphere, and its shape will typically be more complicated, this underlying geometry permits an immediate understanding of our approach. The schematic representation of figure 1 depicts all the points within P in blue and yellow.
In turn, all points outside this area correspond to invalid extended maps with a system component that is not CP. The dynamics induced by (4) is indicated by arrows. The goal is to find conditions which guarantee that, for a given initial condition Λ(t p ) within P (indicated by a black circle in figure 1), the dynamics is such that its trajectory will never leave P.
As this can typically not be shown directly by means of P, one strives for the identification of a region V, such that V contains Λ(t) for all times t > t p V lies inside P.
The arrows show the dynamics according to (4). The thick trajectory corresponds to the dynamics that starts at Λ(0) (black cross) and passes through the point Λ(t p ) (black circle). Equipotential lines of the function F(Λ) = Λ, R Λ are depicted in green (dotted). The line that contains Λ(t p ) is depicted with long dashing and confines the set V (yellow). If V is contained in P (as depicted), and F(Λ(t)) decreases monotonically, CP dynamics is verified for t > t p .
This region V shall be confined by the equipotential line of a function F(Λ) ≡ Λ, R Λ (figure 1 depicts some equipotential lines of F in green (dotted/dashed)) that is determined by R. More precisely, V is given by all extended maps Γ = i Γ i ⊗ |i) for which F(Γ) < F(Λ(t p )). When the coordinate system is chosen such that the asymptotic extended map Λ(t → ∞) lies in its origin, then F(Λ(t)) can be considered a distance between Λ(t) and the asymptotic map around which the object V is centred.
According to (i), the function F(Λ(t)) decreases monotonically under the dynamics induced by (1), for why F(Λ(t)) can never exceed the value F(Λ(t p )) for t > t p . The geometric implication of (i) is that Λ(t) will never leave V. Condition (ii) and (iii), in turn, manifest a relation between P and V such that P contains V. This implies that and all extended maps inside V are valid.
The initial extended map Λ(0) (indicated by a black cross in figure 1) lies on the boundary of P (due to χ(0) having rank one) for why V and P become tangential in the point Λ(0) when t p becomes infinitesimally small.

Simplification of condition (ii') and (iii)
The construction of R is typically most challenging when CP dynamics shall be proven for t ≥ 0 (i.e. when V and P become tangential in Λ(0)). A particularly convenient case is, however, given when the S is of rank one, which, in geometrical terms, means that the valid region P is confined by two parallel planes. The construction of R = R † that proves complete positivity for t ≥ 0 can in this case be simplified as the condition (iii) can be partially replaced by a set of equality conditions, which effectively decrease the number of degrees of freedom in R that have to be determined. This is particularly helpful for an analytical construction of the latter.
More specifically, the inequality R − S ≥ 0 is satisfied (i.e. (iii) holds true) if R is positive semi-definite, normalized as in (ii'), and satisfies for a basis {K (i) } of the kernel of S. Equation (15) is a set of linear constraints that reduces the number of degrees of freedom in R.
To understand why this reformulation of (iii) is possible, we emphasize that Λ(0) lies on the boundary of P and can thus not be from the kernel of the operator S. That means that Λ(0) and {K (i) } are linearly independent and form a complete set (because S is of rank one) such that every extended map v can be decomposed as v = aΛ(0)+bK with K being from the span of {K (i) } and a, b ∈ C. Without loss of generality, it can be assumed ‡ that a * b ∈ R. With this decomposition, we obtain (16) where (ii') has been employed. From (15), one can, however, deduce which, together with R ≥ 0, implies non-negativity of (16) and thus positive semidefiniteness of R − S.
With the framework derived hitherto, we are now in the position to analyse exemplary equations of motion that permit to construct R by purely analytic means. Subsequently, in section 3 we will expand the framework of constructing R towards numerical techniques that will permit to treat more complex cases as exemplified in section 3.3 and 4.3

Examples
2.4.1. The damped Jaynes-Cummings model An instructive example of a HEOM is given by the damped resonant Jaynes-Cummings model, which describes the dynamics of a two-level system inside a leaky cavity. The electro-magnetic cavity field is characterized by a Lorentzian spectral density function that is centred around the transition frequency of the two-level system. The spectral width of the cavity field is denoted by ζ, whereas γ labels the coupling strength of the two-level system and the cavity field. After tracing out the electromagnetic field, this model permits an exact description of the (generally non-Markovian) dynamics in terms of a HEOM, which reads The system state is denoted by 1 (t), whereas the two auxiliary operators 2 (t) and 3 (t) encode information about the environment and system-environment correlation. The introduction of a general framework by means of which this HEOM has be obtained is given in Appendix A.1 and remarks that explicitly concern the construction of (18a) -(18c) are found in Appendix A.1.3.
To turn this HEOM equation into an equation of motion for the dynamical maps Λ 1 (t) to Λ 3 (t), one writes (18a)-(18c) in the Bloch basis. The Bloch representation of the extended state is then given by a 12-dimensional vector, and the 12×12-dimensional Bloch representation for the generator L, is also the generator for the dynamics of the 12 × 4-dimensional Bloch representation of the extended dynamical map. The latter is initialized by Λ 1 (0) = I and Λ i (0) = O for i ≥ 2, where I and O are the identity and the null map, respectively. This corresponds to the initial condition 2 (0) = 3 (0) = O. It has been algebraically verified that the coefficient matrix χ(t), which characterizes the system map Λ 1 (see (5) with µ i = σ i ), is always of the form where the only time-dependent degree of freedom λ 1 (t) is a coordinate of the system map Λ 1 (t). Thus, the condition χ(t) ≥ 0, i.e. complete positivity of the system map Λ 1 (t), does only depend on λ 1 (t) and all other degrees of freedom of Λ 1 (t) can safely be ignored.
is a two-dimensional real-valued vector and The dynamical variable λ 1 (t) is a coordinate of Λ 1 (t) and the initial condition Λ(0) corresponds to λ(0) = [1,0]. Strict positivity of all non-trivial eigenvalues of χ(t) (two of them vanish constantly) in the first non-vanishing time-step is given when γζ > 0. When this product is positive, we can replace the normalization condition (ii) by (ii'). Due to the structure of χ(t) in (19), its highest and the second-highest elementary symmetric polynomials e 4 (χ(t)) and e 3 (χ(t)) are constantly equal to zero. The complete positivity of the system map Λ 1 (t) thus depends on the highest non-vanishing polynomial e h (χ) = e 2 (χ) = (tr(χ) 2 − tr(χ 2 ))/2 = (1 − λ 4 1 )/4 which is positive if and only if |λ To prove strict positivity of e h (χ(t)) for all times t > 0, we strive for a matrix r (the representation of R in the space of λ(t)) that satisfies l † r + r l ≤ 0 (condition (i)), λ † (0) (r −s) λ(0) = 0 (condition (ii')), and r −s ≥ 0 (condition (iii)). Condition (ii') and (iii) hold true, when r is parametrized by r = diag(1, r 22 )/4 with r 22 ≥ 0. Condition (i), in turn, is satisfied if q ≡ l † r + r l ≤ 0, which is the case if and only if trace and determinant of −q are non-negative. The latter is given as det(−q) = −(γr 22 − 2ζ) 2 /4 for why we have to choose r 22 = 2ζ/γ and with tr(−q) = 2r 22 ζ = 4ζ 2 /γ the conditions q ≤ 0 and r 22 ≥ 0 are equivalent to γ ≥ 0 and ζ ≥ 0. Together with γζ > 0 this proves that positivity of both parameters γ and λ is sufficient for e h (χ(t)) > 0 and implies CP dynamics. Even though eigenvalues can become very small (see figure 2(b) for an example), the conditions clearly assert that no eigenvalue of χ(t) vanishes.

2.4.2.
Reviving coherences: a HEOM with two levels Another model for a nonmonotonic decay of phase coherence can be formulated in terms of the hierarchy equation (see Appendix A.1.1 for an explicit derivation) with the Pauli operator σ z and the dephasing Lindbladian D z = σ z σ z − . For α = 0 and the initial condition 2 (0) = 0, this hierarchy equation describes purely Markovian dynamics with exponentially decaying coherence. For non-vanishing values of α, however, the dynamics deviates from a mere exponential decay and can give rise to genuine non-Markovian revivals. The initial condition 2 (0) governs the infinitesimal behaviour of 1 (t) at time t = 0. A generic choice is 2 (0) = 0 although other choices are possible as we will demonstrate later.
In the Bloch basis, the extended state can be expressed in terms of an eightdimensional vector and the according representation of the generator L reads  With the initial conditions Λ 1 (0) = I, and Λ 2 (0) = O and the structure of L one can show by purely algebraic means that the Bloch representations for the extended dynamical map and depends only on the two scalar parameters λ 1 (t) and λ 2 (t). With this form of the dynamical map Λ 1 (t), the corresponding matrix χ(t) is given by [1,0]. Strict positivity of all non-trivial eigenvalues of χ(t) for infinitesimal times is given when either γ 1 > 0, or γ 1 = 0 and ω, α > 0 hold true. The procedure of proving complete positivity is now analogous to the previous case: As λ 1 (t) and λ 2 (t) are the only relevant dynamical variables, the equation of motion for the extended dynamical map can be reduced to ∂ t λ(t) = l λ(t) with λ(t) = [λ 1 (t), λ 2 (t)], the initial condition λ(0) = [1, 0] and the generator Figure 3. Time evolution of the coherence tr(σ x 1 ) for an initial state 1 (0) = |+ +| with |+) = (|0) + |1))/ √ 2, evolving according to (21a), (21b), and its extension (25). The dynamics of (21a) and (21b) with initial conditions 2 (0) = 0 and˙ 1 (0) = 0 is depicted in blue (solid) and red (long dashing), respectively (γ 1 = γ 2 = 1/2 ω, α = 4). The orange (short dashing) line represents the dynamics induced by introducing the additional level in (25) As already seen for the Jaynes-Cummings model, the operator R is characterized by a 2 × 2-dimensional matrix r that is parametrized by r = diag[1, r 22 ]/4 with r 22 ≥ 0 in order to satisfy the conditions (ii') and (iii). The monotonicity condition (i) holds true if l † r + r l ≤ 0, and a good choice for r 22 (the choice is not unique) is given by r 22 = (αω 2 + γ 1 γ 2 )/(2α 2 ω 2 ). In this case, monotonicity and r 22 ≥ 0 are equivalent to γ 1 ≥ 0, γ 2 ≥ 0 and 2α + γ 1 γ 2 ≥ 0 such that these three conditions together with 2 (0) = O and strict positivity of χ(t) for short times are sufficient for the hierarchy in (21a) and (21b) to induce CP dynamics. An example of a process that satisfies these conditions is given in figure 3.
It is well established that the initially linear decay of coherence, that is depicted in the inset of figure 3, is an artefact of the Markov approximation. In the framework of hierarchy equations, it is straight-forward to remove this artefact by choosing the initial condition 2 (0) = −L 11 1 (0), which corresponds to λ(0) = [1, γ 1 /ω] and yieldṡ Figure 3 shows the difference of the time evolution in the first time-steps.
The examination of complete positivity requires a modification of the matrix r in order to still comply with (ii') due to the change of the initial condition λ(0). To this end, we parametrize r by means of the Cholesky decomposition such that it is positive semi-definite and λ(0) (r − s) λ(0) = 0, which (due to the symmetry of r) reduces the number of degrees of freedom in r to two. The condition r − s ≥ 0 is then incorporated into the parametrization of r as shown in section 2.3, which eliminates another degree of freedom. With this parametrization, the two conditions (ii') and (iii) are satisfied and all but one single degree of freedom of r have been determined. This last degree of freedom can for example be chosen such that det(l † r + r l) = 0. In this case one obtains r = 1 4η with η = 2αω 2 + γ 1 γ 2 . The monotonicity condition (i) and r ≥ 0 are then equivalent to γ 1 + γ 2 ≥ 0 and 2αω 2 + γ 1 γ 2 > 0, which (together with the condition that all nontrivial eigenvalues of χ(t) become strictly positive for short times t) are thus sufficient conditions for CP dynamics.
Such HEOM gives rise to a broader variety of dynamical processes and rekindles the question of complete positivity. In favour of a concise discussion, we reduce the number of free parameters by setting γ 1 = γ 2 = γ 3 ≡ γ.
The extended state (t) as well as the generator L are again expressed in the Bloch basis. In similarity to (22), this permits to express the HEOM ∂ t (t) = L (t) in terms of a vector equation and the transition from states to dynamical maps is made by replacing (t) by the extended dynamical map Λ(t), which is initialized by Λ 1 (0) = I and Λ i (0) = O for i = 2, 3. In the Bloch representation, all dynamical maps Λ i (t) preserve the structure diag[δ i1 , λ i (t), λ i (t), δ i1 ] such that the dynamics of Λ(t) is completely characterized by the three-dimensional real-valued vector λ(t), which is initialized by λ(0) = [1, 0, 0] and evolves according to ∂ t λ(t) = l λ(t) with In fact, the question of complete positivity only depends on the two ratiosα = αω 2 /γ 2 andβ = βω 2 /γ 2 , as one can see in terms of the scaled variablesλ i (t) = λ i (t)ω i /γ i , whose equations of motion is described byl that is obtained from (26) by multiplying all elements l ij with (ω/γ) i−j . Since the structure of χ(t) does not change as compared to the HEOM with two levels, its highest non-trivial elementary symmetric polynomial e h = e 2 is still correctly expressed by  with a, b ≥ 0. Considering the more interesting case ofα,β = 0, a very convenient choice for the three free parameters is given by a = 1/|2α|, b = 1/|2αβ| and c = 0. This matrix satisfies the monotonicity condition l † r + rl ≤ 0 if and only if one of the following two conditions is satisfied: Consequently, any of the two conditions (28a) or (28b) ensures complete positivity when strict positivity of all non-trivial eigenvalues of χ(t) is given for short times.
Although this is only one instance of how the matrix r can be chosen and other choices (e.g. with c = 0) are possible, it can prove complete positivity for a rather large parameter regime as shown in figure 4. An individual numeric optimization over r for each numerical instance of the HEOM can even slightly enlarge that parameter set.
3. Numeric approaches to more complex problems

Transformation of the extended dynamical map
A central element of the procedure discussed in section 2 is to express the highest non-trivial elementary symmetric polynomial e h (χ(t)) in terms of a quadratic function in Λ(t) as shown in (7). For the examples discussed in section 2.4, such a representation has been obtained by anticipating the structure of the matrix χ(t) that characterizes potential solution to the initial value problem by purely algebraic means. Because the initial condition Λ 1 (0) = I and Λ i (0) = O for i ≥ 2 is fixed, it is for example often possible to prove that particular subspaces are never occupied. To this end, one shows that the dimension of the space that is spanned by L k Λ(0) with k = 0, . . . , d and d being the number of degrees of freedom in Λ(t) is strictly smaller than d. When that is the case, then this indicates the existence of preserved quantities, which permit predictions about the form of χ(t) and can thus help to obtain (7) (and to reduce the dimensionality of Λ(t)).
To make the approach more universal and to permit e h (χ(t)) to be of a degree m > 2 in χ(t), one can apply a transformation which renders geometries that are always accessible to our technique. To this end, we define a function with a constant c > 0 such that the strict positivity of G[ is of a higher degree in Λ(t), considering the function G[χ(t)] rather than e h (χ(t)) has the big advantage that the prior can always be expressed in a form that is very similar to (7). This is done by defining a new vector Ξ(t) whose first component reads Ξ 1 (t) = e h (χ(t)) − c. Independently of the components Ξ i (t) with i ≥ 2, this vector permits to express G[χ(t)] as with a matrixS that is of the simple formS = diag[1, 0, . . . , 0]. It is essential that Ξ(t) follows a linear differential equation ∂ t Ξ(t) =L Ξ(t) and as the elements of Ξ(t) are functions of χ(t) (and thus Λ(t)), this differential equation needs to be derived from the HEOM. Already the first component Ξ 1 (t), however, is a non-linear function in Λ(t) and its time-derivative does typically also involve terms that are non-linear in Λ(t). For Ξ(t) to anyway evolve according to a linear equation, one can now consider any non-linear contribution to be a new element of Ξ(t) and keep increasing the size of Ξ(t) until one has recovered a linear equation of motion ∂ t Ξ(t) =L Ξ(t) §. The matrix L is constructed accordingly and eventually yields the desired equation of motion for Ξ(t).
In order to prove strict positivity of G[χ(t)], one follows the same procedure that has been introduced in section 2.1 and seeks a matrixR ≥ 0 that satisfies (i) to (iii) (where L and Λ(t) must be replaced byL and Ξ(t) and the scalar product becomes the § A simple example would be the scalar equation of motionẋ(t) = αx(t) and the non-linear coordinate transformation y(t) = x 2 (t) + x(t). The coordinate y(t) satisfiesẏ(t) = 2αy(t) − αx(t); replacing x(t) in terms of y(t) would result in a non-linear equation of motion, but treating x(t) as independent variable, permits to maintain linear equation of motion at the expense of increasing the number of independent variables. standard scalar product in the Euclidean space). To this end, it is beneficial to choose c = lim t→∞ e h (χ(t)) such that the second term in (30) vanishes asymptotically. If such a matrixR can be found and G[χ(t)] increases in first order, then G[χ(t)] > 0 hold true for t > 0 and the system dynamics is completely positive.

Semi-definite programming
The transformation introduced in section 3.1, or a consideration of more complex quantum systems and HEOMs with several levels, can quickly increase the dimensionality of the HEOM to a level at which the operator R can not be constructed analytically. In such cases, however, the convexity of the problem permits the application of a highly efficient numerical approach. To this end, let us parametrize R such that the normalization condition (ii) (or (ii')) is satisfied and define Q ≡ L † R+R L and B is a semi-definite program [29], which can be solved reliably and efficiently. The constraint B ≥ 0 implies R − S ≥ 0 (i.e. condition (iii)) and Q ≤ vI. If the minimum v m is non-positive, one has found R such that (i) holds true. The condition v m ≤ 0 thus verifies complete positivity.

Example: non-Markovian dynamics due to a finite temperature bath
To gain a better intuition of how the transformation introduced in section 3.1 affects the geometric properties of the problem, we consider a two-level system that is coupled to a bath of finite temperature. In the Markovian case, the situation is described by a master equation in Lindblad form˙ with rates γ ± . In order to also permit non-Markovian dynamics, we introduce the traceless auxiliary operator 1 (t) and extend the Lindblad equation to the HEOM (see Appendix A.2 for a motivation of this equation) with γ p = (γ + + γ − )/2, γ m = (γ + − γ − )/2, and ξ and ω being scalar parameters that characterize the additional levels in HEOM (the special cases ω = 0 and ξ = 0 give rise to Markovian dynamics). With = 1 ⊗ |1) + 2 ⊗ |2) and the operator L that is derived from the HEOM (33a)-(33b), the latter can be expressed as ∂ t (t) = L (t). The equation of motion for the dynamical maps is obtained by replacing (t) with The structure of L and the initial condition give rise to a set of identities and restrictions which finally permit to fully represent Λ(t) by means of a four-dimensional vector λ(t), whose initial condition is λ(0) = [−γ m /γ p , 0, 1, 0], and which evolves according to ∂ t λ(t) = l λ(t) with The dynamical map Λ 1 (t) in the Pauli basis is given as in (5) (with µ i = σ i ) and due to the structural properties of L and Λ(0), the coefficient matrix is always of the form where θ(t) ≡ γ m + γ p λ 1 (t). In contrast to the previous examples, this structure of χ(t) does not permit conclusions about constantly vanishing eigenvalues such that the highest possible elementary symmetric polynomial e h = det must be considered. Equation (35) does, however, give rise to the factorisation e h (χ(t)) = det(χ(t)) = αP 2 1 (t) P 2 (t) with a proportionality factor α = γ + γ − /[4(γ 2 − − γ 2 + )] 4 and P 1 (t) = γ p λ 1 (t) + γ m (36a) Here, the constants are given by p . Considering the case γ + , γ − > 0, the proportionality factor α is positive such that it is sufficient to show strict positivity of P 2 1 (t) and P 2 (t) separately in order to prove the positivity of e h (χ(t)). After verifying strict positivity of χ(t) for short times, this ensures CP dynamics for t ≥ 0.
The factor P 2 1 (t) is always non-negative but to show strict positivity for t > 0, we reformulate the problem such that λ † (0) s (1) λ(0) − λ † (t) s (1) λ(t) > 0 with a symmetric matrix s (1) is sufficient for P 2 1 ( λ(t)) > 0 and strive for the construction of a matrix r (1) that satisfies l † r (1) + r (1) l ≤ 0, λ † (0) (r (1) − s (1) ) λ(0) = 0, and r (1) − s (1) ≥ 0 (see (i), (ii'), and (iii)). The procedure is very similar to the previous examples for why we abstain from a detailed discussion and anticipate that P 1 (t) never returns to zero when −2γ 2 p ≤ ωξ holds true. Showing the second condition P 2 ( λ(t)) = λ † (0) s (2) λ(0) − λ † (t) s (2) λ(t) > 0 is more involved, which can be understood best in geometric terms. A cross-section (for λ 2 = λ 4 = 0) of the set P 2 of vectors λ that satisfy P 2 ( λ) > 0 is depicted in figure 3.3(a) and the black cross corresponds to the initial condition λ(0). As one can see, it is not Figure 5. The sets inside which P 2 > 0 is satisfied are depicted in blue before (a) and after (b) the non-linear transformation from λ(t) to Ξ(t). The black lines show the trajectories of λ(t) and Ξ(t) for the initial conditions that correspond to Λ 1 (0) = I (black crosses). For the geometry depicted in (a), it is impossible to find a matrix r (2) such that the equipotential line of points ξ with ξ † r (2) ξ = const contains λ(0) and lies completely inside P 2 (the set with P 2 > 0). After the non-linear transformation a more benign geometry has been obtained such that an according matrix can be found.
possible to place an ellipsoidal equipotential line λ † r (2) λ = const (that would confine the set V 2 ) inside P 2 such that it contains the initial condition λ(0). In algebraic terms, the two conditions (ii') and (iii) can not be satisfied simultaneously for geometric reasons.
Due to the rather high dimensionality of Ξ(t), such a function can not be found analytically but by means of a semi-definite program a very efficient numeric construction has been carried out.

Quantifying the violation of complete positivity
In many cases, an exact description of the dynamics of a physical system is only obtained for infinitely deep hierarchy equations. A truncation of the HEOM is then often employed in order to obtain approximate system dynamics. Such a truncation, however, can give rise to a slight violation of complete positivity of the induced dynamical system map Λ 1 (t). The latter typically occurs for short times, at which the coefficient matrix χ(t) has very small eigenvalues, and vanishes for all times greater than a critical time t p . During this initial time interval 0 ≤ t ≤ t p complete positivity is violated and even if the violation is numerically negligible, the proposed framework will not be able to verify CP dynamics.
In such cases, it is desirable to estimate the degree of violation of complete positivity. This degree is for instance an indicator of whether the approximation is good enough as strong violation of complete positivity suggests unphysical dynamics and thus an insufficient approximation. We will propose two different ways of quantifying this violation and demonstrate their applicability by means of an important example.

Complete positivity after an initial violation time
At the initial time t = 0, when the system map Λ 1 (0) is equal to the identity, the coefficient matrix χ(0) has just one non-vanishing eigenvalue. A truncation (and thus approximation) of the HEOM therefore often causes an eigenvalue of χ(t) to become negative for short times. The identification of a critical time t p at which the dynamics becomes CP again helps to estimate whether the violation is indeed just a short-time phenomenon, or a structural problem that eventually impedes a physical interpretation of the time-evolution. A good guess for t p can be obtained from a numeric propagation of χ(t) for short times and once a candidate for t p (and the according extended map Λ(t p )) has been identified, the proposed method in section 2.1 for proving complete positivity can be employed to show that the dynamics is CP for all t > t p .
The situation is depicted in figure 6(a), where it is schematically visualized how the extended dynamical map leaves the valid region P for short times in order to re-enter at time t p > 0.

Lower bound on the eigenvalues of χ
An alternative way of quantifying the violation of complete positivity of the system map Λ 1 (t) is based on estimating a lower bound −δ min on the eigenvalues of the matrix χ(t). Rather than requiring the positivity of the highest non-trivial elementary symmetric polynomial e h (χ(t)), we strive to prove e h (χ(t) + δ min I) > 0, which permits negative eigenvalues of χ(t) as long as they remain greater than the new parameter −δ min .
Algebraically this means that one needs to find an extended map Ω = i Ω i ⊗ |i) that is such that (7) is modified to  The valid region P of extended maps that induce CP dynamics is depicted in blue and Λ(t) leaves this region for short times. In (a), the initial condition is thus replaced by Λ(t p ) (red cross). The time interval 0 ≤ t ≤ t p is disregarded and the proposed procedure is employed to determine the set V, that contains Λ(t) for t > t p , such that V lies inside P and complete positivity is shown for t > t p . In (b), the condition χ(t) > 0, that defines P, is replaced by χ(t) > −δ min I which enlarges the set P of states that are considered valid. The tangential point Ω of P and V becomes object to an optimization, which increases the chance to find V such that V ⊆ P. In that case, −δ min is a lower bound for the eigenvalues of χ(t).
Whenever R is found such that (i), (ii"), (iii), and (iv) are satisfied, this proves strict positivity of e h (χ(t + δ min )) for all times t ≥ 0. To prove this statement, (ii") is applied to (37), which yields By means of (iii), this can be further bound to which, according to (iv), is strictly positive at the time t = 0. For all times t > 0, the right-hand side of (39) increases monotonically due to (i), such that e h (χ(t) + δ min I) > 0 and all eigenvalues of χ(t) are greater than −δ min for all times.
In geometric terms, this modification enlarges the region P that is considered valid as depicted in figure 6(b). The points Γ with F(Γ) ≤ F(Ω) (where F is defined as in section 2.2) define the region V that will contain the solution Λ(t) for all times and Ω is the tangential point of V and the (enlarged) set P. It is important to point out that Ω is not uniquely defined and becomes itself object to an optimization. This often permits to find an operator R that satisfies (i), (ii"), (iii), and (iv) and proves −δ min to be a lower bound on the eigenvalues of χ(t).

Example: The spin-Boson model
One of the most prominent examples of a microscopically derived HEOM describes the dynamics of a two-level system that is coupled to a bath of harmonic oscillators. For a high temperature case (βγ 1) with weak system-bath coupling (∆ 2 /γ 2 ω), the HEOM is expressed as [10,12,13] Here, 1 (t) represents the state of a system with a resonance frequency ω, whereas [·, ·] and {·, ·} denote the commutator and the anti-commutator, respectively. As discussed in [10,12,13], ∆ characterizes the system-environment coupling, β corresponds to the inverse bath temperature, and γ is a damping constant. Again (40a)/(40b) can be written as ∂ t (t) = L (t), where = 1 ⊗ |1) + 2 ⊗ |2) and the generator L has been deduced form the HEOM. In order to obtain the equation of motion for the extended dynamical map Λ(t), one replaces (t) by Λ(t) = Λ 1 ⊗ |1) + Λ 2 ⊗ |2). Expanding the eigenvalues i (t) of the matrix χ(t) that characterizes the dynamical map Λ 1 (t) for the system state in time, one finds that one of these eigenvalues is given by 1 (χ(t)) = − 1 144 (β 2 γ 2 ∆ 2 ω 2 ) t 4 + O(t 5 ). For non-vanishing values of β, γ, ∆ and ω, this eigenvalue becomes negative for short times, which results in a violation of complete positivity of Λ 1 (t). Sufficient conditions for complete positivity can hence not be found for this important example. It is, however, possible to estimate the degree of violation of complete positivity: A numerical integration yields that the violation is typically not very strong and limited to rather short times, both of which we prove in the following. 4.3.1. Complete positivity for t > t p By propagating the extended dynamical map over a short period of time, it is possible to determine a time t p at which the coefficient matrix χ(t p ) is positive definite. Based on the procedure laid out in section 4.1, we want to prove that complete positivity is given for all subsequent times.
Due to the structure of L and Λ(0), it is possible to represent the extended dynamical map Λ(t) by six degrees of freedom that form the six-dimensional real-valued vector λ(t), which is initialized by λ(0) = [1, 0, 1, 0, 0, 1]. The latter evolves according to the differential equation ∂ t λ(t) = l λ(t) with The simplification of Λ(t) yields a set of identities for the dynamical map Λ 1 (t), which gives rise to in a certain structure of χ(t) given by where all matrix elements are linear functions in λ(t). The structure of χ(t) does not permit conclusions concerning the existence of constantly vanishing eigenvalues such that the determinant is considered to be the highest non-trivial elementary symmetric polynomial. It is evaluated to e h (χ(t)) = det(χ(t)) ∝ P 1 (t)P 2 (t) with a positive proportionality factor and where b = −β 2 ω 2 + 4. For the considered initial condition λ(0), both factors vanish initially and by algebraic means it can be shown that P 1 (t) = P 2 (t) holds indeed true for all times t ≥ 0. It is thus sufficient to keep track of P 1 (t) only, and to show that it never returns to zero, in order to prove e h (χ(t)) > 0 and thus complete positivity of Λ 1 (t) for t ≥ t p . For geometrical reasons that are similar to those discussed in section 3.3, the transformation that has been introduced in section 3.1 is required in order to render a situation in which our procedure is applicable. To this end, a vector Ξ(t) is constructed based on λ(t) that is such that G = Ξ † (0)S Ξ(0) − Ξ † (t)S Ξ(t) > 0 (with the matrix S = diag[1, 0, . . . , 0]) implies P 1 (t) > 0 for t > t p . The strict positivity of G(t) is proven when a matrixR has been found that satisfiesL †R +RL ≤ 0, Ξ † (t p ) (R −S) Ξ(t p ) = 0, andR −S ≥ 0, which can be carried out by means of a semi-definite program.
As an example, we consider the case γ = 3 ω, ∆ = 2 ω and β = 8/10 ω −1 depicted in figure 7. This case features a clear signature of non-Markovianity, which is measured [30] to N (Λ 1 (t → ∞)) = 0.34. Complete positivity is proven after the critical time t p 0.6/ω. It is noteworthy that the non-Markovian revivals arise after the critical time t p .

Lower bound on the eigenvalues of χ(t)
The second possibility to quantify the violation of complete positivity is to estimate a lower bound −δ min on the eigenvalues of the matrix χ(t) as discussed in section 4.2. To this end, we make a guess for δ min > 0 and prove e h (χ(t) + δ min I) = det(χ(t) + δ min I) > 0 for t ≥ 0. The time-evolution of Λ(t) is not affected by this modification of e h for why the generator l given in (41) and the structure of the matrix χ(t) in (42) still apply.
As mentioned before, the vectors Ω (i) are not uniquely defined and become objects to an optimization. Good candidates are determined with a gradient-based method and by means of a semi-definite program we show that −δ min = −10 −2 is a lower bound on the smallest eigenvalue of χ(t) for the test case of γ = 1 ω, ∆ = 2/10 ω and β = 2/10 ω −1 . This HEOM induces a dynamical process with a non-Markovianity of N (Λ 1 (t → ∞)) = 0.08, whose time evolution is depicted in figure 7.

Conclusion
We have demonstrated how complete positivity can be incorporated into the framework of HEOMs. This enables an abstraction of the equations of motion from a derivation based on microscopic models towards a more phenomenological perspective that opens a new approach to non-Markovianity. A construction of elementary models of open system dynamics in which physical quantities like amplitudes and frequencies of revivals of quantum coherence have a clear root in the underlying equations of motion -just like decay constants are rooted in Markovian master equations -will ultimately ease investigations of open system dynamics in which system size or large statistical sampling [31,32,33,34] requires efficient phenomenological models. Although a theorem that describes CP dynamics as beautifully as the Lindblad theorem in the Markovian case is currently far out of reach for non-Markovian systems, the presented approach might be taken as an initial step towards such a theory.
(t) = n i=1 i (t) ⊗ |i) that solves the HEOM for the initial condition 1 (0) =ˆ and The HEOMs that we target for are of triangular form, i.e. their generators satisfy L ij = ωI for j = i + 1 and with 1/ω being a unit of time and I and O being the identity map and the null map, respectively. With this form of L, the level k of the HEOM can affect the system dynamics earliest in the k th time-step.
We will now successively deduce the yet unknown operators L ij , with j ≤ i by imposing the condition that each operator k (t) must vanish in the first k − 2 timesteps: Whereas the system state 1 (t) does not have to vanish at all, the first auxiliary operator 2 (t) shall vanish at the time t = 0, the second auxiliary operator 3 (t) shall satisfy 3 (0) = O as well as ∂ t 3 (t)| t=0 = O, and this series of conditions will be continued for all k (t).
Let us demonstrate the procedure more explicitly: After having defined a dynamical map Λ 1 (t) for the system state, we definẽ as a function of time, the initial system stateˆ , and the yet unknown operator L 11 . We then require˜ 2 (L 11 ,ˆ , t) to vanish at the initial time t = 0 independently of the initial system stateˆ , which can be expressed by ∂ ∂ˆ ij˜ 2 (L 11 ,ˆ , 0) = 0 ∀i, j .
Solving this linear equation with respect to L 11 determines the first level of the HEOM as well as the first auxiliary operator, which can now be obtained froṁ 1 (t) = L 11 1 (t) + ω 2 (t) and reads In order to obtain the operators L k1 to L kk for each higher level k of the HEOM, we iteratively reapply this procedure. That is, we first define a functioñ as a function of yet unknown operators L k1 to L kk . According to the initial conditions, the first k − 2 time derivatives of˜ k+1 (L k1 , . . . , L kk , t) must vanish at the initial time is required for all v = 0 . . . , k − 2. After the operators L k1 to L kk have been determined such that they satisfy all these conditions, the auxiliary operator vanishes in the first k time-steps as it has been required above. As this procedure is iteratively re-applied, the number of levels of the HEOM increases until (A.7) yields n+1 (t) = 0. This indicates that one has obtained a HEOM with n levels, for which the initial condition j (0) = O for j > 1 gives rise to the targeted system dynamics 1 (t) = Λ 1 (t) 1 (0). Appendix A.1.1. Example: the derivation of the HEOM (21a)/(21b) The method that has been introduced above is now applied in order to motivate how the HEOM in (21a)/(21b) has been obtained. To this end, we aim for a HEOM that induces the system dynamics with f (t) = e −γt cos(ωt) being the function that characterises the reviving coherences. The parametersˆ (i) = tr(ˆ σ i ) with i = 0, x, y, z and σ 0 = I determine the initial statê and need to satisfyˆ (0) = 1 and [ˆ (x) ] 2 + [ˆ (y) ] 2 + [ˆ (z) ] 2 ≤ 1 in order forˆ to have a trace equal to one and to be positive semi-definite.
To find a HEOM that induces the targeted dynamics, we parametrise the yet unknown operator L 11 by where the summation is carried out over i, j = 0, x, y, z. In accordance with (A.2), one defines˜ 2 (L 11 , t) ≡ (˙ 1 (t) − L 11 1 (t))/ω such that (A.3) reads ∂ ∂ˆ (k)˜ 2 (L 11 , 0) = 0 with k = 0, x, y, z (A. 10) and is solved with respect to χ L 11 ij . This yields χ L 11 = diag(−γ, 0, 0, γ)/2 or, equivalently, L 11 = γ D z /2 with the dephasing Lindbladian D z ≡ σ z σ z − . Equation (A.7) determines the first auxiliary operator 2 (t) as with k = 2 and the scalar function g 2 (t) = −e −γt sin(ωt)/2. We can now re-applying the same procedure to the operator 2 (t), in order to determine L 21 and L 22 . First, these operators are parametrised in terms of yet unknown scalars χ L 21 where k = 0, x, y, z, with respect to parameters χ L 21 ij and χ L 22 ij . The conditions in (A.12) determine only 24 out of the 32 degrees of freedom in χ L 21 ij and χ L 22 ij and the remaining degrees of freedom can be chosen such that the forms of L 21 and L 22 become as simple as possible, what eventually yields L 21 = ω D z /2 and L 22 = γ σ z σ z . The second auxiliary operator 3 (t) ≡ (˙ 2 (t) − L 21 1 (t) − L 22 2 (t))/ω vanishes constantly such that the HEOM has only two level.
We have now derived a HEOM that induces the system dynamics defined in (A.8). However, we want to generalise this HEOM by permitting an arbitrary pre-factor for L 21 such that the latter reads L 21 = αω D z with a free parameter α. Furthermore, we permit the damping constant γ to take different values γ 1 and γ 2 on the two different levels of the HEOM. For α = 1/2 and γ 1 = γ 2 = γ one obtains the original HEOM but the additional parameters also permit to cover different processes. The final HEOM that has now been obtained is stated in (21a)/(21b).
Appendix A.1.2. Example: extension of the HEOM (21a)/(21b) by means of (25) We can further generalise the dynamics that has been defined in (A.8) by modifying the function f (t) to f (t) = e −γt (1 + 2α[cos(ωt) − 1]). In the case of α = 1/2 one recovers the same dynamics that has been targeted for in Appendix A.1.1. However, different values of α permit a broader range of dynamical processes and call for a new HEOM.
In order to obtain the latter, we follow the proposed procedure in the very same way that has been carried out in Appendix A.1.1 and find that the first operator in the HEOM reads L 11 = γ D z /2. The auxiliary operator 2 (t) is then determined as in (A.7) with k = 1 and is given through (A.11) with k = 2 and g 2 (t) = −e −γt α sin(ωt). As in (A.5) we then define˜ 3 (L 21 , L 22 , t) all conditions on which (compare (A.12)) are satisfied by the operators (the choice is not unique) L 21 = αω D z and L 22 = γ σ z σ z . In contrast to the previous example, the operator 3 (t) does in the present case not vanish constantly but is given through (A.11) with k = 3 and g 3 (t) = e −γt α(1 − 2α)(1 − cos(ωt)). For that reason, an additional level in the HEOM is required, which is characterised by the yet unknown operators L 3k = ij χ L 3k ij σ i † j with k = 1, 2, 3. By imposing the conditions ∂ ∂ˆ (k)˜ 4 (0) = 0 , ∂ ∂ˆ (k) [∂ t˜ 4 (t)| t=0 ] = 0 , and ∂ ∂ˆ (k) ∂ 2 t˜ 4 (t)| t=0 = 0 (A. 13) with k = 0, x, y, z on˜ 4 (L 31 , L 32 , L 33 , t), which has been obtained from (A.5) with k = 3, one can eventually determine the operators L 31 = O, L 32 = ω(1 − 2α) σ z σ z , and L 33 = γ σ z σ z . Other choices for the operators are possible as (A.13) does not determine them uniquely, but the present choice has been found most convenient. Eventually, the HEOM shall again be generalised by replacing the operator L 32 with L 32 = βω σ z σ z . The first two levels of final HEOM are then identical to (21a) and (21b). The latter is, however, extended by an additional summand ω 3 (t), which evolves according to (25).
Appendix A.1.3. Example: the Jaynes-Cummings model in (18a) -(18c) In contrast to the previous two examples, which were not related to any particular physical system, we will now derive a HEOM that describes the dynamics of the resonant Jaynes-Cummings model [35]. To this end, we consider a two-level system inside a leaky cavity whose coupling to the cavity field is determined by the coupling constant γ. The field is characterised by a Lorentzian spectral density that is peaked at the resonance frequency of the two-level system and whose spectral width is denoted by ζ.
The reduced system dynamics has been obtained analytically and reads [35] 1 (t) = |f (t)| 2ˆ (11) f (t)ˆ (12) f * (t)[ˆ (12) ] * 1 − |f (t)| 2ˆ (11) , (A.14) where the parametersˆ (11) ∈ R andˆ (12) ∈ C determine the initial system state and with α = ζ 2 − 2γζ characterises the system dynamics. Knowing the analytical solution 1 (t), one can apply the procedure that has been described above in order to find a HEOM for the resonant Jaynes-Cummings model. The derivation is carried out in analogy to Appendix A.1.2 and eventually yields the HEOM (18a)-(18c). of 1 (t) would be completely encoded into the time-evolution of 2 (t). However, due to the summand ω 2 (t) in the equation of motion for 1 (t), the auxiliary operator 2 (t) couples back to the system state and affects the latter not at the time t but only in the next time-step t + δt. This contributes a kind of "inertia" to the system dynamics and gives rise to non-Markovian oscillation of quantum coherence.