Collision model for non-Markovian quantum trajectories

TL;DR

This paper introduces a collision model-based algorithm to simulate measurement-conditioned quantum trajectories in non-Markovian systems, accommodating different detection schemes and complex environmental memory effects.

Contribution

It provides a novel algorithm for simulating non-Markovian quantum trajectories using collision models, applicable to various detection methods and environmental memory types.

Findings

Successfully simulates trajectories with delayed feedback

Handles systems with continuous environmental memory

Adapts to photodetection and homodyne detection

Abstract

We present an algorithm to simulate genuine, measurement-conditioned quantum trajectories for a class of non-Markovian systems, using a collision model for the environment. We derive two versions of the algorithm, the first corresponding to photodetection and the second to homodyne detection with a finite local oscillator amplitude. We use the algorithm to simulate trajectories for a system with delayed coherent feedback, as well as a system with a continuous memory.

Equations

$$b_k$$

$$⟺B_n$$

$$H=H_S+H_E+H_I.$$

$$H_E=\sum _{k=0}^{N-1}{\omega }_k{b}_k^{†}{b}_k$$

$$H_I=\sum _{k=0}^{N-1}\left({\kappa }_k{a}^{†}{b}_k+{\kappa }_k^{\ast }{b}_k^{†}a\right),$$

$${\kappa }_k=\frac{1}{\sqrt{L}}\sum _{n=0}^{N-1}{\gamma }_n{e}^{i{\omega }_{k}n\Delta t},$$

$$H_I=\frac{1}{\sqrt{\Delta t}}\sum _{n=0}^{N-1}\left({\gamma }_n{a}^{†}{B}_n+{\gamma }_n^{\ast }{B}_n^{†}a\right).$$

$$U(\Delta t)\approx \Delta U_E\Delta U_{SI},$$

$$\Delta U_E^{†}{B}_n\Delta U_E=\frac{1}{\sqrt{N}}\sum _{k=0}^{N-1}{b}_k{e}^{i{\omega }_k(n-1)\Delta t}=B_{n-1},$$

$$|k_{M-1},\dots ,k_0⟩=|k_{M-1}⟩\cdots |k_0⟩,$$

$$\Delta U_E|k_{N-1},\dots ,k_1,k_0⟩=|k_0,k_{N-1},\dots ,k_1⟩.$$

$${\widetilde{\Delta U}}_E|k_{N-1},\dots ,k_1⟩|\phi ⟩=|0,k_{N-1},\dots ,k_1⟩.$$

$${\widetilde{\Delta U}}_E=\sum _{k_{N-1},\dots ,k_1}|0,k_{N-1},\dots ,k_1⟩⟨k_{N-1},\dots ,k_1|⟨\phi |.$$

$$H_I=\frac{1}{\sqrt{\Delta t}}\left(\gamma a^{†}B+{\gamma }^{\ast }{B}^{†}a\right).$$

$$|\psi (\Delta t)⟩\approx \left(1-\frac{\Delta t}{2}|\gamma {|}^2{a}^{†}a\right)|{\psi }_S⟩|0⟩-i\sqrt{\Delta t}{\gamma }^{\ast }a|{\psi }_S⟩|1⟩,$$

$$N_{\pm }=\frac{1}{2}{\left(B\pm C\right)}^{†}\left(B\pm C\right),$$

$$N_{+}|n\pm ⟩={\lambda }_{n\pm }|n\pm ⟩\text{and}{N}_{-}|n\pm ⟩={\lambda }_{n\mp }|n\pm ⟩,$$

$$\begin{gathered} |\psi (\Delta t)⟩\approx \left(1-\frac{\Delta t}{2}|\gamma {|}^2{a}^{†}a\right)|{\psi }_S⟩|0,\alpha \sqrt{\Delta t}⟩ \\ -i\sqrt{\Delta t}{\gamma }^{\ast }a|{\psi }_S⟩|1,\alpha \sqrt{\Delta t}⟩. \end{gathered}$$

$$⟨0,0⟩\psi (\Delta t)$$

$$⟨1\pm ⟩\psi (\Delta t)$$

