Deterministic multi-mode nonlinear coupling for quantum circuits

The fundamental practical applicability of quantum systems largely depends on the dimension of the Hilbert spaces used for their description. Quantum systems with large dimension are prohibitively difficult to classically simulate but, in turn, their controlled dynamics can provide us with answers we cannot reach classically, be it for the lofty quantum computation [1, 2], or the more approachable quantum simulation [3–5] or metrology [6–8]. Large dimension can be obtained by parallelization—combining several quantum systems and treating them as one. The most straightforward example are the recently applied multi-qubit architectures with photons [9], trapped ions [10], or superconducting systems [11]. A complementary approach uses infinite dimensional Hilbert space of quantum modes propagating in circuits to reduce necessary parallelism. However, to optimally utilize this potential, even the already infinite dimensional harmonic oscillators can benefit from multi-oscillator architecture [12, 13].

The processing power of such circuit networks is largely given by the amount of control we can exert over the systems. The control comes in two broad varieties—local manipulation of individual modes and interaction coupling of different modes. Their availability is given by the physical nature of the harmonic oscillators. For systems in the form of electromagnetic field modes at microwave or optical frequencies, the local operations start as simple time delay elements that introduce phase shifts and progress up to to linear quantum amplifiers [14]. Similarly, the multi-mode interactions start with a basic linear coupler, for light represented by a partially transparent mirror, and again progress up to linear multi-mode interactions capable of generating Gaussian quantum entanglement [14]. These elementary linear operations can be realized taking advantage of existing materials. However, to gain the full benefit of technologies based on quantum systems we need to be able to implement operations that are far from linear.

In principle, an arbitrary, even highly nonlinear, operation on modes of propagating radiation can be decomposed into a sequence of elementary unitary gates [15, 16]. Majority of these gates implement Gaussian operations, both single-mode and two-mode [14]. The Gaussian operations can be considered well accessible in contemporary experiments; their largest disadvantage is that alone they are not sufficient [17]. To generate an arbitrary operation we also need access to non-Gaussian gates that employ nonlinearity of higher order to deterministically transform pure Gaussian states into pure non-Gaussian ones. One of the simplest gates of this kind is the cubic gate which implements an unitary operation with Hamiltonian proportional to the third power of a quadrature operator. An individual cubic gate can be implemented in the measurement induced fashion [18] in which the required nonlinearity is imprinted on the target system through Gaussian interaction with a suitably prepared ancilla and subsequent feed-forward. The biggest obstacle is the availability of the requisite ancillary states. They can be prepared approximately [19–22] but these approximations are imperfect and add unwanted noise during implementing the operation. This extra noise makes repeated application of the gates impractical. For single mode operations the number of gates can be reduced by directly implementing non-Gaussian operations of higher order [23]. Our aim is to provide similar methodology for non-Gaussian interactions between the propagating modes and thus cover all possible complex quantum circuits with the infinite Hilbert spaces.

In this paper we present a general description of a class of non-Gaussian couplings between several propagating field modes together with a measurement induced method for experimental implementation in quantum circuits. The method is a final generalization of the original GKP proposal [18]. The architecture of the N-mode coupling is based on N sequential Gaussian quantum non-demolition interactions with N auxiliary modes prepared in a specific nonlinear state, homodyne detection of the ancillary modes, and nonlinear feed-forward control applied to the modes that remain. We consider a third-order non-Gaussian coupling as an example and discuss several different approaches towards the entangled resource state generation. Ultimately we show that genuinely non-Gaussian coupling between several harmonic oscillators is feasible, which opens up the avenue of both quantum technologies and fundamental research of highly nonlinear quantum mechanics and thermodynamics.

A general unitary coupling between N modes can be characterized by interaction Hamiltonian expressed as a functional ${V}_{{xp}}({\hat{x}}_{1},{\hat{p}}_{1},\,\cdots ,\,\,{\hat{x}}_{N},{\hat{p}}_{N})$, where ${\hat{x}}_{j}$ and ${\hat{p}}_{j}$ are quadrature operators of the individual modes with commutation relations $[{\hat{x}}_{j},{\hat{p}}_{k}]={\rm{i}}{\delta }_{{jk}}$.The required Hermitian nature of the Hamiltonian demands the underlying real function ${V}_{{xp}}(.)$ to be symmetrical with respect to exchanging operators ${\hat{x}}_{j}$ and ${\hat{p}}_{j}$ for all j. Implementation of such general operation is still an open problem. However, the mono-quadrature subclass consisting of operation with interaction Hamiltonians that are real functions of only single type of quadrature operators, ${V}_{x}({\hat{x}}_{1},\,\cdots ,\,\,{\hat{x}}_{N})$ can be implemented by generalizing the original proposal for cubic gates [18]. Hamiltonians of such operations are implicitly Hermitian and no further constrains on the form of the function are required. When used as part of larger circuits over the elementary single mode cubit gates, these more complex operations can significantly decrease the overall number of gates required and thus simplify the realization of quantum circuits [16, 23].

Any operation from this class transforms an unknown quantum state $| \psi \rangle$ into ${{\rm{e}}}^{-{{\rm{i}}{V}}_{x}({\hat{x}}_{1},\cdots ,{\hat{x}}_{N})}| \psi \rangle$. It can be implemented by the measurement induced circuit schematically depicted in figure 1. In the circuit, the N propagating field modes in an unknown state

