Efficient representation of Gaussian states for multi-mode non-Gaussian quantum state engineering via subtraction of arbitrary number of photons

TL;DR

This paper introduces a comprehensive formalism for representing multi-mode Gaussian states and their non-Gaussian derivatives via photon subtraction, enabling analytical characterization of state fidelity, heralding probability, and applications in quantum state engineering.

Contribution

It develops a novel $K$ function formalism for multi-mode Gaussian states, allowing analytical analysis of non-Gaussian states generated by photon subtraction and heralding.

Findings

Heralding probability expressed as a Hafnian, linking to Gaussian boson sampling.

Analytical fidelity and success probability for non-Gaussian state preparation.

Proposed method for near-perfect two-mode coherent-cat Bell state generation.

Abstract

We introduce a complete description of a multi-mode bosonic quantum state in the coherent-state basis (which in this work is denoted as "$K$" function ), which---up to a phase---is the square root of the well-known Husimi "$Q$" representation. We express the $K$ function of any $N$-mode Gaussian state as a function of its covariance matrix and displacement vector, and also that of a general continuous-variable cluster state in terms of the modal squeezing and graph topology of the cluster. This formalism lets us characterize the non Gaussian state left over when one measures a subset of modes of a Gaussian state using photon number resolving detection, the fidelity of the obtained non-Gaussian state with any target state, and the associated heralding probability, all analytically. We show that this probability can be expressed as a Hafnian, re-interpreting the output state of a circuit claimed to demonstrate quantum supremacy termed Gaussian boson sampling. As an example-application of our formalism, we propose a method to prepare a two-mode coherent-cat-basis Bell state with fidelity close to unity and success probability that is fundamentally higher than that of a well-known scheme that splits an approximate single-mode cat state---obtained by photon number subtraction on a squeezed vacuum mode---on a balanced beam splitter. This formalism could enable exploration of efficient generation of cat-basis entangled states, which are known to be useful for quantum error correction against photon loss.

Figures

Equations

$$I=\frac{1}{(2\pi {)}^N}\int d^{2N}{\vec{x}}_{\alpha }|\vec{\alpha }⟩⟨\vec{\alpha }|,$$

$$|{\Psi }_0⟩$$

$$|\Psi ⟩=\int d^{2N}{\vec{x}}_{\alpha }K({\vec{x}}_{\alpha })G({\vec{x}}_{\alpha },{\vec{x}}_{\beta })|\vec{\alpha }⟩,$$

$$K({\vec{x}}_{\alpha })$$

$$G({\vec{x}}_{\alpha },{\vec{x}}_{\beta })$$

$$\mathcal{B}$$

$$\mathcal{D}$$

$$\mathcal{X}$$

$${\Gamma }^{-1}=\left(\begin{array}{cc}A& C\\ C^T& B\end{array}\right).$$

$$|{\Psi }_{-\vec{m}}⟩$$

$$|\Phi ⟩$$

$$\displaystyle P_{\vec{n}}=\frac{1}{\det\mathcal{H}\sqrt{\det\Gamma}\prod\limits_{i=1}^{N}n_{i}!2^{n_{i}}}\big{|}\mathcal{I}_{\vec{n}}\big{|}^{2},$$

$${\mathcal{I}}_{\vec{n}}$$

$$R({\vec{x}}_{\alpha })$$

$$\displaystyle\mathcal{I}_{\vec{n}}=\left\{\begin{array}[]{ll}0&\Sigma=\textrm{odd},\\ \textrm{Hf}\left(F\right)&\Sigma=\textrm{even},\end{array}\right.$$

$$\hat{H}=-\frac{i}{2}\sum _{i,j}^N{G}_{ij}\left({\hat{a}}_i^{†}{\hat{a}}_j^{†}-{\hat{a}}_i{\hat{a}}_j\right),$$

$$\mathcal{B}=\frac{1}{2}I+\frac{1}{2}\tanh r\left(\begin{array}{cc}-G& iG\\ iG& G\end{array}\right).$$

$$|C⟩=\frac{1}{N_{+}}\left(|\gamma ,\gamma ⟩+|-\gamma ,-\gamma ⟩\right),$$

$$G=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),$$

$$P_5$$

$$F_5$$

$$P_{(i)}$$

$$P_{(ii)}$$

$$F_{(i)}$$

$$F_{(ii)}$$

$${\hat{A}}_l=\sqrt{\frac{(1-\tau {)}^l}{l!}}{\tau }^{\hat{n}/2}{\hat{a}}^l$$

$$\hat{\rho }=\sum _{l_1,\dots ,l_n=0}^{\infty }{\hat{A}}_{l_1}\dots {\hat{A}}_{l_N}|{\Psi }_0⟩⟨{\Psi }_0|{\hat{A}}_{l_N}^{†}\dots {\hat{A}}_{l_1}^{†}.$$

$$\hat{\rho }$$

$$\left(\begin{array}{c}{\hat{b}}_1\\ {\hat{b}}_2\end{array}\right)=\left(\begin{array}{cc}\sqrt{\tau }& \sqrt{1-\tau }\\ -\sqrt{1-\tau }& \sqrt{\tau }\end{array}\right)\left(\begin{array}{c}{\hat{a}}_1\\ {\hat{a}}_2\end{array}\right).$$

$$⟨m|-\sqrt{1-\tau }\alpha ⟩|\alpha \sqrt{\tau }⟩=\frac{{\left(-\sqrt{1-\tau }\right)}^m}{\sqrt{m!}}{\alpha }^m{e}^{-(1-\tau )\frac{|\alpha {|}^2}{2}}|\alpha \sqrt{\tau }⟩,$$

$${\hat{P}}_{-m}[|\alpha ⟩]$$

$$c_s$$

$$|{\alpha }_{-m}⟩$$

