Open quantum systems with delayed coherent feedback

TL;DR

This paper introduces a generalized method for simulating open quantum systems with delayed feedback, including multiple delays and cascaded systems, and extends the quantum regression formula to such systems.

Contribution

It provides a new simulation technique for delayed feedback in quantum systems and generalizes the quantum regression formula for these scenarios.

Findings

Simulated systems with multiple delays and cascaded feedback.

Derived a generalized quantum regression formula for delayed feedback.

Proposed a quantum teleportation protocol for simulating delayed feedback.

Abstract

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.

Equations

$$\rho (t)=\sum _{j_1}\text{tr}\left[\rho (0)e_{j_1}^{†}\right]e_{j_1}(t).$$

$$e_{j_1}(t)={\text{tr}}_1\cdots {\text{tr}}_{n-1}\sum _{j_2\cdots j_n}{e}_{j_1\cdots j_n}(\xi ;t^{\prime })[e_{j_2\cdots j_n}^{†}\otimes I],$$

$$e_{j_1\cdots j_n}(\xi ;t^{\prime })={\Phi }^{(n-1)}(\xi ,t^{\prime }){\Phi }^{(n)}(t^{\prime },0)e_{j_1\cdots j_n}.$$

$$\rho (0)=\sum _{\mu \nu }⟨\mu ⟩\rho (0)\nu \left(|\mu ⟩⟨\nu |\right),$$

$$U(t,0)\left(|\mu ⟩⟨\nu |\right)U^{†}(t,0)$$

$$=U(t,\xi )\left[U(\xi ,0)\left(|\mu ⟩⟨\nu |\right)U^{†}(\xi ,0)\right]U^{†}(t,\xi )$$

$$\displaystyle=U(t,\xi)\Bigl{[}\textstyle\sum_{\mu^{\prime}\nu^{\prime}}\braket{\mu^{\prime}}{U(\xi,0)\left(\ket{\mu}\bra{\nu}\right)U^{\dagger}(\xi,0)}{\nu^{\prime}}\ket{\mu^{\prime}}\bra{\nu^{\prime}}\Bigr{]}U^{\dagger}(t,\xi)$$

$$=\sum _{{\mu }^{\prime }{\nu }^{\prime }}U(t,\xi )\left(|{\mu }^{\prime }⟩⟨{\nu }^{\prime }|\right)U^{†}(t,\xi )\text{tr}\left[U(\xi ,0)\left(|\mu ⟩⟨\nu |\right)U^{†}(\xi ,0){\left(|{\mu }^{\prime }⟩⟨{\nu }^{\prime }|\right)}^{†}\right].$$

$$U(t,0)e_{j_1}{U}^{†}(t,0)=\sum _{j_2}U(t,\xi )e_{j_2}{U}^{†}(t,\xi )\text{tr}\left[U(\xi ,0)e_{j_1}{U}^{†}(\xi ,0)e_{j_2}^{†}\right],$$

$$e_{j_{1}}(t)=\tr_{1}\sum_{j_{2}}\left[\left(U(\xi,0)e_{j_{1}}U^{\dagger}(\xi,0)\right)\otimes\left(U(t,\xi)e_{j_{2}}U^{\dagger}(t,\xi)\right)\right]\bigl{[}e_{j_{2}}^{\dagger}\otimes I\bigr{]}\,.$$

$$e_{j_1\cdots j_n}(\xi ;t^{\prime })=\left[U(\xi ,0)\otimes U(2\xi ,\xi )\otimes \cdots \otimes U(t,(n-1)\xi )\right]e_{j_1\cdots j_n}$$

$$\times \left[U^{†}(\xi ,0)\otimes U^{†}(2\xi ,\xi )\otimes \cdots \otimes U^{†}(t,(n-1)\xi )\right],$$

$$A_m\equiv \underset{(m-1)\mathrm{times}}{\underbrace{I\otimes \cdots \otimes I}}\otimes A\otimes \underset{(n-m)\mathrm{times}}{\underbrace{I\otimes \cdots \otimes I}},$$

$$U^{(m)}(t_1,t_0)\equiv T\exp \left[-i{\int }_{t_0}^{t_1}\mathrm{ds}\sum _{l=1}^m{H}_l(s+(l-1)\xi )\right],$$