$$\begin{eqnarray}&&| \psi \rangle ={\int }_{{{\mathbb{R}}}^{N}}\psi ({x}_{1},\,\cdots ,\,{x}_{N})| {x}_{1},\,\cdots ,\,{x}_{N}\rangle {{\rm{d}}{x}}_{1}\cdots {{\rm{d}}{x}}_{N}\end{eqnarray} \tag{ 1 }$$

interact with N auxiliary modes prepared in a specific resource state

$$\begin{eqnarray}&&| { \mathcal A }\rangle ={\int }_{{{\mathbb{R}}}^{N}}{{\rm{e}}}^{-{{\rm{i}}{V}}_{x}({y}_{1},\cdots ,{y}_{N})}| {y}_{1},\,\cdots ,\,{y}_{N}\rangle {{\rm{d}}{y}}_{1}\cdots {{\rm{d}}{y}}_{N}.\end{eqnarray} \tag{ 2 }$$

These interactions are local, coupling only two modes at a time by the non-demolition quadrature sum gates (QSG) [24] characterized by interaction Hamiltonians ${\hat{H}}_{j}={\hat{x}}_{T,j}{\hat{p}}_{A,j}$, where the subscripts T and A denote the target and the ancilla, respectively. The effective strength of the coupling does not affect the principal performance of the gate and we will therefore consider it to be equal to one. After the interactions which transform the joint state of all modes into:

$$\begin{eqnarray}&&{\int }_{{{\mathbb{R}}}^{2N}}\psi ({x}_{1},\,\cdots ,\,{x}_{N}){{\rm{e}}}^{-{{\rm{i}}{V}}_{x}({y}_{1},\cdots ,{y}_{N})}| {x}_{1},\,\cdots ,\,{x}_{N}\rangle | {y}_{1}-{x}_{1},\,\cdots ,\,{y}_{N}-{x}_{N}\rangle {{\rm{d}}{x}}_{1}\cdots {{\rm{d}}{x}}_{N}{{\rm{d}}{y}}_{1}\cdots {{\rm{d}}{y}}_{N}\end{eqnarray} \tag{ 3 }$$

the ${\hat{x}}_{j}$ quadratures of the auxiliary modes are measured, yielding a set of classical values q_j. For each set of such values, the joint state of the target modes can be expressed as:

$$\begin{eqnarray}&&{{\rm{e}}}^{-{\rm{i}}{F}({\hat{x}}_{1},\cdots ,{\hat{x}}_{N};{q}_{1},\cdots ,{q}_{N})}{{\rm{e}}}^{-{{\rm{i}}{V}}_{x}({\hat{x}}_{1},\cdots ,{\hat{x}}_{N})}| \psi \rangle .\end{eqnarray} \tag{ 4 }$$

Up to the unitary operation ${{\rm{e}}}^{-{\rm{i}}{F}({\hat{x}}_{1},\cdots ,{\hat{x}}_{N};{q}_{1},\cdots ,{q}_{N})}$ depending on the values obtained by the measurement, this is the ideal transformed state we would expect to obtain after application of the multi-mode coupling. The measurement dependent operation is represented by effective Hamiltonian

$$\begin{eqnarray}&&\begin{array}{l}F({\hat{x}}_{1},\,\cdots ,\,{\hat{x}}_{N};{q}_{1},\,\cdots ,\,{q}_{N})={V}_{x}({\hat{x}}_{1}+{q}_{1},\,\cdots ,\,{\hat{x}}_{N}+{q}_{N})-{V}_{x}({\hat{x}}_{1},\,\cdots ,\,{\hat{x}}_{N})\\ \quad =\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{k}_{1},\,\cdots ,\,{k}_{N}=0}{\displaystyle \sum {k}_{j}\gt 0}}^{\infty }\displaystyle \frac{{q}_{1}^{{k}_{1}}\cdots {q}_{N}^{{k}_{N}}}{{k}_{1}!\cdots {k}_{N}!}\displaystyle \frac{{\partial }^{{k}_{1}}}{\partial {\hat{x}}_{1}^{{k}_{1}}}\cdots \displaystyle \frac{{\partial }^{{k}_{N}}}{\partial {\hat{x}}_{N}^{{k}_{N}}}{V}_{x}({\hat{x}}_{1},\cdots {\hat{x}}_{N}).\end{array}\end{eqnarray} \tag{ 5 }$$

This undesired term can be compensated in several ways. We can accept only those realizations in which all the measurements yielded values sufficiently close to zero. Such post-selection would result in perfect imprint of the multi-mode nonlinearity of the ancilla on the joint state of the target oscillators, albeit at the cost of exponentially diminishing success rate. Since the low success rate can be an issue in some quantum processing applications, it is worthwhile to find ways to compensate the term deterministically. This is possible for operations whose Hamiltonian can be expressed as a finite polynomial of the individual quadrature operators. This encompasses a broad range of operations which can be used to approximate any single-quadrature coupling for any given set of quantum states with arbitrary precision. For such operations, term (5) can be removed by applying unitary operator

$$\begin{eqnarray}&&{\hat{U}}_{{\rm{cor}}}={{\rm{e}}}^{{\rm{i}}{F}({\hat{x}}_{1},\cdots ,{\hat{x}}_{N};{q}_{1},\cdots ,{q}_{N})}.\end{eqnarray} \tag{ 6 }$$

