# Quantum-jump vs stochastic Schrödinger dynamics for Gaussian states with
quadratic Hamiltonians and linear Lindbladians

Robson Christie∗, Jessica Eastman∗, Roman Schubert∗∗ and Eva-Maria Graefe∗
∗Department of Mathematics, Imperial College London, London SW7 2AZ, United
Kingdom  ∗∗School of Mathematics, University of Bristol, Bristol, BS8 1TW,
United Kingdom<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>
###### Abstract

The dynamics of Gaussian states for open quantum systems described by Lindblad
equations can be solved analytically for systems with quadratic Hamiltonians
and linear Lindbladians, showing the familiar phenomena of dissipation and
decoherence. It is well known that the Lindblad dynamics can be expressed as
an ensemble average over stochastic pure-state dynamics, which can be
interpreted as individual experimental implementations, where the form of the
stochastic dynamics depends on the measurement setup. Here we consider
quantum-jump and stochastic Schrödinger dynamics for initially Gaussian
states. While both unravellings converge to the same Lindblad dynamics when
averaged, the individual dynamics can differ qualitatively. For the stochastic
Schrödinger equation, Gaussian states remain Gaussian during the evolution,
with stochastic differential equations governing the evolution of the phase-
space centre and a deterministic evolution of the covariance matrix. In
contrast to this, individual pure-state dynamics arising from the quantum-jump
evolution do not remain Gaussian in general. Applying results developed in the
non-Hermitian context for Hagedorn wavepackets, we formulate a method to
generate quantum-jump trajectories that is described entirely in terms of the
evolution of an underlying Gaussian state. To illustrate the behaviours of the
different unravellings in comparison to the Lindblad dynamics, we consider two
examples in detail, which can be largely treated analytically, a harmonic
oscillator subject to position measurement and a damped harmonic oscillator.
In both cases, we highlight the differences as well as the similarities of the
stochastic Schrödinger and the quantum-jump dynamics.

  *

##  1 Introduction

A large class of open quantum systems can be described by Lindblad equations
for the dynamics of density operators. In many atomic physics and quantum
optics applications, this is an accurate description of physical phenomena.
Lindblad dynamics can be viewed as an ensemble average over stochastic pure-
state trajectories. The stochastic trajectories may be interpreted as the
state of an individual experimental quantum system conditioned on a
measurement record. The nature of the measurement record determines the types
of trajectories in the ensemble [1, 2, 3], and intriguingly can lead to
strikingly different microscopic dynamics and behaviours, while averaging to
the same Lindblad dynamics. An example of this has been reported in [4] where
it has been demonstrated that depending on the specific unravelling a
dissipative Kerr dynamics can lead to switches between either odd and even
Schrödinger cat states or between coherent states of opposite phase. It has
been observed in examples that individual stochastic trajectories can even
display chaotic dynamics that are not seen in the ensemble average [5].

Two of the most important unravellings are given by stochastic Schrödinger
equations (SSEs) and quantum-jump trajectories, respectively. The SSE dynamics
is a good description of certain homodyne and heterodyne detection schemes in
quantum optics arising from constant weak continuous measurements [6, 7].
Quantum-jump dynamics, on the other hand, may arise in photodetection
experiments and is theoretically described by periods of deterministic
evolution under an effective non-Hermitian Hamiltonian, stochastically
interrupted by discrete measurements [8]. In both SSE and quantum-jump
scenarios, if the measurement channel is not recorded, then the best estimate
for the state is obtained by averaging over all possibilities, resulting in
Lindblad dynamics [6]. There are other interpretations of the unravellings of
the Lindblad equation, viewing them as candidate laws of nature in an effort
to explain the wavefunction collapse in quantum measurement (see e.g. [9, 10,
11, 12] for more details).

In the present paper, we analyse in detail the structural differences of these
unravellings as compared to the full Lindblad dynamics for the case of
quadratic Hamiltonians and linear Lindbladians for initial Gaussian states. In
these cases, the full quantum dynamics in the Lindblad case can be reduced to
simple phase-space dynamics, as has been discussed in [13, 14]. This follows
the spirit of Heller’s and Littlejohn’s approach to closed system dynamics
[15, 16] where for a quadratic Hamiltonian initially Gaussian states remain
Gaussian under time evolution. This idea carries through to the SSE [17, 18,
7]. Since Gaussian states may be parameterised entirely by their centres and
covariances, the dynamics of these systems may be transformed from state
vector dynamics in infinite-dimensional Hilbert space to a handful of
differential equations for classical phase-space observables. This simple
approach does not carry through to the quantum-jump dynamics, which in general
do not preserve Gaussian states for arbitrary linear Lindbladians. Here we
present a different approach that builds on an extension of Hagedorn’s
wavepacket dynamics [19] that has been adapted to non-Hermitian systems in
[20].

In detail, the paper is organised as follows. In section 2 we introduce the
Lindblad equation along with the quantum-jump and SSE unravellings. In section
3 we insert an initially Gaussian state into the SSE and derive a set of
stochastic differential equations which is exact in the case of quadratic
Hamiltonians and linear Lindbladians. The resulting parameter dynamics are
less computationally expensive to simulate than the full state vector
dynamics. In section 4 we briefly review the non-Hermitian dynamics of
Gaussian states and then apply the non-Hermitian Hagedorn wavepacket dynamics
developed in [20] to formulate several semi-analytical algorithms for
simulating quantum-jump trajectories. In section 5 we illustrate our findings
for two examples which we consider in some detail: A harmonic oscillator
Hamiltonian with either a position measurement, modelled by a Lindbladian
$\hat{L}=\hat{x}$, or a radiative decay, modelled by a Lindbladian
$\hat{L}=\hat{a}$. We conclude with a summary and outlook in section 6.

##  2 Unravellings of the Lindblad Equation

The Lindblad equation was initially derived as the most general Markovian
dynamical equation that preserves the trace, Hermiticity and positivity of the
density matrix [21, 22]. Physically, it can be used to describe certain
quantum systems that are weakly coupled to a memoryless environment. In this
spirit, dynamics of Lindblad form can be obtained for the reduced density
matrix by averaging over the effect of a bath of quantum harmonic oscillators
[23]. Any linear and Markovian (local in time) master equation that preserves
the Hermiticity and trace of the density matrix may be expressed in Lindblad
form,

$i\hbar\frac{d}{dt}\hat{\rho}(t)=[\hat{H},\hat{\rho}(t)]+i\sum_{k}[\hat{L}_{k}\hat{\rho}(t)\hat{L}_{k}^{\dagger}-\frac{1}{2}\hat{L}_{k}^{\dagger}\hat{L}_{k}\hat{\rho}(t)-\frac{1}{2}\hat{\rho}(t)\hat{L}_{k}^{\dagger}\hat{L}_{k}].$
(1)

Here $\hat{H}$ is a Hamiltonian and ${\hat{L}_{k}}$ are general Lindblad
operators, the properties of which are system-specific. For simplicity in what
follows we shall confine the discussion to a single Lindblad operator.

Much like the deterministic Fokker-Planck equation for the dynamics of
probability distributions in classical physics admits unravellings in terms of
single trajectories of the stochastic Langevin equation, the deterministic
Lindblad equation for the dynamics of the density operator may be unravelled
in terms of stochastic pure-state trajectories. There are infinitely many such
unravellings that differ from each other in the stochastic driving processes.
The two types of unravellings most commonly considered in the literature are
SSE trajectories [12] driven with continuous Gaussian distributed noise, and
quantum-jump [8] trajectories driven by discrete Poissonian distributed noise.

The SSE we consider here is given by [12, 17]

$\displaystyle\ket{d\psi}$
$\displaystyle=\frac{1}{\hbar}(-{i}\hat{H}-\frac{1}{2}\hat{L}^{\dagger}\hat{L}+\braket{\hat{L}^{\dagger}}\hat{L}-\frac{1}{2}\braket{\hat{L}^{\dagger}}\braket{\hat{L}})\ket{\psi}dt$
(2)
$\displaystyle\phantom{=}+\frac{1}{\sqrt{2\hbar}}(\hat{L}-\braket{\hat{L}})\ket{\psi}(d\xi_{R}+id\xi_{I}),$

where $d\xi_{R}$ and $d\xi_{I}$ are independent
($\mathbb{E}[d\xi_{R}d\xi_{I}]=0$) Itô stochastic processes with mean zero
($\mathbb{E}[d\xi_{R}]=\mathbb{E}[d\xi_{I}]=0$) and normalisation
$d\xi_{R}^{2}=d\xi_{I}^{2}=dt$. The SSE trajectories are driven with a
continuous stochastic process and are often used to model systems undergoing
weak continuous measurement such as heterodyne detection in quantum optics [6]
or quantum Brownian motion [23].

In the quantum-jump description, on the other hand, the pure-state
trajectories deterministically evolve under an effective non-Hermitian
Hamiltonian $\hat{H}-i\hat{\Gamma}$ with
$\hat{\Gamma}=\frac{1}{2}\hat{L}^{\dagger}\hat{L}$, periodically interrupted
by stochastic quantum jumps. These jumps may be used to represent random
discrete measurements of quantum systems such as photodetection from a
microwave cavity [6]. The cumulative effect of these jumps when averaged over
many trajectories induces the “jump term” contribution
$\hat{L}\hat{\rho}(t)\hat{L}^{\dagger}$ in the density operator dynamics.
Concretely, quantum-jump pure-state dynamics can be described by the dynamical
equation

$\displaystyle\ket{d\psi}$
$\displaystyle=\frac{1}{\hbar}\left(-{i}\hat{H}-\frac{1}{2}\hat{L}^{\dagger}\hat{L}-\frac{1}{2}\braket{\hat{L}^{\dagger}\hat{L}}\right)\ket{\psi}(1-dN)dt$
(3)
$\displaystyle\phantom{=}+\left(\frac{\hat{L}}{\sqrt{\braket{\hat{L}^{\dagger}\hat{L}}}}-1\right)\ket{\psi}dN.$

Here $dN$ is a Poisson process, taking the values 0 (no jump) or 1 (jump) with
expectation value $\mathbb{E}[dN]=\braket{\hat{L}^{\dagger}\hat{L}}dt$. The
following algorithm [1] is equivalent to eq. 3:

  1. 1.

Discretise the time interval $(t_{0},t_{N})\rightarrow\\{t_{k}\\}$ with
$\Delta t=t_{k+1}-t_{k}$;

  2. 2.

Pick a random number $R$ from the uniform distribution on the interval
$[0,1]$;

  3. 3.

For all $t_{k}$ with $t_{0}<t_{k}\leq t_{N}$ evaluate the following loop:

⬇

for k=1:N-1

if $R\geq||\psi_{k}||^{2}$ %Jump

$\ket{\psi_{k+1}}=\hat{L}\ket{\psi_{k}}/\braket{\hat{L}^{\dagger}\hat{L}}$

$R$=rand$\sim U(0,1)$

else %No jump

$\ket{\psi_{k+1}}=\ket{\psi_{k}}-\frac{i}{\hbar}(\hat{H}-\frac{i}{2}\hat{L}^{\dagger}\hat{L})\ket{\psi_{k}}\Delta
t$

end

end

  4. 4.

Normalise the entire trajectory;
$\ket{\psi_{k}}\rightarrow\ket{\psi_{k}}/||\psi_{k}||$ for all $k$.

Considering quantum-jump trajectories of a Markovian open system described by
the Lindblad master equation and postselecting only trajectories in which no
jumps have occurred, thus leads to effective non-Hermitian Hamiltonian
dynamics [24, 25, 26]. Quantum dynamics generated by non-Hermitian
Hamiltonians is an active area of research on its own, and we will make use of
some techniques developed in this context [27, 20]. It should be noted that
$\hat{L}^{\dagger}\hat{L}$ is a positive operator and thus non-Hermitian
Hamiltonians, $\hat{H}-\frac{i}{2}\hat{L}^{\dagger}\hat{L}$, arising in the
context of postselection of Lindblad/quantum-jump dynamics may only describe
loss (but not gain).

Despite converging to the same ensemble dynamics, the SSE and quantum-jump
trajectories can differ not just quantitatively but qualitatively. We will
analyse these differences in detail for quadratic Hamiltonians and linear
Lindbladians with initially Gaussian states. While for Lindblad and SSE
dynamics the state remains Gaussian for all times for quadratic Hamiltonians
and linear Lindbladians, this is in general not the case for quantum-jump
dynamics, even though averaging over quantum-jump trajectories recovers the
Gaussian Lindblad results. We will return to this issue after considering the
Gaussian Lindblad and SSE dynamics.

##  3 Gaussian dynamics

Gaussian states are well suited for the analysis of quantum dynamics, since
they are localised in phase space on order $\hbar$ (in appropriate coordinates
they are minimum uncertainty states), and are the only states with a
completely positive Wigner function [28]. In this regard, they may be thought
of as the most “classical” of quantum states [29], and Gaussian approximations
of the full quantum dynamics lead to simple phase-space dynamics. As has been
observed by Schrödinger already in the early days of quantum mechanics,
Gaussian wavepackets remain Gaussian in the dynamics of quantum harmonic
oscillators and follow classical trajectories. For open systems described by
Lindblad equations or Schrödinger dynamics generated by non-Hermitian
Hamiltonians, a similar statement holds, which allows one to reduce the full
quantum Hilbert space dynamics to a simple phase-space dynamics described by a
handful of parameters for Gaussian states for quadratic Hamiltonians and
linear Lindbladians [27, 30, 14]. Gaussian evolution for stochastic
Schrödinger equations with quadratic Hamiltonian and linear Lindbladians in
position representation has been considered in [18, 31, 32], and here too, an
initial Gaussian state remains Gaussian.

