Generalized theory of pseudomodes for exact descriptions of non-Markovian quantum processes

We develop an exact framework to describe the non-Markovian dynamics of an open quantum system interacting with an environment modeled by a generalized spectral density function. The approach relies on mapping the initial system onto an auxiliary configuration, comprising the original open system coupled to a small number of discrete modes, which in turn are each coupled to an independent Markovian reservoir. Based on the connection between the discrete modes and the poles of the spectral density function, we show how expanding the system using the discrete modes allows for the full inclusion non-Markovian effects within an enlarged open system, whose dynamics is governed by an exact Lindblad master equation. Initially we apply this result to obtain a generalization of the pseudomode method [B. M. Garraway, Phys. Rev. A 55, 2290 (1997)] in cases where the spectral density function has a Lorentzian structure. For many other types of spectral density function, we extend our proof to show that an open system dynamics may be modeled physically using discrete modes which admit a non-Hermitian coupling to the system, and for such cases determine the equivalent master equation to no longer be of Lindblad form. For applications involving two discrete modes, we demonstrate how to convert between pathological and Lindblad forms of the master equation using the techniques of the pseudomode method.

We develop an exact framework to describe the non-Markovian dynamics of an open quantum system interacting with an environment modeled by a generalized spectral density function. The approach relies on mapping the initial system onto an auxiliary configuration, comprising the original open system coupled to a small number of discrete modes, which in turn are each coupled to an independent Markovian reservoir. Based on the connection between the discrete modes and the poles of the spectral density function, we show how expanding the system using the discrete modes allows for the full inclusion non-Markovian effects within an enlarged open system, whose dynamics is governed by an exact Lindblad master equation. Initially we apply this result to obtain a generalization of the pseudomode method [B. M. Garraway, Phys. Rev. A 55, 2290(1997] in cases where the spectral density function has a Lorentzian structure. For many other types of spectral density function, we extend our proof to show that an open system dynamics may be modeled physically using discrete modes which admit a non-Hermitian coupling to the system, and for such cases determine the equivalent master equation to no longer be of Lindblad form. For applications involving two discrete modes, we demonstrate how to convert between pathological and Lindblad forms of the master equation using the techniques of the pseudomode method.

I. INTRODUCTION
The theory of open quantum systems, which concerns the interaction between a quantum system of interest and a large macroscopic reservoir or heat bath [1,2], plays a fundamental role in several applications of quantum physics, ranging from quantum information [3], quantum technologies, and decoherence [4,5], through to quantum optics [6], condensed matter [7] and quantum thermodynamics [8]. In many such applications a standard approach is to model the effect of the environment in terms of a Markovian master equation, whose general validity requires the environmental noise (as measured through the reservoir correlation function) to be correlated over a much shorter time interval than the characteristic decay time of the open system. In particular, this condition is known to be well satisfied in quantum optical and mesoscopic systems [9] where the reservoir coupling constants vary slowly with frequency, and the coupling to the system of interest is typically weak.
For many other situations, however, especially those involving environments that are structured-i.e., with long correlation times and frequency dependent coupling constants-the assumption of a large separation of time scales between the system and environment no longer applies, and for these cases the resulting dynamics is non-Markovian [10][11][12]. In recent years there has been renewed attention paid to non-Markovian open quantum systems, not only out of fundamental interest, but * gpleasance1@gmail.com also due to the growing number of practical applications. From one perspective, a wide variety of experimental platforms, including atom-cavity and trapped ion systems [13], solid-state devices [14], and photonic band gap materials [15], have been shown to feature regimes where non-Markovian and strong coupling effects play an significant role in the description of the dynamics. At the same time, the increasing ability to coherently control the non-Markovian dynamics of quantum systems through, e.g., the use of reservoir engineering techniques [16][17][18][19][20][21], has provided new avenues to explore how certain types of environmental noise might be useful for the implementation of quantum technologies: notably, quantum information processing and quantum metrology have been recognized to possibly benefit from non-Markovian noise sources [22][23][24][25].
Within the validity of the Markov and weak coupling (Born) approximations, it is well-known that the quantum master equation describing the reduced dynamics can be generally expressed in Lindblad form [26]. Master equations of this type have long been the focus of both theoretical and experimental research, not only because of their ability to describe essential features of dissipation and decoherence, but also due to the existence of efficient numerical methods for their solution [27,28]. By contrast, while it is possible to derive a generalized form of master equation without the use of such approximations [29], the resulting non-Markovian equations of motion are often far too demanding to solve for an exact description to be feasible. For this reason, a certain class of methods for treating complex open system problems have been developed on the alternative idea of mapping the initial sys-tem onto a simpler, auxiliary configuration, consisting of the original open system coupled to small number of auxiliary bosonic (fermionic) modes, which in turn are coupled to an external Markovian reservoir [30][31][32][33][34][35][36][37][38][39][40][41] (see Fig.  1). In particular, the pseudomode method has effectively utilized such a mapping to describe the non-Markovian dynamics of a two-level system interacting with a bosonic environment [38,39]. In this approach the environment is replaced by a set of auxiliary discrete modes-the pseudomodes-which are identified through evaluating the poles of the spectral density function (i.e., the Fourier transform of the reservoir correlation function) when analytically continued to the lower-half complex frequency plane. By expanding the system over the pseudomodes one, can then derive a Lindblad master equation describing the dynamics induced by the non-Markovian interaction between the pseudomodes and two-level system, in addition to the coupling of the pseudomodes to an external Markovian reservoir. Importantly, while this method is exact, its application is currently restricted to regimes where only one excitation is initially present in the system, as well as to interactions valid within the rotating wave approximation.
Beyond this approach, we note in Ref. [40] that a similar type of mapping has been employed in conjunction with the Fano diagonalization technique to extend the treatment of the pseudomode method to multiple excitation regimes. The method is distinct from [39] in that it relies on 'undressing' the environment operators into a set of auxiliary quasi-modes, whose parameters-including the couplings to the system and overall configuration (i.e., site energies and inter-mode couplings)-are chosen so as to recover the spectral density function of the original environment. However, owing to the general difficulty of determining these parameters exactly, applications of the mapping so far have only focused on specific models where the system of interest is either coupled to a high-Q cavity or photonic band gap reservoir [41]. It is also worth noting that examples of other mappings related to the Fano diagonalization technique have been considered in Refs. [30,31,33,37]. With these approaches the environmental memory is instead accounted for by an ad hoc collection of discrete modes whose parameters are obtained as those which most accurately represent the spectral density function of the original model. Such mappings can be applied very generally, but unlike [39,40] have the disadvantage of not being grounded in exact relations between the physical and auxiliary environments, requiring their accuracy to often be validated against exact numerical methods.
Recently, a proof of an exact mapping of a non-Markovian open system onto a Markovian one for a system interacting with a Gaussian (bosonic) environment was given in [42] (see also Refs. [43,44]). There it was shown that the reduced dynamics of a non-Markovian system can be equivalently described in terms of an exact Lindblad-type master equation for a enlarged Markovian open system (system plus discrete modes), which in the context of the pseudomode description was used to generalize Ref. [39] beyond single excitation regimes for cases where the reservoir spectral density adopts a Lorentzian form. In this paper, we show that this treatment may be extended to physical environments modeled by a generalized spectral density function. In particular, our approach relies on the connection between the poles of the spectral density in the lower-half complex frequency plane and the properties of the discrete modes used represent the memory part of the environment. For certain spectral densities we find that the exact system dynamics may be obtained from a pathological master equation which lacks a suitable physical interpretation for the dynamics of the enlarged system. As we shall see, this arises due to the coherent part of the system-discrete mode evolution being generated by a non-Hermitian interaction Hamiltonian, which may be brought into an appropriate Hermitian form by applying a change of basis to the discrete modes. This paper is organized as follows. After outlining the physical model in Sec. II, we proceed Sec. III to introduce a mapping of the initial problem onto an auxiliary model and subsequently prove the reduced system dynamics to be indistinguishable between the two. In Sec. IV we then derive an exact form of Lindblad master equation for the enlarged Markovian system and present an initial application of this result. In Sec. V, we address for certain cases how to convert between pathological and Lindblad forms of the derived master equation. Finally, a summary and outlook is presented in Sec. VI.

II. DYNAMICAL MODEL
We start by considering a generic microscopic model of an open quantum system (OQS) S interacting with a bosonic environment E, as depicted in Fig. 1(a). The total Hamiltonian of the model is written as where Here, a λ (a † λ ) is the bosonic annihilation (creation) operator for an excitation of frequency ω λ satisfying the usual commutation relation [a λ , a † λ ] = δ λλ , c j (c † j ) is a generic OQS operator associated to the j transition of S involved in the coupling, and g jλ denotes the coupling strength between the ω λ mode of the field and the j transition of the OQS.
In what follows the system Hamiltonian H S is to be left unspecified and may in general have an explicit time dependence. On the other hand, the free evolution of the OQS (i.e., the evolution occurring in the absence of any driving or coupling between internal degrees of freedom) is described by the Hamiltonian with the set of discrete energy levels (eigenenergies) of the system denoted by {|e n j } S ( nj ). Based on this decomposition, if we assume each |e n j to be separated by a constant energy difference ω j = nj +1 − nj ( nj +1 > nj , ∀j, n), the transition (jump) operators c j , c † j appearing in Eq. (3) may be formally defined as such that c j (c † j ) lowers (raises) the internal energy of the OQS by an amount ω j . It is also noted from this same property that the OQS transition operators satisfy the eigenoperator relations [H S,0 , c j ] = −ω j c j ([H S,0 , c † j ] = ω j c j ). Hence, in an interaction picture generated by the unitary transformation the interaction Hamiltonian H I becomes with defining the environmental noise operators. Notice in particular that the absence of terms oscillating at frequencies ±i(ω λ + ω j ) in Eq. (7) implies the use of the rotating wave approximation (RWA). Following Ref. [42], we are interested in examining the time dependent behavior of the OQS in cases involving initially factorizing conditions ρ SE (0) = ρ S (0) ⊗ ρ E (0), where for simplicity the environment is taken at t = 0 to be in the vacuum state Because ρ E (0) is then Gaussian and satisfies Tr E [B j (t)ρ E (0)] = 0, the OQS dynamics described by the reduced density operator will only depend on the second-order moments of the noise operators B j (t), B † j (t). For this model, these are explicitly written in terms of the two-time correlation with The continuum limit of Eq. (12) can now be taken by replacing the sum over the coupling constants g jλ with an integral weighted by the density of states ρ λ of the reservoir modes. Since the only quantities entering into the physical description through Eq. (12) are ρ λ and |g jλ | 2 , we may combine their joint frequency dependence into a single spectral density function D(ω λ ), with so that in turn, the quantities Ω j act to measure the coupling strength of the j-transition of the OQS to the full set of environment modes via the expressions For the remainder of this paper our focus will be on describing the non-Markovian dynamics of OQSs coupled to various types of structured reservoir. Hence, for this purpose we shall model the system-reservoir interaction (7) using a generalized form of spectral density function which may vary strongly over ω with respect to the frequency scales of the OQS. In particular, it will only be assumed that D(ω) is a meromorphic function when analytically continued to the lower half complex ω-plane, and that D(ω) tends to zero at least as fast as ∼ O(1/|ω|) for |ω| → ∞. Under these assumptions, and with all other non-analytic features of the spectral density function removed (e.g., branch cuts), the two-time correlation function (12) may then be evaluated solely in terms of the poles and residues of D(ω) via contour integration methods. In this way, we proceed to write Eq. (12) in terms of the integral where C is a contour defined along the full real line and closed by a semicircular arc in the lower half complex plane. We note that the construction of (16) relies on the use of the RWA, which in the limit ω j → ∞ (or more loosely for ω j Ω j , λ l , |∆ jl |) formally allows an extension of the environment definition to include modes of negative frequency ω < 0. Furthermore, the poles of D(ω) in the lower half complex ω-plane are located at positions z 1 , z 2 , . . . z l , . . . with their corresponding residues denoted by r 1 , r 2 , . . . r l , . . . , while each z l has real and imaginary parts We can now apply the residue theorem to Eq.
where in the following, it will also prove useful to define the coupling constants which are in general complex quantities. For simplicity we shall first restrict ourselves to real couplings g jl . Besides the assumptions already made on the spectral density (13), this imposes no extra limitations on the model given that f j (0) = −iΩ 2 j l r l must always evaluate to the real quantity Ω 2 j [c.f. Eqs. (12)-(15)], and hence (g jl ) 2 has no net imaginary part, i.e., The most general case involving complex couplings g jl will be considered later on in Sec. V.

III. AUXILIARY MODEL
In this section we proceed to introduce the auxiliary model that will allow us to represent the reduced evolution of Eq. (10) within an enlarged open system whose dynamics is Markovian. To this end, let us first consider a mapping of Eqs. (2)- (3) in which the original Hamiltonian (1) is replaced by The mapping modifies the environment E so that the open system S is now coupled to a set of auxiliary discrete modes M , which in turn are each coupled to an independent (Markovian) reservoir with vanishing correlation time (c.f. Fig. 1(b)). Thus, the Hamiltonian of the new environment configuration E reads is the annihilation (creation) operator for an excitation of frequency ω in the reservoir R l , and . The free Hamiltonian of the full set of discrete modes M is written as while their coupling to the system S is described by Moreover, we stress that the parameters of the Hamiltonians (23)- (24) have been chosen in such a way that the discrete modes have the same one-to-one association with the poles of D(ω) as the pseudomodes introduced in Ref. [39]. Our next step is to transform H I to the same interaction picture as Eq. (7) by means of the unitary operator which following the technique of Eq. (21), has simply been obtained by replacing the free Hamiltonian H E in Eq. (6) with H E . The interaction Hamiltonian in this frame of reference reads with detunings ∆ jl = ω j − ξ l from the OQS transition Here, U M R (t, 0) ≡ U M R (t) is a unitary operator describing the time evolution of the free environment oscillators (i.e., the degrees of freedom of M +R with no coupling to the OQS), T is the chronological time-ordering operator [1], and (29) Note that we have also defined the general form of noise operator B j (t), which in the auxiliary model is the counterpart to that given in Eq. (8), Now, fixing the environment M + R to have the same initial conditions as E above-namely, by choosing an initially factorized state ρ SM R (0) = ρ S (0) ⊗ ρ M R (0), and for ρ M R (0) to be given as with

the expectation values of the environmental noises satisfy Tr
Hence, in an analogous way to Eq. (10), the reduced dynamics of will only depend on the second-order moments of the environmental noise operators B j (t), B † j (t). The only non-zero contribution written in terms of these moments With the relevant details in place we now look to prove an equivalence between the OQS dynamics generated by the two forms of interaction in Eqs. (3) and (25). Based on the discussion so far, our proof exploits the fact that given the initial choice of vacuum states (9) and (31), the reduced system dynamics of the two models will be identical as long as the two-time correlation functions of the original and auxiliary environments share the same time dependence [31,42] (see details in Appendix A). Thus, using Eq. (34) as our starting point, it remains for us to solve the Heisenberg equations of motion for the operators b l (t, s), and subsequently use this result to obtain a time dependent expression for b l (t, s)b † i (0) M R . In a Heisenberg picture generated via Eq.
Note that here we have introduced the noise operators which adopt the same definition as the so-called 'input' fields introduced in the Gardiner-Collet description of (Markovian) quantum white noise [2,45]. Accordingly, since the reservoirs R l are taken from Eq. (31) to each initially be in the vacuum states ρ R l (0) = |0 0| R l , the expectation values of the noise operators satisfy Tr R l [a in (l, t)ρ R l (0)] = 0. This then allows us to directly substitute Eq. (35) into (34) to obtain which proves the full equivalency of the two-time correlation functions of E and M + R, given that all other correlation functions equally match due to having a trivial time dependence. Therefore, because the environmental noise operators B j (t) and B j (t) have been shown to be identical correlation-wise, we may conclude that the dynamics of the OQS are indistinguishable between the two models. In other words, the non-Markovian response of the system S is invariant under replacing the original environment E by a finite number of discrete modes M coupled to independent Markovian reservoirs R; from Eqs. (10) and (33), it subsequently follows Finally, we re-state our main assumptions of D(ω) being a meromorphic function in the lower half complex ω-plane, and of the coupling constants g jl being real.

IV. EXACT SOLUTION TO THE PROBLEM
Dealing with the auxiliary model in place of Eq. (1) now enables us to reproduce the exact OQS dynamics without making any form of approximation involving weak-coupling or separation of time scales between the system and environment. To show this explicitly, we will proceed to derive the quantum Langevin equation for the enlarged open system comprising the original system S and the discrete modes M . For convenience we choose to work in an interaction picture with respect to the free Hamiltonian H 0 = H S,0 + H M + H R . In this frame of reference, the time evolution for an arbitrary operator A of the enlarged system is defined and Notice that since the two configurations of environment E and M + R are interchangeable at the level of the OQS, here we have dropped the dash used to label terms of the auxiliary model. Furthermore, the system Hamiltonian of Eq.
where H M R (t) again is given by Eq. (23).
The Heisenberg equation of motion for an arbitrary operator A(t) is written as so that by formally eliminating the reservoir variables a R l (ω, t) = U † (t)a R l (ω)U (t) from Eq. (41), one obtains the quantum Langevin equation At this stage we may equally derive the master equation for the reduced density matrix ρ(t) ≡ Tr R [ρ SM R (t)] by taking the expectation value of both sides of Eq. (42) with respect to an initially factorized density matrix ρ SM R (0) = ρ(0) ⊗ ρ R (0) (ρ R (0) = |0 0| R ). Since the noise terms proportional to a in (l, t) and a † in (l, t) do not contribute from the remaining commutators are easily expanded to obtain The resulting equation now is expressed solely in terms of operators pertaining to the enlarged OQS S +M . Hence, one may use the cyclic trace property where the superoperators describe the local dissipation of each discrete mode occurring at rate λ l . Equation (46) establishes one of the main results of this paper. The master equation is of a standard Lindblad form [26] and represents the general evolution of a reduced system density matrix ρ S (t) = Tr M [ρ(t)] as the projection of a larger quantum Markov process, where all memory effects contained in the non-Markovian dynamics are incorporated into the coupling between the OQS and discrete modes. Importantly, this general property of Eq. (46) also allows an unravelling of the master equation into Markovian pure state trajectories such that a numerically efficient simulation to the problem may be readily obtained using the quantum jump (Monte Carlo wavefunction) method [27] or other related approaches [28]. In this respect our result may be considered as a non-Markovian generalization of quantum jump method. In fact, a number of recent studies have shown similar extensions of this approach to be possible by considering the evolution a non-Markovian OQS embedded within a larger Markov process [46].

A. Application of the result
For the sake of concreteness let us briefly examine an application of our result. As a convenient example we consider a two-level system (TLS) S with excited (ground) state |e (|g ) interacting with a bosonic environment E at zero temperature. Within the RWA, the interaction Hamiltonian reads H I (t) = σ + B(t) + σ − B † (t), where σ + = σ † − = |e g| denote the Pauli matrices such that σ + |g is an eigenstate of the free Hamiltonian H S,0 = ω 0 σ + σ − with eigenvalue ω 0 . The only non-zero correlation function of environment is taken to be of the form of Eq. (18), which as before is associated to real coupling constants g l = Ω 0 √ −ir l .
To now find the exact master equation generated from the TLS-discrete mode interaction, we can apply Eqs. (40) and (46) The result here agrees with that recently obtained by Tamascelli et al. in Ref. [42] within the RWA. Moreover, when the total system contains at most a single excitation, i.e., H S (t) = 0, the above reduces to the same form of master equation derived via the pseudomode method [39], where for this case each discrete mode adopts the role of a pseudomode. Therefore, not only have we recovered the result of Ref. [39] within the restrictions considered above, but we have also generalized it via Eq. (46) to a much greater variety of open system models, including those applying to the description of multiple excitation dynamics.
To further illustrate the connection of our approach with the pseudomode theory, we note that one may obtain the equivalent spectral density to Eq. (18) [and indeed f (τ )] by inverting the general expression for the correlation function in Eq. (16). Since D(ω) is by definition a real positive function, f jk (τ ) must be Hermitian in time, i.e., and so (Re[r l ] = 0) (51) Here the normalization condition in Eq. (14) is satisfied from −i l r l = 1 [c.f. Eq. (20)]. In the case of real couplings g jl , the spectral density function is thus found to reduce to a linear combination of Lorentzians, where each of the poles in the lower half complex plane of D(ω) is connected to one of the N discrete modes in the auxiliary model. Based on this association the discrete modes are then connected to precisely same feature of the spectral density as the pseudomodes. On the other hand, it is worth also pointing out that the current restriction to real couplings g jl equally limits the range of applicability of our theory to spectral density functions with only simple poles. This illustrated by the fact that in the case of complex g jl , the form of spectral density obtained by inverting Eq. (16) changes D(ω) from a sum of Lorentzians to a more general (non-Lorentzian) form of function with possibly higher order poles. Therefore, if we are extend our generalization of the pseudomode method even further it is now clear that we must look beyond the assumption of real couplings g jl .

V. NON-HERMITIAN INTERACTION HAMILTONIAN
Remarkably, as is shown in Appendix A, the considered proof of Eq. (38) does not actually rely on the assumption of g jl being real-rather, the only condition necessary to guarantee the equivalence of ρ S (t) and ρ S (t) is for the correlation functions of the original and auxiliary environments to satisfy Eq. (37). Hence, with the identity (38 still valid, one may proceed to write Eq. (19) in terms of its real and imaginary parts, The interaction Hamiltonian (40) will now comprise of both a Hermitian and non-Hermitian part where each respectively read By following the steps in the previous section we may again derive the equivalent master equation for enlarged system. However, due to non-Hermitian nature of H I (t), we find the resulting master equation to no longer be of Lindblad form, and in turn for the evolution of ρ(t) to in general be nonpositive. Interestingly, we are then for this case able to generate a pathological form of master equation that, while is still capable of describing the correct OQS behavior, lacks a suitable physical interpretation for the dynamics of the enlarged system S +M . In particular, this may present issues for solving Eq. (46) using quantum jump methods, given that the unravelling of ρ(t) into quantum trajectories from the current form of master equation would admit the possibility of jumps occurring between states of the discrete mode system with probabilities exceeding unity, among other unphysical defects. Below we shall focus how the master equation (46) may be brought into Lindbald form by applying a suitable transformation to the discrete modes.

A. Lindblad construction of the master equation
The same issues connected with replacing an environment by a discrete set of bosonic modes which have complex coupling constants to the OQS have also previously been encountered with the pseudomode method [39]. In that context, i.e., with a RWA form of interaction and single excitations, a procedure was established to convert between pathological and Lindblad forms of the master equation by applying an effective change of basis to the pseudo (discrete) mode operators b l (b † l ). Hence, in the spirit of Ref. [39], we enact a similar procedure below by first rewriting Eq. (44) in the Schrödinger picture: Next we introduce a new set of discrete mode operators with U ml an orthogonal (real) matrix and [b m ,b † m ] = δ mm . Inserting the above decomposition into the interaction term of Eq. (56) now leads to so that fromg jm = l U ml g jl , we can in principle fix the elements U ml of an otherwise arbitrary matrix by requiring the new couplingsg jm to be real. We also note with this definition that the couplings are constrained to satisfy the same normalization as Eq. (20), i.e., m (g jm ) 2 = Ω 2 j . For the remaining terms, the effect of Eq. (57) will follow that of a similarity transformation applied to the diagonal matrices z l and λ l , where The non-Hermitian Hamiltonian H eff can subsequently be written in the form with Γ mm = −Im[z mm ], and (62) Now, since Re[z mm ] is symmetric by definition, the master equation resulting from Eqs. (59) and (60) will only be in Lindblad form if we impose the additional constraint that the matrix Γ mm is Hermitian and positivedefinite [26]. Under this restriction, and in an interaction frame generated by the free Hamiltonian H 0 = H S,0 + m Re[z mm ]b † m b m , one can then obtain the following (Lindblad) master equation: whereH is written in terms of the transition frequenciesξ m = Re[z mm ] and detunings∆ jm = ω j −ξ m . Thus, while the OQS dynamics remains unaffected by the transformation, we find that the removal of pathological terms from Eq. (46) generally introduces a non-zero coupling between the discrete modes. Starting from Eq. (55), the task of converting between a non-Lindblad master equation and a Lindblad one now amounts to finding a transformation matrix U ml with the simultaneous requirements forg jm to be real and for the dissipation matrix Γ mm to be positive-definite. Unfortunately for the general case involving N discrete modes, this constitutes a highly nontrivial problem that has no guarantee of a unique solution. However, for the cases in which the OQS couples to only two discrete modes (N = 2), it is in fact possible to obtain analytical expressions forz mm andg jm following directly the techniques of the pseudomode method [39]. Indeed, as determined from Sec. VB of that paper, one may here parameterize the matrixz mm as (m, m = 1, 2) where the analytical expressions for the new discrete mode decay rates Γ m , coupling constants V 12 = V * 21 , and transition frequenciesξ m are given in Appendix C.

B. Further application of the result
To provide an application of the theory outlined above let us finally return to the OQS model introduced in Sec. IV A-namely, of a TLS S interacting with a zero temperature bosonic field E-to consider the case in which the spectral density is instead given by the difference of two Lorentzian functions: (66) This form of spectral density has previously been used to model an environment exhibiting a photonic band gap, where the localized gap in the density of states occurs for λ 2 W 1 = λ 1 W 2 [39]. For the analysis below, we also note that the conditions W 1 λ 1 > W 2 λ 2 (λ 1 > λ 2 ) and λ 2 W 1 ≥ λ 1 W 2 are generally required since D(ω) is defined always to be positive, while W 1,2 must be normalized to W 1 − W 2 = 1 [c.f. Eq. (14)].
Considering that Eq. (66) has poles located at positions z 1,2 = ξ − iλ 1,2 in the lower half complex ω-plane, the two-time correlation function (18) reads with the discrete mode couplings to the TLS given by Then, as g 2 is clearly pure imaginary, one can apply the relevant formulas from Appendix C to obtain the new frequency parameters for the TLS and discrete modes: This implies the TLS to physically couple to the second auxiliary mode, which in turn interacts with the first mode, as well as additional Markovian reservoirs. The reduced Lindblad-type master equation for the enlarged system can subsequently be obtained by inserting Eqs.
(65) and (68) into (63), where and As found previously with Eq. (49), this result similarly generalizes the equivalent pseudomode master equation derived in [39] to now apply to a much greater variety of OQS models, including those not limited to the single excitation regime.

VI. SUMMARY AND OUTLOOK
In conclusion, we have derived a master equation that provides a nonperturbative and non-Markovian description of an OQS dynamics within the RWA. This has been achieved by showing the reduced dynamics of an OQS to be indistinguishable under the effect of two different types of structured environment-one comprising an infinite collection of harmonic oscillators with a frequency-dependent coupling to the system, and an auxiliary one comprising a small number of discrete modes which are each coupled to an independent Markovian reservoir. The equivalence of these two models has subsequently been exploited to construct a general Lindblad master equation. In this way the reduced dynamics of the original problem can be simulated efficiently within the framework of quantum trajectories, assuming a reasonable Hilbert space dimension for the enlarged system. In particular, our approach has shown a full extension of the pseudomode theory to be possible to cases involving multiple excitations, where the strength of the OQS-discrete mode coupling and the degree of separation of timescales between the system and environment poses no restriction on the validity of the result.
The procedure we have introduced to obtain the master equation relies on the connection between the poles of the spectral density contained in the complex ω-plane, their residues, and the properties of the discrete modes used to represent the memory part of the environment. Specifically, the decay rates and couplings of the discrete modes are determined directly from the positions and residues of the poles of the spectral density function located in the lower half complex plane. When the residues of these poles are evaluated to give complex coupling constants g jl , we have shown that our approach can still be applied using a non-Hermitian form of interaction Hamiltonian which reproduces the exact physical dynamics of the original open system. The master equation we obtain in this instance contains pathological terms that in general violate the positivity of the auxiliary system density matrix. For the two-discrete mode case, these issues have been rectified via a change of basis of the discrete mode annihilation (creation) operators to obtain a Lindbladtype master equation.
Finally, a current open problem relating to construction of the master equation is how to convert between the pathological and Lindblad form for when the OQS couples to more than two discrete modes with complex coefficients; as far as we are aware, no generalized form of transformation matrix allowing for the conversion has been determined beyond this case (although see perhaps the related inversion problem explored in Appendix C of Ref. [47]). Therefore, future work in this area could focus on developing a systematic approach to regularizing the master equation for arbitrarily complicated environment structures, which in itself would likely rely on a numerical implementation.

ACKNOWLEDGMENTS
This work is based upon research supported by the South African Research Chair Initiative of the Department of Science and Technology and the National Re-search Foundation of the Republic of South Africa.
Appendix A: Equivalence between the reduced density operators ρS(t) and ρ S (t) Here we outline the proof showing that the reduced system dynamics of ρ S (t) and ρ S (t) will be equivalent as long as the two-time correlation functions (12) and (34) share the same time dependence. To this end, let us first restate our main assumptions of factorizing initial conditions ] is restricted to the vacuum state (9) [ (31)]. We may then proceed by noting that a general solution to the reduced density operator ρ S (t) [ρ S (t)] can be written as an expansion over the noise operators in Eq. (7) [ (27)] [42,44]. Indeed, within an interaction picture generated by the unitary transformation (6), the evolution of the density operator ρ SE (t) is determined by the von Neumann equation which can immediately be solved to obtain If we now expand the Dyson series of Eq. (A2) and trace out the environmental degrees of freedom to obtain ρ S (t) = Tr E [ρ SE (t)], the reduced system density operator can be derived in the form where L(t 1 )L(t 2 ) · · · L(t n ) E (t 1 ≥ t 2 ≥ ...) contains the n th order moments of the environmental noise operators B j (t), B † j (t) (and system coupling operators c j (t), c † j (t)), as well as terms connected to the free system Hamiltonian H S (t). Similarly for ρ S (t), we may also introduce the von Neumann equation and follow an analogous procedure to obtain It can therefore be observed in the general case that the OQS dynamics of the reduced density operators will only be equivalent if first, second, etc., order terms in the two expansions have an equal time dependence. However, taking into account that the states ρ E (0) and ρ M R (0) are Gaussian, Wick's theorem implies the complete set of terms appearing in either expansion to only depend on the first and second order moments of the corresponding noise operators, so that from ], the only difference between ρ S (t) and ρ S (t) will be contained in the expressions for the two-time correlation functions (12) and (34). Now, since the free system dynamics is unchanged between the different environment configurations, as well as the coupling operators c j (t) (c † j (t)), we see that if then the reduced expansions (A3) and (A5) will be indistinguishable in the sense Φ(t, 0) = Φ (t, 0) =⇒ ρ S (t) = ρ S (t).
In particular, we note that the proof makes no requirement for the couplings (19) to be real.
pathological and Lindblad forms of the master equation with equivalent coefficients to Eqs. (55) and (63); here we then utilize this very same procedure to overall obtain a more general result. Considering firstly that the transformation matrix U ml for m, l = 1, 2 can generally be written in the form of a complex rotation about an angle θ 0 , U (θ 0 ) = cos θ 0 sin θ 0 − sin θ 0 cos θ 0 , the similarity transformation in Eq. (59) can be evaluated explicitly to yield z mm = 1 2 z 1 + z 2 + (−1) m ∆z cos 2θ 0 δ mm + ∆z sin 2θ 0 2 δ mm −1 + δ mm +1 , m, m = 1, 2, where ∆z = z 2 − z 1 ≡ |∆z|exp(iθ z ), defines the distance between the two poles of D(ω) in the lower half complex plane. In turn, we may also proceed to parameterize the OQS-discrete mode couplings g jl = Ω j √ −ir l via which by construction satisfies the normalization property (20), i.e., g j 2 = Ω 2 j . Since the relative magnitude of these couplings depends only on a single free parameter θ 1 , we can more conveniently characterize Eq. (C4) in terms a complex ratio µ, where µ = tan θ 1 = g j2 /g j1 . (C5) We will now proceed to write the new couplingsg jm as so that with the help of these basic definitions, the rotation angle θ 0 can in principle be determined according to the physical constraints placed on the quantitiesg jm and Γ mm = −Im[z mm ] (see the details provided in Sec. V A of the main text). This is in part what we shall now discuss below. Following Ref. [39], we first note that the off-diagonal elements of Eq. (C2) are constrained to be real, V 12 = V * 21 ; as such, we may writẽ While this is not a necessary condition for the approach, due to the fact that the resulting matrix Γ mm is diagonal, i.e., Γ mm = Γ m δ mm , its positivity can now be guaranteed by the simpler constraint that only Γ m must be positive: This inequality, which relates the allowed values θ 0 to the parameters of correlation function (18), may or may not be possible to satisfy depending on the choice λ l and ξ l , so that its validity must be checked in the general case (though this should not be too demanding to satisfy based on the examples shown in Ref. [39]). The second constraint comes from the requirement thatg jm must be real, and so Im[θ 0 ] = Im[θ 1 ].
Thus, the problem of determiningz mm andg jm reduces to eliminating θ 0 from Eq. (C8) through the constraints (C9)-(C10). Although this will not be explicitly shown here, it ultimately follows from Ref. [39] that z mm (ξ m , Γ m , V 12 ) andg jm can be expressed solely as functions of the known variables ∆z, µ and θ z . Indeed, for the coupling between discrete modes V 12 , we find