Open quantum systems with delayed coherent feedback

We present an elementary derivation and generalisation of a recently reported method of simulating feedback in open quantum systems. We use our generalised method to simulate systems with multiple delays, as well as cascaded systems with delayed backscatter. In addition, we derive a generalisation of the quantum regression formula that applies to systems with delayed feedback, and show how to use the formula to compute two-time correlation functions of the system as well as output field properties. Finally, we show that delayed coherent feedback can be simulated as a quantum teleportation protocol that requires only Markovian resources, pre-shared entanglement, and time travel. The requirement for time travel can be avoided by using a probabilistic protocol.


Introduction
While the theory of Markovian open quantum systems is well-understood [1,2], simulating the dynamics of non-Markovian open quantum systems is considerably more difficult [3]. A signature non-Markovian open quantum system is a qubit emitting into a discrete feedback reservoir: an environment that 'remembers' the state of the system and feeds this information back coherently -as a quantum field -after one or more discrete time delays. One example of such a system is sketched schematically in Fig. 1. This kind of feedback has previously been studied in work on "atomic" emission in front of a mirror [4][5][6][7], and also in solid state systems with significant propagation delays [8]. Such systems are all the more interesting due to the fact that an environmental memory with a continuous kernel, where the evolution of the system depends most generally on its state at all previous times, may be approximated as a sequence of coherent feedback loops with discrete delays. As such, a tool that proves capable of simulating multiple discrete delays may shed light on the more general problem of dealing with a continuous memory.
Discrete propagation delays also appear in the standard theory of cascaded open quantum systems, which describes the situation where the retarded output from one open quantum system drives a second system [9][10][11]. In the standard treatment of cascaded systems there is no backscatter from the second system, so the coupling between the systems is one-way. In this case the propagation delay between the systems is an arbitrary parameter that can be removed by way of a simple transformation of the time variable, which leads to an irreversible, Markovian coupling between subsystems. The cascaded systems formalism is easily generalised to describe irreversible coupling between the systems in both directions, so long as the propagation delay associated with the coupling is sufficiently small that it may be neglected [12]. We cannot, however, transform away a non-negligible delay in both directions. A treatment of open quantum systems with delayed feedback is therefore necessary to extend the theory of cascaded open quantum systems to encompass cascaded systems with backscatter where there is a nontrivial propagation delay.
There are conditions under which delayed feedback and cascaded systems are well known to be connected. In fact, in a certain limit a coherent feedback loop leads to exactly the same dynamics as a chain of cascaded systems. This equivalence is employed in, for example, work by Menicucci et al. [13]. While it is not the case that cascaded systems and delayed feedback are equivalent in general, it turns out that there is sufficient similarity between these two set-ups that results from the theory of cascaded systems can be exploited to help simulate feedback. Recently, Grimsmo [14] employed tensor network methods to show how to simulate a nonlinear quantum system (such as a two-state 'atom') interacting with a discrete feedback reservoir and indeed, the resulting equations demonstrate a close connection to the master equation for cascaded systems. Pichler and Zoller [15] subsequently published a different technique, based on matrix product states, for simulating quantum circuits in the regime where time delays are significant. We present here an elementary derivation and generalisation of Grimsmo's method, which permits the simulation of multiple delays, as well as cascaded systems with delayed backscatter where the delay may differ in each direction.
Our derivation is based on the fact that the evolution of a generic open quantum system may be decomposed into a nested sequence of evolutions over distinct time intervals-what we refer to as a decomposition into intervals. This decomposition is outlined in Sec. 2. In Sec. 3, we apply this decomposition to open quantum systems interacting with the environment such that the system experiences delayed coherent feedback, with discrete delays of various lengths applying to interactions between different pairs of subsystems. We go on in Sec. 4 to present a series of examples of the use of this algorithm, and demonstrate in Sec. 5 how the algorithm may be extended to enable the calculation of two-time (and more generally multi-time) correlation functions, presenting an example calculation of the second-order photon correlation for the output field of a system interacting with a delayed coherent feedback loop. We then show in Sec. 6 that our algorithm for simulating open quantum systems with delayed coherent feedback may be interpreted as a quantum teleportation protocol. We summarise our main results and S A S B Figure 1. A bipartite open quantum system, depicted here as a pair of ring cavities, with delayed coherent feedback. The propagation delay from subsystem A to subsystem B is greater than that in the reverse direction. Black arrows denote system fields, while grey arrows denote fields propagating in the environment.
discuss their application to more general systems in Sec. 7.