The required feed-forward is of the same nature as the operation we are attempting to implement, only with lower order. If the initial attempted operation has an effective Hamiltonian that can be expressed as finite polynomial of the quadrature operators

$$\begin{eqnarray}&&{V}_{x}({\hat{x}}_{1},\,\cdots ,\,{\hat{x}}_{N})=\displaystyle \sum _{{k}_{1},\,\cdots ,\,{k}_{N}=0}^{{M}_{1},\,\cdots ,\,{M}_{N}}{c}_{{k}_{1},\cdots ,{k}_{N}}{\hat{x}}_{1}^{{k}_{1}}\cdots {\hat{x}}_{N}^{{k}_{N}},\end{eqnarray} \tag{ 7 }$$

with total maximal power equal to ${M}_{1}+\cdots +{M}_{N}=M$, then the corrective operation has an effective Hamiltonian that is different, but of the same form (7) only with total maximal power M₁ + ⋯ + M_N = M − 1. Due to this similarity, the same kind of circuit can be employed for implementing the corrective operations, possibly merging the required feed-forwards for the future steps, similarly as in the single mode scenario [23].

Zoom In Zoom Out Reset image size

It should be noted that the QSGs can be in experiments replaced by a passive linear couplings—beam splitters. Beam splitters transform the $\hat{x}$ representation kets as $| {x}_{j},{y}_{j}\rangle \to | ({x}_{j}+{y}_{j})/\sqrt{2},({y}_{j}-{x}_{j})/\sqrt{2}\rangle$. After the measurement yields value ${q}_{j}=({y}_{j}-{x}_{j})/\sqrt{2}$, the ket can be found in the form $| {x}_{j}\sqrt{2}+{q}_{j}\rangle$ while the corresponding wave functions remain identical. A form equivalent to (4) can now be obtained by squeezing each mode by factor $\sqrt{2}$ and displacing its ${\hat{x}}_{j}$ quadrature by $-{q}_{j}/\sqrt{2}$.

The crucial component of the gate is the joint state of the auxiliary modes (2). The ideal resource state can be approached only asymptotically, which is a feature these states share with the quadrature squeezed states. The exact role of this resource state becomes apparent when we describe the operation in the Heisenberg picture. For the interaction Hamiltonian $V({\hat{x}}_{1},\,\cdots ,\,{\hat{x}}_{N})$ the quadrature operators of the target modes are supposed to evolve as:

$$\begin{eqnarray}&&{\hat{x}}_{j}\to {\hat{x}}_{j},\quad {\hat{p}}_{j}\to {\hat{p}}_{j}-\displaystyle \frac{\partial }{\partial {\hat{x}}_{j}}V({\hat{x}}_{1},\,\cdots ,\,{\hat{x}}_{N}).\end{eqnarray} \tag{ 8 }$$

In the measurement induced scheme, the non-demolition nature of the QSG ensures the desired preservation of the ${\hat{x}}_{T}$ operators, while the ${\hat{p}}_{T}$ operators evolve into $\hat{p}{{\prime} }_{T,j}={\hat{p}}_{T,j}+{\hat{p}}_{A,j}$. The corrective feed-forward now relies on transforming the ${\hat{p}}_{T}$ quadratures by applying the unitary operator (6):

$$\begin{eqnarray}\begin{array}{rcl}\hat{p}{{\prime\prime} }_{T,j} & = & \hat{p}{{\prime} }_{T,j}+\displaystyle \frac{\partial }{\partial {\hat{x}}_{T,j}}F({\hat{x}}_{T,1},\,\cdots ,\,{\hat{x}}_{T,N};{q}_{1},\,\cdots ,\,{q}_{N})\\ & = & {\hat{p}}_{T,j}-\displaystyle \frac{\partial }{\partial {\hat{x}}_{T,j}}V({\hat{x}}_{T,1},\,\cdots ,\,{\hat{x}}_{T,N})+{\hat{p}}_{A,j}+\displaystyle \frac{\partial }{\partial {\hat{x}}_{A,j}}V({\hat{x}}_{A,1},\,\cdots ,\,{\hat{x}}_{A,N}).\end{array}\end{eqnarray} \tag{ 9 }$$

Here ${q}_{j}={\hat{x}}_{A,j}-{\hat{x}}_{T,j}$ are the measured values and we can see that the quadratures are transformed in the manner we desire but the imperfect ancilla contributes by N noise terms, one for each of the modes. Fortunately these noise terms commute and their variances can be, in principle, simultaneously reduced to zero, which is exactly the case of the ideal state (2). Complete removal of the noise terms requires states with infinite energy and thus can be accomplished only asymptotically. This is the long term goal. However, for the short term experimental tests we need to find states which, under feasible and realistic constraints, minimize the variances of these operators and thus produce a high order nonlinear squeezed and entangled states suitable for implementing the operations with a better-then-Gaussian performance.

Let us now discuss the preparation of the nonlinear resource states. Such states are required to reduce the added fluctuations and enable implementations of the multi-mode gates that are observably superior to using only Gaussian resources. The defining property of such states is that the joint variance of the nonlinear quadratures (9)