$$J_{\pm }=\frac{\alpha \mp i{\gamma }^{\ast }a}{\sqrt{2}}.$$

$$P_0$$

$$P_{\pm }$$

$${\gamma }_n=\sqrt{\gamma }\lambda \Delta t{e}^{-\lambda n\Delta t},$$

Full text

I Introduction

A quantum trajectory is a sequence of states of an open quantum system, conditioned on a sequence of measurements of the system’s output. Markovian quantum trajectory theory describes “unravelings” of a Lindblad master equation: each unraveling is a decomposition of the system density matrix corresponding to a specific measurement set-up Carmichael (1993). This approach leads to a family of efficient Monte Carlo algorithms that can be used to solve the dynamics of Markovian open quantum systems. It has been shown that pure-state quantum trajectories for non-Markovian systems do not exist in general Wiseman and Gambetta (2008), but this does not preclude the existence of mixed-state trajectories for such systems. While various Monte Carlo methods for analyzing non-Markovian systems have been proposed Strunz (1996); Diósi and Strunz (1997); Diósi et al. (1998); Strunz et al. (1999); Breuer (2004); Piilo et al. (2008); Megier et al. (2018), there remains—in contrast to the Markovian case—no generally applicable computational tool that produces genuine measurement-conditioned trajectories for systems with environmental memories. One way around this is to model part of the environment along with the system, as any non-Markovian system can be rendered Markovian by simply enlarging the system. This is the approach we follow in this paper.

Collision models, in which the dynamics of open quantum systems are modeled through repeated unitary interactions between the system and constituent subsystems of the environment Campbell et al. (2018), are emerging as a general tool for the understanding of non-Markovian systems Kretschmer et al. (2016); Ciccarello (2017). In recent years, models of this kind have been applied to the problem of coherent feedback—a signature example of an environmental memory—by Grimsmo (2015) and separately by Pichler and Zoller (2016). Collision models are particularly convenient as a computational tool: by simulating a portion of the environment along with the system, these models are able to treat memory effects while avoiding some of the challenges associated with analytical approaches. A connection between models of this kind and quantum trajectory theory was pointed out by Brun (2002), who used a collision model (although this terminology was not in use at the time) as a simple illustration of the essential physics of Markovian, pure-state quantum trajectories. Quantum trajectory theory has also been used in conjunction with collision models by Cresser to analyze the micromaser Cresser (2006, 2019), and by Dąbrowska and coauthors to model systems where non-Markovianity arises due to correlations in the input field Dąbrowska et al. (2017, 2019); Dąbrowska (2019). In this paper we show how to use a collision model to simulate quantum trajectories of a class of non-Markovian open systems with environmental memories. The algorithm we derive provides a fairly efficient and general method of simulating the dynamics of these systems.

In Sec. II we introduce a collision model for an open quantum system, which in general can be non-Markovian. In Sec. III we describe a quantum trajectory algorithm within this model, focusing on trajectories for two different detection schemes: photodetection (Sec. III.1) and homodyne detection (Sec. III.2). In Sec. IV we apply this algorithm to a specific example of a non-Markovian quantum system, namely a system with delayed coherent feedback, and also illustrate how the approach can be generalized to systems with continuous memories. We conclude in Sec. V. Natural units in which $\hbar=c=1$ are used throughout.

II Model

Our starting point is a simple, prototypical open quantum system that interacts with a one-dimensional, freely propagating electromagnetic field. This set-up has been studied elsewhere Grimsmo (2015); Pichler and Zoller (2016), but we recapitulate the basic model here. The system is allowed to couple to the environment at multiple spatial locations, which in general can make the system’s dynamics non-Markovian, as we discuss in detail in Sec. IV. Sampling the environment at discrete intervals leads to a collision model that approximately describes the evolution of the combined system.

We model the environment on an interval $[-L,0]$. We denote the frequency domain annihilation operators for an environment mode by $b_{k}$, while the corresponding time-domain operators are denoted by $B_{n}$. We sample the time domain at equal intervals $\Delta t=L/N$, so that these operators are related by the discrete Fourier transform:

[TABLE]