Decomposing the evolution of an open quantum system into intervals
The evolution of an open quantum system may be decomposed into intervals. This decomposition, based on the work of Grimsmo [14], is outlined below. In this section we will justify it for a generic system, before applying the resulting simulation algorithm to the case of a system with delayed coherent feedback in Sec. 3. We begin by summarising the structure of the decomposition implemented as dynamical map. Following Gardiner and Zoller [16] we introduce a complete set, {e j }, of basis operators for the system that are orthogonal with respect to the trace: tr e † j e k = δ jk . The state of the reduced system (with the environment traced out) can be expanded, as usual, in terms of these basis operators: The evolution of the basis operators (and hence the system) is divided into intervals of length ξ, with each interval represented by a formally separate Hilbert space. To divide the system's evolution into n ≡ t/ξ intervals, we require n copies of the system's Hilbert space. The total evolution time t will not necessarily be an integer multiple of the interval length, so we define the auxiliary time variable t ≡ t − (n − 1)ξ, as illustrated in Fig. 2. The evolved basis operator e j 1 (t) is obtained using the mapping formula e j 1 (t) = tr 1 · · · tr n−1 where we have defined the product basis operators e j 1 ···jn ≡ e j 1 ⊗ · · · ⊗ e jn , and where e j 1 ···jn (ξ; t ) is a basis operator of the entire fictitious 'chain' of system copies, given by Decomposition of an open quantum system's evolution into intervals. Here, we split the evolution of a system up to time t into n = 3 intervals of length ξ, making use of the auxiliary time variable t = t − (n − 1)ξ. Black arrows denote the evolution of all system copies up to time t , corresponding to the map Φ (n) (t , 0). Gray arrows illustrate the evolution of all systems save the last from t up to ξ, corresponding to Φ (n−1) (ξ, t ).
At this stage, we simply assume the existence of the dynamical map Φ (m) (t 1 , t 0 ) that describes the evolution of the first m system copies in the chain from t 0 up to t 1 ; in Sec. 2.3 we will demonstrate a convenient method of deriving this map, by using an effective representation of the environment's entire timeline. We can see that Eq. (3) describes the same decomposition of the system's timeline as Fig. (2): the dynamical map is used to evolve the chain of system copies in its entirety up to time t , after which all systems up to but not including the last in the chain (the 'present' interval) are evolved all the way up to time ξ. The complete timeline of each basis operator's evolution is reconstructed using Eq. (2), and these basis operators are then used in Eq. (1) to obtain the state of the system. The map Φ (m) (t 1 , t 0 ) is met in its simplest form by considering a closed quantum system and starting out with more familiar notation.