$$\begin{eqnarray}&&\displaystyle \sum _{j=1}^{N}\left\langle {\left\{{\rm{\Delta }}[{\hat{p}}_{A,j}+\displaystyle \frac{\partial }{\partial {\hat{x}}_{A,j}}V({\hat{x}}_{A,1},\cdots ,{\hat{x}}_{A,N})]\right\}}^{2}\right\rangle \end{eqnarray} \tag{ 10 }$$

is minimal. Here subscript A marks the operators as belonging to the auxiliary systems but, since we are going to be dealing solely with the resource states from now on, we will be omitting it in the following text. Let us discuss the techniques for minimization on a concrete example of ancillary states for nonlinear coupling with Hamiltonian $\hat{H}=\kappa {\hat{x}}_{1}{\hat{x}}_{2}^{2}$. This is operation of the third order, which is the lowest order the operation can have and still be non-Gaussian. The operation has the potential of being experimentally tested in the foreseeable future and, from a practical point of view, it can be used for implementing nonlinear non-demolition gates such as those for creating the so called grid state [18]. Since the operation is of the third order, the corrective feed-forward operations (5) will be of the second order and therefore Gaussian [25]. The only source of non-Gaussian nonlinear behavior is therefore the resource state which should exhibit squeezing in nonlinear quadratures

$$\begin{eqnarray}&&\hat{p}{{\prime} }_{1}={\hat{p}}_{1}+\kappa {\hat{x}}_{2}^{2},\quad \hat{p}{{\prime} }_{2}={\hat{p}}_{2}+2\kappa {\hat{x}}_{1}{\hat{x}}_{2}.\end{eqnarray} \tag{ 11 }$$

The desired ancillary states are those that under chosen conditions minimize the sum of the variances of the nonlocal quadrature operators

$$\begin{eqnarray}&&\langle {\left({\rm{\Delta }}\hat{p}{{\prime} }_{1}\right)}^{2}\rangle +\langle {\left({\rm{\Delta }}\hat{p}{{\prime} }_{2}\right)}^{2}\rangle =\left\langle \psi \right|\left({\left({\hat{p}}_{1}+\kappa {\hat{x}}_{2}^{2}\right)}^{2}-{d}_{1}^{2}+{\left({\hat{p}}_{2}+2\kappa {\hat{x}}_{1}{\hat{x}}_{2}\right)}^{2}-{d}_{2}^{2}\right)\left|\psi \right\rangle ,\end{eqnarray} \tag{ 12 }$$

where ${d}_{1}=\left\langle \psi \right|({\hat{p}}_{1}+\kappa {\hat{x}}_{2}^{2})\left|\psi \right\rangle$ and ${d}_{2}=\left\langle \psi \right|({\hat{p}}_{2}+2\kappa {\hat{x}}_{1}{\hat{x}}_{2})\left|\psi \right\rangle$. This moment is the balanced sum of the variances of the nonlinear combination of quadratures (11) and is therefore a nonlinear functional of the states. For the ideal state (2) the term is equal to zero but for any physical state it serve as a source of added noise in the circuit. Both the quadratures are nonlinear combination of quadratures of different modes. As a consequence the resource quantum state should exhibit non-Gaussian quantum entanglement manifesting as nonlocal nonlinear squeezing.

The minimization can be separated into two steps. The ancillary states can be expressed as a non-Gaussian core state prepared as a finite superposition of Fock states, which is subsequently adjusted by an optimized two-mode Gaussian operation

$$\begin{eqnarray}&&| \psi \rangle ={\hat{U}}_{G}\displaystyle \sum _{m=0,n=0}^{M,N}{c}_{{mn}}| m,n\rangle .\end{eqnarray} \tag{ 13 }$$

This is a feasible two-mode extension of the core state principle [26, 27]. With the help of the Bloch–Messiah decomposition [28] the Gaussian operation on two oscillators can be decomposed into two passive linear coupling operations with two squeezing operations sandwiched in between, together with two displacement operations at the output. This amounts to 12 free parameters which need to be considered in the optimization. Fortunately, some of these parameters can be removed due to the symmetry of the problem. Displacement of $\hat{p}$ quadratures can be neglected because it does not affect the variances. Similarly, the optimal displacement of the $\hat{x}$ quadratures is zero because the operators (11) are symmetric with respect to $\hat{x}$. The linear coupling operations can be decomposed into rotation operations and beam splitter operations [29]. The nonlinear quadratures (11) behave quite differently with regards to their $\hat{x}$ and $\hat{p}$ operators. As a consequence, the rotation operations which mix the different single oscillator quadratures do not contribute towards reduction of the variance. The final Gaussian operation can be therefore composed as