$$U^{(n-1)}(\xi ,t^{\prime })U^{(n)}(t^{\prime },0)=U(\xi ,0)\otimes U(2\xi ,\xi )\otimes \cdots \otimes U(t,(n-1)\xi ),$$

$$e_{j_1\cdots j_n}(\xi ;t^{\prime })=U^{(n-1)}(\xi ,t^{\prime })U^{(n)}(t^{\prime },0)e_{j_1\cdots j_n}{U}^{(n)†}(t^{\prime },0)U^{(n-1)†}(\xi ,t^{\prime }).$$

$$H_{SE}(t)=\sum _{\alpha }\left[a_{\alpha }^{†}{B}_{\alpha }(t)+B_{\alpha }^{†}(t)a_{\alpha }\right],$$

$$\rho (0)=\sum _{\mu \nu }⟨\mu ⟩{\rho }_S(0)\nu \left({\rho }_E\otimes |\mu ⟩⟨\nu |\right),$$

$$U(t,0)\left({\rho }_E\otimes |\mu ⟩⟨\nu |\right)U^{†}(t,0)$$

$$=U(t,\xi )\left\{U(\xi ,0)\left({\rho }_E\otimes |\mu ⟩⟨\nu |\right)U^{†}(\xi ,0)\right\}{U}^{†}(t,\xi )$$

$$\displaystyle=\sum_{\mu^{\prime}\nu^{\prime}}U(t,\xi)\left\{\tr_{S}\Bigl{[}U(\xi,0)\left(\rho_{E}\otimes\ket{\mu}\bra{\nu}\right)U^{\dagger}(\xi,0){(\ket{\mu^{\prime}}\bra{\nu^{\prime}})}^{\dagger}\Bigr{]}\ket{\mu^{\prime}}\bra{\nu^{\prime}}\right\}$$

$$\times U^{†}(t,\xi ).$$

$$U(t,0)({\rho }_E\otimes e_{S;j_1})U^{†}(t,0)$$

$$\displaystyle=\sum_{j_{2}}U(t,\xi)\left\{\tr_{S}\Bigl{[}U(\xi,0)(\rho_{E}\otimes e_{S;j_{1}})U^{\dagger}(\xi,0)e_{S;j_{2}}^{\dagger}\Bigr{]}e_{S;j_{2}}\right\}U^{\dagger}(t,\xi)\,.$$

$$\displaystyle U^{\prime}(t_{1},t_{0})\equiv T\exp\Biggl{(}-i\int_{t_{0}}^{t_{1}}\mathop{d{s}}\nolimits\sum_{m}1_{(m-1)\xi\leq s<m\xi}$$

$$\displaystyle\hskip 85.35826pt\times\Biggl{\{}H_{S;m}+\sum_{\alpha}\left[B_{\alpha}(s)\otimes a_{\alpha;m}^{\dagger}+B_{\alpha}^{\dagger}(s)\otimes a_{\alpha;m}\right]\Biggr{\}}\Biggr{)}\,.$$

$$U(t,0)({\rho }_E\otimes e_{S;j_1})U^{†}(t,0)=$$

$${\text{tr}}_{S;1}\cdots {\text{tr}}_{S;n-1}\sum _{j_2\cdots j_n}\left[U^{\prime }(t,0)({\rho }_E\otimes e_{S;j_1\cdots j_n})U^{\prime †}(t,0)\right][I_E\otimes e_{S;j_2\cdots j_n}^{†}\otimes I_S].$$

$${\rho }_S(t)=\sum _{j_1}{\text{tr}}_S\left[\rho (0)e_{S;j_1}^{†}\right]e_{S;j_1}(t),$$

$$e_{S;j_1}(t)={\text{tr}}_{S;1}\cdots {\text{tr}}_{S;n-1}\sum _{j_2\cdots j_n}{e}_{S;j_1\cdots j_n}(\xi ;t^{\prime })[e_{S;j_2\cdots j_n}^{†}\otimes I_S],$$

$$e_{S;j_1\cdots j_n}(\xi ;t^{\prime })\equiv {\text{tr}}_E\left[U^{\prime }(t^{\prime }+(n-1)\xi ,0)({\rho }_E\otimes e_{S;j_1\cdots j_n})U^{\prime †}(t^{\prime }+(n-1)\xi ,0)\right].$$