A closed quantum system
We consider the evolution of a density operator, ρ(t), with initial condition written in matrix notation, where µ|ν = δ µν . The unitary time-evolution operator is denoted where the time ordering operator T orders products of time dependent operators such that their time arguments increase from right to left. We break the evolution from the initial time 0 to time t into two intervals: an initial interval of length ξ and the remainder. The time evolution is generated by propagating the elementary operators (|µ ν|). By interrupting the evolution at ξ and re-expanding the result in matrix notation before proceeding, we get Then returning to our more compact notation, we may write where we allow a single index to stand in for the pair (µ, ν), with e j 1 -a basis operator as introduced above -standing in for (|µ ν|) and e j 2 for (|µ ν |). Finally, we note that the results required from the two separate steps in the time evolution may be computed in formally distinct Hilbert spaces, which we index as space 1 (interval 0 to ξ) and space 2 (interval ξ to t), and then "parked" in appropriate locations in a tensor product of those two spaces. This allows us to re-express Eq. (6) using tensor product notation: The division into n intervals gives an obvious extension of the formula, which turns out to be exactly Eq. (2) with where we have used the fact that operators of different system copies commute to group all the unitary operators together. Equation (8) evolves the basis operators for each of the system copies 1 to (n − 1) through an interval of length ξ, and Eq. (2) maps these final conditions onto the corresponding initial conditions for system copies 2 through n, respectively. The n th system is evolved through an interval of length t , where we recall that t ≡ t + (n − 1)ξ. As such, this process reconstructs the evolved basis operator e j 1 (t) and (by completeness) the state of the 'real' system. In Eq. (8), we have separate unitaries describing the evolution of each system copy. In Sec. 2.3, however, we will consider open systems with delayed feedback, with the feedback simulated by way of a single effective environment shared between all system copies in the chain. The feedback will appear in this formalism as interactions between different system copies, meaning that we will no longer be able to treat the evolution of each interval independently. For that reason we now generalise our notation, creating a combined unitary operator that evolves all systems in the chain simultaneously. To write this generalisation in a compact way, we first need to define some notation which will be used throughout the rest of this paper. We define for any system operator A, where I is the identity operator for a single system copy. We now define an evolution operator acting on all n system copies: where in this context 0 ≤ t 0 ≤ t 1 ≤ ξ, as this unitary will only be used to evolve the chain on the interval [0, ξ]. It is easily verified that which means that the evolved basis operators in Eq. (8) can now be written in the form It follows that in the simple case of a closed system the map introduced in Eq. (3) is

An open quantum system
We now begin the process of generalising to open quantum systems. Firstly we define, in general terms, the system whose evolution we wish to decompose. Suppose the evolution of an open quantum system is described, in a rotating frame, by the Hamiltonian where H S involves only system operators and describes the internal dynamics of the system, and the interaction between system and environment is generated by where Here ω 0 is some fiducial frequency that may be freely chosen. The index α labels different subsystems, while j labels modes of the environment. System and environment operators commute at equal times, and the environment will be assumed to be an assemblage of harmonic oscillators with It is also useful to define the dissipation (memory) kernel of the reservoir: We assume that the initial combined state of the system and environment is separable and given by ρ(0) = ρ E ⊗ ρ S (0), where ρ S (0) and ρ E are the initial states of the system and environment respectively.
The derivation for an open quantum system proceeds much as in Sec. 2.1. We start by making the trivial expansion where in this case the operators (|µ ν|) form a basis for the system. We once again break the evolution from 0 to t into two intervals, interrupting the evolution at ξ and re-expanding the system-but not the environment: Just as we did above in the case of the open quantum system, we replace the pair of indices (µ, ν) with a single index by introducing basis operators {e S;j } for the system, which gives We now introduce a new unitary operator where a α;m is an operator of the m th system copy, expressed using the notation introduced in Eq. (9), and where 1 (m−1)ξ≤s<mξ is an indicator function, defined as 1 when (m − 1)ξ ≤ s < mξ and 0 otherwise. This unitary allows us to re-write Eq. (16) using tensor products. Generalising at the same time to n intervals, we find Finally, we trace out the environment to obtain the state of the reduced system, where e S;j 1 (t) = tr S;1 · · · tr S;n−1 and where, recalling that t ≡ t + (n − 1)ξ, we have defined Equations (19) and (20) are, up to minor notational differences, identical to Eqs. (1) and (2). In Sec. 2.3 we will derive an effective representation of the system-reservoir dynamics that will allow us to re-write Eq. (21) in terms of a dynamical map, in the form of Eq. (3).