$$\frac{1}{{\pi }^N}\int d^{2N}\vec{\alpha }|\vec{\alpha }⟩⟨\vec{\alpha }|=\frac{1}{(2\pi {)}^N}\int d^N{\vec{q}}_{\alpha }{d}^N{\vec{p}}_{\alpha }|\vec{\alpha }⟩⟨\vec{\alpha }|=\frac{1}{(2\pi {)}^N}\int d^{2N}{\vec{x}}_{\alpha }|\vec{\alpha }⟩⟨\vec{\alpha }|=I$$

$$|\Psi ⟩=\frac{1}{(2\pi {)}^N}\int d^{2N}{\vec{x}}_{\alpha }⟨\vec{\alpha }|\Psi ⟩|\vec{\alpha }⟩=\int d^{2N}{\vec{x}}_{\alpha }K({\vec{x}}_{\alpha })|\vec{\alpha }⟩,$$

Full text

I Introduction

Gaussian states of bosonic modes—quantum states of light that can be prepared using quadrature squeezed light and passive linear optics—form an important set of quantum states whose elegant mathematical description Ferraro et al. (2005) and feasibility of experimental production Yoshikawa et al. (2016a) make Gaussian quantum information processing a major success Weedbrook et al. (2012). However, it is well known that Gaussian states and Gaussian measurements (homodyne and heterodyne detection) do not constitute a universal set, i.e., resources that would allow universal quantum computation Lloyd and Braunstein (1999). Moreover, various important protocols for quantum enhanced information processing cannot be performed when restricted to Gaussian states, Gaussian unitaries, and Gaussian measurements alone. Such no-go theorems have appeared for universal quantum computing Bartlett et al. (2002), entanglement distillation Eisert et al. (2002); Fiurášek (2002); Giedke and Ignacio Cirac (2002), optimal cloning of coherent states Cerf et al. (2005), optimal discrimination of coherent states Takeoka and Sasaki (2008); Tsujino et al. (2011); Wittmann et al. (2010), receivers for optical communications Namiki et al. (2014a), quantum error correction Niset et al. (2009), quantum-enhanced sensing Gagatsos et al. (2016), and quantum repeaters Namiki et al. (2014b).

Therefore, having access to non-Gaussian states becomes imperative in pretty much any application of quantum enhanced photonic information processing. Introducing non-Gaussianity into an optical system can be challenging. For example, large $\chi^{(3)}$ non-linearities are very difficult to be implemented at optical frequencies, and obtaining a strong-enough non-Gaussian interaction through a $\chi^{(2)}$ medium with a depleted pump Wolinsky and Carmichael (1988) is hard. An alternative way to inject non-Gaussianity is to utilize detection-induced, often probabilistic, methods such as photon number subtraction. Theoretical, numerical, and experimental studies Dakna et al. (1997); Averchenko et al. (2016); Takahashi et al. (2008); Tualle-Brouri et al. (2009); Marek et al. (2008); Ra et al. (2017); Dufour et al. (2017); Barnett et al. (2018); Glancy and de Vasconcelos (2008); Ra et al. (2019) have shown that photon subtraction on a single-mode Gaussian (squeezed vacuum) state yields approximations of coherent cat-states and have validated the non-Gaussian character of photon-subtracted multi-mode states. Further, photon subtraction has been shown to enhance entanglement Opatrný et al. (2000); Kitagawa et al. (2006); Navarrete-Benlloch et al. (2012), and the fidelity of continuous-variable teleportation as was originally shown Opatrný et al. (2000) and also later studied  Seshadreesan et al. (2015).

Evaluating the state obtained after subtracting $m$ photons from a state $|\psi\rangle$, i.e., $\hat{a}^{m}|\psi\rangle$ using the photon number (Fock) basis $\{|n\rangle\}$, and even methods using the Husimi $Q$ representation of the state $|\psi\rangle$ lead to onerous calculations. This is because one has to calculate expressions such as $\hat{a}^{m}|n\rangle$ and $\hat{a}^{\dagger m}|\alpha\rangle$ within difficult-to-handle summations and integrals, where $\hat{a}$ is the modal photon annihilation operator. Similar difficulties apply when using the Wigner representation. Furthermore, taking into account the deviations of the photon-subtracted state from $\hat{a}^{m}|\psi\rangle$ pursuant to actual experimental methods of implementing such operation using a beamsplitter and photon number resolving (PNR) detectors, creates additional complexities. Despite photon number subtraction being a very promising tool for non-Gaussian state engineering, this analytical difficulty has come in the way of theoretical progress in the field.

In this work, we remedy the above situation by expressing the state on the coherent basis Bargmann (1961); Klauder and Sudarshan (1968). Specifically, we utilize the positive $P_{+}$ representation of a quantum state, which is essentially expressing a general density operator in the coherent-state over-complete basis Drummond and Gardiner (1980). This representation always exists unlike the Glauber-Sudarshan $P_{\textrm{GS}}$ function, which is not always defined, especially for squeezed states which are of interest in creation of cat states. The $P_{+}$ representation has been utilized for the numerical and analytical study of Fokker-Planck equations of dynamical systems Drummond and Gardiner (1980); Wolinsky and Carmichael (1988); Zhu and Lu (1989); Schack and Schenzle (1991); Gilchrist et al. (1997); yi Fan and Xiao (1998); Olsen and Bradley (2009), Ising systems Barry and Drummond (2008), and single-mode quantum information analyses Dodonov et al. (1994). Formulas for $P_{+}$ have been given for Gaussian states Hel (1976) but only for the cases where the Glauber-Sudarshan $P_{\textrm{GS}}$ function is well defined.