In what follows we shall provide a brief review of the derivation and result
for Lindblad dynamics, and then extend the idea to stochastic Schrödinger
dynamics in quantum phase space, where we use the Wigner-Weyl formalism that
illuminates the underlying phase-space geometry and allows for a better direct
comparison with the Lindblad dynamics. At the heart of the Wigner-Weyl scheme
is the Weyl transformation, a bijective map that maps observables on Hilbert
spaces to their corresponding Weyl symbols. The Weyl symbol corresponding to
an observable $\hat{O}(\hat{X},\hat{P})$ is a distribution on classical phase
space and given by

${O}(x,p)=\int_{-\infty}^{\infty}ds\bra{x-\frac{s}{2}}\hat{O}(\hat{X},\hat{P})\ket{x+\frac{s}{2}}e^{\nicefrac{{ips}}{{\hbar}}}.$
(4)

The Weyl symbol corresponding to the density operator is known as the Wigner
function. We will consider Gaussian states with Wigner functions of the form

$\displaystyle W(z)=\frac{\sqrt{\det G}}{\pi\hbar}e^{-\frac{1}{\hbar}\delta
z\cdot G\delta z}$ with $\displaystyle\delta z=\begin{pmatrix}x-\braket{x}\\\
p-\braket{p}\end{pmatrix}=z-\tilde{z},$ (5)

with a real symmetric matrix $G$, that encodes the phase-space covariance
matrix of the system as

$\Sigma_{ij}=\Delta(z_{i}z_{j})^{2}=\frac{\hbar}{2}G^{-1}_{ij}.$ (6)

This describes a pure state if and only if $\det G=1$ [33]. Inserting an
ansatz of the form (5) with time-dependent (and in the SSE case stochastic)
parameters $G$ and $\tilde{z}$ into the evolution equation for the Wigner
function yields dynamical equations for the parameters.

###  3.1 Gaussian dynamics of the Lindblad equation

The Lindblad equation in Weyl representation takes the form

$\displaystyle\frac{dW}{dt}$ $\displaystyle=\frac{1}{\hbar}\left(-i(H\star
W-W\star H)+L\star W\star\bar{L}-\frac{1}{2}\bar{L}\star L\star
W-\frac{1}{2}W\star\bar{L}\star L\right),$ (7)

with the Moyal (star) product of Weyl symbols given by

$(A\star
B)(q,p)=A(q,p)e^{\frac{i\hbar}{2}\overleftarrow{\nabla}\cdot\Omega\cdot\overrightarrow{\nabla}}B(q,p).$
(8)

Since we are considering only quadratic Hamiltonians and linear Lindbladians,
the Moyal products in eq. 7 can be fully expanded to yield [14]

$\frac{dW}{dt}=-i\nabla\bar{L}\cdot\Omega\nabla LW+\nabla H\cdot\Omega\nabla
W+\operatorname{Im}(L\nabla\bar{L})\cdot\Omega\nabla
W-\frac{\hbar}{2}\operatorname{Re}(\nabla L\cdot\Omega
W^{\prime\prime}\Omega\nabla\bar{L}),$ (9)

where

$\Omega=\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}$ (10)

is the symplectic matrix. An initial Gaussian state remains Gaussian for all
times [14]. It is useful to rewrite $H$ and $L$ as polynomials in $\delta z$

$\displaystyle H(z)$ $\displaystyle=H(\tilde{z})+\nabla
H|_{z=\tilde{z}}\cdot\delta z+\frac{1}{2}\delta z\cdot
H^{\prime\prime}|_{z=\tilde{z}}\delta z$ (11) $\displaystyle L(z)$
$\displaystyle=L(\tilde{z})+\nabla L|_{z=\tilde{z}}\cdot\delta z,$

where we have used a Taylor series (which is exact here) to identify the
linear and quadratic coefficients. We substitute these expansions into eq. 9
to obtain

$\frac{dW}{dt}=[-\frac{2}{\hbar}(\nabla H+H^{\prime\prime}\delta z)\cdot\Omega
G\delta z-\frac{2}{\hbar}\operatorname{Im}(L\nabla\bar{L})\cdot\Omega G\delta
z-\frac{2}{\hbar}\operatorname{Im}(\nabla L\nabla\bar{L})\cdot\Omega G\delta
z\\\ -i\nabla\bar{L}\cdot\Omega\nabla L+\operatorname{Re}(\nabla L\cdot\Omega
G\Omega\nabla\bar{L})-\frac{2}{\hbar}\delta z\cdot
G\Omega\operatorname{Re}(\nabla L\nabla\bar{L}^{T})\Omega G\delta z]W.$ (12)

On the other hand, taking the time derivative of eq. 5 we find

$\displaystyle\frac{dW}{dt}$
$\displaystyle=\left(\frac{1}{2}\Tr(G^{-1}\dot{G})+\frac{2}{\hbar}\dot{\tilde{z}}\cdot
G\delta z-\frac{1}{\hbar}\delta z\cdot\dot{G}\delta z\right)W.$ (13)

Comparing coefficients of $\delta z$ in (12) and (13) we obtain the dynamical
equations

$\displaystyle\frac{d\tilde{z}}{dt}$ $\displaystyle=\Omega\nabla
H+\Omega\text{Im}(L\nabla\bar{L})$ (14) $\displaystyle\frac{dG}{dt}$
$\displaystyle=\left(H^{\prime\prime}+\text{Im}(\nabla
L\nabla\bar{L}^{T})\right)\Omega
G-G\Omega\left(H^{\prime\prime}-\text{Im}(\nabla L\nabla\bar{L}^{T})\right)$
$\displaystyle\quad+2G\Omega\text{Re}(\nabla L\nabla\bar{L}^{T})\Omega G$

for the dynamics of the Gaussian parameters. The first order differential
equation describing the central motion is linear and can be trivially
integrated. The dynamical equation for the covariance matrix G decouples from
the central motion as the Hamiltonian and Lindbladian dependent terms become
constant. The general solution can be obtained as (see [34] and references
therein)

$G(t)=\left[2D(t)+(R^{T}(-t)G^{-1}(0)R(-t))^{-1}\right]^{-1},$ (15)

where

$R(t)={\rm e}^{\left(\Omega H^{\prime\prime}+\text{Im}(\nabla
L\nabla\bar{L}^{T})\Omega\right)t}$ (16)

and

$D(t)=\int_{0}^{t}R(s)\text{Re}(\nabla L\nabla\bar{L}^{T})R^{T}(s){\rm d}s.$
(17)

For Hermitian Lindbladians, the Lindblad term in the central dynamics vanishes
leaving only Hamiltonian dynamics, as expected for these purely decohering
systems. In general non-Hermitian Lindbladians lead to both decoherence (that
may be characterised by the evolution of $G$ in our case) and dissipation that
leads to non-Hamiltonian dynamics of the centre of the Gaussian.

###  3.2 Gaussian Stochastic Schrödinger dynamics

Let us now use the same approach to derive parameter dynamics for the SSE. We
begin by writing the SSE in projector form as

$\displaystyle d(\ket{\psi}\bra{\psi})$
$\displaystyle=\ket{d\psi}\bra{\psi}+\ket{\psi}\bra{d\psi}+\ket{d\psi}\bra{d\psi}$
$\displaystyle=\frac{1}{\hbar}(-{i}\hat{H}-\frac{1}{2}\hat{L}^{\dagger}\hat{L}+\braket{\hat{L}^{\dagger}}\hat{L}-\frac{1}{2}\braket{\hat{L}^{\dagger}}\braket{\hat{L}})\ket{\psi}\bra{\psi}dt$
(18)
$\displaystyle\quad+\frac{1}{\sqrt{2\hbar}}(\hat{L}\hat{\rho}+\hat{\rho}\hat{L}^{\dagger}-\braket{\hat{L}+\hat{L}^{\dagger}}\rho)d\xi_{R}+\frac{i}{\sqrt{2\hbar}}(\hat{L}\hat{\rho}-\hat{\rho}\hat{L}^{\dagger}-\braket{\hat{L}-\hat{L}^{\dagger}}\rho)d\xi_{I}.$

The deterministic ($dt$) part of equation (3.2) is the same as that of the
Lindblad eq. 9 which can evolve pure states into mixed states. In the SSE
case, however, the state remains pure for all times, since the stochastic
terms conspire to conserve the purity of the state.

We translate section 3.2 into the Wigner-Weyl representation to obtain

$dW=-\frac{i}{\hbar}[H\star W-W\star H+i(L\star
W\star\bar{L}-\tfrac{1}{2}\bar{L}\star L\star W-\tfrac{1}{2}W\star\bar{L}\star
L)]dt\\\
+\tfrac{1}{\sqrt{2\hbar}}[L\\!\star\\!W\\!+W\\!\star\\!\bar{L}-\\!\braket{\hat{L}+\hat{L}^{\dagger}}W]d\xi_{R}+\tfrac{i}{\sqrt{2\hbar}}[L\star
W\\!-W\star\bar{L}-\\!\braket{\hat{L}-\hat{L}^{\dagger}}W]d\xi_{I}.$ (19)

Since the deterministic part is the same as that of the Lindblad system eq. 9,
we need only calculate the stochastic terms. We have

$\frac{1}{\sqrt{2\hbar}}(L\star
W+W\star\bar{L}-\braket{\hat{L}+\hat{L}^{\dagger}}W)=\sqrt{\frac{2}{\hbar}}(L^{R}-\braket{L^{R}})W+\sqrt{\frac{\hbar}{2}}\\{W,L^{I}\\},$
(20)

and

$\frac{i}{\sqrt{2\hbar}}(L\star
W-W\star\bar{L}-\braket{\hat{L}-\hat{L}^{\dagger}}W)=-\sqrt{\frac{2}{\hbar}}(L^{I}-\braket{L^{I}})+\sqrt{\frac{\hbar}{2}}\\{W,L^{R}\\},$
(21)

where $\\{,\\}$ denotes the Poisson bracket and $L^{R}/L^{I}$ denote the
real/imaginary parts of $L$. Using the expressions (20) and (21) equation (19)
becomes

$dW=[\\{H,W\\}+\text{Im}(L\\{\bar{L},W\\})-i\\{\bar{L},L\\}W+\tfrac{\hbar}{2}\text{Re}\\{L,\\{\bar{L},W\\}\\}]dt\\\
\quad+(\sqrt{\tfrac{2}{\hbar}}(L^{R}-\\!\braket{L^{R}})W\\!+\\!\sqrt{\tfrac{\hbar}{2}}\\{W,L^{I}\\})d\xi_{R}\\!-\\!(\sqrt{\tfrac{2}{\hbar}}(L^{I}-\\!\braket{L^{I}})-\sqrt{\tfrac{\hbar}{2}}\\{W,L^{R}\\})d\xi_{I}.$
(22)

As in the Lindblad case in section 3.1 we replace the Weyl symbols with their
finite Taylor series eq. 11. Using

$\displaystyle\begin{split}\sqrt{\frac{2}{\hbar}}(L^{R}-\braket{L^{R}})W+\sqrt{\frac{\hbar}{2}}\\{W,L^{I}\\}&=\sqrt{\frac{2}{\hbar}}(\nabla
L^{R}\cdot\delta z+\nabla L^{I}\cdot\Omega G\delta z)W,\\\
-\sqrt{\frac{2}{\hbar}}(L^{I}-\braket{L^{I}})+\sqrt{\frac{\hbar}{2}}\\{W,L^{R}\\}&=-\sqrt{\frac{2}{\hbar}}(\nabla
L^{I}\cdot\delta z-\nabla L^{R}\cdot\Omega G\delta z)W,\end{split}$ (23)

we find

$\displaystyle dW$ $\displaystyle=\big{[}\tfrac{2}{\hbar}\left(\delta z\cdot
G\Omega(H^{\prime\prime}+\text{Im}(\nabla\bar{L}\nabla{L}^{T}))\delta z-\delta
z\cdot G\Omega\text{Re}(\nabla\bar{L}\nabla L^{T})\Omega G\delta z\right.$
(24) $\displaystyle\quad\left.-(\nabla H+\text{Im}(L\nabla\bar{L}))\cdot\Omega
G\delta z-i\tfrac{\hbar}{2}\nabla\bar{L}\cdot\Omega\nabla
L+\tfrac{\hbar}{2}\nabla L\cdot\Omega G\Omega\nabla\bar{L}\right)dt$
$\displaystyle\quad+\sqrt{\tfrac{2}{\hbar}}(\nabla L^{R}\cdot\delta z+\nabla
L^{I}\cdot\Omega G\delta z)d\xi_{R}-\sqrt{\tfrac{2}{\hbar}}(\nabla
L^{I}\cdot\delta z-\nabla L^{R}\cdot\Omega G\delta z)d\xi_{I}\big{]}W.$

On the other hand, using the Itô chain rule [35]

$dW=\frac{\partial W}{\partial
y_{i}}dy_{i}+\frac{1}{2}\frac{\partial^{2}W}{\partial y_{i}\partial
y_{j}}dy_{i}dy_{j},$ (25)

where

$dy_{j}=\mu_{j}dt+\sigma_{j}^{R}d\xi_{R}+\sigma_{j}^{I}d\xi_{I},$ (26)

we have

$\displaystyle dW$ $\displaystyle=\frac{\partial
W}{\partial\tilde{z}_{k}}d\tilde{z}_{k}+\frac{\partial W}{\partial
A_{kl}}dA_{kl}$ (27)
$\displaystyle\quad+\frac{\partial^{2}W}{\partial\tilde{z}_{k}\partial
A_{mn}}d\tilde{z}_{k}dA_{mn}+\frac{1}{2}\frac{\partial^{2}W}{\partial\tilde{z}_{k}\partial\tilde{z}_{l}}d\tilde{z}_{k}d\tilde{z}_{l}+\frac{1}{2}\frac{\partial^{2}W}{\partial
A_{kl}\partial A_{mn}}dA_{kl}dA_{mn}.$