where $\omega_{k}=2\pi k/L$. Because we are using natural units, the time domain operator $B_{n}$ can also be thought of as representing the environment field between positions $-n\Delta t$ and $-(n+1)\Delta t$.

We now introduce a system with creation and annihilation operators $a^{\dagger}$ and $a$. The Hamiltonian for the combined system can be split up as

[TABLE]

Here $H_{S}$ describes the internal dynamics of the system, while

[TABLE]

generates the free evolution of the environment. The interaction between system and environment is described, in the rotating-wave approximation, by

[TABLE]

where

[TABLE]

with $\gamma_{n}$ being a position-dependent coupling strength. Substituting $\kappa_{k}$ from Eq. (5) into Eq. (4) and reexpressing in terms of the time domain operators $B_{n}$ gives

[TABLE]

We use this form of the interaction Hamiltonian in what follows.

The evolution of the system is described by the unitary operator $U(t)=e^{-iHt}$. For small $\Delta t$ we can make the approximation

[TABLE]

where $\Delta U_{E}=e^{-iH_{E}\Delta t}$ and $\Delta U_{SI}=e^{-i(H_{S}+H_{I})\Delta t}$. Equation (7) leads to an algorithm for approximating the evolution of the system and environment. First, the state is evolved with $H_{S}+H_{I}$ through a time $\Delta t$, which is equivalent to applying $\Delta U_{SI}$. The free evolution of the environment is then accounted for by applying $\Delta U_{E}$. Using $e^{iH_{E}t}b_{k}e^{-iH_{E}t}=b_{k}e^{-i\omega_{k}t}$, we can see that

[TABLE]

showing that excitations in the environment propagate in the direction of decreasing $n$.

To make the effect of $\Delta U_{E}$ more explicit, we fix a basis for the environment by defining the number states

[TABLE]

where $0<M\leq N$, and where $\ket{k_{n}}$ is a number state of a single environment oscillator. We can express any state of the environment, which we recall comprises $N$ oscillators, as a linear combination of the states $\ket{k_{N-1},\ldots,k_{0}}$. Equation (8) then leads to

[TABLE]

We thus arrive at a collision model describing our system. In this picture, the system interacts stroboscopically with an environment comprising a chain of harmonic oscillators whose annihilation operators are $B_{n}$. In the first stage of the evolution, the system interacts with these environment oscillators according to the Hamiltonian $H_{S}+H_{I}$ for a time interval $\Delta t$. In the second stage, the unitary operator $\Delta U_{E}$ is applied; as Eq. (10) shows, this shifts any excitations in the $n$th oscillator into the $(n-1)$st for $n>0$, while excitations in the zeroth oscillator end up in the $(N-1)$st. These stages are iterated to produce $\Delta t$-periodic samples of the system state. Two time steps of this evolution are illustrated, for a Markovian system, in Fig. 1.

III Trajectory algorithm

Due to the periodic boundary condition implicit in Eq. (II), the model presented in the previous section will eventually begin to display unphysical behavior, as excitations in the zeroth oscillator “loop around” to the $(N-1)$st. We can avoid this by simulating a measurement of the zeroth environment oscillator to disentangle this oscillator from the system and the rest of the environment. The state of the zeroth oscillator can then be ignored in future time steps, thereby creating an absorbing boundary condition at the zeroth oscillator and allowing the simulation to continue indefinitely.

There are two ways in which this measurement can be represented. We can take the partial trace over the zeroth subsystem, representing a probability-weighted sum over the possible measurement outcomes, which would yield a density-matrix representation of the state of the system and environment. Alternatively, we can perform a Monte Carlo simulation of the application of the Born rule. Here we focus on the latter approach, which leads to a quantum trajectory unraveling of the system–environment combined state. Trajectories for the system alone can then be obtained by tracing out the environment, which will in general be a mixed-state trajectory. These trajectories for the system alone provide an unraveling of the system density matrix, in the same sense as the conventional algorithm does for Markovian systems, and are contextual: the specific decomposition of the density matrix obtained from the algorithm depends on the measurement set-up being simulated.

