Low dimensional manifolds for exact representation of open quantum systems

TL;DR

This paper introduces methods to efficiently approximate open quantum systems by identifying low-dimensional manifolds of quasi-classical states, optimizing transformations to minimize quantum relative entropy, and enabling simpler numerical simulations.

Contribution

It develops a family of analytical and computational techniques for deriving optimal unitary transformations based on Lie group representations to simplify open quantum system descriptions.

Findings

Transforms minimize quantum relative entropy to quasi-classical states

Coupled equations enable efficient numerical simulation

Quantifies non-classicality and compares system complexities

Abstract

Weakly nonlinear degrees of freedom in dissipative quantum systems tend to localize near manifolds of quasi-classical states. We present a family of analytical and computational methods for deriving optimal unitary model transformations based on representations of finite dimensional Lie groups. The transformations are optimal in that they minimize the quantum relative entropy distance between a given state and the quasi-classical manifold. This naturally splits the description of quantum states into quasi-classical coordinates that specify the nearest quasi-classical state and a transformed quantum state that can be represented in fewer basis levels. We derive coupled equations of motion for the coordinates and the transformed state and demonstrate how this can be exploited for efficient numerical simulation. Our optimization objective naturally quantifies the non-classicality of states occurring in some given open system dynamics. This allows us to compare the intrinsic complexity of different open quantum systems.

Equations

$${\chi }_{\beta ,\theta }=U_{\theta }\underset{{\chi }_{\beta }}{\underbrace{Z(\beta {)}^{-1}{e}^{-\beta N}}}{U}_{\theta }^{†}$$

$$d\alpha$$

$$|d{\varphi }_t⟩$$

$$-\eta {⟨a⟩}_{\varphi }{a}^{†}|{\varphi }_t⟩dt+\eta {⟨a⟩}_{\varphi }^{\ast }a|{\varphi }_t⟩dt$$

$$-\sqrt{\kappa }{\overline{n}}_{\mathrm{th}}{a}^{†}|{\varphi }_t⟩dW+\sqrt{\kappa }({\overline{n}}_{\mathrm{th}}+1)a|{\varphi }_t⟩d{W}^{\ast }$$

$$d\alpha \to d\alpha -\sqrt{\kappa }ϵdt$$

$$H$$

$$\mathbf{L}$$

$$d{\alpha }_{\det }$$

$$\dot{Q}\approx -\left(\frac{\kappa }{2}-\chi +\frac{\beta }{2}{Q}^2\right)Q.$$

$$|{\psi }_t⟩=U_t|{\varphi }_t⟩\iff |{\varphi }_t⟩=U_t^{†}|{\psi }_t⟩.$$

$${\rho }_t=U_t{\sigma }_t{U}_t^{†}\iff {\sigma }_t=U_t^{†}{\rho }_t{U}_t.$$

$$F_j^{>}(\theta ):=i{U}_{\theta }^{†}\frac{\partial U_{\theta }}{\partial {\theta }^j}.$$

$${\theta }_t:={\theta }_{\ast }^{(\mathcal{J})}({\rho }_t)={\mathrm{argmin}}_{\theta }\mathcal{J}(\theta ;{\rho }_t).$$

$${\mathcal{J}}_M(\theta ;\rho ):=\underset{{⟨M⟩}_{U_{\theta }^{†}\rho U_{\theta }}}{\underbrace{{⟨U_{\theta }M{U}_{\theta }^{†}⟩}_{\rho }}}.$$

$$M(\beta )=\sum _{k=1}^e{\beta }^k{M}_k,{\beta }^k>0,k=1,2,\dots ,e.$$

$$N=\sum _{k=0}^{d}k|k⟩⟨k|,d\in \mathbb{N}\cup \{\infty \},$$

$$N(\beta ):=\sum _{k=1}^e{\beta }^k{N}_k,{\beta }^k>0,k=1,2,\dots ,e.$$

$$[N(\beta ),A_{l,\pm }]=\pm {\beta }^k{A}_{k,\pm }.$$

$$\left[\sum _{k=1}^e{\beta }^k{a}_k^{†}{a}_k,a_l^{†}{a}_m\right]$$

$$\left[\sum _{k=1}^e{\beta }^k{a}_k^{†}{a}_k,a_l^{†}{a}_m^{†}\right]$$

$$\left[\sum _{k=1}^e{\beta }^k{a}_k^{†}{a}_k,a_l{a}_m\right]$$

$$\mathcal{J}(\theta +\delta \theta ;\rho )$$

$$+\frac{1}{2}\sum _{j,k=1}^n{h}_{jk}(\theta ;\rho )\delta {\theta }^j\delta {\theta }^k+O\left(\delta {\theta }^3\right)$$

$$y_j(\theta ;\rho )$$

$$h_{jk}(\theta ;\rho )=h_{kj}(\theta ;\rho )$$

$$Y_j^{>}(\theta )$$

$$H_{jk}^{>}(\theta )$$

$${\mathbb{P}}_{\rho }\left[N_{\theta }>N_0\right]\le \frac{{⟨N_{\theta }⟩}_{\rho }}{N_0}.$$

$${\mathbb{P}}_{\rho }\left[N_{\theta }\ge N_0\right]={\mathbb{P}}_{\rho }\left[f(N_{\theta })\ge f(N_0)\right]\le \frac{{⟨f(N_{\theta })⟩}_{\rho }}{f(N_0)}.$$

$${\mathbb{P}}_{\rho }\left[N_{\theta }>N_0\right]$$

$$=e^{-\lambda N_0+{\mathcal{J}}_N^{\text{CGF}}}.$$

$${\mathcal{J}}_N^{\text{CGF}}(\theta ;\rho ,\lambda ):=\ln \left[\underset{{⟨e^{\lambda N}⟩}_{U_{\theta }^{†}\rho U_{\theta }}}{\underbrace{{⟨e^{\lambda N_{\theta }}⟩}_{\rho }}}\right],$$

$${\mathcal{J}}_N^{\text{CGF}}(\theta ;\rho ,\lambda )\approx \lambda {⟨N_{\theta }⟩}_{\rho }+\frac{{\lambda }^2}{2}\mathrm{var}{\left(N_{\theta }\right)}_{\rho }+O({\lambda }^3),$$

$$\frac{⟨n+1|U_{\theta }^{†}\rho U_{\theta }|n+1⟩}{⟨n|U_{\theta }^{†}\rho U_{\theta }|n⟩}\le \alpha ,$$

Full text

I Introduction

A given quantum mechanical system can be described in more than one way. Our choice of description is usually motivated by the insight it provides, its economy, its accuracy and, when dealing with sufficient complexity, the efficiency with which it can be numerically simulated and analyzed.

Closed quantum systems evolve unitarily, and if their Hamilton operator admits a sufficient set of individually tunable control operators one can – at least in principle – realize arbitrary unitary operations on the system’s Hilbert space Schirmer et al. (2001). This enables powerful quantum computing and quantum simulation schemes that derive their advantage over classical computers from the exponential scaling of the Hilbert space dimension with system size.

We may then ask, what are the implications for open quantum systems which exhibit dissipative dynamics? As dissipation increases it becomes increasingly difficult to use them for unitary quantum computing, but there also exist applications in quantum engineering that explicitly require coupling to input and output fields, ranging from quantum limited signal amplification Caves (1982); Mollow and Glauber (1967), via quantum key distribution to autonomously correcting quantum memories Kerckhoff et al. (2010). An important class of such applications can be described in the language of quantum feedback networks Gough and James (2008, 2009) or quantum input-output models Combes et al. (2017).

In general, our ability to design quantum systems for specific tasks is severely limited by the state space complexity. This is true of closed and open systems, but as fewer guarantees exist on what dynamics are achievable with open quantum systems than in the closed system case we are even more reliant on efficient numerical schemes. Fortunately, the dynamics of open quantum systems tend to exhibit phase space localization Percival (1999), which implies that there exist nonlinear sets within the system’s Hilbert space that act as attractors for the quantum dynamics.

A very appealing feature of quantum network models based exclusively on nonlinear oscillators is that they allow continuously tuning between the coherent qubit regime Mabuchi (2012) and the semi-classical weakly nonlinear limit. Several near-term applications exist for nonlinear oscillator networks: in the weakly nonlinear regime these include frameworks for photonic logic Mabuchi (2011); Santori et al. (2014), photonic Ising machines Wang et al. (2013); Inagaki et al. (2016) and for all-optical machine learning Tezak and Mabuchi (2015). In the strongly nonlinear regime novel quantum error correcting schemes have been proposed Mirrahimi et al. (2014) and implemented Leghtas et al. (2015). For all of these systems it is an extremeley interesting question to ask how their dynamics and capabilities change as the ratio of dissipation to non-linearity is increased or decreased, but no existing simulation method has allowed continuous interpolation between these regimes. Our approach remediates this by providing a framework that relies on exact quantum model transformations to exploit semi-classical localization and obtain more efficient system representations.

The key to our method is formulating the problem of finding efficient state parametrizations as an optimization problem. Specifically, we employ smoothly parametrized unitary transformations to represent states in co-moving adaptive bases. The description of, e.g., a pure quantum state $\left|\psi_{t}\right\rangle$ thus splits into the transformation coordinates $\theta_{t}$ and a residual quantum state $\left|\phi_{t}\right\rangle=U_{\theta_{t}}^{\dagger}\left|\psi_{t}\right\rangle$ that is represented in a localized basis. The diagram in Figure 1 intends to visualize this. We quantify localization in terms of complexity functionals that are shown to have precise information geometric meaning.

With this, we establish a self-consistent analytic framework that not only allows for numerical simulation but also serves to derive coupled dynamical equations for the semi-classical group manifold coordinates and the residual quantum state. We provide expressions for different diffusive Schrödinger equations as well as (stochastic) master equations. Although stochasticity is not at all a necessary requirement of our approach, we expect our method to have the greatest impact for systems that are continuously weakly observed by the environment and this naturally leads us to consider stochastic dynamics.

There already exist methods that rely on localization to analytically derive semi-classical stochastic dynamics Carter (1995) from quantum quasi-probability distributions Wigner (1932); Husimi (1940); Sudarshan (1963). These, however, require approximations that can lead to significant discrepancies with simulations based on a full quantum state description Kinsler and Drummond (1991); Santori et al. (2014). Furthermore, there does not exist a general approach to incrementally increase the accuracy of these methods.