Here we have introduced the matrix $A$ with $G=\frac{1}{2}(A+A^{T})$, to
circumvent complications arising from the symmetry condition on $G$. Here and
for the rest of this section, we implicitly sum over repeated indices
(Einstein summation convention).

As we shall see, equation (27) simplifies drastically as many of the terms
vanish. In particular, we find that

$\displaystyle\sigma_{A}^{R}=0=\sigma_{A}^{I},$ (28)

and thus the second order derivatives evolving $A_{jk}$ in equation (27)
vanish. This can be seen as follows. Let us focus our attention on the term
$\frac{1}{2}\frac{\partial^{2}W}{\partial A_{kl}\partial
A_{mn}}dA_{kl}dA_{mn}$. This term would generate a term $O(\delta z^{4})$

$\frac{1}{2}\frac{\partial^{2}W}{\partial A_{kl}\partial
A_{mn}}dA_{kl}dA_{mn}=\frac{1}{2\hbar^{2}}((\delta z\cdot\sigma_{A}^{R}\delta
z)^{2}+(\delta z\cdot\sigma_{A}^{I}\delta z)^{2})dt\\\ +\text{terms of lower
order in $\delta z$},$ (29)

which would not be matched by any term on the right hand side of equation
(24), from which we infer (28). Hence equation (27) simplifies to

$\displaystyle dW=\frac{\partial
W}{\partial\tilde{z}_{k}}d\tilde{z}_{k}+\frac{\partial W}{\partial
A_{kl}}dA_{kl}+\frac{1}{2}\frac{\partial^{2}W}{\partial\tilde{z}_{k}\partial\tilde{z}_{l}}d\tilde{z}_{k}d\tilde{z}_{k},$
(30)

and we have explicitly

$dW=[\frac{2}{\hbar}(\mu_{z}dt+\sigma_{z}^{R}d\xi_{R}+\sigma_{z}^{I}d\xi_{I})\cdot
G\delta z+\frac{1}{2}\text{Tr}(G^{-1}\mu_{A})dt-\frac{1}{\hbar}\delta
z\cdot\mu_{A}\delta zdt\\\ -\frac{1}{\hbar}(\sigma_{z}^{R}\cdot
G\sigma_{z}^{R}+\sigma_{z}^{I}\cdot
G\sigma_{z}^{I})dt+\frac{2}{\hbar^{2}}(\delta z\cdot
G\sigma_{z}^{R}\sigma_{z}^{R^{T}}G\delta z+\delta z\cdot
G\sigma_{z}^{I}\sigma_{z}^{I^{T}}G\delta z)dt]W.$ (31)

Using $dG=\frac{1}{2}(dA+dA^{T})$ and equating eqs (31) and (24) we obtain the
(stochastic) dynamical equations

$\displaystyle d\tilde{z}$ $\displaystyle=\left(\Omega\nabla
H+\Omega\text{Im}(L\nabla\bar{L})\right)dt$
$\displaystyle\quad+\sqrt{\frac{\hbar}{2}}(G^{-1}\nabla L^{R}-\Omega\nabla
L^{I})d\xi_{R}-\sqrt{\frac{\hbar}{2}}(G^{-1}\nabla L^{I}+\Omega\nabla
L^{R})d\xi_{I}$ (32) $\displaystyle\frac{dG}{dt}$ $\displaystyle=-G\Omega
H^{\prime\prime}+H^{\prime\prime}\Omega G+\text{Re}(\nabla\bar{L}\nabla
L^{T})+G\Omega\text{Re}(\nabla\bar{L}\nabla L^{T})\Omega G$ (33)

for the Gaussian parameters.

We notice that the deterministic part of the dynamics of the centre
$\tilde{z}$ is the same as that of the Lindblad equation, however, the SSE
dynamics have an additional stochastic component as expected. This stochastic
component contains covariance dependent terms and unlike in the Lindblad case,
we can no longer simulate the centre trajectories without calculating the
covariance dynamics. For quadratic systems the evolution of the covariance
matrix $G$ is deterministic and independent of the motion of the centre but
different from that of the Lindblad evolution. This difference is not
surprising taking into account that the $G$ matrix in the SSE describes the
covariances of the individual pure-state trajectories, while the $G$ matrix of
the Lindblad evolution describes that of the total density matrix arising from
the ensemble average. In fact, the dynamics of $G$ for the SSE are the same as
those arising from deterministic non-Hermitian Hamiltonian dynamics, which we
shall briefly review in the next section as the first phase of quantum-jump
dynamics. The dynamical equation (33) for $G$ can be explicitly solved by [27]

$G(t)=\left(-\Omega\operatorname{Re}(S(t)\Omega
G(0)+\Omega\operatorname{Im}(S(t))\right)\left(\operatorname{Im}(S(t))\Omega
G(0)+\operatorname{Re}(S(t))\right)^{-1},$ (34)

where $S(t)$ is given by

$S(t)=e^{\Omega K^{\prime\prime}t},$ (35)

With $K^{\prime\prime}=H^{\prime\prime}-i\operatorname{Re}(\nabla\bar{L}\nabla
L^{T})$. A short calculation also confirms that $\frac{d}{dt}\det G=0$,
implying that the Gaussian state remains pure as expected. Having obtained the
evolution equations for Gaussian dynamics according to the full Lindblad and
SSE dynamics, we shall now turn towards the corresponding quantum-jump
dynamics.

##  4 Quantum-Jump dynamics

Quantum-jump dynamics do not preserve Gaussian states for arbitrary linear
Lindbladians. As an example, consider the Lindblad operator
$\hat{L}=\hat{a}^{\dagger}$. The first jump maps a state $\ket{\psi_{g}}$ to
$\hat{a}^{\dagger}\ket{\psi_{g}}$ and thus transforms a Gaussian state into a
non-Gaussian one. We will show in what follows, that it is still possible to
calculate the quantum-jump dynamics building on the propagation of Gaussian
states, leading again to just a handful of time-dependent parameters. For this
purpose, we adapt a method introduced for non-Hermitian dynamics in [20].

In the quantum jump unravelling of the Lindblad dynamics, we propagate the
initial state $\ket{\psi_{0}}$ with the time evolution generated by the non-
Hermitian Hamiltonian $\hat{H}-\frac{i}{2}\hat{L}^{\dagger}\hat{L}$, and
intersperse it with jumps at discrete times $t_{j}$,

$\ket{\psi}\mapsto\frac{\hat{L}\ket{\psi}}{\sqrt{\braket{\psi}{\hat{L}^{\dagger}\hat{L}}{\psi}}}\,\,.$
(36)

The non-Hermitian nature of time evolution
$U(t)=e^{-\frac{i}{\hbar}(\hat{H}-\frac{i}{2}\hat{L}^{\dagger}\hat{L})}$ will
lead to a decreasing norm of the propagated state, and the quantum jumps reset
the norm to $1$ due to the inclusion of the denominator. After $k$ quantum
jumps at times $t_{1},t_{2},\cdots,t_{k}$, i.e, for $t\in[t_{k},t_{k+1})$, the
state will be of the form

$\ket{\psi(t)}=\frac{1}{\lVert\psi(t_{k})\lVert}U(t-t_{k})LU(t_{k}-t_{k-1})LU(t_{k-1}-t_{k-2})L\cdots
U(t_{2}-t_{1})LU(t_{1})\ket{\psi_{0}}$ (37)

To compute this expression we will introduce a basis which is moving with the
state $\hat{U}(t)\psi_{0}$, the so called Hagedorn basis, which we will recall
in the next subsection. We then show how to compute $\hat{U}(t)$ in the moving
bases for the case that the non-Hermitian Hamiltonian
$\hat{H}-\frac{i}{2}\hat{L}^{\dagger}\hat{L}$ is of no higher than quadratic
order in $\hat{p}$ and $\hat{q}$. Combining these results will allow us to
evaluate (37) explicitly by an iterative algorithm.

###  4.1 Hagedorn Basis

It is well known that a coherent state creates an orthonormal basis by
applying powers of creation operators to it. In [19] Hagedorn introduced a
parametrization of coherent states and their associated raising and lowering
operators which is particularly well adapted to the study of the time
evolution of wavepackets. This has been used in numerical analysis of the time
dependent Schrödinger equation in [36, 37] and it has been adapted to non-
Hermitian evolution in [20]. We will now recall some of the notions we use in
the following, for simplicity we will restrict ourselves to one degree of
freedom systems.

A coherent state centred at the origin is characterised by its annihilation
operator which is defined in terms of a complex vector $a\in\mathbb{C}^{2}$ as

$\hat{A}(a):=\frac{i}{\sqrt{2\hbar}}a\cdot\Omega\hat{z}\,\,.$ (38)

The corresponding coherent state $\ket{0,a}$ is defined up to normalisation by
$\hat{A}(a)\ket{0,a}=0$, and one can show that in position representation it
is given by

$\psi(a,x)=(\pi\hbar)^{-1/4}(a_{q})^{-1/2}e^{\frac{i}{2\hbar}\frac{a_{p}}{a_{q}}x^{2}}\,\,,\quad\text{where}\quad
a=(a_{q},a_{p})\,\,.$ (39)

An important role is played by the Hermitian form defined as

$h_{\Omega}(a,b):=\frac{1}{2i}a^{\dagger}\Omega b\,\,.$ (40)

The state (39) is normalizable if $h_{\Omega}(a,a)>0$ and has norm $1$ if
$h_{\Omega}(a,a)=1$. Notice as well that $h_{\Omega}(\bar{a},a)=0$ and
$h_{\Omega}(\bar{a},\bar{a})=-h_{\Omega}(a,a)$, and with these relations it
follows that if $h_{\Omega}(a,a)=1$ then we have for any complex vector $b$
the expansion

$b=h_{\Omega}(a,b)a-h_{\Omega}(\bar{a},b)\bar{a}\,\,.$ (41)

That is, $a,\bar{a}$ form a basis of $\mathbb{C}^{2}$, similar to an
orthonormal basis.

The creation operator is the adjoint of $\hat{A}(a)$ and we have

$\hat{A}^{\dagger}(a)=-\hat{A}(\bar{a})\,\,,$ (42)

and the corresponding orthonormal basis is defined as

$\ket{n,a}:=\frac{1}{\sqrt{n!}}\big{[}\hat{A}^{\dagger}(a)\big{]}^{n}\ket{0,a}\,\,,\quad
n=0,1,2,3,\cdots\,\,.$ (43)

We can move this basis by applying the phase space translation operators
$\hat{T}(z):=e^{-\frac{i}{\hbar}z\cdot\Omega\hat{z}}$, $z\in\mathbb{R}^{2}$,
to the creation and annihilation operators,
$\hat{A}(a,z):=\hat{T}(z)\hat{A}(a)\hat{T}^{\dagger}(z)$,
$\hat{A}^{\dagger}(a,z):=\hat{T}(z)\hat{A}^{\dagger}(a)\hat{T}^{\dagger}(z)$
and to the basis states

$\ket{n,a,z}:=\hat{T}(z)\ket{n,a}=\frac{1}{\sqrt{n!}}\big{[}\hat{A}^{\dagger}(a,z)\big{]}^{n}\ket{0,a,z}\,\,.$
(44)

We have as well the explicit representation

$\hat{A}(a,z)=\frac{i}{\sqrt{2\hbar}}a\cdot\Omega(\hat{z}-z)\,\,,\quad\hat{A}^{\dagger}(a,z)=-\hat{A}(\bar{a},z),,$
(45)

and we note that

$\hat{A}(a,z)\ket{n,a,z}=\sqrt{n}\ket{n-1,a,z}\,\,\quad\hat{A}^{\dagger}(a,z)\ket{n,a,z}=\sqrt{n+1}\ket{n+1,a,z}\,\,.$
(46)

In the following we will encounter operators of the form
$\hat{A}(b,\chi)=\frac{i}{\sqrt{2\hbar}}b\cdot\Omega(\hat{z}-\chi)$ for some
$b,\chi\in\mathbb{C}^{2}$. It is convenient to write them as linear
combinations of $\hat{A}(a,z)$ and $\hat{A}^{\dagger}(a,z)$ to determine their
action on the basis vectors $\ket{n,a,z}$. To that end, we note that

$\hat{A}(b,\chi)-\hat{A}(b,z)=\frac{i}{\sqrt{2\hbar}}b\cdot\Omega(z-\chi)=\sqrt{2/\hbar}\,h_{\Omega}(\bar{z}-\bar{\chi},b)\,\,,$
(47)

and with the expansion (41) and (45) we then find

$\hat{A}(b,\chi)=h_{\Omega}(a,b)\hat{A}(a,z)+h_{\Omega}(\bar{a},b)\hat{A}^{\dagger}(a,z)+\sqrt{2/\hbar}\,h_{\Omega}(\bar{z}-\bar{\chi},b)\,\,.$
(48)

###  4.2 Non-Hermitian evolution

The propagation of a coherent state under a time evolution generated by a non-
Hermitian operator has been studied in [38] and the special case of quadratic
Hamiltonians was analysed in [27] and in the context of Hagedorn wave-packets
in [20]. In order to apply these results we need first to rewrite the
Hamiltonian $\hat{H}-\frac{i}{2}\hat{L}^{\dagger}\hat{L}$ slightly; since we
assume that $L$ is linear, we have $(\bar{L}\star
L)(z)=\bar{L}(z)L(z)+\frac{i\hbar}{2}\\{\bar{L},L\\}$ and the term
$\\{\bar{L},L\\}=\nabla\bar{L}\cdot\Omega\nabla L$ is constant. This gives us
$\hat{L}^{\dagger}\hat{L}=\widehat{\bar{L}L}+\frac{i\hbar}{2}\\{\bar{L},L\\}$
and so we obtain

$\hat{U}(t)=e^{-\frac{i}{\hbar}(\hat{H}-\frac{i}{2}\hat{L}^{\dagger}\hat{L})t}=e^{-\frac{i}{4}\nabla\bar{L}\cdot\Omega\nabla
Lt}e^{-\frac{i}{\hbar}\hat{K}t}\quad\text{where}\quad
K(z)=H(z)-\frac{i}{2}|L(z)|^{2}\,\,.$ (49)

In [38], the dynamics of a Gaussian wavepacket under a non-Hermitian
Hamiltonian were derived following a similar procedure to the one we have
outlined for the Lindblad and SSE cases. Substituting the effective non-
Hermitian Hamiltonian $K(z)$ (49) into the results from [38] yields the
parameter dynamics

$\displaystyle\frac{d\tilde{z}}{dt}$ $\displaystyle=\Omega\nabla
H-G^{-1}\text{Re}(\bar{L}\nabla L)$ (50) $\displaystyle\frac{dG}{dt}$
$\displaystyle=-G\Omega H^{\prime\prime}+H^{\prime\prime}\Omega
G+\text{Re}(\nabla\bar{L}\nabla L^{T})+G\Omega\text{Re}(\nabla\bar{L}\nabla
L^{T})\Omega G.$ (51)

Where the evolution equation for $G$ is the same as the one for the SSE case.
As expected this fulfils $\frac{\rmd}{\rmd t}\det G=0$ for $\det G=1$, and an
initially pure state remains pure.

As has been discussed in [14] the dissipative part of the central motion of
the non-Hermitian dynamics can appear either quite different or very similar
to that of the Lindblad case, depending on the structure of the Lindblad
operator. For a Lindblad operator that is an analytic function of $\hat{a}$ or
$\hat{a}^{\dagger}$, for example, the dissipative term in the central dynamics
in equation (14), given by $\Omega\text{Im}(L\nabla\bar{L})$ can be rewritten
as $-\text{Re}(\bar{L}\nabla L)$, which is very similar to the non-Hermitian
dissipation, with the difference that the latter is modulated by the changing
covariance metric $G$. An example for which Lindblad and non-Hermitian central
dynamics are very different, are Hermitian Lindbladians, for which the
dissipative term in the Lindblad dynamics vanishes entirely. The quantum-jump
evolution turns the non-Hermitian behaviour into the Lindblad one, by
averaging over different quantum jumps, that in general do not leave an
initially Gaussian state Gaussian.

In practice it is often useful to use instead of (50) and (51) the complex
classical dynamics created by $K(z)$. If we write $K(z)=\frac{1}{2}z\cdot
K_{2}z+k_{1}\cdot\Omega z+k_{0}$, where $K_{2}$ is a symmetric complex matrix,
$k_{1}$ is a complex vector, and $k_{0}$ is a constant, then the corresponding
solution to Hamilton’s equations in phase space is given by

$\Phi(t,z)=S(t)z+v(t)\,\,,$ (52)

with

$S(t)=e^{t\Omega K_{2}}\,\,\quad\text{and}\quad
v(t)=\int_{0}^{t}S(t-s)k_{1}\,ds\,\,,$ (53)

where $S(t)$ is complex and symplectic, i.e., $S^{T}\Omega S=\Omega$. It turns
out that $\hat{U}(t)\ket{0,a_{0},z_{0}}$ can be described entirely in term of
$\Phi(t),S(t)$ and its action on $a_{0}$ and $z_{0}$, [27, 20], to that end
let us first define

$N(t):=\frac{1}{\sqrt{h_{\Omega}(S(t)a_{0},S(t)a_{0})}}\quad\text{and}\quad
a_{t}:=N(t)S(t)a_{0}\,\,,$ (54)

so that $a_{t}$ is normalised again. Now we introduce
$J_{t}:=-\operatorname{Re}(a_{t}a_{t}^{\dagger})\Omega$, where
$a_{t}^{\dagger}$ denotes the transposed and the complex conjugate so that
$a_{t}a_{t}^{\dagger}$ is a $2\times 2$ matrix, then the solution to (50) is
given by

$z_{t}:=\operatorname{Re}\Phi(t,z)+J\operatorname{Im}\Phi(t,z)\in\mathbb{R}^{2}\,\,,$
(55)

and this vector is real-valued, in contrast to the complex centre
$z_{t}^{(\mathbb{C})}=\Phi(t,z)$. There is as well a corresponding expression
for $G$ in (51)

$G(t)=\Omega
J_{t}=\Omega^{T}\operatorname{Re}(a_{t}a_{t}^{\dagger})\Omega\,\,.$ (56)

Using these definitions we can write the propagated coherent state as

$\hat{U}(t)\ket{0,a_{0},z_{0}}=e^{\frac{i}{\hbar}\alpha(t)}\sqrt{N(t)}\ket{0,a_{t},z_{t}}$
(57)

where

$\alpha(t)=\int_{0}^{t}\dot{q}_{s}p_{s}-\frac{1}{2}z_{s}Kz_{s}\,{\rm
d}s-\frac{\hbar}{4}\nabla\bar{L}\cdot\Omega\nabla Lt$ (58)

with $z_{s}=(q_{s},p_{s})$ and we have incorporated the factor containing
$\nabla\bar{L}\cdot\Omega\nabla L$ from (49) into $\alpha$. The state
$\ket{0,a_{t},z_{t}}$ is now normalised, so all the information about the
decay of the norm of $\hat{U}(t)\ket{0,a,z}$ is contained in

$\sqrt{N(t)}|e^{\frac{i}{\hbar}\alpha(t)}|=\sqrt{N(t)}e^{-\frac{1}{2\hbar}\int_{0}^{t}|L(z_{s})|^{2}{\rm
d}s+\frac{1}{4}\nabla\bar{L}\cdot\Omega\nabla Lt}\,\,.$ (59)

In the following we will use the moving basis
$\\{\ket{n,a_{t},z_{t}}\,:\,n\in\mathbb{N}_{0}\\}$ associated with $a_{t}$ and
$z_{t}$ by (44), which is an orthonormal basis centred around the moving state
$U(t)\ket{0,a_{0},z_{0}}$.

We want to emphasise that if $\hat{U}(t)$ is non-unitary then
$\hat{U}(t)\ket{n,a_{0},z_{0}}$ is in general not proportional to
$\ket{n,a_{t},z_{t}}$ but will acquire contributions from lower order excited
states. In order to compute the expansion of $\hat{U}(t)\ket{n,a_{0},z_{0}}$
into the basis $\ket{n,a_{t},z_{t}}$ at $t$ we use (117) from A,
$\hat{U}(t)\hat{A}^{\dagger}(a_{0},z_{0})\hat{U}(-t)=-\hat{A}(S(t)\bar{a}_{0},S(t)z_{0})$

$\begin{split}\hat{U}(t)\ket{n,a_{0},z_{0})}&=\frac{1}{\sqrt{n!}}\hat{U}(t)\big{[}\hat{A}^{\dagger}(a_{0},z_{0})\big{]}^{n}\ket{0,a_{0},z_{0}}\\\
&=\frac{1}{\sqrt{n!}}\big{[}-\hat{A}(S(t)\bar{a}_{0},S(t)z_{0})\big{]}^{n}\hat{U}(t)\ket{0,a_{0},z_{0}}\,\,.\end{split}$
(60)