Effective dissipation kernel for an open quantum system
In Sec. 2.2 we demonstrated a decomposition of the evolution of an open quantum system into intervals, but we have not yet derived a map describing the evolution of the reduced system. We now introduce a new system-environment interaction that will allow us to derive a reduced system map. We refer to this as an effective system-environment interaction, because it will be defined in such a way that its effects on the system reproduce the evolution of the real system once the environment has been traced out. Provided that the resulting reduced system map turns out to be divisible, the algorithm described in Sec. 2.1 may be extended to open quantum systems. As we will see in Sec. 3, when dealing with delayed coherent feedback this effective interaction leads to a considerable simplification, because it means the feedback may be modelled as an interaction between intervals/system copies.
To find this effective system-environment interaction we define a new Hamiltoniañ H m (t ) ≡ H S;m +H SE;m (t ) for a single system copy with interaction term We also define the corresponding unitary evolution operator Here, 'effective' quantities that differ from their analogues from Sec. 2.2 are denoted by tildes. Note that, consistent with the notation introduced in Eq. (9), the Hamiltonian (22) acts only on the m th system copy. In an analogue of Eq. (10), we define a corresponding effective unitary for the chain of fictitious system copies, Equation (24) is similar to Eq. (10); however, note the difference in the time argument of the Hamiltonian under the integral. Using this unitary, we define a reduced system map: where χ is some operator of the reduced system copies, and ρẼ is the initial state of the effective environment. We require that this reduced system map reproduce Eq. (21) when used in Eq. (3), in the sense that we should have as well as with e S;j 1 ···jn (ξ; t ) given by Eq. (21). As discussed in Appendix A, Eq. (26) is satisfied whenever the map Φ (m) (t 1 , t 0 ) is divisible [17,18] for all m, meaning that In the following, we will show that the effective environment may be structured such that Eq. (27) also holds. The first restriction we impose on the effective environment operatorsB α;m (t ) is the stipulation that B α;m (t 2 ),B † β;m (t 1 ) = 0 when t 2 ≥ t 1 and m < m .
From this commutation relation, we easily find that a similar relation holds for the effective Hamiltonian: H m (t 2 ),H m (t 1 ) = 0 when t 2 ≥ t 1 and m < m .
Equation (29) thus means that a system copy appearing earlier in the chain of decomposed intervals is independent of (commutes with) any subsequent system copy in the chain when the former is considered at a later time than the latter. This ensures that a given system copy cannot be affected by any system copy in its relative 'future' -a natural requirement of causality for the real system. Put another way, causality in the chain of system copies runs from m = 1 towards m = n. As such, it follows from Eq. (30) that, analogous to Eq. (11), we can writẽ which allows us to re-write Eq. (27) as  (21). We do not have to specify this operator explicitly: if we restrict attention to vacuum reservoirs -by assuming that both ρ E and ρẼ represent the vacuum states of the respective reservoirsthe interaction of the system copies with the environment may be fully characterised by the dissipation kernel. As such, our task becomes to derive the effective dissipation kernel ] such that Eqs. (21) and Eq. (32) are equivalent. After a few simple manipulations of these two equations, as detailed in Appendix B, we find that this condition is satisfied by the effective dissipation kernel with t 2 ≥ t 1 , where the case m < m follows from Eq. (29). Note that the dissipation kernel is, by construction, a Hermitian function:F αβ;mm (t 1 − t 2 ) =F βα;m m (t 2 − t 1 ). In Eq. (33) we focus on the case t 2 ≥ t 1 , as this is the only case relevant to the dynamics. Henceforth we will work entirely with the effective dissipation kernel (33), and omit the tildes denoting effective quantities for the sake of simplicity.
We have thus shown that Eq. (27) is satisfied when the interaction of the system copies with the effective environment is described by the dissipation kernel (33). If the resulting map, given by Eq. (25), is divisible then Eq. (26) is also satisfied. Combining Eqs. (26) and (27), we obtain where e S;j 1 ···jn (ξ; t ) is given by Eq. (21). Eq. (34) is, up to minor notational differences, Eq. (3).