Even prior to Percival’s analytic results on localization, Schack, Brun and Percival developed the MQSD simulation method Schack et al. (1995); Steimle et al. (1999) which allowed simulating stochastic Schrödinger equations in displaced Fock bases. It can be seen as a special case of our approach, but it does not produce analytical equations of motion for the transformation coordinates. The MQSD simulation method works well for strongly dissipative oscillator systems that exhibit nearly Gaussian states with low amounts of squeezing or multi-mode correlations, but it does not generalize to other degrees of freedom. Furthermore, MQSD is based on a heuristic that has lacked a rigorous interpretation prior to our framework.

There have also been several model-agnostic attempts to reduce the computational complexity of simulating open quantum systems: in Tezak et al. (2012) a photonic SR-latch in the weakly quantum regime was approximated by a finite state-space Markov chain that could be re-expressed as an effective quantum network model, and in Shi and Nurdin (2016) a conceptually very appealing approach was proposed based on principal component analysis of simulated quantum states. One important shortcoming of both of these methods, however, is that they require at least one simulation of the full system model which can be very expensive and in many cases impossible. Furthermore, they do not preserve important qualitative features of a quantum model, such as the commutators of system observables.

The contents of this paper are as follows: First, we motivate the technical details presented in subsequent sections by presenting numerical results of our method applied to two separate systems: (a) the Kerr-nonlinearity based SR-latch as first described in Mabuchi (2011) and discussed as an example in Tezak et al. (2012) and Shi and Nurdin (2016), and (b) the degenerate parametric oscillator (DPO) model which is a highly interesting system both in the semi-classical regime where it exhibits a pitchfork bifurcation and in the strongly nonlinear quantum regime, where the semi-classical fix-points limit to low decoherence cat-states that can be used to encode a qubit Mirrahimi et al. (2014).

Next, we formally introduce the underlying ‘state compression’ framework which relies on positive ‘penalty’-operators to quantify the complexity of quantum states. We show that their spectra can be directly related to achievable bounds on numerical truncation errors, and how this can be exploited for adapting our approach to situations in which Gaussian states poorly approximate the quantum dynamics.

We then describe how our method can be understood from an information theoretic point of view by relating our complexity functional to the quantum relative entropy. This allows to characterize the attractors near which the dynamics of a given system localize in terms of generalized Gibbs states generated by our specific choice of penalty operator which is smoothly transformed by our group.

We further present a small library of possible groups and penalty operators and finish by providing an outline of how this method can be extended to higher-dimensional Lie groups for which it is infeasible to derive the transformation differential in closed form.

There are many exciting future directions and valuable applications of our research, some of which we mention in this paper’s conclusion.

II Applications

We first present some numerical results of our method. The technical details will be described in the following sections. Our results were obtained using a custom software package QMANIFOLD Tezak (2016) that not only facilitates analytic derivations but also allows to directly carry out numerical simulations. It uses QuTiP’s Johansson et al. (2012) datastructures and interfaces with our existing package QNET Tezak and Goerz (2012) to automate various tedious symbolic calculations related to deriving the adjoint group representation and thus the derivation of the right generators $F_{j}^{>}(\theta_{t})$. Given a particular transformation $U_{\theta}$ and complexity functional (cf. Section III) it can compute the minimum complexity state $\left|\phi_{t}\right\rangle$ and coordinates $\theta_{t}$ for any input state $\left|\psi_{t}\right\rangle.$ Furthermore, given a dynamical open system model parametrized by its internal Hamiltonian $H$ and some dissipation operators $\mathbf{L}$ it can directly derive and simulate stochastic complex quantum diffusion dynamics.

II.1 Toy Model: An empty cavity

Consider a simple open cavity model described by a Hamiltonian $H=\hbar\omega a^{\dagger}a$ and a single dissipation operator $L=\sqrt{\kappa}a.$ Using the group of coherent displacements and linear excitation minimization (cf. Section VI.1) of the canonical counting operator $N=a^{\dagger}a$ we apply the method of gradient coupled fiducial state dynamics with the fiducial state given by a displaced Gibbs state

[TABLE]

with moving basis excitation $\,\left\langle N\right\rangle_{\chi_{\beta}}=\overline{n}_{\rm th}$. Summarizing the manifold coordinates in a single complex amplitude $\alpha=\frac{Q+iP}{\sqrt{2}}$, this coordinate and the reduced complexity state $\left|\phi_{t}\right\rangle$ have the joint equations of motion (cf Section V)

[TABLE]

where $\left|d\phi_{t}\right\rangle$ does not necessarily conserve the norm of $\left|\phi_{t}\right\rangle$ and where we have dropped a term $\propto\left|\phi_{t}\right\rangle dt$ that does not affect the evolution of any expectation values.

These dynamics can be simulated for any choice of $\overline{n}_{\rm th}$ which can also be turned into a dynamic variable (cf. Section V.6). For any finite choice of $\eta>0$ the mode operator expectation will fluctuate around zero $\,\left\langle a\right\rangle_{\phi_{t}}\approx 0$ and in the limit $\eta\to\infty$ we recover the hard constraint $\left.\,\left\langle a\right\rangle_{\phi_{t}}\right|_{\eta\to\infty}=0$. We see that $\alpha$ only couples to the input noise process if the moving basis excitation is non-zero. On the other hand, for $n_{\rm th}=0\Leftrightarrow\beta\to\infty$ we find that $\left|\phi_{t}\right\rangle=\left|0\right\rangle$ is a stable fix-point of the dynamics in accordance with our intuition about passive linear quantum systems. We could have equally well derived a (stochastic) master equation or a homodyne SSE using the formulas from Section V.

We can easily add a coherent displacement to this model by modifying $L\to L+\epsilon,H\to H+\frac{\sqrt{\kappa}}{2i}(\epsilon a^{\dagger}-\epsilon^{\ast}a)$. A straightforward calculation reveals that this linear displacement is fully absorbed into the dynamics of $\alpha$ such that relative to Equations (2),(3) we have

[TABLE]

while $\left|d\phi_{t}\right\rangle$ is left unchanged.

II.2 Kerr-cavity based NAND latch

Mabuchi recently proposed Mabuchi (2011) designing photonic logic gates by using nonlinear resonators in an interferometric feedback configuration. To achieve maximum power efficiency, devices could operate in a semi-classical regime where several tens of photons in a Kerr-resonator would cause the frequency to shift by one linewidth. A particularly interesting model is that of a photonic NAND-latch Mabuchi (2011) symmetrically constructed from two Kerr-resonators with mutual coherent feedback as visualized in Figure 2. The circuit is designed such that it has two metastable states with either the first or the second resonator in a high photon number state while the other resonator has low internal photon number.

A full quantum trajectory simulation of a Stochastic Schrödinger equation with $D=d^{2}=75^{2}=5625$ levels for this model with the same parameters as in Mabuchi (2011) and Tezak et al. (2012) is feasible on a typical workstation.

In Figure 3 we present the Fock level occupation probability $\mathbb{P}(n_{j}=n):=\,\left\langle\Pi^{(j)}_{n}\right\rangle$ for each oscillator mode $j=1,2$ and each Fock-level $n=0,1,\dots,d-1$ where $\Pi^{(j)}_{n}$ is the projection operator onto the $n$-th Fock state of oscillator $j$. We visualize this for both the original, static basis and for a coherently displaced representation, obtained by solving the linear excitation minimization problem for $N=a_{1}^{\dagger}a_{1}+a_{2}^{\dagger}a_{2}$, i.e., the total number of excitations present in the displaced frame. We find that we can represent the dynamics of the system in a basis set that is between at least four and ten times smaller depending on whether we dynamically update the basis size during simulation without sacrificing any precision in the simulation. If we further allow for modest simulation error, we can achieve up to an additional order of magnitude in basis reduction. Note that dimensional reduction for the same model was studied in Tezak et al. (2012); Shi and Nurdin (2016) and very similar physical models were studied in Santori et al. (2014) but we claim that ours is the only approach that allows to achieve strong dimensional reduction with a controlled error and without having to simulate the full system first.

II.3 The degenerate parametric oscillator

Another simple nonlinear extension of the empty cavity is given by the degenerate parametric oscillator (DPO). DPOs can exhibit very rich dynamics and have been long been employed for amplification Akhmanov et al. (1965) and generation of light at tunable wavelengths Giordmaine and Miller (1965) as well as more recently coherent optical Ising machines Utsunomiya et al. (2011). It can be physically realized by a resonant signal mode coupled to a strongly driven pump mode through a nonlinear parametric interaction that mediates the conversion of signal photon pairs to pump photons and vice versa. In the strongly nonlinear limit and for a low quality factor of the pump mode, the pump mode can be adiabatically eliminated yielding an open system model in a rotating frame parametrized by

[TABLE]

Using the same Ansatz as above but with $\overline{n}_{\rm th}=\eta=0$ for simplicity we find for the deterministic part of the coordinate dynamics

[TABLE]

For positive pump phase $\chi>0$ the dynamics are primarily captured by the evolution of the signal mode’s real quadrature $Q=(\alpha+\alpha^{\ast})/\sqrt{2}$ while the orthogonal $P$ quadrature is suppressed. We can then fix $P=0$ and consider the one dimensional differential equation valid in the classical limit

[TABLE]

Equation (8) is identical to the normal form of a pitchfork bifurcation up to some re-scaling. We visualize the bifurcation diagram in Figure 4.

The bifurcation exists for any non-zero two-photon loss rate $\beta>0,$ however the magnitude of $\beta$ strongly affects how non-classical (which in this context we take to mean non-Gaussian) the state of the signal mode becomes. When the system is pumped only slightly above the threshold, random switching or tunneling between the two equilibria is possible. We present such a trajectory in Figure 5.

For a constant steady state mode amplitude the switching rate strongly depends on the ratio of linear to two-photon loss. Additionally, in the strongly nonlinear case $\beta\geq\kappa,$ the system can spontaneously evolve into cat-like states that feature exhibit significant simultaneous overlap with coherent states centered at either equilibrium.