$$\begin{eqnarray}&&{\hat{U}}_{G}={\hat{U}}_{{BS}}({\theta }_{2}){\hat{S}}_{1}({\lambda }_{1}){\hat{S}}_{2}({\lambda }_{2}){\hat{U}}_{{BS}}({\theta }_{1}),\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&{\hat{U}}_{{BS}}^{\dagger }(\theta )\left(\begin{array}{c}{\hat{x}}_{1}\\ {\hat{x}}_{2}\\ {\hat{p}}_{1}\\ {\hat{p}}_{2}\end{array}\right){\hat{U}}_{{BS}}(\theta )=\left(\begin{array}{cccc}\cos \theta & \sin \theta & 0 & 0\\ -\sin \theta & \cos \theta & 0 & 0\\ 0 & 0 & \cos \theta & \sin \theta \\ 0 & 0 & -\sin \theta & \cos \theta \end{array}\right)\left(\begin{array}{c}{\hat{x}}_{1}\\ {\hat{x}}_{2}\\ {\hat{p}}_{1}\\ {\hat{p}}_{2}\end{array}\right),\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&{\hat{S}}_{i}^{\dagger }({\lambda }_{i})\left(\begin{array}{c}{\hat{x}}_{i}\\ {\hat{p}}_{i}\end{array}\right){\hat{S}}_{i}({\lambda }_{i})=\left(\begin{array}{c}{\hat{x}}_{i}/{\lambda }_{i}\\ {\hat{p}}_{i}{\lambda }_{i}\end{array}\right)\end{eqnarray} \tag{ 16 }$$

represent the composing linear coupling and squeezing operations, respectively. The optimal Gaussian processing therefore has only four relevant parameters, beam splitter parameters $\theta \in (-\pi ,\pi ]$, and squeezing parameters ${\lambda }_{1},{\lambda }_{2}\in (0,\infty )$. Note that even though we neglect arbitrary phase shifts, by admitting the full 2π range of the beam splitter parameters we allow for phase flips. When limited to Gaussian resources the input will be the vacuum state ${\left|\psi \right\rangle }_{12}={\left|0\right\rangle }_{1}{\left|0\right\rangle }_{2}$ and we will denote the minimal variance (12) that can be achieved with such resources by V_G. This variance can be found numerically. By introducing non-Gaussian input state we aim to reduce this variance. We will quantify the reduction by using the relative variance ${R}_{V}=\tfrac{{V}_{{NG}}}{{V}_{G}}$, which is defined as the ratio between the variances of the non-Gaussian and the optimal Gaussian ancillary states.

Let us start by checking how easily can the nonlinearity required for the gate be created or, in other words, what is the easiest setup that can produce suitable approximation of nonlinear state (2). For this we can consider the simplest possible non-Gaussian core states, which are the factorized states of the two oscillators, and reduce the Gaussian processing only to single passive linear coupling. This reduces the difficulty of the numerical optimization as well as the potential experimental implementation. Thus obtained resource states of the form

$$\begin{eqnarray}&&| {\psi }_{M,N}\rangle ={\hat{U}}_{{BS}}(t)\displaystyle \sum _{m=0}^{M}{c}_{m}| m\rangle \otimes \displaystyle \sum _{n=0}^{N}c{{\prime} }_{n}| n\rangle \end{eqnarray} \tag{ 17 }$$

can be prepared in currently available experimental circuits [20, 30] and therefore offer the possibility of experimental verification. The results of the numerical optimization, which is in detail described in appendix A, are summarized in figure 2. The relative variances of the auxiliary states are plotted in relation to the value of the nonlinear parameter κ. We can see that nonlocal nonlinear squeezing which surpasses the optimal Gaussian states with arbitrary squeezing can be constructed even with the limited resources we consider. Already a qubit-like superposition of the two lowest Fock states in one oscillator and vacuum in the other, $| \psi {\rangle }_{\mathrm{1,0}}=(0.8\left|0\right\rangle +i0.58\left|1\right\rangle )\otimes \left|0\right\rangle$ shows this improvement with R_V = 0.94 for the gate strength amplitude κ = 0.46. At the same time, more resources in the form of higher dimension of single or both oscillator states lead to significant improvement. For example, two factorized qubit-like state ${U}_{{BS}}\left(1/\sqrt{2}\right)(0.82\left|0\right\rangle +i0.57\left|1\right\rangle )\otimes (0.82\left|0\right\rangle +i0.57\left|1\right\rangle )$ can provide relative variance of R_G = 0.84 for κ = 0.38. The optimal states for other parameters can be found in appendix B. The improvement following from the higher dimension is quantitative as well as qualitative. For higher dimensions, the relative variances V_R are lower and the non-Gaussian resources surpass the Gaussian benchmark over larger ranges of the nonlinear parameter κ. The reason why the non-Gaussian resources are not universally better than Gaussian is that the figure of merit (12) is composed of moments of quadrature operators. These operators have eigenstates, which can be approached by Gaussian preparation. For example, when κ is approaching zero, the optimal resource state needs very low variances of the $\hat{p}$ quadratures, which are unobtainable by discrete superpositions of low number Fock states, but easily achieved by squeezed vacuum states.

Zoom In Zoom Out Reset image size

This issue can be resolved by considering Gaussian processing (14) as part of the nonlinear state preparation (13). Interestingly, the numerical analysis reveals that the second linear coupling is not even necessary and parameter of the second beam splitter can be se to θ₂ = 0. The squeezing operations are crucial, though. Their main role lies in adjusting the nonlinear strength κ of the produced state. This means that a non-Gaussian core state that is optimal for one value of κ is optimal for all of them. This gives rise to the following invariants for each non-Gaussian core:

$$\begin{eqnarray}&&\displaystyle \frac{\kappa }{{\lambda }_{1}^{2}{\lambda }_{2}},\quad \displaystyle \frac{\kappa }{{\lambda }_{1}{\lambda }_{2}^{2}},\quad \displaystyle \frac{{\lambda }_{1}}{{\lambda }_{2}},\quad \displaystyle \frac{\kappa }{{\lambda }_{1}^{3}},\quad \displaystyle \frac{\kappa }{{\lambda }_{2}^{3}}.\end{eqnarray} \tag{ 18 }$$

The invariants are only necessary conditions the parameters need to satisfy; not all combinations of λ₁ and λ₂ lead to an optimal output state, but those that do obey (18). This property of the Gaussian processing means that all the discussion about the resource states can be reduced to the discussion of the non-Gaussian cores.

For that we compared two categories of states:

$$\begin{eqnarray}\begin{array}{rcl}| {\psi }_{M,N}\rangle & = & {\hat{U}}_{G}\displaystyle \sum _{m=0}^{M}{c}_{m}| m\rangle \otimes \displaystyle \sum _{n=0}^{N}c{{\prime} }_{n}| n\rangle \\ | {\psi }_{M\times N}\rangle & = & {\hat{U}}_{G}\displaystyle \sum _{m=0}^{M}\displaystyle \sum _{n=0}^{N}{c}_{{nm}}| m\rangle | n\rangle .\end{array}\end{eqnarray} \tag{ 19 }$$

The first one, denoted as M, N, consists of factorized states of the two oscillators, while the second one, marked by M × N, refers to inseparable entangled states. The comparison of these classes is summarized in figure 3. We can see that the entangled state is more powerful than a factorized state with the same dimensions, which in turn is more powerful than the completely asymmetric case with one of the states being the vacuum state. Interestingly, even though the nonlinear coupling with Hamiltonian $\hat{H}=\kappa {\hat{x}}_{1}{\hat{x}}_{2}^{2}$ is not symmetric, the optimal non-Gaussian resource states $| {\psi }_{N,N}\rangle$ are.

Zoom In Zoom Out Reset image size

The numerical methods for finding the proper non-Gaussian core states need to find the minimum of a function with a large value of variables. For example, for a state $| {\psi }_{N,M}\rangle$ it is 2(M + N) for the non-Gaussian core state alone. The optimal state $| {\psi }_{N\times M}\rangle$ requires even more, 2(NM + N + M). As we look for stronger nonlinear squeezing and the dimension of the problem increases, the numerical optimization starts becoming challenging. One possibility is observing symmetries in the states. The observed symmetry of states $| {\psi }_{N,N}\rangle$, for example, reduces the number of free parameters into half. The observation that the coefficients of both separable and entangled non-Gaussian cores are, up to a π/2 phase shift, always real does that as well. Another observation we feel worthy of mentioning is that the constituents of the found factorized non-Gaussian core states have very large overlap with the non-Gaussian resource states $| {\gamma }_{N}\rangle$ for the single mode cubic gate with Hamiltonian $\hat{H}=\kappa {\hat{x}}^{3}$ [22]. The overlap also translates into the nonlinear properties. This can be seen in figure 3 where the variance of cubic resource states is very close to the variance of resource states optimized with the two-mode interaction in mind. Not only the gap vanishes when the dimension N approaches infinity, but for the state $| {\gamma }_{N}\rangle \otimes | {\gamma }_{N}\rangle$ the total added variance approaches zero (see appendix C for details). For the case of only a single cubic ancilla, $| {\gamma }_{N}\rangle \otimes | 0\rangle$, the minimal achievable variance was found to be R_V = 0.70, but in this case the optimality of this scenario is not guaranteed. However, the pair of perfect single mode cubic ancillas is optimal and together with Gaussian coupling they are sufficient for obtaining the perfect two-mode resource state. The nonlocal aspect of entangled core states only expedites the process. This is very promising outcome, because the cubic phase can be produced by other platforms, like quantum optomechanics [31, 32], and then, they become resource for two-mode cubic couplings as well. In classical levitating optomechanics, such the cubic interaction has been already experimentally achieved [33, 34].

We have presented a general method for measurement induced implementation of multi-mode non-Gaussian interactions that couple only single quadratures of propagating modes of microwaves or light. While the non-Gaussian couplings can be decomposed into a sequence of elementary gates [15], the dedicated circuits require significantly reduced amount of non-Gaussian resources, similarly as in [23]. For example, the basic non-Gaussian coupling with Hamiltonian $\hat{H}={\hat{x}}_{1}{\hat{x}}_{2}^{2}$ studied in this paper needs to be decomposed into a sequence requiring three ideal cubic gates [35]. The direct implementation of the couplings presented here can therefore reduce the overall resources required. It can be also used for immediate proof-of-principle experimental tests which, we hope, will in turn stimulate the theoretical investigation of quantum nonlinear interactions.

The proposed method is based on multi-mode extension of the original GKP scheme [18]. It relies on a non-demolition non-Gaussian measurement [22] of each of the individual modes followed by a suitable joint nonlinear feed-forward operation. The method can be also translated to a loop architecture of nonlinear gates on temporal propagating modes [36]. The crucial element of the implementation is the resource state. For N-mode operations, the resource state is an entangled state simultaneously squeezed in N nonlinear nonlocal quadrature operators. The ideal state is a joint eigenstate of these operators but, similarly to single-quadrature eigenstates, it can be approached only asymptotically.

