Deterministic multi-mode nonlinear coupling for quantum circuits

We present a general technique for deterministic implementation of a multi-mode nonlinear coupling between several propagating microwave or optical modes in quantum circuits. The measurement induced technique combines specifically prepared resource states together with feasible feed-forward operations. We explore several ways of generating the suitable resource states and discuss their difference on an illustrative example of cubic coupling between two modes. We also show that the required entangled states with requisite nonlinear properties can be already generated in the present day experiments.


I. INTRODUCTION
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][4][5] or metrology [6][7][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 multimode 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 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 in 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 feedforward. The biggest obstacle is the availability of the requisite ancillary states. They can be prepared approximately [19][20][21][22] but these approximations are imperfect 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 nondemolition interactions with N auxiliary modes prepared in a specific nonlinear state, homodyne detection of the ancillary modes, and nonlinear feedforward 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 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.

II. MULTI-MODE NON-GAUSSIAN COUPLING
A general unitary coupling between N modes can be characterized by interaction Hamiltonian expressed as a functional V xp (x 1 ,p 1 , · · · ,x N ,p N ), wherex j andp j are quadrature operators of the individual modes with commutation relations [x j ,p k ] = iδ jk .The required Hermitian nature of the Hamiltonian demands the underlying real function V xp (.) to be symmetrical with respect to exchanging operatorsx j andp 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 (x 1 , · · · ,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 |ψ into e −iVx(x1,··· ,x N ) . It can be implemented by the measurement induced circuit schematically depicted in Fig. 1. In the circuit, the N propagating field modes in an unknown state interact with N auxiliary modes prepared in a specific resource state These interactions are local, coupling only two modes at a time by the non-demolition quadrature sum gates (QSG) [24] characterized by interaction HamiltoniansĤ j =x T,jpA,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: thex 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: Up to the unitary operation e −iF (x1,··· ,x N ;q1,··· ,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 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. When the multi-mode interaction Hamiltonian The input modes interact with the auxiliary modes A1, · · · , AN through a sequence of quadrature sum gates (QSG). The auxiliary modes are then measured by homodyne detectors (HD), yielding values q1, · · · , qN , which are then used for the corrective feed-forwardF (q1, · · · , qN ), which is a shorthand notation for (6), in order to produce the output.
can be expressed as a finite polynomial of the individual quadrature operators, the term can be removed by applying unitary operatorÛ 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, with total maximal power equal to M 1 + · · · + 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 1 + · · · + 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]. It should be noted that for the quadrature sum gates can be in experiments replaced by a passive linear couplings -beam splitters. Beam splitters transform thex representation kets as |x j , y j → |(x j + y j )/ √ 2, (y j − x j )/ √ 2 . After the measurement yields value q j = (y j − x j )/ √ 2, the ket can be found in the form |x j √ 2 + q j while the corresponding wave functions remain identical. A form equivalent to (4) can now be obtained by squeezing each mode by factor √ 2 and displacing itsx j quadrature by −q j / √ 2. The crucial component of the gate is the joint state of the auxillary 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 (x 1 , · · · ,x N ) the quadrature operators of the target modes are supposed to evolve as: In the measurement induced scheme, the non-demolition nature of the QSG ensures the desired preservation of thex T operators, while thep T operators evolve intop T,j =p T,j +p A,j . The corrective feed-forward now relies on transforming thep T quadratures by applying the unitary operator (6): Here q j =x A,j −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.

III. RESOURCE STATES FOR THE TWO-MODE CUBIC GATE
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 the classical approaches. The defining property of such states is that the joint variance of the nonlinear quadratures (9 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 H = κx 1x 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 a photon number resolving measurement [25]. Since the operation is of the third order, the corrective feed-forward operations (5) will be of the second order and therefore Gaussian. The only source of non-Gaussian nonlinear behavior is therefore the resource state which should exhibit squeezing in nonlinear quadratureŝ The desired ancillary states are those that, under chosen conditions minimize the sum of the variances of the the nonlocal quadrature operators where d 1 = ψ| (p 1 + κx 2 2 ) |ψ and d 2 = ψ| (p 2 + 2κx 1x2 ) |ψ . 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, 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 ofp quadratures can be neglected because it does not affect the variances. Similarly, the optimal displacement of thex quadratures is zero because the operators (11) are symmetric with respect tox. 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 theirx andp 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Û whereÛ † represent the composing linear coupling and squeezing operations, respectively. The optimal Gaussian processing therefore has only four relevant parameters, beam splitter parameters θ ∈ (−π, π], and squeezing parameters λ 1 , λ 2 ∈ (0, ∞). 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 |ψ 12 = |0 1 |0 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 = V N G 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 can be prepared in currently available experimental circuits [30,31] 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 Fig. 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, |ψ 1,0 = (0.8 |0 + i0.58 |1 ) ⊗ |0 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 1/ √ 2 (0.82 |0 + i0.57 |1 ) ⊗ (0.82 |0 + i0.57 |1 ) 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 κ.
While the simplified resource states (17) are interesting for proof-of-principle tests, we are ultimately interested in the best performance available. We can approach it by considering the full Gaussian processing (14) as part of the 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 θ 2 = 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 two primary invariants for each non-Gaussian core: which consequently lead to secondary invariant quantities: λ1 λ2 , κ . The invariants are only necessary conditions the parameters need to satisfy; not all combinations of λ 1 and λ 2 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.
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 Fig. 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Ĥ = κx 1x 2 2 is not symmetric, the optimal non-Gaussian resource states |ψ N,N are.
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 |ψ N,M it is 2(M + N ) for the non-Gaussian core state alone. The optimal state |ψ N ×M requires even more, 2(N M + 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 |ψ N,N , 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 |γ N for the single mode cubic gate with HamiltonianĤ = κx 3 [22]. The overlap also translates into the nonlinear properties. This can be seen in Fig. 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 |γ N ⊗ |γ N the total added variance approaches zero (see appendix C for details). For the case of only a single cubic ancilla, |γ N ⊗ |0 , the minimal achievable variance was found to be V R = 0.7022, 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 it is 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 [32,33], 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 [34,35].

IV. SUMMARY AND DISCUSSION
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. The full ramifications and potential applications of these nonlinear quantum interactions are still not completely understood, mainly because the difficulty accompanying both the theoretical analysis the potential experimental implementation. We believe that our method can be used for first proof-of-principle experimental tests which will in turn stimulate the theoretical investigation of this sizeable yet unexplored area.
The proposed method is based on multi-mode extension of the original GKP scheme [18]. It relies on a nondemolition 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Ĥ = κx 1x 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 and |1 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 in the Gaussian two-mode couplings to prepare required entangled resource states for 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,31,37], preparation of multi-mode state has been devised only for specific classes of states [38] 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.
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 e • When limited to two photons in both inputs. In this section we show that using prefect cubic phase states with Gaussian resources one can approach two mode cubic gate in principle.