Now we expand $\hat{A}(S(t)\bar{a}_{0},S(t)z_{0})$ in terms of the
annihilation and creation operators of our moving basis at time $t$ using (48)
which gives

$-\hat{A}(S(t)\bar{a}_{0},S(t)z_{0})=h_{+}\hat{A}^{\dagger}(a_{t},z_{t})+h_{-}\hat{A}(a_{t},z_{t})+h_{0}$
(61)

where $h_{+}=-h_{\Omega}(\bar{a}_{t},S(t)\bar{a}_{0})$,
$h_{-}=-h_{\Omega}(a_{t},S(t)\bar{a}_{0})$ and
$h_{0}=-\sqrt{2/\hbar}h_{\Omega}(z_{t}-\bar{S}(t)z_{0},S(t)\bar{a}_{0})$. It
will be useful to introduce

$M(t):=\frac{h_{\Omega}(S(t)a_{0},S(t)\bar{a}_{0})}{h_{\Omega}(S(t)a_{0},S(t)a_{0})}$
(62)

and then we have $h_{+}=N(t)$, $h_{-}=-M(t)/N(t)$ and

$h_{0}(t)=\sqrt{\frac{2}{\hbar}}\
\big{[}h_{\Omega}(a_{0},z_{0})(N(t)^{2}-1)+h_{\Omega}(\bar{a}_{0},z_{0})M(t)\big{]}\,\,,$
(63)

valid if $h_{\Omega}(a_{t},a_{t})=1$.

In (117) in B we show that

$\frac{1}{\sqrt{n!}}\big{[}h_{+}\hat{A}^{\dagger}(a_{t},z_{t})+h_{-}\hat{A}(a_{t},z_{t})+h_{0}\big{]}^{n}\ket{0,a_{t},z_{t}}=\sum_{m=0}^{n}B_{nm}(t)\ket{m,a_{t},z_{t}}$
(64)

with

$B_{nm}(t)=N(t)^{m}\sum_{k=0}^{[\frac{n-m}{2}]}\sqrt{\frac{n!}{m!}}\frac{(-M(t))^{k}h_{0}^{n-m-2k}}{2^{k}(n-m-2k)!k!}\,\,.$
(65)

where $[\frac{n-m}{2}]$ denotes the floor of $\frac{n-m}{2}$, and hence we get
an explicit expression for the matrix elements

$\bra{m,a_{t},z_{t}}\hat{U}(t)\ket{n,a_{0},z_{0}}=\begin{cases}0&m>n\\\
e^{\frac{i}{\hbar}\alpha(t)}\sqrt{N(t)}\,B_{mn}(t)&m\leq n\end{cases}\,\,.$
(66)

For the first few states, we find explicitly

$\hat{U}(t)\ket{1,a_{0},z_{0}}=e^{\frac{i}{\hbar}\alpha(t)}\sqrt{N(t)}\big{[}N(t)\ket{1,a_{t},z_{t}}+h_{0}\ket{0,a_{t},z_{t}}\big{]}\,\,,$
(67)
$\begin{split}\hat{U}(t)\ket{2,a_{0},z_{0}}&=e^{\frac{i}{\hbar}\alpha(t)}\sqrt{N(t)}\bigg{[}N(t)^{2}\ket{2,a_{t},z_{t}}+\sqrt{2}\,h_{0}N(t)\ket{1,a_{t},z_{t}}\\\
&\hskip
142.26378pt+\frac{1}{\sqrt{2}}\big{(}h_{0}^{2}-2M(t)\big{)}\ket{0,a_{t},z_{t}}\bigg{]}\,\,.\end{split}$
(68)

We can generalise (66) to the situation that the initial state is at a time
$t_{0}\neq 0$. To that end, it is useful to introduce a notation for $N(t)$
and $M(t)$ which makes the dependence on the vector $a$ explicit,

$N(t,a_{0}):=\big{[}h_{\Omega}(S(t)a_{0},S(t)a_{0})\big{]}^{-1/2}\,\,,\quad
M(t,a_{0}):=N^{2}(t,a)h_{\Omega}(S(t)a_{0},S(t)\bar{a}_{0})$ (69)

and we assume that $a$ satisfies $h_{\Omega}(a,a)=1$. It is not hard to show
that $N(t,a)$ satisfies

$N(t_{2}-t_{0},a_{t_{0}})=N(t_{2}-t_{1},a_{t_{1}})N(t_{1}-t_{0},a_{t_{0}})$
(70)

for $t_{2}\geq t_{1}\geq t_{0}\geq 0$ and where $a_{t}:=N(t,a_{0})S(t)a_{0}$.
We then use (63) to extend this notation to $h_{0}$,

$h_{0}(t,a_{0},z_{0}):=\sqrt{\frac{2}{\hbar}}\
\big{[}h_{\Omega}(a_{0},z_{0}))(N(t,a_{0})^{2}-1)+h_{\Omega}(\bar{a}_{0},z_{0})M(t,a_{0})\big{]}\,\,.$
(71)

With these notations we have for
$\mathbf{U}_{mn}(t_{2},t_{1}):=\bra{m,a_{t_{2}},z_{t_{2}}}\hat{U}(t)\ket{n,a_{t_{1}},z_{t_{1}}}$

$\mathbf{U}_{mn}(t_{2},t_{1})=\begin{cases}0&m>n\\\
e^{\frac{i}{\hbar}[\alpha(t_{2})-\alpha(t_{1})]}\sqrt{N(t_{2}-t_{1},a_{t_{1}})}\,B_{mn}(t_{2}-t_{1},a_{t_{1}},z_{t_{1}})&m\leq
n\end{cases}\,\,,$ (72)

where

$B_{mn}(t,a_{0},z_{0}):=N(t,a_{0})^{m}\sum_{k=0}^{[\frac{n-m}{2}]}\sqrt{\frac{n!}{m!}}\frac{(-M(t,a_{0}))^{k}h_{0}(t,a_{0},z_{0})^{n-m-2k}}{2^{k}(n-m-2k)!k!}\,\,.$
(73)

###  4.3 Quantum jumps in a Hagedorn basis