Delayed coherent feedback
We want to use the decomposition into intervals described in the previous section to simulate a coherent environmental feedback loop with discrete delays. With that in mind we specify the dissipation kernel with τ αβj ≥ 0, which describes exactly such a system. We assume that the delays τ αβj are commensurable, with the interval length ξ chosen to match the greatest common divisor of the delays so that k αβj = τ αβj /ξ is an integer for all α, β, and j.
Our task now is to derive the map Φ (m) (t , t 0 ), t 0 ≤ t ≤ ξ, describing the reduced system dynamics associated with the dissipation kernel Eq. (35). Substituting the dissipation kernel into Eq. (33) yields Recalling that, in the context of Eq. (3), t ≤ ξ, and noting that both m − m and k αβj are non-negative integers, we can see that only those terms in Eq. (36) for which k αβj = m − m contribute to the evolution of system operators on each interval. As such, it is easily seen that the dissipation kernel (36) describes n identical systems coupled to a common reservoir in which the output of each propagates in the direction of increasing m.
As is well-known, this situation is described by the theory of cascaded open quantum systems. It is therefore straightforward to show that, when the initial state of the combined system is separable and the environment is in its vacuum state, Eq. (22) leads to the following Liouvillian generator of the dynamics: The correspondence shown here between delayed feedback and cascaded systems is analogous to the well-known method of solving a classical delay-differential equation by re-casting it as a multivariate Markov process [19]. The fact that the dissipation kernel (35) describing delayed coherent feedback maps to the generator (37) for cascaded open quantum systems is what makes the algorithm presented in Sec. 2 so useful in this particular case: we already know how to solve the dynamics obtained by splitting the evolution into intervals.
Having found the Liouvillian (37), we are in a position to write down the map Φ (m) (t , t 0 ) that first appeared in Eq. (25): The density matrix of the real system is then obtained by substituting this map into Eq. (3), and subsequently using Eqs. (2) and (1). The Liouvillian generator (37) may be written in Lindblad form, and as such the map (38) is divisible [17,20] and, as discussed in 2.3, satisfies Eq. (3).
In the case of a single system with a single delay, this algorithm is equivalent to that derived by Grimsmo [14]. The presentation given above is, however, more general, in the sense that it admits multiple subsystems and multiple delays. We show below how to use this formalism to describe feedback with more than one delay, as well as cascaded systems with delayed backscatter in the case where the delay differs in either direction.
The algorithm reported here is not computationally efficient for long times. Because a copy of the system Hilbert space is required for each and every ξ-interval we wish to simulate, both the number of basis operators that must be evolved and the dimension of these operators increases exponentially in the number of intervals. This means that it is not practical to simulate beyond a few ξ. In particular, while we have established that the algorithm can handle multiple commensurable delays, the exponential scaling of memory requirements will pose difficulties if any individual k αβj is large. For this reason, in the following examples we have restricted ourselves to situations where all k αβj are small integers. Furthermore, the restriction to small numbers of intervals means that in many situations the steady state is not accessible using this algorithm. The advantage of our presentation is that it highlights the connection between networks of cascaded systems and delayed coherent feedback, generalising the earlier work of Grimsmo. Although the Liouvillian (37) is exactly that of an array of cascaded systems, delayed feedback is not equivalent to such an array at the level of physical systems. The mapping rule (2), enabled by the decomposition into basis operators, is required to correctly account for correlations between the different time intervals when these intervals are represented as separate systems; it does so simply by mapping the final state of each fictitious system copy onto the initial state of the next system copy in the chain. A completely analogous algorithm can be used to simulate classical delay-differential equations.