In the limit of vanishing linear loss $\kappa/\beta\to 0,$ the system has a decoherence-free sub-manifold spanned by the two equilibrium amplitude coherent states. In Mirrahimi et al. (2014) Mirrahimi et al. outline a scheme to encode quantum information in such a system. A detailed study of the switching dynamics was carried out in Kinsler and Drummond (1991).

In Figure 6 we compare how each basis level contributes to a whole trajectory of states when represented in the original fixed basis to one obtained via linear excitation minimization (cf. Section VI.1.1) using either a coherently displaced basis or a displaced and squeezed basis. We see that the adaptive schemes perform well in the case of strong linear dissipation but not so well in the case of strong two-photon loss.

We can understand this better by inspecting typical states that occur in each evolution. In Figure 7 we present snapshots of the signal mode’s Wigner function. For strong linear dissipation, the Wigner function of the signal mode typically appears quite Gaussian in shape, whereas in the strong two-photon loss case we see significant non-Gaussian features both in the transition states and when the mode is at one of the equilibria.

The bad performance of the excitation minimization functional in the non-Gaussian case is much improved by the cumulant generating function (CGF) minimization approach (cf. Section III.5) . In Figure 8 we compare the efficiency of the fixed basis with a coherently displaced basis where the coordinates are determined either by excitation minimization or by CGF minimization. We find that the CGF minimization (for $\lambda=3/2$) outperforms both the fixed basis and the excitation minimization method (which is equivalent to the QSD package’s approach). Here we have not even exploited the additional advantages that a displaced and squeezed basis combined with the CGF approach may yield.

Further details of the CGF minimized state representation can be found in Appendix D.

III Quantum state compression

In this section we discuss different options for quantifying the efficiency of representing a given state in a particular basis and we introduce a corresponding optimization problem. We will always assume that our adaptive Hilbert Space bases are related to the original fixed basis by means of a unitary transformation $U\in\mathcal{G}$ belonging to a connected Lie group that is locally generated by a finite dimensional Lie algebra $\mathfrak{g}$.

If our quantum state in the original, fixed basis is $\left|\psi_{t}\right\rangle$ we assume that it can be related to a reduced complexity state $\left|\phi_{t}\right\rangle$ via

[TABLE]

The generalization to mixed states is obvious, the fixed basis state $\rho_{t}$ and the reduced complexity state $\sigma_{t}$ are mutually related via

[TABLE]

In the following we will work with mixed states for full generality. Special results applying only to pure states will be discussed as they arise.

By itself, equations (9) and (10) give an over-parameterization of the original state. In the next section we outline how to remove this redundancy by deriving additional constraints on $\sigma_{t}.$

We will usually assume an explicit, smooth parametrization $U_{t}\equiv U_{\theta_{t}}$ by a set of coordinates $\theta_{t}\in D\subset\mathbb{R}^{n}$. Explicit parametrizations allow us to obtain analytic insight into the coupled dynamics of the reduced complexity quantum state and the corresponding group element. We will sometimes drop the explicit time index $t$ and write $U,\phi,\psi,\rho,\sigma,\theta$ when it is clear that they are to be taken at the same time coordinate. This, however, depends on our ability to derive explicitly 111And in a form that can be efficiently evaluated. the partial derivatives of the transformation with respect to the coordinates

[TABLE]

Note that these right generators $\{F_{j}^{>}(\theta),\;j=1,2,\dots,n\}$ are necessarily Hermitian elements of the group’s Lie algebra $\mathfrak{g}$.

In Appendix A we outline how to derive the right generators for some types of parametrizations. In general, however, this requires the ability to explicitly compute exponentiated matrices of the adjoint representation of $\mathfrak{g}$. For high-dimensional Lie algebras this can be fairly challenging. In Section VII we formulate a version of our method that does not require an explicit parametrization of the transformation.

Although our method allows for arbitrary unitary representations of Lie groups some group manifolds cannot be fully parametrized by a single coordinate patch. This can lead to additional technical difficulties which we will usually avoid by limiting ourselves to a single convex coordinate patch $\theta\in D\subset\mathbb{R}^{n}$ that includes the origin [math] which we always assume to map to the identity $U_{0}=\mathbbm{1}$. This poses no serious limitation in most cases of interest. Finally, note that any such parametrization is not unique. We can smoothly re-parametrize the coordinates and then derive the transformed generators via the chain rule. Our use of upper indices for the coordinates and lower indices for the generators is thus motivated by their covariant and contravariant transformation under such re-parametrizations.

III.1 The complexity functional

Consider first a canonical example: for a single bosonic degree of freedom with lowering operator $a$, and a transformation group given by coherent displacements $U_{\theta}=D(\alpha)=e^{\alpha a^{\dagger}-\alpha^{\ast}a},$ an intuitive constraint would be to demand that $\,\left\langle a\right\rangle_{\sigma_{t}}=0,$ or equivalently $\,\left\langle a\right\rangle_{\rho_{t}}=\theta.$ This fully fixes the coordinates and removes the redundancy. This is precisely the constraint which the QSD software package implements. As we demonstrate below, this approach is equivalent to finding the coherent displacement coordinates of the nearest – in the sense of minimal quantum relative entropy distance – symmetric variance Gaussian oscillator state. This approach works well for nearly coherent states $\rho_{t}$ but as we saw in Section II.3 it can actually increase the complexity when $\rho_{t}$ has significant non-Gaussian features.

Below we re-derive this as the result of an optimization problem, which enables us to generalize and improve the approach. To formulate the problem, we introduce a complexity functional $\mathcal{J}(\theta;\rho)$ which, given a state $\rho$, attains a unique global minimum on the space of coordinates $\theta$, we can then fix our coordinates at all times to be

[TABLE]

III.2 Expectation minimization

The simplest choice of complexity functional is obtained by evaluating the expectation of a lower bounded operator $M$ that penalizes population of undesired basis levels in the transformed basis.

[TABLE]

For a bosonic degree of freedom, the penalty operator $M$ could simply be the canonical number operator $M=a^{\dagger}a$. As we will see below, when the transformation is given by coherent displacements $U_{\theta}:=D(\theta^{1}+i\theta^{2})$ this particular choice leads to the QSD scheme $\theta_{\ast}=\left({\,\rm{Re}}\left[\,\left\langle a\right\rangle_{\rho}\right],{\,\rm{Im}}\left[\,\left\langle a\right\rangle_{\rho}\right]\right).$

More generally, however, any lower bounded operator defines a partial order relation for its eigenspaces via the ordering of its eigenvalues. According to such an order, a low complexity state would be characterized by being confined to a subspace spanned by basis states of low order. For composite systems comprising multiple degrees of freedom or even just competing measures of complexity for a single degree of freedom, we can define composite penalty operators by taking linear combinations of $e$ mutually commuting, positive penalty operators $\{M_{k}\;|\;M_{k}\geq 0,[M_{k},M_{l}]=0,\;k,l=1,2,\dots,e\}$ with positive weights $\{\beta_{k}\geq 0,\;k=1,2,\dots,e\}$, i.e.,

[TABLE]

We discuss how to choose the weights $\beta$ below. We will generally refer to this class of optimization problems as expectation minimization.

III.3 Counting Operators

For single degrees of freedom there often exists a special class of such positive operators that has an evenly spaced and non-degenerate spectrum $\lambda_{k}=\lambda_{0}+\lambda k.$ Examples of this are the number operator for a bosonic degree of freedom and any projection to a single spatial axis of the angular momentum operator such as $J_{z}$ for a system with conserved total angular momentum $\vec{J}^{2}=\hbar^{2}J(J+1).$ Any such operator can be normalized such that $\lambda_{0}=0$ and $\lambda=1$ by rescaling and translation $M\to(M-\lambda_{0}\mathbbm{1})/\lambda$. We will refer to such normalized operators as single degree counting operators

[TABLE]

where $\{\left|k\right\rangle;k=0,1,\dots,d\}$ is the orthonormal eigenbasis of $N$. Such counting operators may also admit certain raising and lowering operators $A_{\pm}=A_{\mp}^{\dagger}$ satisfying $[N,A_{\pm}]=\pm A_{\pm}.$ This implies that $A_{\pm}\left|k\right\rangle\propto\left|k\pm 1\right\rangle$. It is then possible to show that under the additional condition that $[A_{-},A_{+}]=\alpha_{1}\mathbbm{1}+\alpha_{N}N$ for some real numbers $\alpha_{1},\alpha_{N}$ this reduces exactly to the above mentioned examples, i.e., in the infinite dimensional case $d=\infty$ the above conditions imply that our problem is equivalent to a harmonic oscillator with the typical raising and lowering operators $A_{+}\propto a^{\dagger}$, while the finite dimensional case $d=2J+1$ is equivalent to an angular momentum space of fixed integral or half-integral $J$ with the angular momentum ladder operators $A_{\pm}\propto J_{\pm}$.

As above we can combine individual counting operators $\{N_{k},k=1,2,\dots,e\}$ to define composite counting operators

[TABLE]

In most cases of interest to us the resulting operators have a unique ground state $N(\beta)\left|\Omega\right\rangle=0$, and they always have a very simple spectrum. An important special case is realized when the individual counting operators $N_{k}$ count excitations of different physical subsystems and each has its own pair of raising and lowering operators $[N_{k},A_{l,\pm}]=\pm\delta_{kl}A_{l,\pm}.$ It is easy to see that these must still be raising and lowering operators of the composite counting operator

[TABLE]

In some cases there exist additional generalized raising and lowering operators for the composite counting operator beyond the raising and lowering operators associated with individual subsystems corresponding to particle exchange or correlated particle creation. As an example, for a collection of harmonic oscillators with raising and lowering operators $\{a_{1},a_{1}^{\dagger},\dots,a_{e},a_{e}^{\dagger}\}$ we could define $N(\beta)=\sum_{k=1}^{e}\beta^{k}a^{\dagger}_{k}a_{k}$ as a composite counting operator and we would then find that the quadratic operators $\{a_{j}^{\dagger}a_{k},\;j,k=1,2,\dots,e\}$ which induce unitary mixing of multiple oscillator modes or $\{a_{j}^{\dagger}a_{k}^{\dagger},\;j,k=1,2,\dots,e\}$ which – along with their hermitian conjugates – induce multi-mode squeezing are also generalized raising and lowering operators:

[TABLE]