We have investigated techniques for finding the required state for the example of two-mode cubic gate with Hamiltonian $\hat{H}=\kappa {\hat{x}}_{1}{\hat{x}}_{2}^{2}$. First of all, we have found that for some values of the nonlinear parameter κ the entangled resource outperforming purely Gaussian approximation of the coupling can be obtained simply by creating a suitable superposition of Fock states $| 0\rangle$ and $| 1\rangle$ and splitting them on a beam splitter. The squeezing of the nonlocal quadratures can be then significantly improved either by using higher Fock states, or by interfering a pair of such states on the beam splitter. Moreover, singe-mode cubic phase states can be merged by the Gaussian two-mode couplings to prepare required entangled resource states for the two-mode-cubic phase gate. This opens possibility for other platforms, like nonlinear quantum optomechanics, to deliver the required resources for nonlinear optical circuits. All these techniques are feasible with the available experimental techniques and the nonlinearity and its properties can therefore be readily experimentally tested.

Universal ancillary states can be obtained by preparing a specific single or multi-mode finite superposition of Fock states and altering it by suitable Gaussian coupling which includes squeezing operations. The squeezing operations can be used for tuning the value of parameter κ so the underlying Fock superposition state, the non-Gaussian core of the resource state, is universal for all κ and depends only on the selected dimension. Its exact form needs to be found by numerical methods and the nonlinear squeezing significantly depends on the complexity: entangled superposition of Fock states leads to a better resource state than factorized one. This is an intriguing example of genuine non-Gaussian entangled state significantly outperforming state which is also non-Gaussian but with entanglement created only by Gaussian operations. This opens up an interesting avenue of research: while the conditional techniques for preparation of single mode superpositions of Fock states are well known [20, 30, 37], preparation of multi-mode state has been devised only for specific classes of states [38, 39] and the general case remains and open problem even in the conditional regime. The deterministic preparation of such states then arises as the ultimate challenging goal of future investigations.

We were supported by Project No. GB14-36681G of the Czech Science Foundation. PM acknowledges LTC17086 of INTER-EXCELLENCE program of Czech Ministry of Education. RF has received national funding from the MEYS of the Czech Republic (Project No. 8C18003) under Grant Agreement No. 731473 within the QUANTERA ERA-NET cofund in quantum technologies implemented within the European Unions Horizon 2020 Programme (Project TheBlinQC).

The value that we are trying to minimize is:

$$\begin{eqnarray}&&\left\langle \psi \right|A(a,{a}^{\dagger })\left|\psi \right\rangle -{\left(\left\langle \psi \right|B(a,{a}^{\dagger })\left|\psi \right\rangle \right)}^{2}\quad {\rm{s}}.{\rm{t}}.\left\langle \psi | \psi \right\rangle =1.\end{eqnarray} \tag{ A1 }$$

Here A and B are known operators which are represented as matrices in this particular vector space that we are limited to. The vector that minimizes this value can be expanded in an orthogonal basis

$$\begin{eqnarray}&&\left|\psi \right\rangle ={c}_{0}\left|0\right\rangle +{c}_{1}\left|1\right\rangle +...{c}_{f}\left|f\right\rangle .\end{eqnarray} \tag{ A2 }$$

Plugging it into equation (A1) one can derive the following quartic polynomial that has to be minimized over a complex unit sphere:

$$\begin{eqnarray}&&\displaystyle \sum _{{nm}}{A}_{{nm}}{c}_{n}^{* }{c}_{m}-{\left(\displaystyle \sum _{{nm}}{B}_{{nm}}{c}_{n}^{* }{c}_{m}\right)}^{2}.\end{eqnarray} \tag{ A3 }$$

For calculating the variances following matrix elements has to be calculated:

$$\begin{eqnarray}&&\left\langle \psi \right|A(a,{a}^{\dagger })\left|\psi \right\rangle ,\quad \left\langle \psi \right|B(a,{a}^{\dagger })\left|\psi \right\rangle .\end{eqnarray} \tag{ A4 }$$

We have calculated the matrix elements in the following way:

$$\begin{eqnarray}&&{A}_{{mn}}=\left\langle m\right|A(a,{a}^{\dagger })\left|n\right\rangle ,\quad {B}_{{mn}}=\left\langle m\right|B(a,{a}^{\dagger })\left|n\right\rangle .\end{eqnarray} \tag{ A5 }$$

Here m and n are the Fock states, i.e. ${a}^{\dagger }a\left|n\right\rangle =n\left|n\right\rangle$ However note that since we are dealing with complex numbers the number of parameters is two times the dimension of the vector space. (Due to normalization condition it is one less than two times the dimension to be more precise.) The most convenient way to deal with such system is to expand the coefficients into real and imaginary parts and minimize the equation in real vector space.

Eventually this method relies on minimization of a quartic polynomial over a unit sphere that contains large number of variables which is computationally demanding. This is actually a NP-hard problem [40]. We have used a particular built in minimization algorithm from Mathematica software. More specifically we have employed a quasi newton type BFGS algorithm through using tens of thousands of initial points for the search.

In this appendix we present some optimal input states and the corresponding Gaussian transformation. Note how the two-mode states are symmetrical under the exchange of inputs. Also the coefficients are all real up to a ${{\rm{e}}}^{\tfrac{{\rm{i}}\pi }{2}({n}_{1}+{n}_{2})}$ operation. For setups that contain only passive elements:

• When limited to single photon in the first input for κ = 0.46
  $$\begin{eqnarray}&&{U}_{{BS}}(0.86)(0.8\left|0\right\rangle +i0.58\left|1\right\rangle )\otimes \left|0\right\rangle .\end{eqnarray} \tag{ B1 }$$
• When limited to two photons in the first input For κ = 0.38
  $$\begin{eqnarray}&&{U}_{{BS}}(0.87)(0.47\left|0\right\rangle +i0.78\left|1\right\rangle -0.4\left|2\right\rangle )\otimes \left|0\right\rangle .\end{eqnarray} \tag{ B2 }$$
• When limited to one photons in both inputs for κ = 0.38
  $$\begin{eqnarray}&&{U}_{{BS}}\left(\displaystyle \frac{\pi }{4}\right)(0.82\left|0\right\rangle +i0.57\left|1\right\rangle )\otimes (0.82\left|0\right\rangle +i0.57\left|1\right\rangle ).\end{eqnarray} \tag{ B3 }$$
• When limited to two photons in both inputs for κ = 0.29
  $$\begin{eqnarray}&&{U}_{{BS}}\left(\displaystyle \frac{\pi }{4}\right)(0.51\left|0\right\rangle +i0.76\left|1\right\rangle -0.38\left|2\right\rangle )\otimes (0.51\left|0\right\rangle +i0.76\left|1\right\rangle -0.38\left|2\right\rangle ).\end{eqnarray} \tag{ B4 }$$

For setups that employ the full Gaussian processing with passive coupling and single mode squeezing operations:

• When limited to single photon in the first input
  $$\begin{eqnarray}&&{U}_{{BS}}(-0.16){S}_{1}(1.04){S}_{2}(1.47){U}_{{BS}}(1.07)(0.8\left|0\right\rangle +i0.59\left|1\right\rangle )\otimes \left|0\right\rangle .\end{eqnarray} \tag{ B5 }$$
• When limited to two photons in the first input
  $$\begin{eqnarray}&&{U}_{{BS}}(-0.22){S}_{1}(1.04){S}_{2}(1.47){U}_{{BS}}(1.14)(0.47\left|0\right\rangle +i0.77\left|1\right\rangle -0.41\left|2\right\rangle )\otimes \left|0\right\rangle .\end{eqnarray} \tag{ B6 }$$
• When limited to one photons in both inputs
  $$\begin{eqnarray}&&{S}_{1}(1.06){S}_{2}(1.58){U}_{{BS}}\left(\displaystyle \frac{\pi }{4}\right)(0.7\left|00\right\rangle +i0.42\left|01\right\rangle +i0.42\left|10\right\rangle -0.38\left|11\right\rangle ).\end{eqnarray} \tag{ B7 }$$
• When limited to two photons in both inputs
  $$\begin{eqnarray}&&\begin{array}{l}{S}_{1}(1.07){S}_{2}(1.79){U}_{{BS}}\left(\displaystyle \frac{\pi }{4}\right)(0.33\left|00\right\rangle +i0.36\left|01\right\rangle -0.11\left|02\right\rangle \\ \quad +i0.36\left|10\right\rangle -0.59\left|11\right\rangle -i0.3\left|12\right\rangle -0.11\left|20\right\rangle -i0.3\left|21\right\rangle +0.21\left|22\right\rangle ).\end{array}\end{eqnarray} \tag{ B8 }$$

In this section we show that using prefect cubic phase states with Gaussian resources one can approach two-mode cubic gate in principle.

We first generate two single mode cubic phase gate ancillae and then couple it using a fifty–fifty beam splitter:

$$\begin{eqnarray}&&{U}_{{BS}}\int {{\rm{e}}}^{{{\rm{i}}{t}{X}}_{1}^{3}}\left|{x}_{1}\right\rangle {{\rm{d}}{x}}_{1}\int {{\rm{e}}}^{{{\rm{i}}{t}{X}}_{2}^{3}}\left|{x}_{2}\right\rangle {{\rm{d}}{x}}_{2}=\int {{\rm{e}}}^{{\rm{i}}{t}^{\prime} {\left({X}_{1}-{X}_{2}\right)}^{3}}{{\rm{e}}}^{{\rm{i}}{t}^{\prime} {\left({X}_{1}+{X}_{2}\right)}^{3}}\left|{x}_{1}\right\rangle \left|{x}_{2}\right\rangle {{\rm{d}}{x}}_{1}{{\rm{d}}{x}}_{2}.\end{eqnarray} \tag{ C1 }$$

After that we can implement squeezing operations in modes 1 and 2:

$$\begin{eqnarray}&&\int {{\rm{e}}}^{{\rm{i}}{t}^{\prime} (2{X}_{1}^{3}{\lambda }_{1}^{3}+6{X}_{1}{X}_{2}^{2}{\lambda }_{1}{\lambda }_{2}^{2})}\left|{x}_{1}\right\rangle \left|{x}_{2}\right\rangle {{\rm{d}}{x}}_{1}{{\rm{d}}{x}}_{2}\end{eqnarray} \tag{ C2 }$$

If we choose the parameters so that λ₁λ₂² = 1 and ${\lambda }_{1}\to 0$, only the two-mode term remains relevant and the ideal resource state is obtained.