We first define the $P_{+}$ function of an $N$-mode pure Gaussian state, which we call the $K$ function. It is a unique representation of any pure state, and can be interpreted as a square root of the Husimi $Q$ function up to a phase. The latter mathemnatical obesrvation, allow us to derive clean, closed form, and easy to use formulas for the $P_{+}$ representaion (called $K$ in this work), for any Gaussian state. We begin with developing a closed-form expression of the $K$ function of a general $N$-mode Gaussian state. This lets us analytically characterize non-Gaussian states created by photon number detection and/or photon number subtraction on a subset of modes of any $N$-mode Gaussian state in an analytic integral form. We show that this reproduces—in a rather simple set of steps—the theory behind Gaussian boson sampling, where it was argued that sampling from the photon number distribution of a random $N$-mode entangled Gaussian state is a classically-hard computational task as was proven in Lund et al. (2014) and also subsequntly studied Hamilton et al. (2017); Quesada (2018). As a new example application of our formalism, we consider the problem of engineering coherent cat basis entangled cluster states. We propose a method to prepare a two-mode cat-basis Bell state by subtracting photons from both modes of a Gaussian two-mode entangled squeezed state. We show that the fidelity versus success probability trade-off of our method is higher than that of the conventional method—that of splitting an approximate single-mode cat state, obtained by photon number subtraction on a squeezed vacuum mode, on a balanced beamsplitter. The above analysis would be extremely cumbersome (and not scalable to a larger entangled state) if done in the traditional way in the photon number basis. We expect generalization of the above, to enable exploration of efficient generation of cat-basis cluster states, which have recently emerged as a very powerful resource for quantum error correction against photon losses, with applications both to photonic quantum repeaters as well as superconducting quantum computing Michael et al. (2016); Li et al. (2017); Albert et al. (2018).

II The $K$ function of a pure Gaussian state

We work in units of $\hbar=1$, where $N$-mode vacuum state’s covariance matrix (CM) is $V_{0}=I/2$, with $I$ being the $N$-mode identity operator. Coherent states of $N$ modes $|\vec{\alpha}\rangle$ are not mutually orthogonal. Yet they form an over-complete basis. In other words, they resolve the identity operator, viz.,

[TABLE]

where ${\vec{x}_{\alpha}^{T}=(\vec{q}_{\alpha}^{T},\vec{p}_{\alpha}^{T})}$, and the volume element $d^{2N}\vec{x}_{\alpha}=dq_{\alpha_{1}}\ldots dq_{\alpha_{N}}dp_{\alpha_{1}}\ldots dp_{\alpha_{N}}$. We take $\alpha_{i}=(q_{\alpha_{i}}+ip_{\alpha_{i}})/\sqrt{2}$. Using Eq. (1), we can express any $N$-mode pure state $|\Psi_{0}\rangle$ as

[TABLE]

where we call $K(\vec{x}_{\alpha})=1/(2\pi)^{N}\langle\vec{\alpha}|\Psi_{0}\rangle$ the $K$-function of the state $|\Psi_{0}\rangle$. When compared to the $Q$ function $Q(\vec{x}_{\alpha})=1/(2\pi)^{N}|\langle\vec{\alpha}|\Psi_{0}\rangle|^{2}$, the $K$-function resembles something that could be called the square root of the $Q$ function. However, one has to be careful as $\langle\vec{\alpha}|\Psi_{0}\rangle$ is a complex number and its square root will contain a phase that if omitted will produce wrong results since it depends on $\vec{x}_{\alpha}$.

Let us assume that $|\Psi_{0}\rangle$ is a zero-mean Gaussian state, such that $Q(\vec{x}_{\alpha})=1/(2\pi)^{N}\langle\vec{\alpha}|\Psi_{0}\rangle\langle\Psi_{0}|\vec{\alpha}\rangle$ is a Gaussian function. In order to calculate the $K$ function, one must break up the Gaussian $Q$ function’s exponent into two conjugate parts, yielding a Gaussian $K$ function. This step becomes easier if instead of working with Cartesian coordinates $(\vec{q}_{\alpha},\vec{p}_{\alpha})$ we move to complex coordinates $(\vec{\alpha},\vec{\alpha}^{*})$ with a $\pi/4$ phase space rotation. After we finish the calculation we rotate back to Cartesian coordinates.

Let us now consider a general $N$-mode Gaussian pure state $|\Psi\rangle=D(\vec{\beta})|\Psi_{0}\rangle$, where $D(\vec{\beta})$ is the displacement operator. With $|\Psi_{0}\rangle$ expressed in its $K$-function form (2), it is straightforward to evaluate the $K$-function of $|\Psi\rangle$ since $D(\vec{\beta})\left|\vec{\alpha}\right>=\exp\left(\vec{\beta}\vec{\alpha^{*}}-\vec{\beta^{*}}\vec{\alpha}\right)|{{\vec{\alpha}}+{\vec{\beta}}}\rangle$.

Using the above method, we show that any $N$-mode pure Gaussian state with CM $V$ and displacement vector $\vec{x}^{T}_{\beta}=(\vec{q}^{T}_{\beta},\vec{p}^{T}_{\beta})$ 111we work in the $qqpp$ representation, i.e., the upper left (lower right) block of the CM concerns position (momentum), while the off-diagonal blocks hold information of correlations thereof. can be written as follows (see App. Sec. 2 and 3 for the complete derivation),

[TABLE]

where,

[TABLE]

with $\Gamma=V+I/2$, and

[TABLE]

where $A=A^{T},B=B^{T},\ C$ are defined as the blocks of the CM $\Gamma$ defined as follows 222Since the CM $V$ is symmetric, $\Gamma$ and $\Gamma^{-1}$ will be symmetric.:

[TABLE]

III Photon subtraction from a general multi-mode Gaussian state.