These are useful for reducing representation complexity associated with multi-degree correlations. When combining several individual counting operators, it may not be a priori known how to pick good weights $\beta:=(\beta^{1},\dots,\beta^{e})^{T}$. In Section IV.1 we will see that $\rho$ itself contains all the information required to choose $\beta$.

We will generally refer to this class of optimization problems, i.e., minimizing the expectation of single or composite counting operators, as linear excitation minimization. The expectation value of counting operators and polynomials of counting operators can generally be expressed as a finite linear combination of operator expectations with coefficients depending on the coordinates. In some cases this allows to solve for the optimal coordinates directly and in closed form.

III.4 Iterative complexity reduction

While some choices of complexity functionals lead to solutions in closed form, we must usually resort to numerical optimization. For compact Lie groups represented on finite dimensional Hilbert spaces there exist some provably powerful gradient flow methods as described in Li et al. (2008); Schulte-Herbrueggen et al. (2010) but the non-compact and (at least formally) infinite dimensional case that is often of interest to us is more complicated and less well understood. Nonetheless, for specific examples of groups and parametrizations we are able to prove the convexity of some linear and even nonlinear expectation minimization schemes which allows us to employ existing schemes such as generalized Newton’s method Mahony and Manton (2002) for obtaining optimal coordinates. Here we derive explicit expressions for the gradient and Hessian in terms of operator moments resulting from a second order expansion of the functional

[TABLE]

If the complexity functional is strictly convex then the Hessian $h(\theta;\rho)=(h_{jk})_{j,k=1}^{n}$ is positive definite everywhere and an appropriate variant of Newton’s method can be applied to find the optimal coordinates which are implicitly defined by requiring the gradient to vanish $y_{j}(\theta;\rho)=0,\;j=1,2,\dots,n.$

For the case of expectation minimization with penalty operator $M$, the explicit expressions for the gradient and Hessian are also given by simple expectation values

[TABLE]

where

[TABLE]

Although it is not immediately obvious from Equation (25), the symmetry of the Hessian operators $H_{jk}^{>}(\theta)=H_{kj}^{>}(\theta)$ follows straightforwardly from $\frac{\partial^{2}U_{\theta}}{\partial\theta^{j}\partial\theta^{k}}=\frac{\partial^{2}U_{\theta}}{\partial\theta^{k}\partial\theta^{j}}$.

When $\rho_{t}$ evolves with time we can now either solve the minimization problem (12) at each time and use this to obtain the coordinates $\theta_{t}$ or alternatively derive explicit (stochastic) differential equations for the coordinates. While the former will allow us to adapt our scheme to arbitrary stochastic dynamics – jump equations and diffusive dynamics – the latter method can provide us with more insight into the dynamics and open up interesting opportunities for designing control schemes.

III.5 Nonlinear excitation minimization

In this section we explain the importance of the level spacing of the penalty operator spectrum. Assume that we are starting with a counting operator $N$. Since $N$ is positive, we can bound the probability of highly excited states using Markov’s inequality

[TABLE]

Unfortunately, this bound decays only as $O\left(N_{0}^{-1}\right).$

A useful alternative is then to exploit Markov’s extended inequality by applying a monotonically increasing map to $N\to f(N).$ In this case we must have

[TABLE]

If $f(n)$ increases super-linearly with $n$, then the resulting penalty operator $f(N)$ has an increasing spacing of eigenvalues versus $N$ and thus the bound on high excitations will decrease faster than $O\left(N_{0}^{-1}\right)$.

A very strong bound can be achieved with the exponential map $f(n)=\exp(\lambda n)$ parametrized by $\lambda>0.$ This gives

[TABLE]

where we have naturally been lead to introduce the cumulant generating function (CGF) of the penalty operator

[TABLE]

which, for very small $0<\lambda\ll 1$ reduces to

[TABLE]

showing that for small $\lambda$ this optimization problem is equivalent to directly minimizing the expectation of the penalty operator while for increasing $\lambda$ also penalizing large variance.

Note that when $N$ is unbounded there may exist normalizable states ${\rho}$ for which $\,\left\langle e^{\lambda N_{\theta}}\right\rangle_{{\rho}}$ diverges for any $\lambda>0,$ but this is generally true for $\,\left\langle N_{\theta}\right\rangle_{\rho}$ itself. We will assume here without proof that such states do not actually arise in the dynamical evolution of open quantum systems. For a constantly spaced, unbounded spectrum if there exist constants $N_{0}\in\mathbb{N}$ and $\alpha\in[0,1)$ such that $\forall n\geq N_{0}$ we find

[TABLE]

then $\,\left\langle e^{\lambda N_{\theta}}\right\rangle_{{\rho}}$ exists for all $\lambda<\ln 1/\alpha.$ The inequality in (28)-(29) is an example of a Chernoff bound. Since the bound is satisfied for any $\lambda$ we can minimize the right hand side over $\lambda$ to achieve the most restrictive bound yielding

[TABLE]

Given $N_{0}$ the optimal $\lambda_{\ast}$ leading to the lowest bound is implicitly defined via

[TABLE]

The negated left hand side in equation (32) above is proportional to the digits of relative accuracy obtained when truncating the basis at the level $N_{0}.$ The above demonstrates that it is related to the cumulant generating function via a Legendre transformation:

[TABLE]

We thus see that for a fixed $N_{0}$ our complexity functional guarantees a minimal accuracy with which a given quantum state can be represented in a $d$-dimensional subspace of the overall state space, where $d$ is the number of eigenvalues $\lambda_{k}$ of $N$ such that $\lambda_{k}<N_{0}$. This is visualized in Figure 9. Conversely, assuming that $\mathcal{A}(N_{0})$ is one-to-one, for any desired accuracy there exists a specific truncation level $N_{0}$ at which the Chernoff bound is tightest, making it the most efficient truncation level to achieve a certified accuracy.

We will refer to this family of nonlinear excitation minimization schemes as CGF minimization.

IV Information theoretic interpretation

In this section we demonstrate that expectation minimization can also be interpreted as minimizing the quantum relative entropy between $\rho$ and a transformed canonical Gibbs state $\chi_{\theta,\beta}:=U_{\theta}\chi_{\beta}U_{\theta}^{\dagger}$ where $\chi_{\beta}:=Z(\beta)^{-1}e^{-\beta M}$ and $Z(\beta)={\rm Tr}\left(e^{-\beta M}\right),$ which is itself the canonical Gibbs state associated with the transformed penalty operator

[TABLE]

To see this, note that the quantum relative entropy between $\rho$ and $\chi_{\theta,\beta}$ is given by

[TABLE]

Here the first and second terms $-H(\rho)+\ln Z(\beta)$ do not depend on $\theta$ and therefore minimizing the quantum relative entropy over all possible coordinates $\theta$ is equivalent to excitation minimization as discussed above. We may also minimize (38) over $\beta$ to derive an optimal Gibbs weight $\beta_{\ast}$ for a given quantum state $\rho$ and the associated optimal coordinates $\theta_{\ast}.$ By the construction of the Gibbs state, this is equivalent to $\rho$ and $\chi_{\theta,\beta_{\ast}}$ having equal expectation values of the penalty operator, i.e.,

[TABLE]

This then affords a very nice interpretation of expectation minimization: It finds the transformation coordinates that minimize the quantum relative entropy between the actual state and a manifold of generalized Gibbs states. These form the set of maximum Von-Neumann entropy states constrained to have a specific expectation of the penalty operator. In this sense, given a penalty operator $M$, they form the the least biased manifold of states to compare a given state with.

When $\beta=\beta_{\ast}$ the minimum relative entropy is given by the difference in entropies of the actual state and the nearest thermal state

[TABLE]

Note that for a pure state $\rho_{t}=\left|\psi_{t}\right\rangle\left\langle\psi_{t}\right|$ the second term vanishes, whereas the positivity of $S(\rho||\chi)\geq 0$ (with equality if and only if $\rho=\chi$) implies that $\chi_{\theta_{\ast},\beta_{\ast}}$ really is the maximum (over all states) Von-Neumann entropy state with minimum (over all parameters) expected penalty $\,\left\langle M_{\theta_{\ast}}\right\rangle_{\rho}=\,\left\langle M_{\theta_{\ast}}\right\rangle_{\chi_{\theta_{\ast},\beta_{\ast}}}$. By its definition, the Von-Neumann entropy only depends on the spectrum of the density matrix, which is invariant under unitary transforms. This allows us to compute it in either the original or the moving basis $H(\rho_{t})=H(\sigma_{t})$.

In the specific case of collective oscillator states with $M=N(\beta)=\beta^{1}N_{1}+\dots+\beta^{e}N_{e}$, i.e., multi-mode counting operators, the minimal quantum relative entropy between a given state and a generalized Gaussian state has been employed as a measure of non-Gaussianity Genoni et al. (2008); Genoni and Paris (2010). If we take $U_{\theta}$ to allow for all Gaussianity preserving transformation (coherent displacements, unitary and more generally symplectic transformations) then excitation minimization corresponds exactly to finding the nearest Gaussian state. Often, however, we may restrict the transformation to a sub-manifold of the full set of Gaussian transformations to avoid a large decrease of the sparsity of the operators generating the moving basis dynamics. More generally, for any composite penalty operator $M=M(\beta)=\beta^{1}M_{1}+\dots+\beta^{e}M_{e}$ the quantum relative entropy splits into a sum

[TABLE]

where $Z(\beta)={\rm Tr}\left(\exp(-\sum_{k=1}^{e}\beta^{k}M_{k})\right)$ and $M_{k,\theta}=U_{\theta}M_{k}U_{\theta}^{\dagger}$. This information geometric interpretation then affords us a method for optimally composing penalty operators: We demand that the relative entropy be minimized not only over the transformation coordinates $\theta$ but also over the complexity weights $\beta$.

A special situation arises when the penalty operators act on different degrees of freedom, in which case we can decompose the partition function

[TABLE]

Additionally, in certain cases the group transformation may be decomposable into single degree transformations

[TABLE]

such that each factor $U_{k,\theta_{k}}$ is a particular unitary operator parametrized by its own disjoint set of coordinates $\theta_{k}=(\theta_{k}^{1},\dots,\theta_{k}^{n_{k}})$ and acting non-trivially only on the degree of freedom labeled $k$. In this case we have

[TABLE]