Examples
We turn now to our examples. Figures 3 and 4 illustrate four categories of system, listed below. In each case, we display a schematic depiction of the system in question, as well as a sketch of the corresponding array of cascaded systems; finally, we display the results of simulations (performed using a program based on QuTiP [21,22]) with selected parameters.
While the dissipation kernel Eq. (35) is quite general and describes multiple reservoirs, for simplicity we will initially focus on a single system coupled to a single reservoir. Firstly, we reproduce for comparison the now well-understood case of a driven qubit emitting into a feedback loop with a discrete delay [14]. The feedback loop is described by the dissipation kernel where φ = ω 0 τ . The internal system Hamiltonian, describing coherent driving with Rabi frequency Ω/2, is Ω(σ − + σ + ). This system is depicted schematically in Fig. 3(a), and the corresponding cascade of system copies is shown in Fig. 3(b). Our simulation results are shown in Fig. 3(c). Secondly, we examine a qubit, driven as above, that instead couples to the reservoir at N locations with equal spacing τ . This is described by the coupling constant κ j = N −1 n=0 γ/N e iω j nτ ; the corresponding dissipation kernel is Note that this dissipation kernel can also be thought of as describing N separate reservoirs, each with a single feedback loop: from the perspective of the system alone, the dynamics are identical. Figures 3(f) and 3(i) show sample simulations with the dissipation kernel (40), with finite N . In the limit N → ∞, this dissipation kernel describes the situation depicted in Fig. 3(j), in which the system couples to a reservoir that loops back on itself without any irreversible dissipation. This can be thought of as a system coupled to a (very long) multi-mode cavity which acts as a one-dimensional waveguide that feeds anything that gets into the waveguide back to the system first after one round trip, then also after two round trips, three, and so on. Thirdly, we turn attention to the case of more than one reservoir: we consider cascaded qubits with backscatter, that is a pair of driven qubits arranged such that the output from each subsystem drives the other subsystem with a propagation delay τ in both directions. This system has previously been examined by Pichler and Zoller [15] using a different technique. The reservoir is characterised by the dissipation kernels: The internal system Hamiltonian is given by H = Ω(σ −,A + σ +,A ) + Ω(e iφ σ −,B + e −iφ σ +,B ). Finally, we generalise the model of cascaded qubits with backscatter to consider delays that differ in each direction. More specifically, we consider the case where the delay from subsystem B to A is twice the delay from A to B. The reservoir is characterised by the kernels (41) and (42), along with We have here supposed that the phase advance in each direction is the same. The internal system Hamiltonian is as in the previous example. These examples are presented only to illustrate the applicability of the derived algorithm. Of course, there is much more that could be said about the physics of any one of them, or variations on the driven qubit setup.

Two-time correlation functions
The algorithm can be further developed to allow computation of multi-time correlation functions. We consider here two-time correlation functions, though the method is easily S (a) (g)-(i) Three feedback loops, with round-trip times τ , 2τ , and 3τ . (j)-(l) An "infinite loop", as described by the dissipation kernel (40) in the limit N → ∞, where the system is coupled to a multi-mode cavity, which acts as a one-dimensional waveguide that feeds the system output back after one round trip, then after two round trips, three, and so on. (a) generalised. Denoting by ρ S+E (t) the combined state of the system and reservoir, we trivially find and where we recall that U (t, t 0 ) is the unitary generated by H S +H SE (t), with interaction part given by Eq. (13). Equations (45) and (46) are sometimes referred to as the quantum regression formula. Once again, we divide the integration time t into intervals of length ξ. As one might intuitively expect, the operators C and A are applied in the l 1 th fictitious system copy at time t 1 , with t n = t n − (l n − 1)ξ with l n = t n /ξ . That is to say, we modify Eq. (3) to read: We similarly modify Eq. (2): e j 1 ;CA,l 1 (t 2 , t 1 ) = tr 1 · · · tr n−1 j 2 ···jn e j 1 ···jn;CA,l 1 (ξ; t 2 , t 1 )[e † j 2 ···jn ⊗ I] . We can use Eqs. (47) and (48) to calculate properties of the environment. As an example, consider the second-order photon correlation function is the field in the α th subenvironment outside the feedback loop, measured at the location of the output from the system. It is easy to show using Heisenberg picture methods [24][25][26] that which allows us to express g (2) α (t 1 , t 2 − t 1 ) in terms of system operators. A sample calculation is presented in Fig. 5, which shows the second order correlation function (in the transient regime) for the output field from a driven qubit emitting into a feedback reservoir with a single discrete delay, as described by the dissipation kernel (39), for various different t 1 .