Our algorithm, illustrated in Fig. 2, involves three steps, which are iterated to obtain a discrete-time quantum trajectory for the combined system:

Apply $\Delta U_{SI}$ by evolving the state of the combined system with the Hamiltonian $(H_{S}+H_{I})$ through a time $\Delta t$. 2. 2.

Make a simulated measurement of some observable of the zeroth environment oscillator. 3. 3.

Apply a truncated form of $\Delta U_{E}$ to the postmeasurement state.

Step 1 is easily performed by using any standard differential equation solver. We used ZVODE from ODEPACK Hindmarsh (1983), interfaced through SciPy Jones et al. (01). We consider steps 2 and 3 in more detail below.

Step 2 of our algorithm is to make a simulated measurement of some observable of the zeroth environment oscillator using the Born rule. The specific observable to be measured determines what kind of quantum trajectory we obtain. Two examples are considered in Secs. III.1 and III.2, but for the time being we keep things generic. Consider an observable of the zeroth oscillator with discrete eigenstates $\{\ket{q_{j}}\}$, corresponding to eigenvalues $\{q_{j}\}$. The probability of each measurement result is easily calculated as $P(q_{j})=\lvert\braket{q_{j}}{\psi}\rvert^{2}$, where $\ket{\psi}$ is the state of the system and environment. A standard weighted pseudorandom choice is used to determine the measurement outcome at each step, and the zeroth oscillator is then projected into some fiducial state $\ket{\varphi}$: $\ket{\psi}\to\ket{\varphi}\braket{q_{j}}{\psi}$. The choice of fiducial state is influenced by our chosen basis, which we discuss below.

Step 3 of our algorithm is to apply a truncated form of $\Delta U_{E}$, denoted $\widetilde{\Delta U}_{E}$, which is defined analogously to Eq. (10). The states $\ket{k_{N-1},\ldots,k_{1}}\ket{\varphi}$ form a basis for the system after it has undergone the simulated measurement process described above, with the zeroth oscillator projected into the fiducial state $\ket{\varphi}$. We therefore define $\widetilde{\Delta U}_{E}$ by its action on these basis states:

[TABLE]

The effect of applying $\widetilde{\Delta U}_{E}$ is therefore to shift excitations in the $n$th oscillator into the $(n-1)$st, while “resetting” the $(N-1)$st to its vacuum state. We can easily construct a suitable $\widetilde{\Delta U}_{E}$:

[TABLE]

Alternatively, if $\ket{\varphi}$ is chosen to be the vacuum state, the unitary $\Delta U_{E}$ in fact already satisfies Eq. (11).

Finally, we need to consider the truncation of the environment Hilbert space. Our environment comprises $N$ harmonic oscillators, and we have chosen to represent these in a number-state basis. Our chosen truncation is to constrain the values of $k_{n}$ that appear in Eq. (9). In particular, we choose $k_{n}\leq 1$, so that each environment oscillator effectively becomes a qubit, and require the total number of excitations in the environment not to exceed some number: $\sum_{n}k_{n}\leq K_{\text{max}}$. As such, the basis states of the environment are all the combinations of $0,1,\ldots,K_{\text{max}}$ excitations. The simulations presented in this paper were performed with $K_{\text{max}}=2$ except where otherwise specified.

We have yet to specify exactly what is measured in step 2 of our algorithm. In the following two subsections, we consider two different measurement schemes, and show in each case that our collision model reduces to the corresponding Markovian quantum trajectory algorithm in the appropriate limit.

III.1 Photodetection trajectories

Perhaps the simplest measurement scheme is to simulate detection of the excitation number in the zeroth oscillator. We show here that for a Markovian system this corresponds to the well-known quantum trajectory theory of direct photodetection.

For a Markovian system it is sufficient to consider a single environment oscillator, which means that the interaction Hamiltonian (6) simplifies to

[TABLE]

As we have chosen to truncate the environment Hilbert space such that each oscillator contains at most one excitation, $B$ is the lowering operator for a qubit. Coherent evolution through $\Delta t$ followed by a measurement performed on the environment qubit represents a single time step of our collision model. Suppose the initial state of the combined system is $\ket{\psi(0)}=\ket{\psi_{S}}\ket{0}$. For small $\Delta t$, we can make the approximation