which implies that the optimization can be carried out independently for the coordinates $\theta_{k}$ for each $k=1,2,\dots,e$ and also prior to optimizing over the complexity weights $\beta^{k}.$

IV.1 Quadratic expansion of the relative entropy

Any choice of parametrized quantum state family combined with a quantitative measure of distances between quantum states induces a local distance metric for the parameters $\theta$ and $\beta$. For small coordinate displacements $\delta\theta$ and $\delta\beta$ the quantum relative entropy can be expanded to second order which yields

[TABLE]

where the coefficients are

[TABLE]

Near a minimum the linear contributions vanish and in the vicinity of such minima we thus have a quadratic form in the coordinate displacements which can be interpreted as a local distance measure. We will use this expansion later on when we derive effective dynamics for the complexity weights. When $\rho$ is itself a member of the parametrized manifold $\rho=\chi_{\theta_{\ast},\beta_{\ast}}$, then for $(\theta,\beta)=(\theta_{\ast}+\delta\theta,\beta_{\ast}+\delta\beta)$ the linear contribution to the above expansion clearly vanishes and the quadratic coefficients are given by the Kubo-Mori quantum Fisher information metric Petz and Ghinea (2011). The operator derivatives appearing in these definitions can be transformed to the moving bases where they can be expressed in terms of the transformation generators as

[TABLE]

The derivatives of the partition function will generally depend on the particular choice of penalty operators, but we point out that the gradient of the relative entropy with respect to $\beta$ vanishes iff

[TABLE]

As in the case of a single penalty operator these conditions are equivalent to fixing each $\beta^{k}$ such that the expected penalty agrees between the Gibbs state and the actual state $\rho$:

[TABLE]

This then answers our previous question from Section III.2 of how to optimally choose the weights $\beta$. The joint minimization of the quantum relative entropy over both $\beta$ and $\theta$ can also be understood as fitting a generalized Gibbs state parametrized by $(\beta,\theta)$ to $\rho.$

V Dynamics in a moving basis

In this section we describe how to derive coupled equations of motion for the transformation coordinates and the low complexity quantum state such that the coordinates rapidly converge either exactly or approximately to an optimal value.

V.1 Pure state dynamics

We will start by discussing the case of pure state evolution, induced by a (stochastic) Schrödinger equation

[TABLE]

The generator for a closed system and deterministic dynamics is simply $dG=\hbar^{-1}Hdt,$ whereas for an open system evolving according to an unnormalized stochastic Schrödinger equation (SSE) we might have

[TABLE]

where the $\{L_{k},k=1,2,\dots,f\}$ are the Lindblad collapse operators. This particular complex diffusive stochastic unraveling of the open system dynamics can be interpreted as a system whose output modes are all measured using heterodyne detection with perfect fidelity. The associated complex heterodyne measurement processes are given by $dM_{t}^{k}=\,\left\langle L_{k}\right\rangle_{\psi_{t}}dt+dW_{k,t}$, with complex Wiener processes satisfying

[TABLE]

where

[TABLE]

Instead of the complex valued measurement process we could also consider single quadrature homodyne measurements. This yields

[TABLE]

Here the measurement processes are real valued $dM_{t}^{{\rm(hom)},k}=\,\left\langle L_{k}+L_{k}^{\dagger}\right\rangle_{\psi_{t}}dt+dY_{k,t}$ with real Wiener increments satisfying

[TABLE]

where analogously we have introduced

[TABLE]

In either case, it is essential to represent the SSE in the Stratonovich picture as this enables us to carry out stochastic projections Handel and Mabuchi (2005). Contrary to the usual notation we will always take $XdY$ to indicate a Stratonovich stochastic differential as opposed to an Ito differential.

We could also consider situations where some output channels are measured via homodyne and some via heterodyne measurements and we can even generalize to dynamics with discrete jumps such as the quantum jump trajectories encountered when modeling direct photon detection.

As outlined above the state vector $\left|\psi_{t}\right\rangle$ in the fixed basis is related to the reduced complexity state vector $\left|\phi_{t}\right\rangle$ via

[TABLE]

It then follows that the transformed state evolves according to a modified SSE

[TABLE]

with effective generator

[TABLE]

We see that the transformed state has dynamics generated not only by the transformed SSE generator $dG_{t,\theta_{t}}$ but also by the explicit time dependence of the unitary mapping $-\sum_{j=1}^{n}F_{j}^{>}(\theta_{t})d\theta_{t}^{j}$. Note that the coordinate dynamics will in general be stochastic. If we were working in the Ito picture, $dK_{t,\theta_{t}}$ would have to include terms induced by second order stochastic differentials $d\theta^{j}d\theta^{k}$ (and even $d\theta^{k}dM^{j\ast}$) for which the Ito table is a priori unknown. Thus, working in the Stratonovich picture is essential as it allows us to manipulate the processes with the product and chain rules of ordinary calculus. As of yet, we have not specified $d\theta_{t}^{j}$. In Section V.3 we will derive a family of dynamics that constrain the moving basis state exactly or approximately to a minimum of a complexity functional as discussed in Section III.

For pure states the expectation of any hermitian observable $X=X^{\dagger}$ evaluated according to the moving basis state $\phi_{t}$ evolves as

[TABLE]

In deriving this, it was taken into account that (68) does not preserve the norm. Note that when $X$ explicitly depends on time or $\theta_{t}$ one needs to add additional partial derivatives accordingly. We have used Percival’s notation Percival (1999) for the quantum correlation

[TABLE]

which defines a semi-definite inner product on the space of operators (cf. Appendix B).

V.2 Mixed state dynamics

Here we discuss the case of mixed quantum states and deterministic or stochastic Lindblad master equations which are necessary when some of the output channels are unobserved or some measurements have non-ideal fidelity. Note that usually imperfect measurement can be modeled by splitting the output channel on a beamsplitter with transmissivity equal to the fidelity and then observing the transmitted output with perfect fidelity and the reflected output not at all.

The generalization to stochastic master equations is straightforward. In general the fixed basis state $\rho_{t}$ may evolve according to

[TABLE]

where $dG_{t}$ is defined exactly as in the pure state case and the additional Lindblad collapse operators $\{c_{1},c_{2},\dots,c_{f^{\prime}}\}$ correspond to unobserved output channels. The moving basis dynamics follow from a straightforward extension of the pure state case

[TABLE]

where $c_{j,\theta_{t}}:=U_{\theta_{t}}^{\dagger}c_{j}U_{\theta_{t}}.$ In the mixed state case, the expectations of Hermitian observables evolve as

[TABLE]

Here we have overloaded the notation for the quantum correlation Percival (1999) for mixed states in the natural way, i.e.,

[TABLE]

We thus find almost the same result as for pure states except for an additional contribution ${\,\rm{Re}}\left[\,\left\langle[c_{j,\theta_{t}}^{\dagger},X]c_{j,\theta_{t}}\right\rangle_{\sigma_{t}}\right]$ for each unobserved output channel.

V.3 Gradient flow coupled dynamics

Assume now that we have fixed a complexity functional with Hessian $(h_{jk}(\theta_{t}))$ and gradient $(y_{j}(\theta_{t}))$ and that at a given time $t$ we are starting at optimal coordinates, i.e., we have already solved for $\theta_{t}$ such that the complexity gradient $y_{j}(\theta_{t})=0,\;j=1,2,\dots,n.$ Then we can determine the coordinate increments $d\theta_{t}^{j}$ by requiring $y_{j}(\theta_{t+dt})=y_{j}(\theta_{t})+dy_{j}(\theta_{t})=0$. We may then derive the coordinate dynamics by computing the differential change of the gradient coefficients $dy_{j}(\theta_{t})$ as a function of $\left|d\psi_{t}\right\rangle$ and $d\theta_{t}$ solve for $d\theta_{t}$ such that $dy_{j}(\theta_{t})=0$.

More generally, if we assume that we are not starting exactly at optimal coordinates but close to the optimum, then we can instead choose a decay parameter $\eta>0$ and solve for $d\theta$ such that

[TABLE]

This reduces to the above case when we are already at the optimum coordinates, but for a good choice of $\eta$ it leads to increased robustness to slight deviations as they are exponentially damped over time.

Inserting $X=Y_{j}^{>}(\theta)$ into (75) while accounting for its explicit dependence on $\theta$ we obtain

[TABLE]

where we have defined the bias flow

[TABLE]

The second contribution to the bias flow only arises for mixed state dynamics with unobserved output channels.

Combining Equations (76) and (78) we find for $j=1,2,\dots,n$

[TABLE]

Assuming a strictly convex complexity functional and thus a positive definite Hessian $h(\theta_{t})=(h_{jk}(\theta_{t}))_{j,k=1}^{n}$ this relationship can be solved for the coordinate differentials

[TABLE]

for $k=1,2,\dots,n$. We point out that the bias flow is linear in $dG_{t,\theta_{t}}$ and has a single contribution per $c_{j,\theta_{t}}$ which simplifies its derivation.

We will briefly review the independent contributions of Hamiltonian and dissipative deterministic contributions as well as stochastic terms in $dG_{t,\theta_{t}}$. For a Hamiltonian and thus hermitian contribution $\hbar^{-1}Hdt=\hbar^{-1}H^{\dagger}dt$ we find

[TABLE]

For an anti-hermitian dissipative drift term $Vdt=\frac{1}{2i}L^{\dagger}Ldt$ we find

[TABLE]

Finally, each complex diffusion term $dQ=idM^{\ast}L$ contributes

[TABLE]

Considered as a function of $\sigma_{t}$ and $\theta_{t}$ this can be used to derive the corresponding Ito SDEs which, however, are in general quite complicated.

V.4 Gradient coupled fiducial state dynamics

If the system state remains localized near a semi-classical manifold of generalized Gibbs states $U_{\theta}\Omega U_{\theta}^{\dagger},$ then it may be advantageous to evaluate the Hessian $(h_{jk}(\theta_{t}))$ and the bias flow $(dq_{j})$ not in the actual moving basis state $\sigma_{t}$ but in a reference state $\Omega=Z(\beta)^{-1}e^{-\sum_{l=1}^{e}\beta^{l}M_{l}}$ instead. In this case, Equation (80) becomes

[TABLE]