Subtraction of $m$ photons from a single-mode quantum state $|\psi\rangle$ can be implemented by transmitting $|\psi\rangle$ through a beam splitter of transmissivity $\tau$ (chosen to be close to $1$) while detecting the low-transmissivity output of the beam splitter with a PNR detector. If the detector registers $m$ photons, the transmitted state projects to $\mathcal{P}_{-{m}}\left[|\psi\rangle\right]$, which is an approximation of the $m$-photon subtracted state $\hat{a}^{m}|\psi\rangle$. Since $\hat{a}^{m}$ is not a unitary, photon subtraction only succeeds probabilistically.

Let us consider subtracting a vector $\vec{m}=(m_{1},\ldots,m_{N})$ photons from an $N$-mode pure Gaussian state $|\Psi\rangle$ using an array of beam splitters of transmissivities $\tau_{i}$, and PNR detectors. The post-subtraction state will be denoted $\mathcal{P}_{-\vec{m}}[|\Psi\rangle]$, implying $m_{i}$ photons were subtracted from the $i$-th mode, $i=1,2,\ldots,N$. Using the $K$ function of $|\Psi\rangle$ (3), we see that $\mathcal{P}_{-\vec{m}}$ acts only on the coherent states (see App. Sec. 1), i.e., $\mathcal{P}_{-\vec{m}}[|\vec{\alpha}\rangle]$, which assumes a simple form, $\mathcal{P}_{-\vec{m}}[|\vec{\alpha}\rangle]=\prod_{i=1}^{N}c_{i}|\alpha_{i}\sqrt{\tau_{i}}\rangle,\ c_{i}=[\left(-\sqrt{1-\tau_{i}}\right)^{m_{i}}]/[\sqrt{m_{i}!}]\alpha^{m_{i}}e^{-(1-\tau_{i})|\alpha_{i}|^{2}/2}$.

The photon subtracted state $|\Psi_{-\vec{m}}\rangle$ is given as:

[TABLE]

where $P=\langle\Psi_{-\vec{m}}|\Psi_{-\vec{m}}\rangle$ is the probability of success of the $N$-mode vector photon subtraction. $P$ is a $4N$-dimensional integral with the elementary volume $d^{2N}\vec{x}_{\alpha}d^{2N}\vec{x}_{\gamma}$ ($\vec{x}_{\gamma}$ are the coordinates of $\langle\Psi_{-\vec{m}}|$), with a Gaussian kernel, and polynomial terms $(q_{\alpha_{i}}+ip_{\alpha_{i}})^{m_{i}}(q_{\gamma_{i}}-ip_{\gamma_{i}})^{m_{i}}$. This kind of integrals can be analytically calculated (see App. Sec. 4).

If one wishes to use photon subtraction to produce a desired non-Gaussian multimode entangled state $|C\rangle$ (for example a cat-basis Bell state that we consider later), one can evaluate analytically the fidelity $F=|\langle C|\Psi_{-\vec{m}}\rangle|^{2}$ between the desired state $|C\rangle$ and the actual state obtained $|\Psi_{-\vec{m}}\rangle$ if $|C\rangle$ is expressed in its $K$ function form. For cat states for example, which are superpositions of coherent states $|\vec{\gamma}\rangle$, the fidelity calculation will require us to calculate the amplitude $|\langle\vec{\gamma}|\Psi_{-\vec{m}}\rangle|$, which again is a $4N$-dimensional integral, with Gaussian kernel and polynomial terms $(q_{\alpha_{i}}+ip_{\alpha_{i}})^{m_{i}}$, which can be analytically calculated (see App. Sec.5).

For the rest of this paper we will restrict our attention to zero-mean states, to keep the exposition simple. Including non-zero means is a trivial extension. Further, we will assume that all the beamsplitters employed for photon subtraction on an $N$-mode Gaussian state have the same transmissivity, $\tau$.

IV Gaussian boson sampling and non-Gaussian state engineering

Consider a pure $N$-mode Gaussian state $|\Psi\rangle$, the first $M<N$ modes of which are detected using PNR detectors, obtaining the outcome ${\vec{n}}=(n_{1},\ldots,n_{M})$. It is simple to show that the resulting state $|\Phi\rangle$ on the unmeasured modes is given by (see App. Sec. 6),

[TABLE]

where we used $|\alpha\rangle=\exp(-|\alpha|^{2}/2)\sum_{n}\alpha^{n}/(\sqrt{n!})|n\rangle$. The probability $P_{M}$ of detecting the photon number pattern $\vec{n}$ and hence heralding the state $|\Phi\rangle$, can be calculated by setting $\langle\Phi|\Phi\rangle=1$.

Gaussian boson sampling is the special case of $M=N$, where all $N$ modes are detected Hamilton et al. (2017); Quesada (2018). The success probability of detecting a photon-number pattern $\vec{n}$, $P_{\vec{n}}=|\langle\vec{n}|\Psi\rangle|^{2}=|\langle n_{1}\ldots n_{N}|\Psi\rangle|^{2}$ can be evaluated using our formalism, and shown to be (see App. Sec. 7),

[TABLE]

where,

[TABLE]

and $\mathcal{H}=\mathcal{B}+I/2$. Since $\mathcal{H}=\mathcal{H}^{T}$ and its real part is positive definite (see App. Sec. 8), Eq. (14) is a proper Gaussian distribution. Therefore, Eq. (13) is the mean value $\langle f_{1}^{n_{1}}\ldots f_{N}^{n_{N}}\rangle$, where $f_{i}=q_{\alpha_{i}}+ip_{\alpha_{i}}$, under the distribution of Eq. (14). Using Wick’s theorem Zvonkin (1997); Luque and Thibon (2002) we can write it as,

[TABLE]

where $\Sigma=\sum_{i=1}^{N}n_{i}$ and $\textrm{Hf}\left(F\right)$ is the Hafnian of the matrix $F$ with elements $F_{ij}=\langle f_{i}f_{j}\rangle,\ 1\leq i,j\leq\Sigma$.

V Photon subtraction from multi-mode squeezed cluster states