To implement the quantum jumps generated by the Lindblad operator $\hat{L}$ at
time $t$ we need to represent $\hat{L}$ in the basis at time $t$. If we
parametrize $\hat{L}$ as
$\hat{L}=\frac{i}{\sqrt{2\hbar}}l\cdot\Omega(\hat{z}-\chi)=\hat{A}(l,\chi)$
then (48) gives immediately

$\hat{L}=h_{\Omega}(a_{t},l)\hat{A}(a_{t},z_{t})+h_{\Omega}(\bar{a}_{t},l)\hat{A}^{\dagger}(a_{t},z_{t})+\sqrt{2/\hbar}\,h_{\Omega}(z_{t}-\bar{\chi},l)\,,$
(74)

hence the corresponding matrix of $\hat{L}$ in the basis at time $t$ is given
by

$\mathbf{L}_{nm}(t)=h_{\Omega}(a_{t},l)\sqrt{n}\,\delta_{n-1,m}+h_{\Omega}(\bar{a}_{t},l)\sqrt{n+1}\delta_{n+1,m}+\sqrt{2/\hbar}\,h_{\Omega}(z_{t}-\bar{\chi},l)\delta_{n,m}\,\,.$
(75)

With the explicit expression (72) and (75) we can now use the following scheme
to compute the quantum-jump dynamics.

  1. 1.

Choose a time interval $[t_{0},t_{end}]$ and initial parameters $z_{0}$,
$a_{0}$ defining the Hagedorn basis.

  2. 2.

Expand the initial state $\ket{\psi_{0}}$ in the basis $\ket{n,a_{0},z_{0}}$,
this defines a normalised vector $\mathbf{c}(t_{0})$.

  3. 3.

Pick a random number $R$ from the uniform distribution on the interval $[0,1]$

  4. 4.

Solve

$\mathbf{c}^{\dagger}(t_{0})\mathbf{U}^{\dagger}(t_{J},t_{0})\mathbf{U}(t_{J},t_{0})\mathbf{c}(t_{0})-R=0$

for $t_{J}$

  5. 5.

Calculate new coefficients

$\mathbf{c}^{\prime}(t_{J})\to\frac{\mathbf{L}(t_{J})\mathbf{c}(t_{J})}{\sqrt{\mathbf{c}(t_{J})^{\dagger}\mathbf{L}^{\dagger}(t_{J})\mathbf{L}(t_{J})\mathbf{c}(t_{J})}}$

  6. 6.

While the chosen end time $t_{end}\geq t_{J}$, set $t_{0}=t_{J}$ and repeat
steps (iii)-(vi)

  7. 7.

Evaluate the end state

$\ket{\psi(t)}=\sum_{n}c_{n}\ket{n,a_{t},z_{t}}$

at $t=t_{end}$

The above scheme is especially efficient if the initial state $\ket{\psi_{0}}$
is a Gaussian as our coefficient vector $\mathbf{c}(t_{0})$ will only have one
non-zero element. We also have that at each jump only increases the number of
non-zero elements of $\mathbf{c}(t)$ by one in the Hagedorn basis, since any
linear $\hat{L}$ may be written as a linear combination of a
$\hat{A}(a_{t},z_{t})$, $\hat{A}^{\dagger}(a_{t},z_{t})$ and the identity by
eq. 74. Further efficiency gains over standard methods are made by dynamically
changing the basis size and using a root finding algorithm to solve step (v)
since we do not need to calculate the propagator at many points. This scheme
is most useful for studying the long-term behaviour in systems with long times
in between jumps (weakly coupled Lindbladians).

We also will give an alternative scheme that calculates the time evolved
states $\hat{U}(t)\ket{n,a_{0},z_{0}}$ at all times and uses them as a non-
orthogonal basis. This allows one to calculate all deterministic Hamiltonian
evolution a priori greatly reducing the computational cost to calculate
repeated trajectories. We start by calculating the norms and overlaps of the
time evolved states using eq. 73

$\displaystyle O_{mn}(t)$
$\displaystyle=\bra{m,a_{0},z_{0}}\hat{U}^{\dagger}(t)\hat{U}(t)\ket{n,a_{0},z_{0}}$
(76)
$\displaystyle=e^{-\frac{2}{\hbar}\operatorname{Im}(\alpha(t))}N(t)\sum_{i,j=0}B^{\dagger}_{im}(t)\braket{i,a_{t},z_{t}}{j,a_{t},z_{t}}B_{nj}(t)$
$\displaystyle=e^{-\frac{2}{\hbar}\operatorname{Im}(\alpha(t))}N(t)(\overline{\mathbf{B}}\mathbf{B}^{T})_{mn}(t).$

We also find a matrix representation of $\hat{L}$ in the non-orthogonal basis.
We use a formalism found in [39] and find a set of states
$\hat{\widetilde{U}}(t)\ket{m,a,z}$ dual to $\hat{U}(t)\ket{n,a_{0},z_{0}}$
such that

$\bra{m,a_{0},z_{0}}\hat{\widetilde{U}}\,^{\dagger}(t)\hat{U}(t)\ket{n,a_{0},z_{0}}=\delta_{mn}.$
(77)

We thus have $\hat{\widetilde{U}}(t)=(\hat{U}^{-1})^{\dagger}(t)$.

$\hat{\widetilde{U}}(t)\ket{n,a_{0},z_{0}}=\sum_{m=0}\frac{e^{\frac{i}{\hbar}\bar{\alpha}(t)}\widetilde{B}_{nm}(t)}{\sqrt{N(t)}}\ket{m,a_{t},z_{t}}$
(78)

with

$\widetilde{B}_{nm}(t)=N(t)^{-m}\sum_{k=0}^{[\frac{m-n}{2}]}\sqrt{\frac{m!}{n!}}\frac{(-1)^{m+n}\overline{M}(t)^{k}h_{0}^{m-n-2k}}{2^{k}(m-n-2k)!k!}.$
(79)

Using this basis the matrix elements of $\hat{L}$ are given simply by

$\displaystyle\mathbb{L}_{mn}(t)$
$\displaystyle=\bra{m,a_{0},z_{0}}\hat{\widetilde{U}}\,^{\dagger}(t)\hat{L}\hat{U}(t)\ket{n,a_{0},z_{0}}$
(80)
$\displaystyle=\sum_{i,j=0}\widetilde{B}^{\dagger}_{im}(t)\bra{i,a_{t},z_{t}}\hat{L}\ket{j,a_{t},z_{t}}B_{nj}(t)$
$\displaystyle=\sum_{i,j=0}\widetilde{B}^{\dagger}_{im}(t)\mathbf{L}_{ij}B_{nj}(t)$

Substituting eq. 75 into the above we obtain

$\mathbb{L}_{mn}(t)=\sum_{i,j=0}\overline{\widetilde{B}}_{mi}(t)\big{(}\sqrt{j+1}h_{\Omega}(\bar{a}_{t},l)\delta_{i,j+1}\\\
+\sqrt{i}h_{\Omega}(a_{t},l)\delta_{i,j-1}+\sqrt{2/\hbar}h_{\Omega}(z_{t}-\bar{\chi},l)\delta_{i,j}\big{)}B^{T}_{jn}(t).$
(81)

Using this matrix representation of the Lindblad operator we may calculate the
effect of a jump by acting directly on the coefficient vector $\mathbf{c}$.
The jump maps

$\ket{\psi}=\sum_{n}c_{n}\hat{U}(t)\ket{n,a_{0},z_{0}}\mapsto\ket{\psi^{\prime}}=\sum_{n}c^{\prime}_{n}\hat{U}(t)\ket{n,a_{0},z_{0}},\quad\text{with}\quad\mathbf{c}^{\prime}=\frac{\mathbb{L}\mathbf{c}}{\sqrt{\mathbf{c}^{\dagger}\mathbb{L}^{\dagger}\mathbf{O}\mathbb{L}\mathbf{c}}}.$
(82)

The coefficients after the jump are then used as the new input coefficients
for the next stretch of non-Hermitian evolution until the next jump. This can
be implemented using the following algorithm, although we may propagate any
arbitrary state initial state some simplifications occur when using an
initially Gaussian state.

  1. 1.

Discretise the time interval $(t_{0},t_{N})\rightarrow\\{t_{k}\\}$ with
$\Delta t=t_{k+1}-t_{k}$;

  2. 2.

Choose $l(t_{0})$ and $z_{0}$ such that the Gaussian we wish to propagate is
given by eq. 39. With this choice of parameters we have $c_{0}(0)=1$ and
$c_{k\neq 0}(0)=0$

  3. 3.

Calculate $B_{kl}(t_{n}),\widetilde{B}_{kl}(t_{n}),O_{kl}(t_{n})$ and
$L_{kl}(t_{n})$ for all $t_{n}$.

  4. 4.

Pick a random number $R$ from the uniform distribution on the interval
$[0,1]$;

  5. 5.

For all $t_{k}$ with $t_{0}<t_{k}\leq t_{N}$ evaluate the following loop:

⬇

for k=1:N-1

if $R\geq\textbf{c}(t_{n})^{\dagger}\textbf{O}(t_{n})\textbf{c}(t_{n})$

$\textbf{c}(t_{n+1})=\mathbb{L}(t_{n})\textbf{c}(t_{n})/\sqrt{\textbf{c}^{\dagger}(t_{n})\mathbb{L}^{\dagger}(t_{n})\textbf{O}(t_{n})\mathbb{L}(t_{n})\textbf{c}(t_{n})}$

$R$=rand

else

$\textbf{c}(t_{n+1})=\textbf{c}(t_{n})$

end

end

  6. 6.

Normalise the trajectory; for all $t_{n}$

$\mathbf{c}(t_{n})\to\frac{\mathbf{c}(t_{n})}{\sqrt{\mathbf{c}^{\dagger}(t_{n})\mathbf{O}(t_{n})\mathbf{c}(t_{n})}}$

.

  7. 7.

Construct state trajectory
$\ket{\psi(t_{n})}=\sum_{k=0}c_{k}(t_{n})\hat{U}(t_{n})\ket{k,a_{t_{0}},z_{t_{0}}}$
for all $t_{n}$

An advantage of using a Hagedorn basis to simulate quantum jumps over the
standard method described in section 2 is that a much smaller basis size may
be used to accurately propagate an initially Gaussian state as each jump only
increases the number of non-zero elements of the state vector $\mathbf{c}$ by
one. In addition one can see from eq. 65 that the state
$\hat{U}(t)\ket{n,a_{0},z_{0}}$ is damped by a factor $N(t)^{n}$ and as
$N(t)\leq 1$ and strictly decreasing, the higher states have diminishing
contribution to the overall dynamics. For a single trajectory, it may still be
more efficient to avoid the extra effort of propagating the whole basis set in
time and calculating the overlap matrices, however, this may be compensated
for if the scheme is used to generate a large number of quantum-jump
trajectories since the time-dependent basis is the same for every realisation.
In the examples discussed, we will use the above algorithm to implement
quantum-jump trajectories numerically as we wish to obtain observable
expectation values over the full trajectory.

In summary, the quantum-jump dynamics for a quadratic Hamiltonian and linear
Lindbladian while in general non-Gaussian, can be understood and simulated
almost entirely on the grounds of the single dynamical quantity $S(t)$
describing the linearised complexified flow intercepted with discrete quantum
jumps. In what follows we shall explore the resulting dynamics for two
instructive examples.

##  5 Examples

To illustrate the results above, let us consider the Lindblad dynamics and the
two unravellings for a harmonic oscillator Hamiltonian
$\hat{H}=\frac{\omega}{2}\left(\hat{p}^{2}+\hat{x}^{2}\right)$ with two
different Lindblad operators, one Hermitian and one non-Hermitian.

###  5.1 Example 1: Position measurement

Let us first consider a quantum harmonic oscillator with a Hermitian Lindblad
operator

$\hat{L}=\sqrt{\gamma}\hat{x},$ (83)

which can be thought of as modelling a position measurement. Since the
Lindbladian is Hermitian (purely decohering), in the Lindblad dynamics it
yields no contribution to the dynamics of the expectation values $z_{t}$,
which simply follow the familiar harmonic oscillator trajectories. This is
different for individual SSE and quantum-jump trajectories.

The dynamical equations for the Gaussian parameters in the Lindblad dynamics
eq. 14 simplify to

$\displaystyle\frac{d\tilde{z}}{dt}$ $\displaystyle=\omega\Omega\tilde{z},$
(84) $\displaystyle\frac{dG}{dt}$ $\displaystyle=\omega(\Omega
G-G\Omega)+2G\Omega\Gamma\Omega G,$ (85)

where we have defined

$\Gamma=\operatorname{Re}(\nabla
L\nabla\bar{L}^{T})=\begin{pmatrix}\gamma&&0\\\ 0&&0\end{pmatrix}.$ (86)

While the central dynamics is that of the unitary harmonic oscillator, the
dynamics of the covariances encoded by $G(t)$, are influenced by the position
measurement. Let us consider the simple example of an initially squeezed state
with, $G(0)=\left(\begin{smallmatrix}\zeta&&0\\\
0&&\nicefrac{{1}}{{\zeta}}\end{smallmatrix}\right)$. Solving eq. 85 and
substituting the result into eq. 6 we find the physical variances as