[TABLE]

where we have neglected any internal dynamics of the system for simplicity. The excitation of the environment qubit $B^{\dagger}B$ is then measured: the result $1$ is identified with a “click” in a photodetector, while the result [math] is identified with no detection. From Eq. (14), we can immediately see that the probability of such a click is $\lvert\gamma\rvert^{2}\braket{\psi_{S}}{a^{\dagger}a}{\psi_{S}}\Delta t$, and that the normalized system state after such a detection will be $a\ket{\psi_{S}}/\sqrt{\braket{\psi_{S}}{a^{\dagger}a}{\psi_{S}}}$, up to a phase factor. The state conditioned on no detection is $\left(1-\frac{\Delta t}{2}\lvert\gamma\rvert^{2}a^{\dagger}a\right)\ket{\psi_{S}}$, ensuring conservation of probability to first order in $\Delta t$. These are exactly the results of the conventional quantum trajectory theory of direct photodetection Carmichael (1993); Dalibard et al. (1992).

III.2 Homodyne trajectories

We now wish to extend our trajectory treatment to homodyne detection. In this case, the system and environment are augmented with an ancillary system, the local oscillator. This is a harmonic oscillator prepared, at each time step, in a coherent state $\ket{\alpha\sqrt{\Delta t}}$. The Hamiltonian describing the coherent evolution of the combined system, in the Markovian case, is again Eq. (13). Note that the local oscillator does not evolve under this Hamiltonian because it is already prepared in the desired state. Coherent evolution through $\Delta t$, followed by a joint measurement on the zeroth environment qubit and local oscillator, represents a single step of the stroboscopic evolution of a collision model.

Balanced homodyne detection consists in simultaneously measuring two environment observables, $N_{+}$ and $N_{-}$, where

[TABLE]

and where $B$ and $C$ are annihilation operators for the zeroth environment qubit and the local oscillator respectively. The eigenstates of these operators are $\ket{n\pm}=(\ket{0,n}\pm\ket{1,n-1})/\sqrt{2}$ with $n\geq 1$, as well as the vacuum. We find

[TABLE]

where $\lambda_{\pm}=(n\pm\sqrt{n})/2$. If $\ket{\alpha}$ is a coherent state and $\ket{\varphi}$ is some generic state of the qubit, we have the inner product $\braket{n\pm}{\varphi,\alpha}=e^{-\lvert\alpha\rvert^{2}/2}\alpha^{n-1}(\alpha\braket{0}{\varphi}\pm\sqrt{n}\braket{1}{\varphi})/\sqrt{2n!}$.

We now aim to show that this measurement scheme is equivalent to homodyne detection, in the relevant parameter regime. The argument proceeds much as it did for photodetection. The initial state is now $\ket{\psi(0)}=\ket{\psi_{S}}\ket{0,\alpha\sqrt{\Delta t}}$, meaning that Eq. (14) becomes

[TABLE]

Homodyne detection corresponds to a measurement in the $\ket{n\pm}$ basis, and provided $\Delta t$ is small enough the probability of observing states with $n>1$ will be negligible. We note that $N_{\pm}\ket{1\pm}=\ket{1\pm}$, while $N_{\pm}\ket{1\mp}=0$. As such, we have three eigenstates to consider: $\ket{0,0}$, the vacuum state of the environment, corresponding to no detection; $\ket{1+}$, corresponding to a click in the $N_{+}$ detector; and $\ket{1-}$, corresponding to a click in the $N_{-}$ detector. To first order in $\Delta t$, the corresponding inner products are

[TABLE]

where we have defined

[TABLE]

These operators are, up to a phase convention, the jump operators that appear in the conventional trajectory algorithm for homodyne detection with finite local oscillator amplitude Carmichael (2008). Finally, we arrive at the probabilities

[TABLE]

We note in particular that Eq. (22) is exactly what we would expect from the conventional trajectory algorithm. From this, we can conclude that the measurement scheme described above does indeed correspond to a trajectory simulation of homodyne detection, in the limit of small $\Delta t$.