for $k=1,2,\dots,n$. The advantage of this is that $[h^{(\Omega)}(\theta_{t})^{-1}]^{kj}$ and $dq_{j}^{(\Omega)}$ are purely functions of $\theta$ and other scalar model parameters. Thus, (87) can be understood as semi-classical equations of motion for the coordinates coupled to the true quantum state via the gradient $\eta\;y_{j}(\theta_{t})\;dt$.

V.5 Information projected fiducial state dynamics

It turns out that for an expectation minimization functional the quantum state independent contributions to the dynamics $d\theta^{k,\rm sc}:=\sum_{j=1}^{n}[m^{(\Omega)}(\theta_{t})^{-1}]^{kj}dq_{j}^{(\Omega)}$ are fully equivalent to the manifold projection method proposed in Mabuchi (2008) except that the metric on the tangent space is taken to be the Kubo-Mori metric Bengtsson and Zyczkowski (2006) that naturally arises as a generalized quantum Fisher information metric associated with differential increases in quantum relative entropy Petz and Ghinea (2011). In fact, in this setting one may assume the complexity weights $\beta$ to be time dependent and derive coupled dynamics for them, as well. We discuss a special case of this in the following.

V.6 Gibbs manifold projection

Here we briefly sketch how to use our framework to obtain approximate low dimensional dynamics that arise when the state $\rho_{t}$ is constrained to the Gibbs-manifold parametrized by the group transformation coordinates $\theta_{t}$ and the complexity weights $\beta_{t}$ at all times. For a Gibbs state $\chi_{\theta,\beta}=Z(\beta)^{-1}\exp\left(-\sum_{k=1}^{e}\beta^{k}M_{k,\theta}\right)$ both the quantum relative entropy as well as its linear variation $\sum_{k=1}^{e}z_{k}\delta\beta^{k}+\sum_{k=1}^{n}y_{k}\delta\theta^{k}$ (cf. Equations (48) and (50)) with respect to $\theta$ and $\beta$ vanish because it is already on the Gibbs manifold we project to. We can now apply the full procedure outlined in Mabuchi (2008) to also project the stochastic dynamics of $\rho_{t}$ onto the manifold. Generally this amounts to choosing $d\theta$ and $d\beta$ such that $dy=(dy_{1},\dots,dy_{n})$ and $(dz_{1},\dots,dz_{e})$ are orthogonal to the tangent space spanned by $d\theta$ and $d\beta$ with respect to the Kubo-Mori associated with the second order variation of the quantum relative entropy.

In the interest of brevity we simply present the result of that derivation. The coordinate increments $d\theta$ and $d\beta$ are the implicit solutions of the following linear system

[TABLE]

with $h_{jl},y_{kj}$ and $g_{km}$ defined as in Equations (49) - (52) for $\rho\equiv\chi_{\theta,\beta}$. The flow vectors $dq=(dq_{1},\dots,dq_{n})$ and $dm=(dm_{1},\dots,dm_{e})$ are given by

[TABLE]

both of which can be evaluated in the moving basis via the diagonal Gibbs state $\chi_{\beta}$ for which all off-diagonal operator moments conveniently vanish. In many cases of interest the above expressions can be evaluated analytically as all operator moments become explicit functions of only $\beta$ and $\theta$, yielding coupled, low dimensional ODEs or SDEs that approximately describe the original quantum dynamics. We note that the transformation coordinates $\theta$ can often be directly associated with classical or semi-classical quantities such as generalized canonical position and momentum variables, but the complexity weights $\beta$ have a purely statistical interpretation as they encapsulate the uncertainty originally encoded in the exact quantum state.

VI Examples of Manifolds

Here we present some examples of transformation groups, penalty operators and the corresponding generators and sensitivity variables. A given transformation $U_{\theta}$ is characterized by its differential form $U_{\theta_{t}}^{\dagger}dU_{\theta_{t}}=-i\sum_{j=1}^{n}F_{j}^{>}(\theta_{t})d\theta_{t}^{j}$, its adjoint action on its own Lie-Algebra $U_{\theta_{t}}^{\dagger}XU_{\theta_{t}}$ and the sensitivity operators whose expectations form the gradient and Hessian of a particular expectation minimization problem. The sensitivity operators can be derived for any given complexity functional.

VI.1 Coherent displacement

The simplest example is that of coherent displacements $\theta=(Q,P)$. Working in this representation we have

[TABLE]

where the generators and right generators are given by

[TABLE]

Coherent displacements of different modes commute and therefore the extension to multiple oscillator modes follows trivially.

VI.1.1 Linear excitation minimization

For the canonical penalty operator given simply by the photon number operator $M=N:=a^{\dagger}a$ the gradient and Hessian operators are

[TABLE]

We see that the Hessian is constant and strictly positive definite. Since the domain of the coordinates is $\mathbb{R}^{2}$ and thus clearly convex and closed, we always have a unique optimum.

VI.1.2 Nonlinear excitation minimization

Assuming a nonlinearly transformed penalty operator $M=f(N)$ the gradient operators are

[TABLE]

where we have introduced

[TABLE]

which allows us to write down two useful (and equivalent) rules

[TABLE]

The Hessian operators are then given by

[TABLE]

Here the second order differences have been introduced:

[TABLE]

To analyze under what conditions the Hessian may be positive, observe first that the Hessian operator matrix $H^{>}:=(H_{jk}^{>}(\theta_{t}))$ factors as

[TABLE]

where we have defined

[TABLE]

via

[TABLE]

$V$ may be interpreted as a linear transformation from a real valued quadrature representation of the displacement $\delta\theta=(\delta Q,\delta P)^{T}$ to complex amplitudes $(\alpha,\alpha^{\ast})^{T}=V\delta\theta.$ Since $V$ contains only scalars, the matrices $h,\tilde{h}$ of element-wise expectations, i.e.,

[TABLE]

are related by a similarity transform $h=V^{\dagger}\tilde{h}V$ and therefore have the same eigenvalues. Equipped with the above definitions we are thus able to formulate and prove the following

Theorem: A sufficient condition for the positive semi-definiteness of the Hessian $h$ is that the following operators are non-negative

[TABLE]

where $\Gamma$ is defined as above and

[TABLE]

Proof: See Appendix C.

This theorem is helpful because the operators appearing in its conditions are all diagonal in the eigenbasis of $N$, which implies that we can verify the positivity simply by evaluating the inequalities for each eigenstate of $N$. Note furthermore that the non-negativity of $\Gamma$ is also a necessary condition and that when $\Gamma-\Xi>0$ is strictly positive, the Hessian is strictly positive.

Sometimes we only wish to determine whether the Hessian is positive when the states $\sigma=\Pi\sigma\Pi$ are restricted to a subspace of the full state space characterized by a projection operator $\Pi$. In this case the operator inequalities must only hold for the projected operators $\Pi\Sigma\Pi,\Gamma\Pi,\Xi\Pi.$

Applying the above theorem to some interesting test cases we find that $f(n)=n^{2}$ leads to a convex objective, as this leads to $\Delta_{f}^{(2)}(N-1)\geq 0,\Gamma=1+4N,\Xi=1+2n.$

For $f(n)=e^{\lambda n}$ with $\lambda>0$, which is relevant for CGF minimization, we find that $\Gamma-\Xi$ is not positive everywhere, but for small enough $\lambda$ there exists an $n_{\ast}=\left.\left\lceil(3x^{2}-x^{3})/(-1+x)^{3}\right\rceil\right|_{x=e^{\lambda}}$ such that at least on the subspace of number states lower than $n_{\ast}$ the condition holds: $\Pi_{N\leq n_{\ast}}(\Gamma-\Xi)\geq 0.$ This upper bound on $N$ becomes arbitrarily large as $\lambda$ goes to zero.

VI.2 Squeezing and displacement

Our next example will be a mixture of squeezing and displacement albeit each only parametrized by a single variable $\theta=(Q,R)$, which does not allow to realize the most general pure Gaussian state.

[TABLE]

with generators

[TABLE]

For the canonical complexity functional $N=a^{\dagger}a$ the sensitivity operators are

[TABLE]

This still leads to a positive definite Hessian, but it is now dependent on the coordinates and the state. It is straightforward to derive the more general case of arbitrary displacements and squeezing, but the expressions become more tedious. They are implemented in our software package QMANIFOLD Tezak (2016).

VI.3 Spin coherent displacement

Our final example here is for angular momentum states and a single irreducible representation labeled by $J$ such that $\mathbf{J}^{2}=J(J+1)$ and the spectrum of $J_{z}$ is given by $-J,-J+1,\dots,J-1,J.$ The commutator relationships are $[J_{z},J_{\pm}]=\pm J_{\pm}$ and $[J_{+},J_{-}]=2J_{z}.$ The ladder operators can also be represented in terms of the x-y generators $J_{\pm}=J_{x}\pm iJ_{y}.$ We can then define a unitary Lie group explicitly via

[TABLE]

The action of this unitary on a fixed reference state, e.g., $U_{\theta}\left|j_{z}=-J\right\rangle$ can be identified with the spin coherent states as introduced by Radcliffe Radcliffe (2001). For the canonical spin complexity functional $N=J_{z}$ the sensitivity operators are

[TABLE]

where ${\,\rm{Re}}\left[\mu J_{-}\right]=(\mu J_{-}+\mu^{\ast}J_{+})/2$ and ${\,\rm{Im}}\left[\mu J_{-}\right]=(\mu J_{-}-\mu^{\ast}J_{+})/2i$. This does not generally lead to a positive definite Hessian, e.g., consider the $j_{z}=+J$ eigenstate of $J_{z}$ for which the expected Hessian is clearly negative definite. For any given state the eigenvalues of $\left(\,\left\langle H_{jk}^{>}(\theta_{t})\right\rangle_{\sigma}\right)$ are given by

[TABLE]

Assuming that $\,\left\langle J_{z}\right\rangle_{\sigma}<0$, these are non-negative iff

[TABLE]

Two sufficient but not necessary conditions for this are given by

[TABLE]

and a more restrictive version that depends only on $\,\left\langle J_{z}\right\rangle_{\sigma}^{2}.$

[TABLE]

This follows directly from $\,\left\langle J_{k}\right\rangle_{\sigma}^{2}\leq\,\left\langle J_{k}^{2}\right\rangle_{\sigma}$ for $k=x,y,z$ and from $\sum_{k\in\{x,y,z\}}J_{k}^{2}=\mathbf{J}^{2}=J(J+1).$