$\displaystyle\Delta x^{2}(t)$
$\displaystyle=\frac{\hbar}{4}\left(\frac{\zeta^{2}+2\gamma\zeta
t+1}{\zeta}-\frac{\gamma}{\omega}\sin(2\omega
t)-\frac{\zeta^{2}-1}{\zeta}\cos(2\omega t)\right)$ (87) $\displaystyle\Delta
p^{2}(t)$ $\displaystyle=\frac{\hbar}{4}\left(\frac{\zeta^{2}+2\gamma\zeta
t+1}{\zeta}+\frac{\gamma}{\omega}\sin(2\omega
t)+\frac{\zeta^{2}-1}{\zeta}\cos(2\omega t)\right)$ $\displaystyle\Delta
xp(t)$
$\displaystyle=\frac{\hbar}{4}\left(\frac{\gamma}{\omega^{2}\zeta}+\frac{\zeta^{2}\omega-\omega-\gamma\zeta}{\omega^{2}\zeta^{2}}\cos(2\omega
t)\right)$

That is, we observe the typical harmonic oscillations with frequency $2\omega$
in the covariances as they appear in the unitary harmonic oscillator,
accompanied by a linear growth of the position and momentum uncertainties
$\Delta x^{2}$ and $\Delta p^{2}$, associated with the effect of the position
measurement. This behaviour is illustrated in figure 1 which depicts the
expectation values of position and momentum and their uncertainties as a
function of time for an example with $\omega=1$ and $\gamma=0.2$ for an
initially squeeze state centred at $\tilde{z}_{t}=(2,0)^{T}$. The Lindblad
dynamics are depicted as solid black lines.

Figure 1:  Lindblad dynamics (solid black line) compared with single
trajectories of the SSE (dashed red line) and quantum-jump method (dotted blue
line) for the position measurement model eq. 83, with $\omega=1$ and
$\gamma=0.2$. The initial Gaussian is a squeezed state (i.e.
$a_{0}=(\tfrac{1}{\sqrt{2}},i\sqrt{2})^{T}$ or
$G=\left(\begin{smallmatrix}2&0\\\
0&\nicefrac{{1}}{{2}}\end{smallmatrix}\right)$), centered at
$\tilde{z}=\left(2,0\right)^{T}$. We show the time dependence of the position
expectation $\braket{\hat{x}}$ (top left), the momentum expectation
$\braket{\hat{p}}$ (top right), the positional variance $\Delta x^{2}$ (bottom
left) and the momentum variance $\Delta p^{2}$ (bottom right).

The SSE dynamical equations (32) and (33) become

$\displaystyle d\tilde{z}$
$\displaystyle=\omega\Omega\tilde{z}dt+\sqrt{\frac{\hbar\gamma}{2}}\,G^{-1}\begin{pmatrix}1\\\
0\end{pmatrix}d\xi_{R}+\sqrt{\frac{\hbar\gamma}{2}}\begin{pmatrix}0\\\
1\end{pmatrix}d\xi_{I},$ (88) $\displaystyle\frac{dG}{dt}$
$\displaystyle=\omega(\Omega G-G\Omega)+\Gamma+G\Omega\Gamma\Omega G.$

That is, for the central motion we again have the familiar unitary Hamiltonian
flow term, and no damping term, but now there is an additional width-dependent
stochastic noise. The equation for the covariances differs from that in the
Lindbladian case, as expected. As discussed above $G(t)$ can be solved
analytically by equation (34) or (56). The linearised flow $S(t)$ is given by

$S(t)=\begin{pmatrix}\cosh(\sqrt{\frac{\omega}{2}}\Lambda
t)&\sqrt{\frac{\omega}{2\lambda^{2}}}\Lambda^{*}\sinh(\sqrt{\frac{\omega}{2}}\Lambda
t)\\\ \sqrt{\frac{1}{2\omega}}\Lambda\sinh(\sqrt{\frac{\omega}{2}}\Lambda
t)&\cosh(\sqrt{\frac{\omega}{2}}\Lambda t)\end{pmatrix},$ (89)

where we have defined

$\lambda=\sqrt{\gamma^{2}+\omega^{2}}\quad\text{and}\quad\Lambda=\sqrt{\lambda-\omega}+i\sqrt{\lambda+\omega}.$
(90)

That is, we have an oscillatory contribution with frequency
$\sqrt{\frac{\omega(\lambda+\omega)}{2}}$, which reduces to oscillations with
frequency $\omega$ in the limit $\gamma=0$, and an additional exponential
growth with rate $\sqrt{\lambda-\omega}$. The asymptotic behaviour for large
times is given by

$\lim_{t\to\infty}S(t)=e^{\sqrt{\frac{\omega}{2}}\Lambda
t}\begin{pmatrix}1&\sqrt{\frac{\omega}{2\lambda^{2}}}\Lambda^{*}\\\
\sqrt{\frac{1}{2\omega}}\Lambda&1\end{pmatrix}.$ (91)

As a result, independent of the initial value $G(0)$, $\Sigma(t)$ tends to a
fixed point as $t\to\infty$ given by

$\Sigma(t)\to\frac{\hbar}{2\gamma}\begin{pmatrix}\sqrt{2\omega(\lambda-\omega)}&&\lambda-\omega\\\
\lambda-\omega&&\lambda\sqrt{\frac{2(\lambda-\omega)}{\omega}}\end{pmatrix}.$
(92)

This is in stark contrast to the behaviour of the Lindblad covariances, with
their linear growth in $\Delta x^{2}$ and $\Delta p^{2}$. We can see this in
fig. 1 which shows the SSE dynamics as red dashed lines. For the parameter
choices corresponding to the example in fig. 1 ($\omega=1,\gamma=0.2$) we have
$\lambda=\frac{\sqrt{26}}{5}$, and the covariances approach

$\displaystyle\Delta x^{2}(t)$
$\displaystyle\to\frac{5\sqrt{2(\lambda-1)}}{2}\approx 0.4975$ (93)
$\displaystyle\Delta p^{2}(t)$
$\displaystyle\to\frac{5\lambda\sqrt{2(\lambda-1)}}{2}\approx 0.5074$
$\displaystyle\Delta xp(t)$ $\displaystyle\to\frac{5(\lambda-1)}{2}\approx
0.0495,$

that is, the final state is very close to a coherent state, due to the
relatively small value of $\gamma$. For the central dynamics we observe
stochastic fluctuations around the average Lindblad dynamics.

Figure 2: Top row: Wigner functions of the single quantum jump trajectory in
figure 1 at selected times (from left to right: t=3.36 (shortly before the
first jump) and t=3.37 (just after the first jump) and t=6.47 (just after the
second jump)). The white line traces the preceding central motion. Bottom row:
Relative magnitudes of the coefficients of the state in the time evolved basis
($\hat{U}(t)\ket{n,a_{0},z_{0}}$) at the same times as in the top row.

The quantum-jump trajectories, on the other hand, depicted for an example run
as blue dotted lines in figure 1, show very different behaviour. Here up to
the first jump, the centre of the Gaussian state follows the non-Hermitian
dynamics, which in the present case reduce to

$\begin{pmatrix}\dot{x_{t}}\\\
\dot{p_{t}}\end{pmatrix}=\omega\begin{pmatrix}p_{t}\\\
-x_{t}\end{pmatrix}-\frac{2\gamma x_{t}}{\hbar}\begin{pmatrix}\Delta
x^{2}(t)\\\ \Delta xp(t)\end{pmatrix},$ (94)

where $G(t)$ evolves dynamically as in the SSE case. That is, there is an
additional position dependent damping term in the evolution, modulated by the
covariances of the state. This damping in comparison to the Lindblad evolution
is visible in the example depicted in figure 1. We also observe in figure 1
that, as expected, the dynamics of position and momentum variances agree
between the quantum jump and the SSE dynamics up to the first jump. What is
not shown here, but has been numerically verified, is that averaging over many
quantum jump trajectories simulated in the Hagedorn basis does indeed recover
the Lindblad dynamics, the same is true of the SSE parameter dynamics as
expected.

In figure 2 we show the Wigner functions of the quantum-jump trajectory for
the same realisation as in figure 1 at three selected times, where the central
trajectory up to the respective time is depicted as a solid white line. The
damped Gaussian motion is visible in the first figure, just before the first
jump in this realisation. The remaining figures illustrate the effect of the
quantum jumps, resulting in a sudden displacement of the centre as well as the
expected deviations from a Gaussian state. At the jump, the state is acted on
by the position operator and the resulting state is no longer Gaussian. We
clearly observe interference patterns in the Wigner function corresponding to
the excitation of higher harmonic oscillator states. It is remarkable that
averaging over these non-classical excited states results in the same Gaussian
state as the Lindblad equation.

Figure 3: Hagedorn basis parameter evolution for the position measurement
model eq. 83. The figure on the left shows the norm $N(t)$ of the ground
state, and the right figure the evolution of the parameter $M(t)$ with dotted
lines depicting the asymptotic fixed values of $M(t)$.

The coefficients of the state in the moving Hagedorn basis are depicted in the
histograms in the bottom panel of figure 2, for the same times as in the upper
panel. We have chosen the initial ground state of the Hagedorn basis to
coincide with the initial state, and thus there is no contribution from higher
states before the first jump. We observe how each jump leads to a contribution
from the next higher basis states, as expected.

Figure 4: Wigner function for the first four Hagedorn states
$\ket{n,a_{t},z_{t}}$ (n from 0-4 from left to right) for the position
measurement model eq. 83, at times $t=0,5,10$ (top to bottom)

Analytically the evolution of the basis states is solely described by the
linearized flow in equation (89). We may use this expression together with
eqs. 54 and 62 to calculate $N(t)$ and $M(t)$. Note that the dynamics of
$N(t)$ and $M(t)$ do not depend on the initial position of the wavepacket, but
solely on the values of $\omega$ and $\gamma$. Using the asymptotic behaviour
of $S(t)$ we find that $N(t)$ for long times follows a simple exponential
decay

$N(t)\to
N_{\infty}e^{-\sqrt{\frac{\omega(\lambda-\omega)}{2}}t},\quad\text{with}\quad
N_{\infty}=\sqrt{\frac{8\zeta\lambda\omega}{2\zeta\omega(\lambda+\omega)+\left(\zeta^{2}\omega+\lambda\right)\sqrt{2\omega(\lambda+\omega)}}}.$
(95)

$M(t)$ on the other hand tends to the fixed value

$M(t)\to\frac{\lambda-\zeta^{2}\omega-i\zeta\sqrt{2\omega(\lambda-\omega)}}{\lambda+\zeta^{2}\omega+\zeta\sqrt{2\omega(\lambda+\omega)}}.$
(96)

The dynamics of $N(t)$ and $M(t)$ for the example considered here are depicted
in figure 3. We observe the expected decay in the ground state normalisation
and the approach to the fixed point value in $M(t)$.

The Hagedorn raising operator $\hat{A}^{\dagger}(a_{t},z_{t})$ is also
determined by $S(t)$ and in the long time limit it takes the form

$\hat{A}^{\dagger}\to\frac{N_{\infty}\rme^{-i\sqrt{\frac{\omega(\lambda+\omega)}{2}}t}}{\sqrt{2\hbar}}\left(\left(\sqrt{\zeta}+\frac{i}{\sqrt{\zeta\omega}}\Lambda^{*}\right)(\hat{x}-x_{t})-i\left(\frac{1}{\sqrt{\zeta}}-i\sqrt{\frac{\zeta\omega}{2\lambda^{2}}}\Lambda\right)(\hat{p}-p_{t})\right).$
(97)

Which we recognise as squeezed and shifted harmonic oscillator creation
operator with a rotating phase of angular frequency
$\sqrt{\tfrac{\omega(\lambda+\omega)}{2}}$.

Figure 5: Wigner function for the first four time evolved Hagedorn states
$\hat{U}(t)\ket{n,a_{0},z_{0}}$ as defined in eq. 60, for the position
measurement model eq. 83 (n from 0-4 from left to right) at times $t=0,5,10$
(top to bottom)

These quantities determine the time evolution of the Hagedorn basis states
$\ket{n,a_{t},z_{t}}$, the first four of which are depicted for different
times in fig. 4. In fig. 4 we can see how the entire Hagedorn basis is
squeezed and shifted uniformly with $a_{t}$ and $z_{t}$. In contrast in fig. 5
we see how the propagated initial states $\hat{U}(t)\ket{n,a_{0},z_{0}}$
behave differently with each state decaying towards a squeezed state at the
origin, albeit at different rates. Note that the evolution of the Hagedorn
basis states does not depend on the specific realisation of the quantum-jump
trajectory. What differs between different quantum-jump realisations are the
coefficients of the state in this basis. For a given initial state they remain
constant between the jumps and are updated at each jump according to equation
(82). The coefficients for our example just before and just after the first
jump, and just after the second jump are depicted in the bottom row of figure
2.

Figure 6: Lindblad dynamics (left) compared with a single trajectory of the
SSE (middle) and quantum jump (right) for the position measurement model eq.
83, with parameters $\omega=1$ and $\gamma=0.2$. The initial Gaussian is a
squeezed state (i.e. $a_{0}=(\tfrac{1}{\sqrt{2}},i\sqrt{2})^{T}$ or
$G=\left(\begin{smallmatrix}2&0\\\
0&\nicefrac{{1}}{{2}}\end{smallmatrix}\right)$), centered at
$z_{0}=\left(2,0\right)^{T}$. In each case, a snapshot of the Wigner function
at $t=10$ is plotted in phase space, with a white line displaying the
precedent central motion.

To summarise, we observe clear differences between the Lindblad, SSE and
quantum-jump dynamics for the harmonic oscillator with position measurement,
which can be understood to a large degree using the analytical treatment
developed in the previous sections. Figure 6 shows the Wigner functions of the
state at $t=10$ for the three different realisations together with the central
trajectory up until this time. In the right panel, corresponding to the
Lindblad dynamics the state remains Gaussian, and its central motion follows
the usual unitary harmonic oscillator trajectory. The increased uncertainties
in position and momentum lead to the broadening of the Gaussian apparent here.
The SSE dynamics in the central picture, on the other hand, also remains
Gaussian in shape and stays well localised as predicted by the dynamical
behaviour of $G$. The central trajectory performs a Brownian motion around the
harmonic oscillator trajectory. Finally, the quantum-jump trajectory performs
smooth stretches of damped harmonic oscillations interrupted by discrete
jumps, and crucially, the state does not remain Gaussian.