Continuous variable (CV) quantum computing is a field that explores the use of multimode entangled squeezed states for all-photonic quantum computing. Such Gaussian cluster states of thousands of modes have been prepared experimentally Chen et al. (2014); Zhang et al. (2017); Yoshikawa et al. (2016b). It is known however that Gaussian cluster states by themselves are not a resource sufficient for universal quantum processing. Photon number detection being the most practical “de-Gaussification” tool, and given it is known that approximate cat states can be prepared using photon number subtraction from a single-mode squeezed vacuum, we will explore the creation of cat-basis cluster (graph) states by photon number subtraction on Gaussian cluster states.

Let us consider the Gaussian graph state $|G\rangle$ which is the result of the unitary evolution of an $N$-mode vacuum state under the unitary $\hat{U}_{r}=\exp{(-ir\hat{H})}$ whose generating Hamiltonian is,

[TABLE]

where $\hat{a}_{i}$ and $\hat{a}_{i}^{\dagger}$ are the annihilation and creation operators of the $i$-th mode respectively. The state $|G\rangle$ is a squeezed entangled state among its $N$ modes. The information about which modes are entangled is described by the graph (a symmetric matrix) $G$. We assume that the squeezing parameter $r>0$ is the same for all modes 333This assumption can be dropped just by having a symmetric matrix $G$ whose elements are the different values of squeezing parameter.. In the limit $r\to\infty$, $|G\rangle$ is a continuous variable cluster state if $G$ is a full rank matrix  Menicucci et al. (2011). For the same $r$, we will consider a matrix $G$ which is its own inverse, i.e., $G^{2}=I$. Under this assumption on $G$, we show that (see App. Sec. 9),

[TABLE]

To demonstrate the power of our method, as a first example we consider a two-mode squeezed vacuum state (TMSV), from which we subtract five photons per mode (ten in total). We calculate the photon subtracted state $|\Psi_{-5,-5}\rangle$, the probability of success $P_{5}$, and the fidelity $F_{5}=|\langle C|\Psi_{-5,-5}\rangle|^{2}$, where

[TABLE]

with normalization $|N_{+}|^{2}=2[1+e^{-2(q_{\gamma}^{2}+p_{\gamma}^{2})}]$. We compare the state $|\Psi_{-5,-5}\rangle$ with the specific state of Eq. (20), because both states are parity $(-1)^{\hat{n}_{1}+\hat{n}_{2}}$ eigenstates with eigenvalue $+1$. If the $K(\vec{x}_{\alpha})$ function is known, then the state $|\Psi_{-5,-5}\rangle$ is known from Eq. (10) for zero displacement. The only thing required to find the $K(\vec{x}_{\alpha})$ is the matrix $\mathcal{B}$ 444The matrix $\Gamma^{-1}$ can be easily found to be $\Gamma^{-1}=\mathcal{B}+\mathcal{B}^{\dagger}$, from that we calculate $\det\Gamma=\cosh^{4}r$., which is given by Eq. (19) for,

[TABLE]

which describes the graph corresponding to the TMSV, as can also be seen by Eq. (18). The probability $P_{5}$ and the fidelity $F_{5}$ are given by:

[TABLE]

where $p(\mu)=[(1+\mu^{2})(1+24\mu^{2}+76\mu^{2}+24\mu^{6}+\mu^{8})]/[(1-\mu^{2})^{11}]$, $w(z)=1+(5z)/(2)+(5z^{2})/(4)+(5z^{3})/(24)+(5z^{4})/(384)+(z^{5})/(3840)$, $\mu=\tau\tanh r$, and $z=(q_{\gamma}-ip_{\gamma})\mu$. For example for $q_{\gamma}=0.5$, $p_{\gamma}=0$, $\tau=0.01$, and $r=1$ we get $P_{5}=0.025$ and $F_{5}=0.979$. Note that in the above example, the analytical complexity would not have changed if we decided to subtract more (e.g., $10$ photons) from each mode, whereas a traditional Fock basis calculation would become completely intractable.

As a second example we consider two ways to produce the cat-basis Bell state $|C\rangle$: (i) a single-mode squeezed state from which we subtract two photons and the resulting state is known to be an approximation of the cat state $|\delta\rangle+|-\delta\rangle$, which if then split in a 50-50 beam splitter, is known to produce the state $|C\rangle$ with $\delta=\sqrt{2}\gamma$ Ralph et al. (2003). In scenario (ii) we subtract one photon from each of the two modes of a TMSV. In both scenarios two photons are subtracted in total. Also, the beam splitter used in scenario (i) is a Hadamard gate, which if used to mix a position-squeezed state with a momentum-squeezed state we get a TMSV, see Fig. 1 We set $p_{\gamma}=0$ and we calculate the probabilities of success $P_{(i)},\ P_{(ii)}$ and the fidelities $F_{(i)},\ F_{(ii)}$ for scenarios (i) and (ii) to the desired state $|C\rangle$, as:

[TABLE]

Comparative results for these two scenarios are shown in Figs. 3 and  3. To produce a cat-basis Bell state $|C\rangle$ with a small amplitude, scenario (ii) is better than (i) in both fidelity and probability of success. As the amplitude of $|C\rangle$ increases, the situation begins to change: scenario (i) favors high fidelity, at the expense of smaller probability of success compared to scenario (ii). For example, for $q_{\gamma}=0.1,\ p_{\gamma}=0,\ r=0.9,\ \tau=0.4$, $P_{(i)}=0.179,\ P_{(ii)}=0.249,\ F_{(i)}=0.999,\ F_{(ii)}=0.999$ and for $q_{\gamma}=1,\ p_{\gamma}=0,\ r=0.9,\ \tau=0.4$, $P_{(i)}=0.093,\ P_{(ii)}=0.126,\ F_{(i)}=0.990,\ F_{(ii)}=0.806$.