A further simplification is possible: in a balanced homodyne detection scheme, the results of the $N_{-}$ measurement are subtracted from the results of the $N_{+}$ measurement. As such we can, equivalently, measure the operator $Q=N_{+}-N_{-}=C^{\dagger}B+B^{\dagger}C$. We can see that $Q\ket{n\pm}=\pm\sqrt{n}\ket{n\pm}$, along with, of course, $Q\ket{0,0}=0\ket{0,0}$. Provided we only see eigenvalues with $n\leq 1$ we can interpret this in the same way as the conventional homodyne algorithm: we have three possible measurement outcomes ($\pm 1$ or [math]) representing a click in each of the two detectors, or neither. This way of expressing the problem is computationally convenient because we only need to simulate measurement in a single output channel, where previously we had two.

For the equivalence between our collision model and a conventional trajectory simulation of homodyne detection to hold, we need to choose a small enough time step that there is at most one click in each interval. For this to be the case, we must decrease $\Delta t$ quadratically as we increase $\alpha$, such that $\alpha\sqrt{\Delta t}$ remains below some threshold. This is not computationally feasible, so for larger values of $\alpha$ we inevitably end up performing simulations with large enough $\Delta t$ that we observe the eigenvalues corresponding to $n\geq 2$, which will in general not be integers. Nonetheless, numerical comparisons to the conventional jump algorithm, such as those illustrated in Fig. 3, indicate that the correspondence remains approximately valid for a range of $\alpha$.

IV Delayed coherent feedback

In Sec. III we showed that our algorithm reproduces the results of conventional quantum trajectory theory for Markovian systems, when applied to both photodetection and homodyne detection. We now turn our attention to non-Markovian systems. The parameters $\gamma_{n}$ that appear in Eq. (6) describe position-dependent coupling of the system to the environment. Coupling at a single point, as considered in Sec. III, leads to a Markovian open quantum system. If, on the other hand, the system couples to the environment at more than one location, the system can create excitations in the environment that interact again with the system at a later time. Put another way, the environment “remembers” the state of the system and feeds this information back coherently after some time delay. Here we consider perhaps the simplest example of such an environmental memory, wherein the system couples to the environment at exactly two spatial locations, creating a coherent environmental feedback loop with a discrete propagation delay. This kind of feedback has previously been studied in work on “atomic” emission in front of a mirror Dorner and Zoller (2002); Carmele et al. (2013); Tufarelli et al. (2013, 2014) and in solid-state systems with significant propagation delays Guo et al. (2017). Figure 4 illustrates a collision model for this set-up.

To make things more concrete, consider the interaction Hamiltonian (6), with $\gamma_{0}=\gamma e^{i\phi}$, $\gamma_{M}=\gamma$, and $\gamma_{n}=0$ for all other $n$. This coupling describes an environmental feedback loop of length $\tau=M\Delta t$, with a phase advance of $\phi$ in the loop. We choose as an example system a driven qubit, with free Hamiltonian $H_{S}=\Omega(a^{\dagger}+a)$ in a frame rotating at the system resonant frequency, and the system initially in its excited state. Sample trajectories, along with ensemble averages, calculated by using our collision model for this system with both photodetection and homodyne detection, are shown in Fig. 5. Without observing the in-loop field, we cannot tell whether a detected photon was emitted directly from the system or via the loop. This ambiguity means that trajectories for this system are, in general, mixed-state trajectories. Note, however, that this is not a property of coherent feedback as such: a detection delay results in mixed-state trajectories even for a system with Markovian dynamics.

We thus see that our algorithm is particularly well suited for simulating non-Markovian systems. The main computational advantage is that the algorithm is constant-space, requiring only enough computer memory to store the state of the system along with enough of the environment to model the feedback loop in question. This is in contrast to the approach introduced by Grimsmo (2015) and explored further in Ref. Whalen et al. (2017), which requires an additional “copy” of the system Hilbert space after each delay cycle and thus consumes an amount of memory that grows exponentially with the number of delay intervals to be simulated. The approach adopted by Pichler and Zoller (2016), on the other hand, also has constant space complexity. Our algorithm has the further advantage that it allows us to simulate not just discrete feedback loops but other kinds of quantum memory as well. A continuous environmental memory kernel, where the evolution of the system depends most generally on its state at all previous times, may be approximated as a series of discrete feedback loops. This can be accomplished in our algorithm by choosing appropriate $\gamma_{n}$. Take as an example the exponential coupling