###  5.2 Example 2: Damped harmonic oscillator

As a second example, we consider again a harmonic oscillator Hamiltonian and
the non-Hermitian Lindbladian

$\hat{L}=\sqrt{\frac{\gamma}{2}}(\hat{x}+i\hat{p}).$ (98)

For an initially coherent state, the resulting dynamics are rather trivial, as
such a state is an eigenstate of both the Hamiltonian and the Lindbladian. As
a result, for both the quantum-jump and SSE the stochastic terms vanish and
the dynamics of all three descriptions are the same, simply transporting the
initial state along the trajectories of the damped oscillator. Using a
squeezed initial state instead, the motion becomes more interesting and
Lindblad, SSE and quantum-jump trajectories differ.

Explicitly the dynamical equations for the Lindblad evolution eq. 14 become

$\displaystyle\frac{dz_{t}}{dt}$
$\displaystyle=(\omega\Omega-\frac{\gamma}{2})z_{t},$ (99)
$\displaystyle\frac{dG}{dt}$ $\displaystyle=\omega(\Omega
G-G\Omega)+\gamma(G-G^{2}).$

The equation for the central dynamics is simply that of a damped oscillator.
It can be written in terms of a second-order differential equation for
$\tilde{x}$ taking the more familiar form

$\frac{d^{2}x_{t}}{dt^{2}}+\frac{\gamma}{2}\frac{dx_{t}}{dt}+\omega x_{t}=0.$
(100)

Thus, the central dynamics are given by

$\begin{pmatrix}x_{t}\\\ p_{t}\end{pmatrix}=e^{-\frac{\gamma
t}{2}}\begin{pmatrix}\cos(\omega t)&\sin(\omega t)\\\ -\sin(\omega
t)&\cos(\omega t)\end{pmatrix}\begin{pmatrix}x_{0}\\\ p_{0}\end{pmatrix}.$
(101)

In the equation for the covariance matrix $G$ we immediately observe the fixed
point for $G=\mathbb{I}$, corresponding to a coherent state, which is indeed
approached asymptotically by any initial Gaussian state. For an initially
squeezed state with $G(0)=\big{(}\begin{smallmatrix}\zeta&&0\\\
0&&\nicefrac{{1}}{{\zeta}}\end{smallmatrix}\big{)}$ we find the time-dependent
covariances

$\Sigma(t)=\frac{\hbar}{2}\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}+\frac{\hbar
e^{-\gamma
t}}{4\zeta}\begin{pmatrix}(\zeta-1)^{2}-\left(\zeta^{2}-1\right)\cos(2\omega
t)&\left(\zeta^{2}-1\right)\sin(2\omega t)\\\
\left(\zeta^{2}-1\right)\sin(2\omega
t)&(\zeta-1)^{2}+\left(\zeta^{2}-1\right)\cos(2\omega t).\\\ \end{pmatrix}$
(102)

That is, we observe the usual oscillations with frequency $2\omega$ around the
coherent state covariances, which are now damped, and $\Sigma(t)$
asymptotically approaches $\Sigma(t)\to\frac{\hbar}{2}\mathbb{I}$, while
$\tilde{z}\to(0,0)^{T}$, resulting in a coherent state at the origin. This
behaviour is clearly visible in figure 7, in which the expectation values and
uncertainties of position and momentum in the Lindblad case are depicted as
black solid lines for an initial state with $\zeta=2$ centred at
$\tilde{z}_{0}=(2,0)^{T}$.

Figure 7:  Lindblad dynamics (solid black line) compared with single
trajectories of the SSE (dashed red line) and quantum-jump method (dotted blue
line) for the damped oscillator model eq. 98 with $\omega=1$ and $\gamma=0.2$.
The initial Gaussian is a squeezed state (i.e.
$a_{0}=(\tfrac{1}{\sqrt{2}},i\sqrt{2})^{T}$ or
$G=\left(\begin{smallmatrix}2&0\\\
0&\nicefrac{{1}}{{2}}\end{smallmatrix}\right)$), centered at
$z_{0}=\left(2,0\right)^{T}$. We show the time dependence of the position
expectation $\braket{\hat{x}}$ (top left), the momentum expectation
$\braket{\hat{p}}$ (top right), the positional variance $\Delta x^{2}$ (bottom
left) and the momentum variance $\Delta p^{2}$ (bottom right).

The dynamical eqs. 32 and 33 in the SSE case become

$\displaystyle\frac{dz_{t}}{dt}$
$\displaystyle=(\omega\Omega-\frac{\gamma}{2})z_{t}dt+\frac{\sqrt{\hbar\gamma}}{2}(G^{-1}-\mathbb{I})\begin{pmatrix}1\\\
0\end{pmatrix}d\xi_{R}-\frac{\sqrt{\hbar\gamma}}{2}(G^{-1}-\mathbb{I})\begin{pmatrix}0\\\
1\end{pmatrix}d\xi_{I},$ (103) $\displaystyle\frac{dG}{dt}$
$\displaystyle=\omega(\Omega G-G\Omega)+\frac{\gamma}{2}(\mathbb{I}-G^{2}).$

The deterministic part of the SSE and the Lindblad central dynamics agree as
usual. The additional drift term is proportional to $G^{-1}-\mathbb{I}$, and
vanishes for the coherent state $G=\mathbb{I}$. Although the evolution of $G$
differs from the Lindblad case, it too asymptotically approaches the fixed
point $G=\mathbb{I}$, independent of the initial conditions. Thus, the
stochastic contributions to the central trajectory become smaller over time,
and the dynamics drive the state towards a coherent state at the origin, just
as in the Lindblad case.

The complexified linearised flow of the damped oscillator model eq. 98 is
explicitly given by

$S(t)=\begin{pmatrix}\cosh\left(\frac{\gamma+2i\omega}{2}t\right)&-i\sinh\left(\frac{\gamma+2i\omega}{2}t\right)\\\
i\sinh\left(\frac{\gamma+2i\omega}{2}t\right)&\cosh\left(\frac{\gamma+2i\omega}{2}t\right)\\\
\end{pmatrix},$ (104)

For an initial squeezed state with $G(0)=\big{(}\begin{smallmatrix}\zeta&&0\\\
0&&\nicefrac{{1}}{{\zeta}}\end{smallmatrix}\big{)}$ this yields the time
dependent covariances in the SSE dynamics

$\Sigma(t)=\frac{\hbar}{2f(t)}\bigg{(}\begin{smallmatrix}(\zeta^{2}+1)\cosh(\gamma
t)+2\zeta\sinh(\gamma t)-(\zeta^{2}-1)\cos(2\omega
t)&(\zeta^{2}-1)\sin(2\omega t)\\\ (\zeta^{2}-1)\sin(2\omega
t)&(\zeta^{2}+1)\cosh(\gamma t)+2\zeta\sinh(\gamma
t)+(\zeta^{2}-1)\cos(2\omega t)\\\ \end{smallmatrix}\bigg{)},$ (105)

with

$f(t)=\left(\zeta^{2}+1\right)\sinh(\gamma t)+2\zeta\cosh(\gamma t).$ (106)

For long times $t\to\infty$ this behaves as

$\Sigma(t)=\frac{\hbar}{2}\begin{pmatrix}1&0\\\
0&1\end{pmatrix}+\frac{\hbar(\zeta-1)}{2(1+\zeta)}e^{-\gamma
t}\begin{pmatrix}-\cos(2\omega t)&\sin(2\omega t)\\\ \sin(2\omega
t)&\cos(2\omega t)\\\ \end{pmatrix},$ (107)

with $\Sigma(t)$ approaching $\frac{\hbar}{2}\mathbb{I}$ just as in the
Lindblad case. The example in figure 7 nicely demonstrates how the SSE and
Lindblad dynamics approach the same limit in different ways. We also observe
the reduced amplitude of the noise in the central dynamics, as the matrix $G$
approaches the identity.

The quantum jump trajectory in this example again initially follows the
Gaussian dynamics, here given by

$\frac{dz_{t}}{dt}=\left(\omega\Omega-\frac{\gamma}{\hbar}\Sigma(t)\right)z_{t},$
(108)

where $G(t)$ evolves dynamically as in the SSE case. This central dynamics is
also a type of damped harmonic oscillator, however, it differs from the
Lindblad case, both in the presence of a damping term in both position and
momentum, and in the modulation of the damping induced by the time-dependent
covariance matrix. In the example in figure 7, however, we observe that this
is a quantitative rather than qualitative difference in this case.

Figure 8: Top row: Wigner functions of the single quantum jump trajectory in
figure 7 at selected times (from left to right: t=0.84 (shortly before the
first jump) and t=0.86 (just after the first jump) and t=1.48 (just after the
second jump)). The white line traces the preceding central motion. Bottom row:
Relative magnitudes of the coefficients of the state in the time evolved basis
($\hat{U}(t)\ket{n,a_{0},z_{0}}$) at the same times as in the top row.

Figure 9: Hagedorn basis parameter evolution for the damped oscillator model
eq. 98. The figure on the left shows the norm $N(t)$ of the ground state, and
the right figure the evolution of the parameter $M(t)$ with dotted lines
depicting the asymptotic fixed values of $M(t)$.

Using eqs. 104, 54 and 62 we can derive the Hagedorn basis parameters

$\displaystyle N(t)$
$\displaystyle=\sqrt{\frac{2\zeta}{\left(\zeta^{2}+1\right)\sinh\left(\gamma
t\right)+2\zeta\cosh\left(\gamma t\right)}},$ (109) $\displaystyle M(t)$
$\displaystyle=-\frac{\left(\zeta^{2}-1\right)\sinh(\gamma
t)}{\left(\zeta^{2}+1\right)\sinh(\gamma t)+2\zeta\cosh(\gamma t)}.$

Much like the position example in the long time limit $N(t)$ tends to a simple
exponential decay whilst $M(t)$ tends to a fixed value

$\displaystyle N(t)$
$\displaystyle\to\frac{2\sqrt{\xi}}{1+\xi}e^{-\frac{\gamma t}{2}}$ (110)
$\displaystyle M(t)$ $\displaystyle\to\frac{1-\gamma}{1+\gamma}.$

Figure 8 illustrates the effect of the first two jumps in the particular
realisation of figure 7, showing the Wigner distributions shortly before the
first jump and shortly after the first and second jumps, as well as the
corresponding coefficients in the Hagedorn basis.

Figure 10: Lindblad dynamics compared with a single quantum jump and SSE
trajectory for the damped oscillator model eq. 98. In each case, a snapshot of
the Wigner function at $t=10$ is plotted in phase space, with a white line
displaying the precedent central motion. The left plot corresponds to the
Lindblad dynamics, the middle one to the SSE and the right plot to the quantum
jump dynamics.

Figure 10 depicts the Wigner functions of the Lindblad, SSE and quantum-jump
dynamics at time $t=10$, as well as the corresponding central trajectory up to
this time. In comparison with the previous example, we now see dissipation for
all three dynamics. The effect of the diminishing stochastic contribution in
the SSE dynamics, specific to the Lindbladian considered here, is also visible
in the figure. While the final quantum-jump state is still non-Gaussian, every
individual run asymptotically approaches a coherent state in the centre, in
contrast to the dynamics resulting from the first example.

##  6 Summary and outlook

We have investigated the dynamics of initially Gaussian states in open quantum
systems described in a Lindblad formalism with a quadratic Hamiltonian and
linear Lindbladian in comparison to the two popular unravellings of the
dynamics given by SSE and quantum-jump trajectories. In the SSE case, the
state remains Gaussian for all times, where the central dynamics has a
possibly damped deterministic contribution accompanied by a stochastic term
that depends on the covariances of the state. The covariances themselves
follow a deterministic time evolution, that is independent of the central
trajectory, and that coincides with the dynamical equation found for quantum
evolution generated by an effective non-Hermitian Hamiltonian. For the
quantum-jump approach, initial Gaussian states do in general not stay Gaussian
over time. Applying results from [20] we have formulated a method to describe
quantum-jump trajectories utilising a family of solutions to the non-Hermitian
Schrödinger equation that depends only on the dynamics of a $2\times 2$
complex matrix known as the linearised flow. We have studied the similarities
and differences of the dynamics resulting in the SSE, the quantum jump and the
Lindblad descriptions for two important examples.

##  7 Acknowledgements

We would like to thank Bradley Longstaff for his time checking calculations of
the SSE Gaussian limit.

We acknowledge support from the Royal Society (Grants. No. URF\R\201034 and
RGF\EA\180169) and from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (grant agreement No
758453).

##  Appendix A Quadratic Propagation

In this appendix we will discuss how to derive an expression for time evolved
operators of the form
$e^{-\frac{i}{\hbar}t\hat{K}}\hat{A}(b,\chi)\rme^{\frac{i}{\hbar}t\hat{K}}$
where $\hat{K}$ is a possibly non-Hermitian operator of order at most two in
$\hat{p}$ and $\hat{q}$. We will assume that the corresponding classical
Hamiltonian is of the form

$K(z)=\frac{1}{2}z\cdot K_{2}z+k_{1}\cdot\Omega z+k_{0}$ (111)

where $K_{2}$ is a symmetric matrix, $k_{1}$ is a vector and $k_{0}$ a scalar.
The solution $z(t)$ to Hamilton’s equations generated by $K(z)$ with initial
conditions $z(0)=z$ is given by