The cascade algorithm as a quantum teleportation protocol
In this section we show that the mapping rule (2) may be interpreted as a quantum teleportation protocol [27]. For simplicity, we consider the evolution up to a time τ < t ≤ 2τ , for which we require two system copies in the fictitious chain of our algorithm; the generalisation to n systems is straightforward.
Alice has two systems S 0 and S 1 , while Bob has a single system S 2 . All three systems are isomorphic. Systems S 0 and S 2 are prepared in the maximally-entangled state D −1/2 µ |µ ⊗ |µ , with D the subsystem dimension, while S 1 is prepared in the initial system state ρ(0). Thus, the initial state of the combined system can be written as An orthonormal basis of maximally entangled states is given by {|ψ (pq) } with where µ ⊕ q ≡ (µ + q) mod D.
To recapitulate the conventional quantum teleportation protocol, suppose Alice performs a measurement on her systems S 0 and S 1 in this basis. The resulting state of Bob's system S 2 , conditioned on Alice measuring the outcome pq, is If Bob then applies the unitary he is left with a copy of ρ(0) in his system. Now suppose Alice is also in possession of a time machine, so that she can send the outcome pq of her measurement to Bob such that it reaches him before the measurement is performed. Bob applies the necessary unitary transformation to his system before the experiment begins, resulting in the combined system state To make contact with the cascade algorithm, the systems S 1 and S 2 -belonging to Alice and Bob respectively -are now acted upon by the map Φ (1) The result is of course Alice then performs the same measurement as before, sending the result back in time so that Bob can perform the correct unitary ahead of time. This procedure leaves Bob with the state ρ (pq) 2 If we now make the same notational change as in Sec. 2.1, with a single index standing in for the pair (µ, ν) and e j standing in for (|µ ν|), Eq. (57) becomes which is equivalent to Eqs. (1), (2) and (3) together in the case n = 2. The need for time travel is of course unphysical, but the result can be reproduced in a probabilistic fashion: Bob applies no unitary (or, equivalently, U (00) = I) and the protocol succeeds whenever Alice obtains the outcome pq = 00.

Conclusion and outlook
We have derived a generalisation of a technique, derived previously by one of us, for simulating feedback in open quantum systems. Our derivation uses only elementary methods and is based on decomposing the time evolution of a general open quantum system into intervals represented as separate system copies-that is to say, this decomposition is not limited to systems exhibiting delayed coherent feedback. The resulting simulation method admits multiple subsystems with multiple delays in cases where those delays are commensurable. We used our generalised method to simulate systems with multiple delays, including cascaded systems with delayed backscatter. In addition, we presented a generalisation of the quantum regression formula that applies to systems with delayed feedback, and demonstrated how to use this formula to compute two-time correlation functions of the system and output field properties. Finally, we showed that delayed coherent feedback can be simulated through either an exotic quantum teleportation protocol requiring time travel, or through a probabilistic teleportation protocol.
We conclude with some general remarks on the relation between the techniques presented above and non-Markovian open quantum systems in general. Consider, for example, a single open quantum system that interacts with a feedback reservoir with a generic memory kernel f (t). This memory kernel may be approximated by requiring that integrals over it become left Riemann sums: for some chosen h. Note that Eq. (59) takes the same form as Eq. (35) provided f (t) is Hermitian. The error in this approximation is O(h 2 ), and the exact memory kernel is of course recovered in the limit h → 0. As such, provided we have the computational resources to consider sufficiently small h, any memory kernel -even a continuous one -may be approximated as a series of discrete delayed feedback loops. Because of this, continuous coherent feedback can be viewed as an infinite chain of cascaded system copies, subject once again to the inter-system mapping formula (2). The method owes its conceptual generality to the ability of this mapping rule to insert (or teleport, as shown in Sec. 6) a history into the system's evolution after the fact.
where we have defined