[TABLE]

with $\lambda>0$. Substituting this coupling into Eq. (5) and taking the continuum limit ($L\to\infty$) results in a Lorentzian spectral density $J(\omega)={(2\pi)}^{-1}\gamma\lambda^{2}/(\lambda^{2}+\omega^{2})$. This spectral density maps to a system coupled to a damped harmonic oscillator initially in the vacuum state Garraway (1997); Lorenzo et al. (2017); if the system is a qubit then this corresponds to the Jaynes–Cummings model Jaynes and Cummings (1963). As a proof of principle Fig. 6 shows sample trajectories obtained with an exponential coupling and the environment truncated at a finite “length”, and compares the ensemble average to results obtained using the Markovian master equation.

V Conclusion

In this paper we introduced an algorithm to simulate quantum trajectories for non-Markovian systems, by using a collision model to represent the environment and its interaction with the system. The algorithm produces trajectory unravelings of the system density matrix. As is the case in the well-known Markovian theory of quantum trajectories, these unravelings are contextual in the sense that they depend on the measurement set-up. We provided two versions of the algorithm, corresponding to photodetection and homodyne detection respectively, and illustrated the algorithm’s application to two signature examples of non-Markovian dynamics: a coherent feedback loop with a discrete delay, and the Lorentzian spectral density that arises in the Jaynes–Cummings model.

Our simulation method has computational advantages over some alternative approaches. In particular, while the algorithm requires us to simulate a portion of the environment—the “memory”—in addition to the system, the Hilbert-space dimension of this environmental memory remains constant irrespective of the simulation time, meaning that the algorithm has constant space complexity. The algorithm acts on a pure state of the system and environment, which consumes less computer memory than a density-matrix representation—an advantage our approach shares with stochastic approaches derived from Markovian quantum trajectory theory.

Because our approach is derived directly from a fundamental model of a system interacting with an assemblage of harmonic oscillators, it has in principle very wide applicability when used as a brute force method, given enough computing power. An efficient implementation, however, requires a sufficiently compact representation of the part of the environment containing the memory. The systems studied in this paper have finite-length memories that are only ever populated with one or two photons. Systems with very long memories, or systems whose timescales require a very small step size $\Delta t$, may pose computational difficulties for this method because the number of environment subsystems would be large. The same is true of environments with more than one spatial dimension. Performance issues could also arise with systems that scatter many photons into the environment, necessitating a large $K_{\text{max}}$. In these cases, alternative methods of “compressing” the environment state may be required. One possibility is to adapt a matrix product state approach, as used for example by Pichler and Zoller (2016), for use in quantum trajectory simulations.

Quantum trajectory theory is of interest beyond its use as a numerical tool. Genuine quantum trajectories provide an accurate description of the evolution of a system, conditioned on a sequence of observations of the system’s output. The contextuality of a specific trajectory unraveling, depending as it does on the measurement scheme chosen by the observer, has been described as “subjective reality” in the context of quantum measurement theory Wiseman (1996). The algorithm presented in this paper generates genuine, measurement-conditioned quantum trajectories for a fairly large class of non-Markovian open quantum systems, namely those where the environment can be represented as a collision model with a position-dependent coupling. This connection to the theory of quantum measurement opens up the possibility of analyzing non-Markovian systems from a new perspective.

Acknowledgements.

The simulations presented in this paper were performed with the aid of the Python library QuTiP Johansson et al. (2013). The author would like to acknowledge Howard Carmichael for his guidance and many fruitful suggestions, and, through him, the support of the Dodd-Walls Centre for Photonics and Quantum Technologies. The author also gratefully acknowledges valuable discussions with Jim Cresser.