$\Phi(t,z)=S(t)z+v(t)\,\,,\quad\text{with}\quad
v(t)=\int_{0}^{t}S(t-s)k_{1}\,ds\,\,,$ (112)

where $S(t)=e^{t\Omega K_{2}}$ is symplectic, i.e., $S^{T}\Omega S=\Omega$.

Now we consider a general operator $\hat{P}$ and
$\hat{P}_{t}:=\rme^{-\frac{i}{\hbar}t\hat{K}}\hat{P}\rme^{\frac{i}{\hbar}t\hat{K}}$,
which satisfies the differential equation

$i\hbar\partial_{t}\hat{P}_{t}=[\hat{K},\hat{P}_{t}]\,\,.$ (113)

In the Weyl representation this equation takes the form

$\partial_{t}P_{t}(z)=\\{K(z),P_{t}(z)\\}\,\,,$ (114)

and here we used that $K(z)$ is at most quadratic, otherwise there are higher
order terms on the right hand side from the Moyal bracket. This equation is up
to a sign equal to the Liouville equation in classical mechanics which is
solved in terms of the classical flow, and hence we find

$P_{t}(z)=P(\Phi(-t,z))\,\,.$ (115)

This is a special case of a much more general result which is known as
Egorov’s Theorem, the only new aspect here is that we allowed a complex
Hamilton function $K(z)$.

If we apply this result to the case that
$P=A(b,\chi)=\frac{i}{\sqrt{2\hbar}}b\cdot\Omega(z-\chi)$ we find with
$\Phi(-t,z)-\chi=S(-t)(z-\Phi(t,\chi))$ and $\Omega S(-t)=S^{T}(t)\Omega$,
which follows from $S^{T}\Omega S=\Omega$ and $S(-t)=S(t)^{-1}$, that

$P_{t}(z)=A(S(t)b,\Phi(t,\chi))$ (116)

or

$\hat{U}(t)\hat{A}(b,\chi)\hat{U}(-t)=\hat{A}\big{(}S(t)b,\Phi(t,\chi)\big{)}\,\,.$
(117)

##  Appendix B A non-commutative Binomial expansion

In this appendix we present the technical details of the expansion (64), we
assume we have a pair of creation and annihilation operators
$\hat{A}^{\dagger},\hat{A}$, and a corresponding orthonormal basis
$\ket{n}=\frac{1}{\sqrt{n!}}\,[\hat{A}^{\dagger}]^{n}\ket{0}$ such that
$\hat{A}\ket{n}=\sqrt{n}\,\ket{n-1}$ and
$\hat{A}^{\dagger}\ket{n}=\sqrt{n+1}\ket{n+1}$, and in particular
$\hat{A}\ket{0}=0$.

We want to find the expansion coefficients $B_{mn}$ in

$\frac{1}{\sqrt{n!}}(h_{+}\hat{A}^{\dagger}+h_{-}\hat{A}+h_{0})^{n}\ket{0}=\sum_{m=0}^{n}B_{mn}\ket{m}\,\,,$
(118)

where $h_{+},h_{-},h_{0}$ are complex constants. The idea is to use that
$(h_{+}\hat{A}^{\dagger}+h_{-}\hat{A}+h_{0})^{n}$ is proportional to the $n$th
therm in the expansion of

$e^{s(h_{+}\hat{A}^{\dagger}+h_{-}\hat{A}+h_{0})}\ket{0}=\sum_{n=0}^{\infty}\frac{s^{n}}{n!}\,(h_{+}\hat{A}^{\dagger}+h_{-}\hat{A}+h_{0})^{n}\ket{0}\,\,,$
(119)

and to use the Baker Campbell Haussdorf formula

$e^{s(h_{+}\hat{A}^{\dagger}+h_{-}\hat{A}+h_{0})}=e^{sh_{+}\hat{A}^{\dagger}}e^{s(h_{-}\hat{A}+h_{0})}e^{-\frac{s^{2}}{2}[h_{+}\hat{A}^{\dagger},h_{-}\hat{A}+h_{0}]}\,\,.$
(120)

With $[\hat{A}^{\dagger},\hat{A}]=-1$ we obtain for the commutator
$[h_{+}\hat{A}^{\dagger},h_{-}\hat{A}+h_{0}]=-h_{+}h_{-}$ and expanding the
exponentials and equating the terms of order $s^{n}$ on both sides gives

$(h_{+}\hat{A}^{\dagger}+h_{-}\hat{A}+h_{0})^{n}=n!\sum_{m+l+2k=n}\frac{1}{m!l!k!2^{k}}(h_{+}h_{-})^{k}(h_{+}\hat{A}^{\dagger})^{m}(h_{-}\hat{A}+h_{0})^{l}\,\,.$
(121)

We can now use that $(h_{-}\hat{A}+h_{0})^{l}\ket{0}=h_{0}^{l}\ket{0}$ and
$(h_{+}A^{\dagger})^{m}\ket{0}=\sqrt{m!}\,h_{+}^{m}\ket{m}$ to find

$\frac{1}{\sqrt{n!}}(h_{+}\hat{A}^{\dagger}+h_{-}\hat{A}+h_{0})^{n}\ket{0}=\sum_{m=0}^{n}\sum_{l+2k=n-m}\frac{\sqrt{n!}(h_{+}h_{-})^{k}h_{+}^{m}h_{0}^{l}}{\sqrt{m!}l!k!2^{k}}\,\,\ket{m}$
(122)

and the inner sum can be rewritten as

$B_{mn}=h_{+}^{m}h_{0}^{n-m}\sum_{k=0}^{[\frac{n-m}{2}]}\frac{\sqrt{n!}}{\sqrt{m!}(n-m-2k)!k!2^{k}}\bigg{(}\frac{h_{+}h_{-}}{h_{0}^{2}}\bigg{)}^{k}\,\,.$
(123)

## References

## References

  * [1] H. M. Wiseman and G. J. Milburn, _Quantum measurement and control_. Cambridge university press, 2009.
  * [2] V. P. Belavkin, “Nondemolition measurements, nonlinear filtering and dynamic programming of quantum stochastic processes,” in _Modeling and Control of Systems_. Springer, 1989, pp. 245–265.
  * [3] H. Carmichael, _An open systems approach to quantum optics: lectures presented at the Université Libre de Bruxelles, October 28 to November 4, 1991_. Springer Science & Business Media, 2009, vol. 18.
  * [4] N. Bartolo, F. Minganti, J. Lolli, and C. Ciuti, “Homodyne versus photon-counting quantum trajectories for dissipative kerr resonators with two-photon driving,” _The European Physical Journal Special Topics_ , vol. 226, no. 12, pp. 2705–2713, 2017.
  * [5] J. K. Eastman, J. J. Hope, and A. R. Carvalho, “Tuning quantum measurements to control chaos,” _Scientific reports_ , vol. 7, no. 1, pp. 1–10, 2017.
  * [6] H. Wiseman, “Quantum trajectories and quantum measurement theory,” _Quantum and Semiclassical Optics: Journal of the European Optical Society Part B_ , vol. 8, no. 1, p. 205, 1996.
  * [7] W. Verstraelen and M. Wouters, “Gaussian quantum trajectories for the variational simulation of open quantum-optical systems,” _Applied Sciences_ , vol. 8, no. 9, p. 1427, 2018.
  * [8] M. B. Plenio and P. L. Knight, “The quantum-jump approach to dissipative dynamics in quantum optics,” _Reviews of Modern Physics_ , vol. 70, no. 1, p. 101, 1998.
  * [9] G. C. Ghirardi, A. Rimini, and T. Weber, “Unified dynamics for microscopic and macroscopic systems,” _Physical review D_ , vol. 34, no. 2, p. 470, 1986\.
  * [10] A. Bassi, K. Lochan, S. Satin, T. P. Singh, and H. Ulbricht, “Models of wave-function collapse, underlying theories, and experimental tests,” _Reviews of Modern Physics_ , vol. 85, no. 2, p. 471, 2013.
  * [11] D. C. Brody and L. P. Hughston, “Efficient simulation of quantum state reduction,” _Journal of Mathematical Physics_ , vol. 43, no. 11, pp. 5254–5261, 2002.
  * [12] N. Gisin and I. C. Percival, “The quantum-state diffusion model applied to open systems,” _Journal of Physics A: Mathematical and General_ , vol. 25, no. 21, p. 5677, 1992.
  * [13] O. Brodier and A. Ozorio de Almeida, “Symplectic evolution of wigner functions in markovian open systems,” _Phys. Rev. E_ , vol. 69, p. 016204, 2004.
  * [14] E. Graefe, B. Longstaff, T. Plastow, and R. Schubert, “Lindblad dynamics of Gaussian states and their superpositions in the semiclassical limit,” _Journal of Physics A: Mathematical and Theoretical_ , vol. 51, no. 36, p. 365203, 2018.
  * [15] E. J. Heller, “Time-dependent approach to semiclassical dynamics,” _The Journal of Chemical Physics_ , vol. 62, no. 4, pp. 1544–1555, 1975.
  * [16] R. G. Littlejohn, “The semiclassical evolution of wave packets,” _Physics reports_ , vol. 138, no. 4-5, pp. 193–291, 1986.
  * [17] W. T. Strunz and I. C. Percival, “Classical mechanics from quantum state diffusion-a phase-space approach,” _Journal of Physics A: Mathematical and General_ , vol. 31, no. 7, p. 1801, 1998.
  * [18] J. Halliwell and A. Zoupas, “Quantum state diffusion, density matrix diagonalization, and decoherent histories: A model,” _Physical Review D_ , vol. 52, no. 12, p. 7294, 1995.
  * [19] G. A. Hagedorn, “Raising and lowering operators for semiclassical wave packets,” _Annals of Physics_ , vol. 269, no. 1, pp. 77–104, 1998.
  * [20] C. Lasser, R. Schubert, and S. Troppmann, “Non-hermitian propagation of hagedorn wavepackets,” _Journal of Mathematical Physics_ , vol. 59, no. 8, p. 082102, 2018.
  * [21] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, “Completely positive dynamical semigroups of n-level systems,” _Journal of Mathematical Physics_ , vol. 17, no. 5, pp. 821–825, 1976.
  * [22] G. Lindblad, “On the generators of quantum dynamical semigroups,” _Communications in Mathematical Physics_ , vol. 48, no. 2, pp. 119–130, 1976\.
  * [23] H.-P. Breuer, F. Petruccione _et al._ , _The theory of open quantum systems_. Oxford University Press on Demand, 2002.
  * [24] A. J. Daley, “Quantum trajectories and open many-body quantum systems,” _Advances in Physics_ , vol. 63, no. 2, pp. 77–149, 2014.
  * [25] F. Klauck, L. Teuber, M. Ornigotti, M. Heinrich, S. Scheel, and A. Szameit, “Observation of pt-symmetric quantum interference,” _Nature Photonics_ , vol. 13, no. 12, pp. 883–887, 2019.
  * [26] E.-M. Graefe, “PT symmetry dips into two-photon interference,” _Nature Photonics_ , vol. 13, no. 12, pp. 822–823, 2019.
  * [27] E.-M. Graefe and R. Schubert, “Complexified coherent states and quantum evolution with non-hermitian hamiltonians,” _Journal of Physics A: Mathematical and Theoretical_ , vol. 45, no. 24, p. 244033, 2012.
  * [28] R. L. Hudson, “When is the Wigner quasi-probability density non-negative?” _Reports on Mathematical Physics_ , vol. 6, no. 2, pp. 249–252, 1974.
  * [29] A. Kenfack and K. Życzkowski, “Negativity of the Wigner function as an indicator of non-classicality,” _Journal of Optics B: Quantum and Semiclassical Optics_ , vol. 6, no. 10, p. 396, 2004.
  * [30] P. de M Rios and A. O. de Almeida, “On the propagation of semiclassical wigner functions,” _Journal of Physics A: Mathematical and General_ , vol. 35, no. 11, p. 2609, 2002.
  * [31] W. T. Strunz, “Path integral, semiclassical and stochastic propagators for markovian open quantum systems,” _Journal of Physics A: Mathematical and General_ , vol. 30, no. 11, p. 4053, 1997.
  * [32] W. T. Strunz, L. Diósi, N. Gisin, and T. Yu, “Quantum trajectories for brownian motion,” _Physical Review Letters_ , vol. 83, no. 24, p. 4909, 1999\.
  * [33] G. Adesso, S. Ragy, and A. R. Lee, “Continuous variable quantum information: Gaussian states and beyond,” _Open Systems & Information Dynamics_, vol. 21, no. 01n02, p. 1440001, 2014.
  * [34] T. Plastow, “Semiclassical methods for investigating open quantum systems and decoherence,” Ph.D. dissertation, University of Bristol, 2020.
  * [35] B. Oksendal, _Stochastic differential equations: an introduction with applications_. Springer Science & Business Media, 2013.
  * [36] C. Lubich, _From quantum to classical molecular dynamics: reduced models and numerical analysis_ , ser. Zurich Lectures in Advanced Mathematics. European Mathematical Society (EMS), Zürich, 2008.
  * [37] C. Lasser and C. Lubich, “Computing quantum dynamics in the semiclassical regime,” _Acta Numer._ , vol. 29, pp. 229–401, 2020.
  * [38] E.-M. Graefe and R. Schubert, “Wave-packet evolution in non-Hermitian quantum systems,” _Physical Review A_ , vol. 83, no. 6, p. 060101, 2011.
  * [39] D. C. Brody, “Biorthogonal quantum mechanics,” _Journal of Physics A: Mathematical and Theoretical_ , vol. 47, no. 3, p. 035305, 2013.