$$\tilde{H}_{SE;m}(t^{\prime})\equiv\sum_{\alpha}\bigl{[}\tilde{B}_{\alpha;m}(t^{\prime})\otimes a_{\alpha;m}^{\dagger}+{\tilde{B}_{\alpha;m}^{\dagger}}(t^{\prime})\otimes a_{\alpha;m}\bigr{]}\,.$$

$${\tilde{U}}_m(t_1,t_0)\equiv T\exp \left[-i{\int }_{t_0}^{t_1}\mathrm{ds}{\tilde{H}}_m(s)\right].$$

$${\tilde{U}}^{(m)}(t_1,t_0)\equiv T\exp \left[-i{\int }_{t_0}^{t_1}\mathrm{ds}\sum _{l=1}^m{\tilde{H}}_l(s)\right].$$

$${\Phi }^{(m)}(t_1,t_0)\chi \equiv {\text{tr}}_E\left\{{\tilde{U}}^{(m)}(t_1,t_0)[{\rho }_{\tilde{E}}\otimes \chi ]{\tilde{U}}^{(m)†}(t_1,t_0)\right\},$$

Full text

1 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 discuss their application to more general systems in Sec. 7.

2 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}^{\dagger}e_{k}=\delta_{jk}$. The state of the reduced system (with the environment traced out) can be expanded, as usual, in terms of these basis operators:

[TABLE]

The evolution of the basis operators (and hence the system) is divided into intervals of length $\xi$, with each interval represented by a formally separate Hilbert space. To divide the system’s evolution into $n\equiv\lceil t/\xi\rceil$ 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^{\prime}\equiv t-(n-1)\xi$, as illustrated in Fig. 2.

The evolved basis operator $e_{j_{1}}(t)$ is obtained using the mapping formula

[TABLE]

where we have defined the product basis operators $e_{j_{1}\cdots j_{n}}\equiv e_{j_{1}}\otimes\cdots\otimes e_{j_{n}}$, and where $e_{j_{1}\cdots j_{n}}(\xi;t^{\prime})$ is a basis operator of the entire fictitious ‘chain’ of system copies, given by

[TABLE]