It is of similar ease to find expressions for $P_{5},\ F_{5},\ P_{(i)},\ P_{(ii)},\ F_{(i)},\ F_{(ii)}$ for $p_{\gamma}\neq 0$ (generality is not lost by assuming real amplitude).

VI Mixed Gaussian states

A mixed Gaussian state $\hat{\rho}$ can be written as $\hat{\rho}=\hat{U}\hat{\rho}_{\textrm{th}}\hat{U}^{\dagger}$, where $\hat{\rho}_{\textrm{th}}$ is a thermal state and $\hat{U}$ is a Gaussian unitary. Using the Glauber-Sudarshan $P_{\textrm{GS}}$ function of the thermal state $P_{\textrm{GS,th}}\equiv P_{\textrm{th}}(\vec{x}_{\alpha})$, we have $\hat{\rho}=\int d^{2N}\vec{x}_{\alpha}P_{\textrm{th}}(\vec{x}_{\alpha})|\Psi\rangle\langle\Psi|$ where $|\Psi\rangle=U|\vec{\alpha}\rangle$. The state $|\Psi\rangle$ can be expressed using Eq. (3) and therefore $\hat{\rho}$ is expressed in the coherent-state basis as two integrals coming from $|\Psi\rangle$ and $\langle\Psi|$ are convoluted into a third integral over $\vec{x}_{\alpha}$ with $P_{\textrm{th}}(\vec{x}_{\alpha})$.

Concerning mixed Gaussian states, things become even easier if an initial pure Gaussian sate $|\Psi_{0}\rangle$ goes through a pure loss channel. We remind the reader that under a pure loss channel, every mode of the state $|\Psi_{0}\rangle$ is coupled with $|0\rangle$ (the environment) via a beam splitter of transmittance $\tau_{i}$, where $i=1,\ldots,N$ counts the modes, i.e., the loss does not have to be uniform across the $N$ modes. Then the environment’s output is traced out. The single-mode pure loss channel is described by the Kraus operators Ivan et al. (2011),

[TABLE]

and the final state is,

[TABLE]

Here we observe that if $|\Psi_{0}\rangle$ is expressed on the coherent basis, the operators $\tau^{\hat{n}/2}\hat{a}^{l}$ in Eq. (28), will act on coherent states resulting to managable expressions. For further simplicity we assume the same transimittance rate $\tau$ per mode (even thouhgh this assumption can be easily dropped). The final state will be,

[TABLE]

an expresion which can be useful, for example, in an analysis of a Gaussian boson sampling with pure loss scheme.

VII Conclusions

We have derived a general representation of Gaussian states in the coherent-state basis, and showed that it opens the door to analytical and thorough investigations of non-Gaussian states prepared via photon subtraction and partial PNR detection of Gaussian states. We showed a simplified analysis of Gaussian boson sampling as a special case of our formalism. As a specific example application of our formalism, we considered cat-basis cluster creation by multi-mode photon subtraction on entangled Gaussian states. We showed that by subtracting photons simultaneously from both modes of a two-mode squeezed vacuum state, a coherent cat basis Bell state can be produced with higher fidelity and probability of success, compared to the well-known method of first creating a cat state via photon number subtraction of a single-mode squeezed vacuum, followed by linear-optical manipulation. The question on whether more general coherent cat basis graph states—known to be an excellent resource for quantum error correction against photon loss—can be systematically engineered from Gaussian cluster states and photon subtraction, is left open for future work. We anticipate that our formalism will prove a powerful tool for non-Gaussian cluster state engineering Arzani et al. (2018); Walschaers et al. (2018); Sabapathy et al. (2018), which is a subject of intense interest in designing scalable solutions for all-photonic quantum computing and other forms of quantum-enhanced photonic information processing such as all-photonic quantum repeaters where photonic cluster states replace the role of quantum memories Azuma et al. (2015); Pant et al. (2017), and optical-domain quantum machine learning via receivers powered with cluster states Zhuang and Zhang (2019).

While preparing this paper, it came to our attention [PrivatecommunicationwithXanadu]private that similar phase space methods have been developed Su et al. (2019a, b) practically concurrently.

Acknowledgements.

CNG was supported by the Army Research Office (ARO) STIR program, contract number W911NF-18-1-0377. SG acknowledges Xanadu Quantum Technologies for supporting multiple useful discussions on this topic. The authors acknowledge Daiqin Su, Krishna Kumar Sabapathy, Hari Krovi, Raf Alexander, and Kaushik Seshadreesan for valuable discussions.

Appendix

1 Photon subtraction from a coherent state using a beam splitter

Subtraction of $m$ photons from a mode of a state $|\Psi\rangle$ can be implemented with a beam splitter of transmittance $\tau$. The beam splitter couples the mode that the photon subtraction will take place with vacuum. Then, if the photon number resolution measurement (PNRM) registers $m$ photons, the resulting state is $|\Psi_{-m}\rangle$ as shown in Fig. A1. Since a measurement is involved, this procedure is probabilistic and heralded. Because of the probabilistic nature of photon subtraction the final state needs to be normalized. The absolute square of the normalization is the probability of finding $m$ photons in the PNRM. This probability is also called the probability of success.

Since we expand $|\Psi\rangle$ on coherent basis, when subtracting photons from some mode of $|\Psi\rangle$, the beam splitter will couple a coherent state with vacuum. If $\hat{a}_{1},\ \hat{a}_{2}$ and $\hat{b}_{1},\ \hat{b}_{2}$ are the input and output annihilation operators respectively we have,

[TABLE]

Therefore, if the global, two-mode input state is $|\alpha,0\rangle$ the final state is $|\sqrt{\tau}\alpha,-\sqrt{1-\tau}\alpha\rangle$. The conditional state on the upper output port, upon finding $m$ photons in the PNRM, is