VII Coordinate free method

The methods we have introduced in the preceding section are quite appealing in that they can allow us to simulate significantly larger open quantum systems than is possible with static bases. They can also yield analytic insight into semi-classical dynamics and provide more intuition. With the analytic results presented in Section VI one can construct commuting product transformations for fairly complex systems, but they are unlikely to work well when considering non-factoring group manifolds of coordinate dimensions beyond $O(10)$ as it becomes very difficult to obtain a fully exponentiated parametrized transform from which we can derive the right generators $\{F_{j}^{>}(\theta_{t})\}$. These, however, are essential to the derivation of most other important quantities. Moreover, having an explicit coordinate based representation of the unitary $U_{\theta}$ is useful for transforming states from the static to the moving basis. There are, however, examples of groups that can at least in principle be employed without ever having to derive an explicit parametrization of $U_{\theta}$. Here we sketch out how this can be achieved in principle.

Consider a Lie group and a state space and dynamics for which

the model operators that generate the dynamics $H,L$ and $c$ can be represented as polynomials of elements of a group’s finite dimensional Lie algebra, and 2. 2.

a full Hilbert space basis can be generated by repeated action of some raising operators $\{A_{k,+}\}$ that are inside the Lie algebra starting from a unique reference state $\left|\Omega\right\rangle,$ which is itself fully characterized as the unique zero-eigenvalue eigenstate of an operator $M\left|\Omega\right\rangle=0$, where $M$ is also either inside the group’s Lie algebra or a polynomial of Lie algebra elements.

In this case we may use the adjoint representation of the group transform itself as the parametrization, i.e., if $\{Y_{1},Y_{2},\dots,Y_{q}\}\subset\mathfrak{g}$ is a hermitian basis of the Lie algebra, then if a group element $g$ is unitarily represented on the Hilbert space as $U(g)$ there also exists a corresponding element in the adjoint representation of the group $S(g)$ such that $U(g)Y_{k}U(g)^{\dagger}=\sum_{j=1}^{q}S_{k}^{j}(g)Y_{j}.$ This adjoint representation will typically require only a polynomial number of parameters in the system size. In some cases where a part of the Lie algebra can be expressed as polynomials of some other basis elements further reduction is possible, e.g., consider two bosonic modes $a,b$ and a Lie algebra that also contains $a^{\dagger}b$ or $ab$. In this case, if we know how $a$ and $b$ transform under $U$ then we also know how $a^{\dagger}b$ or $ab$ transform.

If we start at $t=0$ with the identity $g=e$ and $S_{k}^{j}(e)=\delta_{k}^{j},$ then at each time step we can consider arbitrary elements from the Lie algebra and write $U(g)^{\dagger}\delta U(g)=-i\sum_{j=1}^{q}\delta\mu^{j}Y_{j}.$ This implies that

[TABLE]

Here we have used the convention for the structure constants given in Equation 147 and the final two lines allow us to read off the differential change of the adjoint transformation

[TABLE]

Armed with this, we can now reformulate the complexity reduction problem. Assume that we wish to find a transfom $U_{\ast}$ such that it minimizes a function $\mathcal{J}_{M}(U)=\,\left\langle UMU^{\dagger}\right\rangle_{\rho}$, i.e.,

[TABLE]

Expanding this to second order in a a perturbation $\delta U=-i\delta\mu^{j}Y_{j}$ we find

[TABLE]

Here we have used the notation $[Y_{\{k},[M,Y_{j\}}]]$ to indicate that the expression is symmetrized over $k$ and $j$, i.e., $[Y_{\{k},[M,Y_{j\}}]]\equiv\frac{1}{2}\left\{[Y_{k},[M,Y_{j}]]+[Y_{j},[M,Y_{k}]]\right\}$. As above, the gradient and Hessian coefficients can be more efficiently evaluated in the current lowest complexity state $\sigma=U^{\dagger}\rho U.$ Assuming that $(h_{jk})$ is positive semi-definite at all steps, we can use Newton’s method or a similar technique to identify good coordinate updates $\delta\mu$ and update our adjoint representation $S\to S+\delta S$ as specified in (138) as well as the transformed state $\sigma=U^{\dagger}\rho U$ according to

[TABLE]

As above, we may also combine this complexity reduction with dynamics of $\rho$ or $\sigma$, respectively. The results of our previous sections carry over with little modification. Although the lack of a fixed parametrization of $U$ may cause us to worry about what the specific meaning of $\sigma$ really is, our initial assumptions guarantee that the adjoint representation matrix $S$ uniquely fixes how our Hamiltonian and all other relevant operators are transformed to the low complexity basis and furthermore guarantees that a representation in terms of the original basis states can at least in principle be recovered in the new basis by transforming $M_{U}:=UMU^{\dagger}$ using $S$ and then solving for the unique eigenvector $\left|\Omega_{U}\right\rangle=U\left|\Omega\right\rangle$ with zero eigenvalue of $M_{U}.$ All other basis states can then be recovered by acting on $\left|\Omega_{U}\right\rangle$ with the transformed raising operators $\{UA_{k,+}U^{\dagger}\}$. This prescription should work straightforwardly for more complicated groups such as the $SU(n)$ symmetry groups of multiple bosonic or fermionic modes or the symplectic group $Sp(2n)$ for an ensemble of oscillators or distinguishable particles.

VIII Summary of our methods

In this work we have presented analytical and computational methods that allow to simulate closed and open quantum systems assuming deterministic – unitary or ensemble averaged – or stochastic dynamics affecting either pure or mixed states.

Given a particular dynamical quantum system to be simulated, we suggest the following steps to apply our methods. We assume as the most general case of mixed states and stochastic evolution as described in Equation (V.2) defined by a generator $dG$ and a set of unobserved dissipation operators $\{c_{j},j=1,2,\dots,f^{\prime}\}$.

Based on intuition or prior simulations, identify the degrees of freedom that are most likely to localize. 2. 2.

Define a joint penalty operator $M=\sum_{k=1}^{e}\beta^{k}M_{k}\geq 0$ and choose a set of transformation generators $\{X_{j},j=1,2,\dots,n\}$ that are elements of a finite dimensional Lie algebra and use the methods outlined in the Supplementary Material to obtain an exponentiated form of the parametrization $U_{\theta}$ and consequently the right generators $\{F_{j}^{>}(\theta_{t})=iU^{\dagger}_{\theta_{t}}\frac{\partial U_{\theta_{t}}}{\partial\theta^{j}},j=1,2,\dots,n^{\prime}\leq n\}$. This generally requires explicit exponentiation of matrices in the adjoint representation and can be done using a computer algebra system. In Section VI we have provided a number of examples that can be combined to handle composite systems. 3. 3.

Derive the explicitly transformed generators $dK_{t,\theta_{t}}$ as defined in Equation (69) and the transformed dissipation operators for the unobserved channels $c_{j,\theta_{t}}=U_{\theta_{t}}^{\dagger}c_{j}U_{\theta_{t}}.$ This can most easily be done by writing these operators as polynomials of the transformation generators and then applying the explicitly obtained adjoint representation to transform them. 4. 4.

Simulate the joint equations of motion for the reduced complexity state $\sigma_{t}$ and the coordinates $\theta_{t}$ defined as

[TABLE]

where the expected gradient $y_{t}$ and Hessian $h_{t}$ are defined in Equation (22) and where the bias flow $dq_{t}$ has been defined in Equation (78). The gain parameter $\eta\geq 0$ can be freely chosen. A reasonable heuristic to picking a good value is to linearize the equation of motion for as a function of $\theta_{t}$ and choose $\eta\geq\min\left\{-{\,\rm{Re}}\left[\lambda\right],\lambda\in{\rm spec}(J)\right\}$ where $J$ is the Jacobian of that linearization. The special case of pure state simulations is handled equivalently. 5. 5.

The resulting simulated trajectory for the reduced complexity state $\sigma_{t}$ can be used to test whether further improvements can be achieved by allowing more general transformations. On the other hand, if an inspection of the coordinate trajectory $\theta_{t}$ reveals that some parameters are nearly constant or more generally the coordinate trajectory is itself confined to some lower dimensional set then it may be possible to remove some degrees of freedom from the transformation via an embedding transformation $\tilde{\theta}\mapsto\theta$ where $\tilde{\theta}$ has fewer coordinates. In general it can be useful to iterate the above steps a few times to identify a good parametrization.

IX Conclusion and Outlook

As stated at the outset, we have presented a family of methods that allow one to investigate the inherent complexity of quantum states by attempting to reduce the total number of variables required for their description. We provided an information theoretic interpretation of our method and several alternate prescriptions for deriving exact coupled dynamics of the semi-classical group coordinates $\theta$ and the reduced complexity quantum state $\sigma$. Our method includes earlier work by Schack, Brun and Percival Schack et al. (1995) as a special case but provides a substantially larger, analytic framework that can also be applied in cases where their method does not work. Further, we have shown how our method connects to earlier model-reduction work on nonlinear projection of quantum models Mabuchi (2008). Our approach is extremely flexible and suited for arbitrary Lie groups as long as their associated Lie algebras have finite dimension. In practice the explicit coordinate representations will work best for low dimensional Lie-algebras, and we have also outlined how to implement a method that does not rely on an explicit coordinate parameterization.

This work has many promising future directions, some key examples are the application of our simulation scheme to quantum measurement and control problems. In particular, our method is very nicely suited to model noisy dispersive or high power qubit readout of superconducting quantum circuits Clerk et al. (2010) and it could be used for a more principled approach to the quantum feedback model for autonomous state preparation considered in Andersen et al. (2016). It could also provide more rigorous simulations of the quantum effects in coherent optical Ising machines or coherent machine learning devices than existing techniques based on quasi-probabilities Takata et al. (2015); Maruo et al. (2016); Tezak and Mabuchi (2015).

By completely projecting the reduced complexity state to the nearest Gibbs state, our method can also be used to derive further reduced order projected models such as the Maxwell-Bloch type model considered in Mabuchi (2008) or a semi-classical coupled mode theory for nonlinear resonators similar to the Wigner method proposed in Santori et al. (2014). It would be very appealing to work out in more detail how our approach could improve our understanding of quantum feedback networks Gough and James (2008, 2009) by conceptually separating both the node-systems and the interconnecting fields into quasi-classical and quantum components. A very promising framework for this avenue is given by the recent work of John Gough Gough (2017) on controlled flows.