At this stage, we simply assume the existence of the dynamical map $\Phi^{(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^{\prime}$, 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 $\xi$. 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 $\Phi^{(m)}(t_{1},t_{0})$ is met in its simplest form by considering a closed quantum system and starting out with more familiar notation.

2.1 A closed quantum system

We consider the evolution of a density operator, $\rho(t)$, with initial condition written in matrix notation,

[TABLE]

where $\braket{\mu}{\nu}=\delta_{\mu\nu}$. The unitary time-evolution operator is denoted $U(t,t_{0})\equiv T\exp[-i\int_{t_{0}}^{t}\mathop{d{s}}\nolimits H(s)]$, 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 [math] to time $t$ into two intervals: an initial interval of length $\xi$ and the remainder. The time evolution is generated by propagating the elementary operators $\left(\ket{\mu}\bra{\nu}\right)$. By interrupting the evolution at $\xi$ and re-expanding the result in matrix notation before proceeding, we get

[TABLE]

Then returning to our more compact notation, we may write

[TABLE]

where we allow a single index to stand in for the pair $(\mu,\nu)$, with $e_{j_{1}}$ – a basis operator as introduced above – standing in for $\left(\ket{\mu}\bra{\nu}\right)$ and $e_{j_{2}}$ for $\left(\ket{\mu^{\prime}}\bra{\nu^{\prime}}\right)$. 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 $\xi$) and space 2 (interval $\xi$ 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:

[TABLE]

The division into $n$ intervals gives an obvious extension of the formula, which turns out to be exactly Eq. (2) with

[TABLE]

where we have used the fact that operators of different system copies commute to group all the unitary operators together. Equation (2.1) evolves the basis operators for each of the system copies $1$ to $(n-1)$ through an interval of length $\xi$, 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^{\prime}$, where we recall that $t\equiv t^{\prime}+(n-1)\xi$. As such, this process reconstructs the evolved basis operator $e_{j_{1}}(t)$ and (by completeness) the state of the ‘real’ system.

In Eq. (2.1), 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

[TABLE]

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:

[TABLE]

where in this context $0\leq t_{0}\leq t_{1}\leq\xi$, as this unitary will only be used to evolve the chain on the interval $[0,\xi]$. It is easily verified that

[TABLE]

which means that the evolved basis operators in Eq. (2.1) can now be written in the form

[TABLE]

It follows that in the simple case of a closed system the map introduced in Eq. (3) is given by $\Phi^{(m)}(t_{1},t_{0})e_{j_{1}\cdots j_{n}}=U^{(m)}(t_{1},t_{0})e_{j_{1}\cdots j_{n}}U^{(m)\dagger}(t_{1},t_{0})$.

2.2 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 $H(t)=H_{S}+H_{SE}(t)$, 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

[TABLE]

where $B_{\alpha}(t)=\sum_{j}\kappa_{\alpha j}e^{i(\omega_{0}-\omega_{j})t}b_{j}$, with couplings $\kappa_{\alpha j}$. Here $\omega_{0}$ is some fiducial frequency that may be freely chosen. The index $\alpha$ 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 $\mathinner{[{b_{j}},{b_{j}^{\dagger}}]}=\delta_{jj^{\prime}}$. It is also useful to define the dissipation (memory) kernel of the reservoir: $F_{\alpha\beta}(t_{2}-t_{1})\equiv\mathinner{[{B_{\alpha}(t_{2})},{B_{\beta}^{\dagger}(t_{1})}]}$. We assume that the initial combined state of the system and environment is separable and given by $\rho(0)=\rho_{E}\otimes\rho_{S}(0)$, where $\rho_{S}(0)$ and $\rho_{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

[TABLE]

where in this case the operators $(\ket{\mu}\bra{\nu})$ form a basis for the system. We once again break the evolution from [math] to $t$ into two intervals, interrupting the evolution at $\xi$ and re-expanding the system—but not the environment:

[TABLE]

Just as we did above in the case of the open quantum system, we replace the pair of indices $(\mu,\nu)$ with a single index by introducing basis operators $\{e_{S;j}\}$ for the system, which gives

[TABLE]

We now introduce a new unitary operator

[TABLE]

where $a_{\alpha;m}$ is an operator of the $m$th system copy, expressed using the notation introduced in Eq. (9), and where $1_{(m-1)\xi\leq s<m\xi}$ is an indicator function, defined as $1$ when $(m-1)\xi\leq s<m\xi$ and [math] otherwise. This unitary allows us to re-write Eq. (2.2) using tensor products. Generalising at the same time to $n$ intervals, we find

[TABLE]

Finally, we trace out the environment to obtain the state of the reduced system, $\rho_{S}(t)\equiv\tr_{E}\rho(t)$:

[TABLE]

where

[TABLE]

and where, recalling that $t\equiv t^{\prime}+(n-1)\xi$, we have defined

[TABLE]

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).

2.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 Hamiltonian $\tilde{H}_{m}(t^{\prime})\equiv H_{S;m}+\tilde{H}_{SE;m}(t^{\prime})$ for a single system copy with interaction term

[TABLE]

We also define the corresponding unitary evolution operator

[TABLE]

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,

[TABLE]

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:

[TABLE]

where $\chi$ is some operator of the reduced system copies, and $\rho_{\tilde{E}}$ 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

[TABLE]

as well as

[TABLE]

with $e_{S;j_{1}\cdots j_{n}}(\xi;t^{\prime})$ given by Eq. (21). As discussed in A, Eq. (2.3) is satisfied whenever the map $\Phi^{(m)}(t_{1},t_{0})$ is divisible [17, 18] for all $m$, meaning that

[TABLE]

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 operators $\tilde{B}_{\alpha;m}(t^{\prime})$ is the stipulation that

[TABLE]

From this commutation relation, we easily find that a similar relation holds for the effective Hamiltonian:

[TABLE]

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 write

[TABLE]

which allows us to re-write Eq. (27) as

[TABLE]

This equation is the open systems version of Eq. (2.1). It remains to choose the environment operator $\tilde{B}_{\alpha;m}(t^{\prime})$ such that Eq. (32) does indeed reproduce Eq. (21). We do not have to specify this operator explicitly: if we restrict attention to vacuum reservoirs – by assuming that both $\rho_{E}$ and $\rho_{\tilde{E}}$ represent the vacuum states of the respective reservoirs – the 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 $\tilde{F}_{\alpha\beta;mm^{\prime}}(t_{2}-t_{1})\equiv\mathinner{[{\tilde{B}_{\alpha;m}(t_{2})},{\tilde{B}_{\beta;m^{\prime}}^{\dagger}(t_{1})}]}$ such that Eqs. (21) and Eq. (32) are equivalent. After a few simple manipulations of these two equations, as detailed in B, we find that this condition is satisfied by the effective dissipation kernel

[TABLE]

with $t_{2}\geq t_{1}$, where the case $m<m^{\prime}$ follows from Eq. (29). Note that the dissipation kernel is, by construction, a Hermitian function: $\tilde{F}_{\alpha\beta;mm^{\prime}}(t_{1}-t_{2})=\tilde{F}_{\beta\alpha;m^{\prime}m}(t_{2}-t_{1})$. In Eq. (33) we focus on the case $t_{2}\geq 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. (2.3) is also satisfied. Combining Eqs. (2.3) and (27), we obtain

[TABLE]

where $e_{S;j_{1}\cdots j_{n}}(\xi;t^{\prime})$ is given by Eq. (21). Eq. (34) is, up to minor notational differences, Eq. (3).

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

[TABLE]

with $\tau_{\alpha\beta j}\geq 0$, which describes exactly such a system. We assume that the delays $\tau_{\alpha\beta j}$ are commensurable, with the interval length $\xi$ chosen to match the greatest common divisor of the delays so that $k_{\alpha\beta j}=\tau_{\alpha\beta j}/\xi$ is an integer for all $\alpha$, $\beta$, and $j$.

Our task now is to derive the map $\Phi^{(m)}(t^{\prime},t_{0})$, $t_{0}\leq t^{\prime}\leq\xi$, describing the reduced system dynamics associated with the dissipation kernel Eq. (35). Substituting the dissipation kernel into Eq. (33) yields

[TABLE]

Recalling that, in the context of Eq. (3), $t^{\prime}\leq\xi$, and noting that both $m-m^{\prime}$ and $k_{\alpha\beta j}$ are non-negative integers, we can see that only those terms in Eq. (36) for which $k_{\alpha\beta j}=m-m^{\prime}$ 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:

[TABLE]

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 (3) 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 (3), we are in a position to write down the map $\Phi^{(m)}(t^{\prime},t_{0})$ that first appeared in Eq. (25):

[TABLE]

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 (3) 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 $\xi$-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 $\xi$. 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_{\alpha\beta j}$ is large. For this reason, in the following examples we have restricted ourselves to situations where all $k_{\alpha\beta 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 (3) 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.

4 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

[TABLE]

where $\phi=\omega_{0}\tau$. The internal system Hamiltonian, describing coherent driving with Rabi frequency $\Omega/2$, is $\Omega(\sigma_{-}+\sigma_{+})$. This system is depicted schematically in Fig. LABEL:sub@fig:delay, and the corresponding cascade of system copies is shown in Fig. LABEL:sub@fig:delay_cascade. Our simulation results are shown in Fig. LABEL:sub@fig:delay_results.

Secondly, we examine a qubit, driven as above, that instead couples to the reservoir at $N$ locations with equal spacing $\tau$. This is described by the coupling constant $\kappa_{j}=\sum_{n=0}^{N-1}\sqrt{\gamma/N}e^{i\omega_{j}n\tau}$; the corresponding dissipation kernel is

[TABLE]

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 LABEL:sub@fig:2delay_results and LABEL:sub@fig:3delay_results show sample simulations with the dissipation kernel (40), with finite $N$. In the limit $N\to\infty$, this dissipation kernel describes the situation depicted in Fig. LABEL:sub@fig:loop, 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 $\tau$ 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:

[TABLE]

The internal system Hamiltonian is given by $H=\Omega(\sigma_{-,A}+\sigma_{+,A})+\Omega(e^{i\phi}\sigma_{-,B}+e^{-i\phi}\sigma_{+,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

[TABLE]

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.

5 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 generalised. Denoting by $\rho_{S+E}(t)$ the combined state of the system and reservoir, we trivially find

[TABLE]

where

[TABLE]

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 $\xi$. As one might intuitively expect, the operators $C$ and $A$ are applied in the $l_{1}$th fictitious system copy at time $t_{1}^{\prime}$, with $t_{n}^{\prime}=t_{n}-(l_{n}-1)\xi$ with $l_{n}=\lceil t_{n}/\xi\rceil$. That is to say, we modify Eq. (3) to read:

[TABLE]

We similarly modify Eq. (2):

[TABLE]

Once again it is a simple exercise to show that Eqs. (5) and (48) reconstruct the operator (46). These equations therefore represent a generalisation of the quantum regression formula to the case of open quantum systems with delayed coherent feedback.

We can use Eqs. (5) and (48) to calculate properties of the environment. As an example, consider the second-order photon correlation function

[TABLE]

$t_{2}\geq t_{1}$, where $E_{\alpha}(t)$ is the field in the $\alpha$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

[TABLE]

which allows us to express $g_{\alpha}^{(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}$.

6 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 $\tau<t\leq 2\tau$, 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}\sum_{\mu}\ket{\mu}\otimes\ket{\mu}$, with $D$ the subsystem dimension, while $S_{1}$ is prepared in the initial system state $\rho(0)$. Thus, the initial state of the combined system can be written as

[TABLE]

An orthonormal basis of maximally entangled states is given by $\{\ket{\psi^{(pq)}}\}$ with

[TABLE]

where $\mu\oplus q\equiv(\mu+q)\bmod 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

[TABLE]

If Bob then applies the unitary

[TABLE]

he is left with a copy of $\rho(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

[TABLE]

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 $\Phi^{(1)}(\tau,t^{\prime})\Phi^{(2)}(t^{\prime},0)$, $t^{\prime}=t-\tau$. The result is of course

[TABLE]

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

[TABLE]

If we now make the same notational change as in Sec. 2.1, with a single index standing in for the pair $(\mu,\nu)$ and $e_{j}$ standing in for $\left(\ket{\mu}\bra{\nu}\right)$, Eq. (57) becomes

[TABLE]

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$.

7 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:

[TABLE]

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\to 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.

We would like to thank the anonymous referee for pointing out an error in an earlier version of the manuscript. This work was supported by the Marsden Fund of the Royal Society of New Zealand.

Appendix A Derivation of Eq. (2.3)

We briefly discuss here the reason why Eq. (2.3) is true whenever $\Phi^{(m)}$ is divisible, for all $m$. We can define an auxiliary map $\Phi_{t^{\prime}}$ that satisfies

[TABLE]

where

[TABLE]

As such, we have

[TABLE]

A sum of divisible maps, with positive real coefficients, is also divisible [17, 28, 29, 30]. Because the indicator functions $1_{\cdots}$ take the non-negative values [math] or $1$, if both $\Phi^{(n)}$ and $\Phi^{(n-1)}$ are divisible – as we have assumed – then $\Phi_{t^{\prime}}$ must be divisible as well. As such, we may write

[TABLE]

Eq. (2.3) then follows immediately from the definitions (61) and (62).

Appendix B Derivation of Eq. (33)

Our aim is to show that the environment operator $\tilde{B}_{\alpha;m}(t^{\prime})$ appearing in Eq. (22) may be chosen such that Eq. (32) agrees with Eq. (21). The unitary (2.2) that appears in Eq. (21) may be re-written as

[TABLE]

where we have defined

[TABLE]

Equation (32) may be written in a form analogous to Eq. (21):

[TABLE]

where we have defined

[TABLE]

with

[TABLE]

Observe that the unitaries (B) and (67) differ only in terms of the environment operators (65) and (68). If we restrict attention to vacuum reservoirs (as we do in Sec. 2.3), the right-hand sides of Eqs. (21) and (66) are equal if the two interactions have equivalent dissipation kernels. To be more specific, we need

[TABLE]

for all $t_{2}\geq t_{1}$, which is satisfied provided

[TABLE]

or equivalently

[TABLE]

for all $t_{2}\geq t_{1}$ and $m\geq m^{\prime}$. Combining Eq. (71) with Eq. (29) immediately gives Eq. (33).