[TABLE]

therefore we can write,

[TABLE]

where

[TABLE]

and the probability of success is given by $P=|c_{s}|^{2}$. Subtracting photons from a coherent state yields the same amplitude-damped coherent state $|\alpha\sqrt{\tau}\rangle$ regardless of the PNRM result. Therefore, for applications there is not much meaning in subtracting photons from coherent states. However, it is highly convenient for mathematical manipulation of photon subtraction written on coherent basis.

2 Coherent basis representation of pure Gaussian states without displacement

We define $\vec{x}_{\alpha}=\begin{pmatrix}\vec{q}_{\alpha}\\ \vec{p}_{\alpha}\end{pmatrix}$ and we work with $\hbar=1$. Using the unit resolution on coherent states,

[TABLE]

for any state $|\Psi\rangle$ we can write,

[TABLE]

where we define,

[TABLE]

which up to some constant is the the square root of the $Q(\vec{x}_{\alpha})$ representation,

[TABLE]

therefore we can write,

[TABLE]

such that,

[TABLE]

Equations (A10) and (A11) imply that to find $Q_{1/2}(\vec{x}_{\alpha})$, we have to separate the $Q(\vec{x}_{\alpha})$ representation into a product of two conjugate parts. In that way, if the state $|\Psi\rangle=|\Psi_{0}\rangle$ is Gaussian state with zero displacement, we can express $K(\vec{x}_{\alpha})$ as a function of the states’ covariance matrix (CM). The $Q(\vec{x}_{\alpha})$ representation of a Gaussian state with CM $V$ is,

[TABLE]

where,

[TABLE]

where $I$ is the identity matrix of appropriate dimensions. Any CM is a real, symmetric matrix, and as per Eq. (A13) $\Gamma$ is a real, symmetric matrix. The inverse of a real, symmetric matrix is again real and symmetric, therefore in block form the matrix $\Gamma^{-1}$ is,

[TABLE]

where $A=A^{T}$ and $B=B^{T}$ and $A,\ B,\ C$ real. It is more convenient if we change coordinates in the following manner,

[TABLE]

where,

[TABLE]

Note that $R$ is unitary, i.e., $RR^{\dagger}=I$.

To break Eq. (A12) into two conjugate parts, we must express the term $\vec{x}_{\alpha}^{T}\Gamma^{-1}\vec{x}_{\alpha}$ which appears in its $\exp(.)$ as a summation of two conjugate terms. To this end we express $\Gamma$ in the $\vec{z}_{\alpha}$ basis,

[TABLE]

where,

[TABLE]

is the transformed $\Gamma^{-1}$ in the $\vec{z}_{\alpha}$ basis. From Eqs. (A22) and (A20) we get,

[TABLE]

Therefore we can write,

[TABLE]

Equation (A25) shows that we can readily derive the two conjugate terms where,

[TABLE]

Going back to Cartesian coordinates $\vec{x}_{\alpha}$ we get the matrix $\mathcal{B}$,

[TABLE]

where,

[TABLE]

where we have used Eqs. (A20) and (A26). Therefore given the CM $V$ of a pure Gaussian state $|\Psi_{0}\rangle$, we can find $\Gamma^{-1}$ and from that we can immediately write $\mathcal{B}$ and the expansion on coherent basis is,

[TABLE]

3 Coherent basis representation of pure Gaussian states with displacement

A displaced pure Gaussian state $|\Psi\rangle$ can be derived by applying a displacement $D(\vec{x}_{\beta})$ and multiple-mode squeezing $S(\vec{r})$ (phases can be absorbed into the squeezing operator) Weedbrook et al. (2012) onto a multiple-mode vacuum state $|\vec{0}\rangle$,

[TABLE]

where $|\Psi_{0}\rangle$, is the state for which we worked out its coherent basis expansion in Sec. 2. Therefore we have,

[TABLE]

where in the last step we have used $D(-\vec{\alpha})D(\vec{\beta})=e^{\frac{1}{2}\vec{\beta}\vec{\alpha}^{*}-\frac{1}{2}\vec{\beta}^{*}\vec{\alpha}}D(\vec{\beta}-\vec{\alpha})$, which acts on $\langle\vec{0}|$ and therefore the sign of the displacement should be inversed. In Eq. (A32), $\langle\vec{\alpha}-\vec{\beta}|\Psi_{0}\rangle$ is known from Eq. (A29). Additionally, by defining,

[TABLE]

Eq. (A32) is written,

[TABLE]

From Eqs. (A29) and (A35) we have,

[TABLE]

where,

[TABLE]

with

[TABLE]

4 Probability of success

The photon subtracted state is,

[TABLE]

therefore the probability of success is given by the condition $\langle\Psi_{-\vec{m}}|\Psi_{-\vec{m}}\rangle=1$. Therefore we have,

[TABLE]

By writing,

[TABLE]

Eq. (A40) gives,

[TABLE]

Equation (A42) is a Gaussian integral (represented by the $e^{-\frac{(1-\tau_{i})}{4}\left(|\vec{x}_{\alpha}|^{2}+|\vec{x}_{\gamma}|^{2}\right)-\frac{1}{4}\tau_{i}\vec{x}_{\alpha}^{T}\vec{x}_{\alpha}-\frac{1}{4}\tau_{i}\vec{x}_{\gamma}^{T}\vec{x}_{\gamma}+\frac{1}{2}\tau_{i}\vec{x}_{\gamma}^{T}\mathcal{X}\vec{x}_{\alpha}}$, $K(\vec{x}_{\alpha}),\ K^{*}(\vec{x}_{\gamma})$, kernels) with linear terms (represented by $G(\vec{x}_{\alpha},\vec{x}_{\beta}),\ G^{*}(\vec{x}_{\gamma},\vec{x}_{\beta})$), and polynomial terms $\left(q_{\alpha_{i}}+ip_{\alpha_{i}}\right)^{m_{i}}\left(q_{\gamma_{i}}-ip_{\gamma_{i}}\right)^{m_{i}}$. The way to calculate this analytically and efficiently, is to use the identity,