It is also possible to extend the formalism to quantized fields, in particular it may be very useful to apply it to traveling wave fields inside quantum feedback networks with time delays. In this scenario the reduced complexity state for the bath modes could be modeled using Matrix Product States as in Pichler and Zoller (2016) or by an approximate delay model Tabak and Mabuchi (2016). Recent and exciting work by Sarovar et al. Sarovar et al. (2017) addresses a slightly different question, i.e., that of which quantum models are robustly solvable on an analog quantum simulator, but using similar information geometric techniques as our approach. Ultimately, we expect the question of efficient representation to be closely related to the question of robust representation. Finally, we believe that it is worth investigating whether the inherent complexity of the quantum states a system evolves through can be related to its computational power. A formal framework for classical dynamical input/output systems has been introduced in Dambre et al. (2012).

Appendix A General construction of the coordinate transformation

In this secection we outline how to construct complex transformations and derive the right generators. The simplest construction for the transformation is by chaining single parameter transformations

[TABLE]

Thus far we allow for complex coordinates $\{\eta^{j}\}$ and arbitrary, i.e., not necessarily hermitian, generators $\{X_{j}\}$. We assume that the generators $X_{j}$ are elements of a finite dimensional Lie algebra $\mathfrak{g}$. Given a basis $\{Y_{1},Y_{2},\dots,Y_{q}\}\subset\mathfrak{g}$ for the Lie algebra with structure constants $c_{jk}^{l}$ implicitly defined via

[TABLE]

we represent each transformation generator in this basis as $X_{j}=\sum_{k=1}^{q}R_{j}^{k}Y_{k}.$ Using the structure constants, it is straightforward to compute the conjugation of any basis element by a single parameter transformation Puri (2001) to be

[TABLE]

The matrices $A^{(j)}_{\eta^{j}}$ are typically very sparse and can be exponentiated symbolically using a tool such as Mathematica Wolfram Research (1988) the SymPy package Meurer et al. (2017). The resulting matrices $S^{(j)}_{\eta^{j}}\in\mathbb{C}^{n\times n}$ are elements of the adjoint representation of the transformation group and can be used to directly transform the generators. As they each depend on only a single coordinate, they satisfy $S^{(j)}_{-\eta^{j}}=\left(S^{(j)}_{\eta^{j}}\right)^{-1}$. With this, it is straightforward to see that

[TABLE]

We see that by requiring this particular differential form of $U_{\eta}$, i.e., with all generators on the right hand side, each generator is additionally transformed $X_{j}\to X_{j}^{>}(\eta)$ by all single parameter transformations that appear to its right.

The differential transformation can be equivalently expressed with the differential generators on the left side of $U_{\eta_{t}}$:

[TABLE]

If we now assume that each generator is antihermitian $X_{j}=iF_{j}$ and if we restrict the coordinates to real values $\eta^{j}=\theta^{j}\in\mathbb{R}$, then the resulting transformation is unitary. This is desirable because observables evaluated in the transformed frame $\left\langle\phi_{t}\right|M^{\prime}\left|\phi_{t}\right\rangle=\left\langle\psi_{t}\right|U^{\dagger}_{\theta_{t}}MU_{\theta_{t}}\left|\psi_{t}\right\rangle$ are actually simply the conjugated observables $M^{\prime}=U^{-1}_{\theta_{t}}MU_{\theta_{t}}.$

To transform any operator $M$ from the basis associated with $\left|\psi_{t}\right\rangle$ to the moving basis $\left|\phi_{t}\right\rangle$ one first needs to express it exclusively in terms of functions of the generators $M=f(X_{1},X_{2},\dots,X_{n})$ where $f:\mathbb{C}^{n}\to\mathbb{C}$ is analytic in all variables, typically $f$ is a polynomial. We then have

[TABLE]

We see that having the adjoint representation single parameter transformation matrices $\{S^{(j)}_{\theta^{j}},j=1,2,\dots,n\}$ allows us to do all necessary computations. We remark that representing $M=f(X_{1},X_{2},\dots,X_{n})$ generally does not imply a unique function $f$ as some generators may themselves be polynomials of the other generators.

Finally, note that there exist alternate ways Wilcox (1967) of parameterizing groups and deriving partial derivatives that may come in useful in more complex cases. Our rules for chained single parameter transformations derived here can straightforwardly generalized to chained multi-parameter transformations.

For very complex parameterizations, analytical/symbolic methods may fail but in that case it should still be possible to work in a purely numerical representation that stores and integrates both $\theta$ and elements of the adjoint representation of $U_{\theta}.$

Appendix B Properties of the quantum correlation

The quantum correlation is not a strictly positive definite inner product because the quantum self-correlation of an operator $A$ vanishes in any eigenstate

[TABLE]

The sufficiency “$\Rightarrow$” of this condition is obvious, the necessity “$\Leftarrow$” follows from the Cauchy-Schwarz inequality for the regular Hilbert space inner product.

Restricted to hermitian operators $A^{\dagger}=A,B^{\dagger}=B$ the quantum correlation can be decomposed into its real and imaginary part as

[TABLE]

with the symmetrized covariance function

[TABLE]

Appendix C Proof of our theorem

Since $\Gamma$ is non-negative, the trace of $\tilde{h}$ is non-negative: ${\rm Tr}\left(\tilde{h}\right)=2\,\left\langle\Gamma\right\rangle_{\sigma}\geq 0$. Due to the hermiticity of $\tilde{h}$ its eigenvalues must be real and it thus suffices to prove that its determinant is non-negative. For any non-negative operator $A\geq 0$ we must have $0\leq\,\left\langle(a+a^{\dagger})A(a+a^{\dagger})\right\rangle_{\sigma}$ for any state $\sigma$. Inserting $A=\Delta_{f}^{(2)}(N-1)$ – which is positive by the first condition of our theorem – this yields

[TABLE]

This inequality is satisfied not only for $\sigma$ but also for any unitarily rotated state $\sigma(\phi):=e^{i\phi N}\sigma e^{-i\phi N}$. The second expectation value in (160) is invariant under any such transformation, whereas it imparts a complex phase on the first expectation value.

[TABLE]

We can therefore always find a $\phi_{\ast}$ such that the first expectation is real valued and non-positive $2{\,\rm{Re}}\left[\,\left\langle\underbrace{\Delta^{(2)}_{f}(N)a^{2}}_{\Sigma}\right\rangle_{\sigma(\phi_{\ast})}\right]=-2\left|\,\left\langle\underbrace{\Delta^{(2)}_{f}(N)a^{2}}_{\Sigma}\right\rangle_{\sigma}\right|.$ We have thus shown that

[TABLE]

Under the second condition to our theorem we have $\Gamma-\Xi\geq 0$ which implies $\,\left\langle\Gamma\right\rangle_{\sigma}\geq\,\left\langle\Xi\right\rangle_{\sigma}.$ We thus know that the determinant of $\tilde{h}$ can be bounded from below as

[TABLE]

This proves our claim. ∎

Appendix D Further results for simulating under CGF minimization

Here we present some additional results obtained in simulating the DOPO system in a basis that minimizes the CGF functional. In particular, we investigate the relationship between the manifold coordinates and the expectation of the oscillators’s mode operators. While excitation minimization will always enforce $\,\left\langle a\right\rangle_{\phi_{t}}\equiv 0\Leftrightarrow\,\left\langle a\right\rangle_{\psi_{t}}=\alpha=\frac{Q+iP}{\sqrt{2}}$ CGF optimization generally does not lead to such a linear relationship as can be seen in Figure 10.

We can also see that different regions in phase space lead to different complexity as measured by the CGF (cf Figure 11). This motivates using a simulation method in which even the basis size is adapted to the inherent complexity of the current dynamics.

Finally, we note that even within the displaced basis there appears to be yet lower dimensional attractors for the reduced complexity state $\left|\phi_{t}\right\rangle$. In Figure 12 we have visualized the distribution of the first three moving basis level populations when transforming to the CGF optimal basis.

In Figure 13 we have simulated the gradient coupled fiducial state dynamics for the DPO system above threshold for different values of the coupling gain $\eta.$ All stochastic simulations were carried out with the same random seed and thus the same realization of the innovation process. Figure 13 (c) shows some discrepancies between states reconstructed from simulations with lower $\eta$ and $\eta_{0}=200.$ Specifically we compute $d_{\rm{FS}}^{2}(\psi^{(\eta)},\psi^{(\eta_{0})}),$ where the fixed basis representation states are obtained from the respective moving bases representation states via $\left|\psi^{(\eta)}(t)\right\rangle=U_{\theta^{(\eta)}}\left|\phi^{(\eta)}(t)\right\rangle$ and the Fubini-Study distance is defined as:

[TABLE]

The errors decrease with increasing $\eta$ suggesting that a strong gradient coupling gain $\eta=O(100)$ is preferable. We intend to investigate this more systematically in future work. We note that the errors appear to stay constant over time, which is encouraging. Furthermore, we have simulated our system with different sizes of the basis (cf Figure 13 (d)) and evaluated the error relative to the most accurate simulation. Surprisingly, we find that the truncation error remains roughly constant, i.e., even for this randomly switching system, the low-dimensional approximations to the system state track the actual system state very well.

Acknowledgements.

We wish to thank Dmitri Pavlichin, Ryan Hamerly, Tatsuhiro Onodera, Edwin Ng, Gil Tabak, Michael Goerz, Dodd Gray, Peter McMahon, Jeff Hill, Nicholas Rubin, Mohan Sarovar, Thaddeus Ladd, Hendra Nurdin, John Gough, Aashish Clerk, Patrick Hayden, Surya Ganguli, Mauro Maggioni, Todd Brun and Thomas Schulte-Herbrüggen for insightful discussion and feedback. This work was supported by the Army Research Office under grant W911NF-16-1-0086 and the National Science Foundation under grant PHY-1648807. N.T. additionally acknowledges support through a Stanford Graduate Fellowship and a Math+X fellowship by the Simons Foundation. Also, N. H. A. is grateful for funding from the JCJC INS2I 2016 “QIGR3CF” project and the JCJC INS2I 2017 “QFCCQI” project.