[TABLE]

Using Eq. (A43), we cast Eq. (A42) into a Gaussian integral, i.e.,there is only an exponential and no polynomial terms, with extra lineal terms $\lambda_{i}\left(q_{\alpha_{i}}+ip_{\alpha_{i}}\right)+\mu_{i}\left(q_{\gamma_{i}}-ip_{\gamma_{i}}\right)$ in the exponential. Then one should take the $m_{i}-$order derivatives on the result of the Gaussian integral with respect to $\lambda_{i}$ and $\mu_{i}$ at $\lambda_{i}=\mu_{i}=0$.

5 Fidelity

For any state of the form,

[TABLE]

where $\sum_{i}|c_{i}|^{2}=1$ and $|\vec{\gamma}^{(i)}\rangle=|\gamma^{(i)}_{1}\gamma^{(i)}_{1}\ldots\gamma^{(i)}_{N}\rangle$. Note that a special example of $|\phi\rangle$ is the coherent cat state (CCS) used in the main paper. The fidelity $F=|\langle\phi|\Psi_{-\vec{m}}\rangle|^{2}$ requires the calculation of $\langle\vec{\gamma}|\Psi_{-\vec{m}}\rangle$. From Eq. (A39) we have,

[TABLE]

where the probability of success $P$ should be calculated first as per Sec. 4. We have,

[TABLE]

therefore Eq. (A45) is written as,

[TABLE]

Equation (A47), similarly to $P$ in Sec. 4, is a Gaussian integral with linear terms, and polynomial terms $\left(q_{\alpha_{i}}+ip_{\alpha_{i}}\right)^{m_{i}}$ which can be injected into the exponential of Eq. (A47) by using the identity,

[TABLE]

That way Eq. (A47) will become a Gaussian integral, upon which we take $m_{i}-$order derivatives with respect to $\lambda_{i}$ at $\lambda_{i}=0$.

6 The conditional state and its normalization

We set zero displacements, therefore we work with the $N-$mode Gaussian state $|\Psi_{0}\rangle$. Upon finding a pattern $\{n_{1},\ldots,n_{M}\},\ M<N$ at the photon number resolution measurements (PNRM) at each one of the $M$ modes, the conditional state $|\Phi\rangle$ is,

[TABLE]

The probability $P_{M}$ is given by the normalization $\langle\Phi|\Phi\rangle=1$,

[TABLE]

where we have used $\langle\gamma|\alpha\rangle=\exp(-|\gamma|^{2}/2-|\alpha|^{2}/2+\gamma^{*}\alpha)$, $|\vec{x}_{\alpha,\gamma}|^{2}=\sum_{k=1}^{N}x^{2}_{\alpha_{i},\gamma_{i}}$, and $\mathcal{X}_{kl}$ are the matrix elements of $\mathcal{X}$ of Eq. (A33) for dimensions $(N-M)\times(N-M)$. The same method using ancillary variables $\lambda_{i}$ as in Sec. 4 can be applied to calculate $P_{M}$ of Eq. (A50).

7 The probability distribution $P_{\vec{n}}$

We set zero displacements. The probability of finding a pattern $\{n_{1},\ldots,n_{N}\}$ at each one of all the $N$ modes is,

[TABLE]

From Eq. (A29) and using $\langle n|\alpha\rangle=\exp(-|\alpha|^{2}/2)\alpha^{n*}/(\sqrt{n!})$ we get,

[TABLE]

where $\mathcal{H}=\mathcal{B}+I/2$. As it is shown in Sec. 8, $\mathcal{H}$ is symmetric with positive definite real part. Therefore, the function,

[TABLE]

represents a Gaussian distribution. In that way, Eq. (A52) is written as,

[TABLE]

where $f_{i}=q_{\alpha_{i}}+ip_{\alpha_{i}}$. Mean values of the form $\langle f_{1}^{n_{1}}\ldots f_{N}^{n_{N}}\rangle$ represent Hafnians via Wick’s theorem as argued in the main paper. Equation (A54) yields a complex number result, whose absolute squared is the probability $P_{\vec{n}}$ of Eq. (A51).

8 $\mathcal{H}$ matrix is symmetric and its real part is positive definite

From Eq. (A28) and given that $A^{T}=A$ and $B^{T}=B$ we can readily see that $\mathcal{B}^{T}=\mathcal{B}$. Therefore $\mathcal{H}=\mathcal{B}+I/2$ is symmetric as well. The real part of $\mathcal{H}$ is,

[TABLE]

Since any CM $V$ is positive definite, denoted as $V>0$, then $\Gamma=V+I/2>0\Rightarrow\Gamma^{-1}>0$ since the inverse of a positive definite matrix is also positive definite.

9 $\mathcal{B}$ matrix for multiple modes squeezed states

The Hamiltonian,

[TABLE]

generates the unitary $\hat{U}_{r}=\exp{(-ir\hat{H})}$ which corresponds to the symplectic matrix,

[TABLE]

Therefore the CM is,

[TABLE]

and from Eq. (A28) we get,

[TABLE]

In Eq. (A59) the matrix $G=G^{T}$ is in the argument of $\tanh(.)$ which denotes,

[TABLE]

For a self-inverse matrix $G=G^{-1}$, i.e., $G^{2}=I$, we expand $e^{2rG}$ in Taylor series. That way we get,

[TABLE]

From Eqs. (A60), (A61), and (A62), we have,

[TABLE]

We have that,

[TABLE]

Equations (A63) and (A65) give,

[TABLE]

From Eqs. (A59) and (A66) we find,

[TABLE]