Universal interference-based construction of Gaussian operations in hybrid quantum systems

Abstract

Beam-splitter operations are an indispensable resource for processing quantum information encoded in bosonic modes. In hybrid quantum systems, however, it can be challenging to implement reliable beam-splitters between two distinct modes due to various experimental imperfections. Without beam-splitters, realizing arbitrary Gaussian operations between bosonic modes can become highly non-trivial or even infeasible. In this work, we develop interference-based protocols for engineering Gaussian operations in multi-mode hybrid bosonic systems without requiring beam-splitters. Specifically, for a given generic multi-mode Gaussian unitary coupler, we demonstrate a universal scheme for constructing Gaussian operations on a desired subset of the modes, requiring only multiple uses of the given coupler interleaved with single-mode Gaussian unitaries. Our results provide efficient construction of operations crucial to quantum information science, and are derived from fundamental physical properties of bosonic systems. The proposed scheme is thus widely applicable to existing platforms and couplers, with the exception of certain edge cases. We introduce a systematic approach to identify and treat these edge cases by utilizing an intrinsically invariant structure associated with our interference-based construction.

Introduction

Hybrid quantum systems can exploit the complementary advantages of different physical platforms to accomplish various tasks relevant to quantum information science^1,2. The major challenge is to develop a coherent quantum interface across the different platforms. Most prior investigations have focused on quantum state transfer (i.e. SWAP operations) between different hardware platforms, while information processing tasks are still separated within the individual sub-systems^{3,4,5,6,7,8,9,10,11,12}. To develop a more powerful and hardware-efficient quantum interface for hybrid systems, it is thus desirable to have the capability of processing quantum information directly over multiple physical platforms.

As an essential building block for various information processing tasks, we focus here on the general construction of multi-mode Gaussian operations over the relevant modes of a hybrid quantum system. Existing protocols and theorems for the construction^13,14 and decomposition^14,15 of Gaussian operations crucially require access to exact and reliable beam-splitter operations — a demanding facility usually only afforded by pure-bred optical systems (with no frequency mismatches) in experimental settings¹⁶. By contrast, however, hybrid bosonic systems lack on-demand beam-splitter operations, and thus the existing protocols are often inapplicable. Indeed, bosonic modes hosted on disparate physical platforms (e.g., microwave-optical or microwave-mechanical) have vastly different resonant frequencies, and so cannot easily be coupled without the use of nonlinear mixing processes. This in turn inevitably leads to the system modes coupling to unwanted auxiliary modes — for instance, to stray sidebands caused by the linearization of the intrinsically nonlinear optomechanical or electro-optical interactions^7,8,17,18 — and thus prohibits us from cleanly realizing beam-splitter interactions and more general Gaussian operations.

Therefore, in order to construct arbitrary Gaussian operations in hybrid systems, we would ideally like to have an efficient hardware-aware protocol that functions without needing exact on-demand beam-splitter operations between selected modes. To this end, we consider a theoretical setting in which we instead have access to only one given multi-mode Gaussian Unitary Coupler (GUC) involving all participating modes (i.e. system and auxiliary modes alike), as well as free access to single-mode Gaussian unitary operations. The multi-mode GUC is an irreducible resource for making disparate modes interact, and we allow it to be replicated (i.e. used multiple times). It is, however, immutable due to the difficulty of changing the intrinsic underlying experimental parameters. Meanwhile, the single-mode Gaussian unitaries can be implemented using only phase-shifting and single-mode squeezing¹⁹; this requirement is justified thanks to inspiring recent development of squeezing techniques in hybrid bosonic systems^{16,20,21,22,23,24,25,26,27,28,29,30}. We emphasize that our setup here forgoes the need for infinite squeezing and/or perfect homodyne detection as is required by certain existing hybrid bosonic control schemes³¹.

A similar setup to the one described above was studied in ref. ³², where the authors introduce the notion of using interference for hybrid bosonic mode control. They demonstrate how a sequence of multiple identical copies of a two-mode Gaussian unitary coupler, interspersed with single-mode operations, can completely swap quantum information between the two involved bosonic modes without any additional pre- or post-processing. However, although this result is powerful, it focuses only on quantum transduction (i.e. SWAP gates) rather than on more general Gaussian operations. More importantly, the methods presented in ref. ³² are specific to two-mode systems and cannot be directly generalized. This hinders the applicability of the protocol since, in practice, hybrid devices typically involve many interacting modes.

In this paper, we resolve the aforementioned challenges by developing an interference-based framework for realizing general multi-mode Gaussian operations in hybrid quantum systems. We consider the theoretical setting described above in which we have free access to single-mode Gaussian control operations, but only one immutable multi-mode GUC that can be replicated (i.e., applied repeatedly). Within this setting, we demonstrate how interference may be used to construct arbitrary multi-mode Gaussian operations between a selected subset of the system modes, while simultaneously isolating this interaction from any unwanted auxiliary modes. The basis for our results is an observation that the coupling between a pair of quadratures can be removed via interference — i.e., implementing the given GUC twice, interspersed with single-mode operations. By constructing an inductive multi-pass sequence of this form, we can then successively remove all unwanted coupling terms quadrature-by-quadrature. With slight modification, this ‘mode-decoupling’ scheme can then be applied recursively with finitely-many identical copies of the GUC in order to realize our central goal — a universal framework for constructing arbitrary Gaussian operations in hybrid systems.

The results presented here are a direct consequence of the fundamental commutation relations for bosonic systems, and thus our work is generically applicable to most current hybrid bosonic quantum information platforms. There are, however, certain types of edge cases: i.e. certain initial types of GUC for which additional processing is required to apply our protocol. It turns out that the investigation of these edge cases not only leads to a comprehensive understanding of the power and limitations of interference-based protocols, but also reveals an obscure invariant structure that is intrinsic to Gaussian unitary operations, and that can be identified using an efficient graph algorithm.

Results

Overview of general scheme

Gaussian unitary operations involving linearly-coupled bosonic modes are completely determined by their action on the expectation values of the quadrature operators $\{{\hat{q}}_{k},{\hat{p}}_{k}\}$. Consider an N-mode system described by a vector of quadrature operators $\hat{{{{\bf{x}}}}}:= {({\hat{q}}_{1},{\hat{p}}_{1},\ldots ,{\hat{q}}_{N},{\hat{p}}_{N})}^{T}$. In the Heisenberg picture, any Gaussian unitary transformation ${\hat{U}}_{{{{\boldsymbol{S}}}}}$ mapping ${\hat{{{{\bf{x}}}}}}_{k}\to {\hat{U}}_{{{{\boldsymbol{S}}}}}{\hat{{{{\bf{x}}}}}}_{k}{\hat{U}}_{{{{\boldsymbol{S}}}}}^{{\dagger} }$ can be equivalently characterized by a 2N × 2N real symplectic matrix S mapping $\hat{{{{\bf{x}}}}}\to {{{\boldsymbol{S}}}}\hat{{{{\bf{x}}}}}$¹⁹ (see Methods for details). Without loss of generality, we can work entirely in terms of these symplectic scattering matrices: given two Gaussian unitaries ${\hat{U}}_{{{{\boldsymbol{R}}}}}$ and ${\hat{U}}_{{{{\boldsymbol{S}}}}}$, we have ${\hat{U}}_{{{{\boldsymbol{R}}}}{{{\boldsymbol{S}}}}}={\hat{U}}_{{{{\boldsymbol{R}}}}}{\hat{U}}_{{{{\boldsymbol{S}}}}}$. Therefore, instead of using infinite-dimensional unitary operators, it suffices to track the matrix product of the (2N × 2N)-dimensional symplectic matrices to capture the effects of the entire process.

In this work, we will consider Gaussian interactions between the modes of a hybrid quantum system. As stated in the Introduction, our starting assumption is that we have access to only one given multi-mode Gaussian unitary operation, characterized by its scattering matrix S. We refer to this as a Gaussian Unitary Coupler (GUC), as it couples all modes of the system. In our setting, the GUC is fixed by the system parameters and is thus immutable. For instance, we could consider the GUC to be the ‘bare’ unitary process induced by the multi-mode system Hamiltonian³³. Now, for GUCs typically available in hybrid systems, the associated symplectic matrices usually lack clear structure, and cannot be utilized to implement useful Gaussian controls. For example, in a hybrid system consisting of mutually interacting optical, mechanical, and microwave modes, we cannot obtain a simple beam-splitting operation between any two of the modes due to stray coupling to sidebands and other auxiliary modes^7,8,17,18. However, as noted above, clean and on-demand Gaussian controls are useful technological tools for many quantum information applications. It is therefore an intriguing question as to whether we can convert a complicated Gaussian unitary process (i.e. available GUC) into some desired Gaussian operation, by making use of only the mathematical properties of symplectic matrices.

In pursuit of an answer to this question, we discover the following solution which is also the main message of this work. Suppose we can repeatedly apply the same GUC, and have access to arbitrary single-mode Gaussian unitary controls on every bosonic mode involved in the process; then, a sequential combination of these two types of Gaussian operations can be constructed to produce any other desired Gaussian operation on any subset of the involved modes, using a finite number of steps. For example, with identical copies of a generic GUC as elementary operations, our scheme could be used to generate a beam-splitter or a SWAP operation on any two of the involved modes. This result is most succinctly described using the mathematical language of symplectic matrices. Let S represent the symplectic matrix associated with a given generic GUC, and let L⁽ⁱ⁾ [1 ≤ i ≤ 4^ℓ] be the symplectic matrices associated with 4^ℓ local Gaussian operations (here “local” means each of the L⁽ⁱ⁾ consists of individual single-mode operations). If we carefully engineer the L⁽ⁱ⁾ according to our knowledge of the matrix S, then we claim that any ℓ-mode Gaussian operation on ℓ of the involved bosonic modes can be obtained generically as the symplectic matrix S^eff of the following interference-based sequence:

$${{{{\boldsymbol{S}}}}}^{{{{\rm{eff}}}}}={{{{\boldsymbol{L}}}}}^{({4}^{\ell })}{{{\boldsymbol{S}}}}{{{{\boldsymbol{L}}}}}^{({4}^{\ell }-1)}{{{\boldsymbol{S}}}}\cdots {{{\boldsymbol{S}}}}{{{{\boldsymbol{L}}}}}^{(1)}{{{\boldsymbol{S}}}}.$$

(1)

The specific choice of local operations L⁽ⁱ⁾ needed to realize this result is discussed in the Methods section.

Mode-decoupling protocol

The general-purpose protocol shown in Eq. (1) for constructing Gaussian operations is itself the logical derivative of another intermediate universal protocol. Given four arbitrary but generic symplectic matrices S⁽¹⁾, S⁽²⁾, S⁽³⁾, and S⁽⁴⁾, we can construct an interference-based sequence by interspersing these matrices with carefully chosen local Gaussian operations L⁽¹⁾, L⁽²⁾, and L⁽³⁾ such that the resulting Gaussian operation S⁽⁴⁾L⁽³⁾S³L⁽²⁾S⁽²⁾L⁽¹⁾S⁽¹⁾ operates separately on a selected mode and the rest of the system. We refer to this as ‘mode-decoupling’ since the resulting operation induces no coupling between the selected mode and remaining modes. As shown in Fig. 1, the protocol takes two recursive steps: (i) The construction of sub-sequences T⁽¹⁾ = S⁽²⁾L⁽¹⁾S⁽¹⁾ and T⁽²⁾ = S⁽⁴⁾L⁽³⁾S⁽³⁾ formed by ‘sandwiching’ local operations in between the multi-mode S^(j) matrices. These are chosen to decouple a quadrature of the selected from the system; (ii) The concatenation of the sub-sequences to form another ‘sandwich’ R = T⁽²⁾L⁽²⁾T⁽¹⁾. This step decouples the conjugate quadrature to the selected mode, thus isolating the entire selected mode from the remaining system. During this process, we utilize the defining mathematical properties of symplectic matrices (i.e. the canonical commutation relations) in order to construct the local control operations L⁽ⁱ⁾. This is discussed in detail in the Methods section. As an aside, we highlight here that this universal mode-decoupling protocol does not require the matrices S^(j) to be identical. Consequently, our main result in Eq. (1) could also be realized mathematically using 4^ℓ distinct GUC’s. However, in keeping with constraints discussed in the Introduction for hybrid systems, we limit ourselves to only one available GUC S in our main result.

Fig. 1: Universal interference-based mode-decoupling protocol.

a We consider a sequence of 4 multi-mode symplectic matrices S^(j) (blue) interspersed with local operations ${{{{\boldsymbol{L}}}}}^{(i)}={{{\rm{diag}}}}({{{{\boldsymbol{L}}}}}_{1}^{(i)},\ldots ,{{{{\boldsymbol{L}}}}}_{N}^{(i)})$, where the ${{{{\boldsymbol{L}}}}}_{k}^{(i)}$ (yellow) are single-mode Gaussian unitaries calculated based on the S^(j) matrices. Each solid black arrow, from right to left, represents the evolution of a bosonic mode under the whole sequence of Gaussian operations, where we use ${\hat{a}}_{k}^{{{{\rm{in}}}}}$ and ${\hat{a}}_{k}^{{{{\rm{out}}}}}$ to denote the input and output mode operators in the Heisenberg picture. The full sequence has a double-layer structure, containing two sub-sequences T⁽¹⁾ and T⁽²⁾. b As an example, we demonstrate the decoupling of mode 1. In the first layer, T⁽¹⁾ is constructed using carefully chosen local Gaussian operations in order to remove all coupling terms between one quadrature of mode 1 and the remaining modes k ≠ 1. This results in a matrix T⁽¹⁾ of the specified form, where * denotes an arbitrary matrix element or sub-block. Note T⁽²⁾ has the same structure. c The second recursive layer of the sequence involves sandwiching another set of local operations between T⁽¹⁾ and T⁽²⁾ (purple) in order to further remove the remaining correlations between the selected mode (1) and the rest of the system. The resulting Gaussian operation R has the first mode decoupled (i.e. isolated from the remaining N − 1 modes), and is depicted as two disjoint blocks.

The mode-decoupling process takes 4 symplectic matrices S^(j) to decouple a single mode from the system. This can be straightforwardly generalized: given 4^ℓ arbitrary symplectic matrices and 4^ℓ − 1 local Gaussian unitaries, we can isolate ℓ modes from the N mode system. This is done by repeating the single mode-decoupling sequence recursively to 4^ℓ−1 groups of four S^(j) matrices. This inductive approach is feasible as the above result is independent of N. Finally, as an important clarification, we stress here that mode-decoupling should not be taken to mean directly removing certain entanglement in a quantum state. Instead, in this work, we are focused on operations rather than the quantum states, and we suggest the readers to stick to the Heisenberg picture throughout the text.

On-demand construction of Gaussian operations

We now discuss our central protocol and main results. For simplicity, however, we will only demonstrate the construction of 2-mode Gaussian operations here, and save the construction of more general ℓ-mode Gaussian operations for the Methods section. Suppose that we wish to generate a desired target operation S^⊙ on the first two modes of the N-mode system. We can do so using 16 copies of the generic N-mode GUC S arranged in an interference-based sequence. Our strategy requires a fictitious decomposition of one of the copies of S into ${{{\boldsymbol{S}}}}={{{\boldsymbol{S}}}}^{\prime} \left({{{{\boldsymbol{S}}}}}^{\odot }\oplus {{{{\boldsymbol{I}}}}}_{N-2}\right)$, where I_N−2 is the (N − 2)-mode identity matrix and ${{{\boldsymbol{S}}}}^{\prime}$ is another N-mode operation. This decomposition is purely mathematical. As shown in Fig. 2b, we then apply the mode-decoupling protocol to the interference-type sequence of 15 copies of S and one copy of ${{{\boldsymbol{S}}}}^{\prime}$ (interspersed with 15 local operations L⁽ⁱ⁾) in order to yield a two-mode decoupled intermediate operation as shown in Fig. 2c. After that, one additional set of local Gaussian “recovery” operations ${{{{\boldsymbol{L}}}}}_{k}^{(r)}$ is applied in order to cancel the resulting single-mode operations L₁ and L₂ from the mode-decoupling. This is done by choosing ${{{{\boldsymbol{L}}}}}_{1}^{(r)}={({{{{\boldsymbol{L}}}}}_{1})}^{-1}$ and ${{{{\boldsymbol{L}}}}}_{2}^{(r)}={({{{{\boldsymbol{L}}}}}_{2})}^{-1}$. The resulting sequence in Fig. 2d is then left only with the desired target operation S^⊙ acting on the first two modes. This operation is isolated from the remaining N − 2 modes, which evolve separately according to some arbitrary S^*. Since no specific constraint was imposed on the initial choice of S^⊙, we can thus realize any arbitrary target Gaussian operation on the first 2 modes. Using 4^ℓ copies of S and a multi-mode-decoupling sequence, we can also generalize this to realize arbitrary ℓ-mode Gaussian unitaries. As with the mode-decoupling protocol, our result here works for arbitrary initial choice of GUC S, provided that S is generic — terminology that will be made precise shortly. For such S, both decoupling and the above fictitious decomposition can be carried out as guaranteed by the properties of symplectic matrices.

Fig. 2: Universal construction of a general two-mode Gaussian operation.

a The generic interference-based sequence for constructing an ℓ-mode target Gaussian operation S^⊙ (green), shown for ℓ = 2. This consists of 4^ℓ = 16 copies of the given GUC S (blue) interspersed with local operations ${{{{\boldsymbol{L}}}}}_{k}^{(j)}$ and ${{{{\boldsymbol{L}}}}}_{k}^{(r)}$ (yellow). The first copy of S (i.e. on the input side) in the sequence is fictitiously decomposed into a product of two consecutive symplectic matrices ${{{\boldsymbol{S}}}}={{{\boldsymbol{S}}}}^{\prime} \left({{{{\boldsymbol{S}}}}}^{\odot }\oplus {{{{\boldsymbol{I}}}}}_{N-\ell }\right)$, where I_N−ℓ is the (N − ℓ)-mode identity matrix and ${{{\boldsymbol{S}}}}^{\prime}$ (red) is another N-mode operation. b Recursive layers of the universal sequence. We organize the 16 multi-mode operations (either ${{{\boldsymbol{S}}}}^{\prime}$ or S) into groups of four, and apply the decoupling protocol to each group to yield four matrices R⁽ⁱ⁾ with the first mode decoupled. We then apply the decoupling protocol again using R⁽ⁱ⁾ in order to additionally decouple the second mode. c The recursive step results in a matrix of the form L₁ ⊕ L₂ ⊕ S^*, where L₁, L₂ (brown) are single-mode operations, and S^* (gray) is an arbitrary (N − ℓ)-mode operation on the remaining modes. Finally, we apply local “recovery” operations ${{{{\boldsymbol{L}}}}}_{1}^{(r)}={({{{{\boldsymbol{L}}}}}_{1})}^{-1}$ and ${{{{\boldsymbol{L}}}}}_{2}^{(r)}={({{{{\boldsymbol{L}}}}}_{2})}^{-1}$ in order to cancel out the L_k matrices. d We are then left with the desired effective operation of the form S^⊙ ⊕ S^*.

Dealing with the edge cases

As we discuss in more detail later, Eqs. (14)–(17) provide an explicit construction of the local Gaussian operations for the mode-decoupling protocol. However, given a GUC S and a certain mode m, these equations may become undefined (i.e. not applicable) if both S_k,2m−1 and S_k,2m (or both S_2m−1,k and S_2m,k) are zero for any quadrature k that we hope to engineer. Although this is a rare situation in practical settings, it calls for a more careful look into the applicability of the general scheme. First of all, given the constraints are imposed by the several formulae for calculating the local Gaussian operations, any randomly sampled multi-mode symplectic matrix, which almost always contains no vanishing elements, should be a valid input. If not, this issue may be resolved via randomization and saturation, i.e., simply replacing the given GUC S with ${\widetilde{{{{\boldsymbol{L}}}}}}^{(2)}{{{\boldsymbol{S}}}}{\widetilde{{{{\boldsymbol{L}}}}}}^{(1)}$ where ${\widetilde{{{{\boldsymbol{L}}}}}}^{(1)}$ and ${\widetilde{{{{\boldsymbol{L}}}}}}^{(2)}$ are randomly sampled local Gaussian operations, such that the vanishing elements disappear.

However, there do exist certain exceptional situations, which we refer to as the edge cases, where the vanishing elements cannot be removed in a similar manner as above using randomized interference-based sequences of the form ${\widetilde{{{{\boldsymbol{L}}}}}}^{(k)}{{{\boldsymbol{S}}}}{\widetilde{{{{\boldsymbol{L}}}}}}^{(k-1)}{{{\boldsymbol{S}}}}\cdots {{{\boldsymbol{S}}}}{\widetilde{{{{\boldsymbol{L}}}}}}^{(2)}{{{\boldsymbol{S}}}}{\widetilde{{{{\boldsymbol{L}}}}}}^{(1)}$ — that is, using multiple copies of S interspersed with randomly-sampled local Gaussian operations ${\widetilde{{{{\boldsymbol{L}}}}}}^{(i)}$. The simplest example of such an edge case is the permutation of modes up to additional single-mode Gaussian operations, which can be realized physically using circulators. Clearly no Gaussian operation besides a permutation can be obtained via any interference-based sequences generated from such an edge case; this is because the local Gaussian operations, which themselves may be interpreted as trivial permutations here, cannot introduce more complicated couplings between modes. Therefore, general-purpose Gaussian operations involving more than one mode are not obtainable if the given GUC is a permutation.

It turns out that the permutation structure of a symplectic matrix S plays a crucial role even when S is not just a simple permutation of individual modes. As we will show later, symplectic matrices can always be thought of as permutations acting on groups of modes. In order to formalize this idea, and demonstrate such a grouping in an unambiguous manner, it will prove useful to introduce the notion of a graph representation of a symplectic matrix. Specifically, we define the graph G_S associated with the symplectic matrix S via the following construction: (1) for each bosonic mode ${\hat{a}}_{i}$ in the system, we assign to it a vertex i in the graph; (2) any pair of vertices i, j are linked by an arrow i → j if and only if the sub-block

$$\left(\begin{array}{ll}{S}_{2j-1,2i-1}&{S}_{2j-1,2i}\\ {S}_{2j,2i-1}&{S}_{2j,2i}\end{array}\right)\,\ne \,{{{{\boldsymbol{0}}}}}_{2\times 2},$$

i.e., it is not equal to the 2 × 2 zero matrix. Clearly, a (2N × 2N) symplectic matrix will result in a graph G_S that contains N vertices. Furthermore, we note that such graphs are directional, meaning the presence of an arrow i → j does not imply j → i. We present an example of this graph construction process in Fig. 3, where the symplectic matrix of the form shown in Fig. 3a yields a two-vertex and two-arrow graph as shown in the left panel of Fig. 3c.

Fig. 3: Example edge case and resulting graph contraction algorithm to identify color sets.

a A 16 × 16 edge case symplectic matrix S acting on four modes. The white 2 × 2 sub-blocks are zero rank, while the gray sub-blocks have non-zero rank. The zero blocks prevent the straightforward application of our interference protocols. b Schematic of a possible physical device realizing the matrix S. The four input modes ${\hat{a}}_{k}^{{{{\rm{in}}}}}$ (blue arrows) are first injected into a four-port circulator. Then from the outputs ${\hat{a}}_{k}^{{{{\rm{out}}}}}$ (red arrows), those of mode 2 and 3 are routed into a 50:50 beam-splitter. The resulting outputs, labeled ${\hat{a}}_{2}^{{{{\rm{out}}}}}$ and ${\hat{a}}_{3}^{{{{\rm{out}}}}}$, are the final outputs of the whole scattering process for mode 2 and 3. c Identifying color sets through vertex contraction. We can easily construct the graph G_S corresponding to S. To perform contraction, we choose a vertex ν (e.g., vertex 1 in the left panel) and merge all of its immediately successors {μ_i} (highlighted vertices 2 and 3 on the left), while removing any redundant edges. Successors are vertices directly connected to ν (via edges ν → μ_i highlighted in black). We repeat until all possible contractions are exhausted, resulting in a final graph where each vertex represents a unique color set (e.g., yellow or green in the right panel). d Interference-based sequences permute the color sets. It is easy to check that each panel will yield the same color sets through the vertex contraction process in c. Since a single use of S simply swaps the two color sets (left panel), a double use of S will lead to a trivial permutation removing any interference across the two color sets (middle panel). Then, a triple use of S should once again swap the two color sets (right panel). This demonstrates the invariant grouping behavior of the modes of same colors. Here, the ${\widetilde{{{{\boldsymbol{L}}}}}}^{(i)}$ matrices are randomly-sampled local Gaussian operations.

Once the graph has been set up for a given input S, we can infer the permutation structure by further grouping the vertices together via a set of graph contraction rules listed in the Methods. We refer to this vertex contraction process as “coloring,” and to the resulting groups as color sets. Following this procedure, we find that each color set is uniquely identified by a vertex in the contracted graph; moreover, the arrows connecting these vertices precisely reflect the permutation structure associated with S. As an example, in Fig. 3c, each of the two colored vertices in the contracted graph is a color set containing two bosonic modes (1, 4 and 2, 3 respectively), with the arrows connecting them indicating how they are permuted by the given S.

With the graph-theoretic formulation in hand, we now turn to the main result of this section. Given that a symplectic matrix S induces a permutation on the collection of color sets, any interference-based sequence generated from S will remain a permutation. This implies that interference-based sequences of different lengths form a cyclic group with the generating element being the permutation induced by a single copy of S. As a result, with a sufficiently long randomized interference-based sequence, we can always generate the identity element of this group acting trivially on the color sets. This corresponds to a new symplectic matrix

$${{{\boldsymbol{S}}}}^{\prime} =\mathop{\bigoplus }\limits_{c=0}^{\gamma \left({G}_{{{{\boldsymbol{S}}}}}\right)-1}{{{{\boldsymbol{S}}}}}_{c}^{\prime}.$$

(2)

Here, the color sets are indexed by a set of colors c ∈ {0, 1, …, γ(G_S) − 1}, where γ(G_S) is the total number of colors associated with the graph G_S (see Supplementary Methods); each ${{{{\boldsymbol{S}}}}}_{c}^{\prime}$ represents a fully-randomized symplectic matrix with no vanishing elements. As an example of the idea above, in Fig. 3d, the permutation (denoted as σ) induced by S generates a group {e, σ} with e being the identity element; moreover, e is generated by the sequence ${{{\boldsymbol{S}}}}{\widetilde{{{{\boldsymbol{L}}}}}}^{(1)}{{{\boldsymbol{S}}}}$ which has the form of ${{{\boldsymbol{S}}}}^{\prime}$ above (i.e., acting separately on the two color sets).

In light of the mathematical structure above, we can construct any arbitrary Gaussian operation acting on bosonic modes inside a same color set by applying our general scheme to the corresponding component of ${{{\boldsymbol{S}}}}^{\prime}$; however, universal Gaussian operations cannot be constructed for modes belonging to different color sets, since this is incompatible with the permutation structure.

As a final remark, we note that the graph contraction procedure introduced above for identifying the color sets can be efficiently executed on a classical computer. The overall time complexity of this algorithm depends polynomially on the number of modes involved. More notably, the physical overhead required to randomize and saturate a given Gaussian interaction S in order to make it suitable for the general scheme will not exceed 2N², with N being the total number of involved bosonic modes. For more details on these procedures and on the categorization of edge cases, we refer the reader to the Supplementary Methods.

Discussion

Due to the recursive structure of our scheme, it is difficult to analyze the impact that noise and other experimental imperfections may have in practical implementations of this work. Indeed, the complicated dependence of the form of the local operations on that of the given GUC makes it challenging to track the propagation of physical errors through the nested sequence of operations. While a comprehensive theoretical framework is still elusive, we expect that the impact of simple noise mechanisms (e.g., Gaussian error channels) will be independent of the total number of modes, given the constant overhead of our scheme for a generic GUC. That is, if we focus on constructing two-mode Gaussian operations, the length of the sequence needed to generate any target unitary S^⊙ will be fixed (cf. Fig. 2), requiring at most sixteen copies of S. We also remark that the closer a GUC is to an edge case, the larger the amount of squeezing will be required the construct the local operations, and thus the more errors will be amplified during the whole process; an example of this can be found in the Supplementary Discussion.

Beyond noise considerations, we note that the precision of the GUC may also be an important factor in the practical deployment of our scheme. Specifically, if the elements of the scattering matrix can only be measured to finite precision in an experiment, the resulting imprecise mathematical description of the GUC may cause deviations in our calculation of the local Gaussian operations. Here, however, the recursive pattern of our scheme may help mitigate the issue. The sequence can be decomposed into smaller mode-decoupling sequences consisting of a few local operations, and we can take advantage of the predictable form of these smaller sequences to calibrate and/or optimize the local operations. In doing so, it may be possible to mitigate the effects of an imprecise GUC specification, and achieve the desired target operation.

In comparison to other existing ideas for Gaussian control in hybrid bosonic systems, our current results present several manifest upsides. Unlike ref. ³², we are able to synthesize general multi-mode Gaussian unitaries, rather than being limited to the two-mode situation. Furthermore, our investigation of the edge cases provides a more systematic categorization of symplectic matrices. Besides ref. ³², another relevant result is that presented in ref. ³¹. This scheme is compatible with multi-mode couplers, and can yield a variety of Gaussian operations by tuning only locally accessible parameters. However, it requires overly demanding resources such as infinite squeezing and perfect homodyne measurement; additionally, it is not as yet known whether the protocol in ref. ³¹ can be used to construct any arbitrary Gaussian operation. Such concerns are overcome by our current scheme, which requires only finite squeezing and no homodyne measurements.

Beyond Gaussian operations, our scheme bears similarity to the quantum approximate optimization algorithm (QAOA)³⁴ — in particular, the use of single-mode unitaries modifying a given quantum interaction in order to yield on-demand quantum operations. Due to the correspondence between Gaussian operations and Clifford gates^19,33, our scheme can also be extended to (discrete) qubit-based systems in order to provide a universal method to generate Clifford gates. While QAOA has the advantage of utilizing and generating non-Clifford operations, our scheme offers the benefit of providing deterministic solutions for the local operations in the Clifford case (which would only be approximate with QAOA).

In summary, in this work we demonstrate interference-based protocols for the universal construction of Gaussian operations in a multi-mode hybrid bosonic system. Our results are hardware-aware and highly compatible with a variety of hybrid platforms with complicated interactions between the constituent bosonic modes. We also discovered an invariant structure intrinsic to Gaussian operations which can be useful for characterization and classification of the Gaussian operations. This characteristic structure is discussed in the Supplementary Methods.

Methods

In this section, we present mathematical details of the key steps of our general protocols for isolating bosonic modes and constructing universal Gaussian operations.

Conventions

The conventions and notation used in this work closely follow the standard definitions for continuous-variable quantum information¹⁹. Nevertheless, for the sake of completeness, we review the salient details below.

We consider multi-mode systems comprised of N coupled bosonic modes, which correspond to N pairs of bosonic field operators ${({\hat{a}}_{1},{\hat{a}}_{1}^{{\dagger} },\ldots ,{\hat{a}}_{N},{\hat{a}}_{N}^{{\dagger} })}^{T}\equiv \hat{{{{\bf{a}}}}}$. Here, $[{\hat{a}}_{j},{\hat{a}}_{k}^{{\dagger} }]={\delta }_{jk}$. We can equivalently describe the system using quadrature operators ${\hat{q}}_{k}\equiv ({\hat{a}}_{k}+{\hat{a}}_{k}^{{\dagger} })/\sqrt{2}$ and ${\hat{p}}_{k}\equiv {{{\rm{i}}}}({\hat{a}}_{k}^{{\dagger} }-{\hat{a}}_{k})/\sqrt{2}$, which satisfy the canonical commutation relations. We also define the quadrature vector $\hat{{{{\bf{x}}}}}\equiv {({\hat{q}}_{1},{\hat{p}}_{1},\ldots ,{\hat{q}}_{N},{\hat{p}}_{N})}^{T}$.

In this work, we study Gaussian unitary operations of the form $\hat{U}=\exp (-{{{\rm{i}}}}\hat{H}t)$ where $\hat{H}$ is bilinear in the field operators. In the Heisenberg picture, such operations will realize the transformation $\hat{{{{\bf{a}}}}}\to {\hat{U}}^{{\dagger} }\hat{{{{\bf{a}}}}}\hat{U}$. This is equivalently characterized by the scattering matrix transforming the quadrature operators $\hat{{{{\bf{x}}}}}\to {{{\boldsymbol{S}}}}\hat{{{{\bf{x}}}}}$. In order to respect the canonical commutation relations, this real 2N × 2N matrix S must be symplectic: SΩS^T = Ω, where the symplectic form Ω is block diagonal:

$${{{\boldsymbol{\Omega }}}}{ = \bigoplus }_{i = 1}^{N}{{{\boldsymbol{\omega }}}}={{{\rm{diag}}}}({{{\boldsymbol{\omega }}}},\ldots ,{{{\boldsymbol{\omega }}}}),\quad {{{\rm{with}}}}\,\,\,{{{\boldsymbol{\omega }}}}=\left(\begin{array}{ll}0&1\\ -1&0\end{array}\right).$$

(3)

We refer to single-mode transformations as “local” since they do not induce coupling between the modes; these operations are represented by 2 × 2 symplectic matrices. We also use the label “local” to denote the direct sum of N single-mode operations, e.g., L = diag(L₁, …, L_N). Note: the matrix Ω can be considered a local operation, as defined: it simply corresponds to a π/2 phase shift (ω) on each mode.

We now demonstrate two examples of local symplectic matrices. First, the transformation corresponding to phase-space rotation (i.e., phase shifting) given by $\hat{{{{\mathcal{R}}}}}(\theta )=\exp [-{{{\rm{i}}}}\theta {\hat{a}}^{{\dagger} }\hat{a}]$ is represented in the quadrature basis by the symplectic matrix

$${{{\boldsymbol{R}}}}(\theta )=\left(\begin{array}{ll}\cos \theta &\sin \theta \\ -\sin \theta &\cos \theta \end{array}\right).$$

(4)

For single-mode squeezing $\hat{{{{\mathcal{Z}}}}}(r)=\exp [r({\hat{a}}^{2}-{\hat{a}}^{{\dagger} 2})/2]$, the associated symplectic matrix representation is given in the quadrature basis by

$${{{\boldsymbol{Z}}}}(r)=\left(\begin{array}{ll}{e}^{-r}&0\\ 0&{e}^{r}\end{array}\right).$$

(5)

Any 2 × 2 (local) symplectic matrix can be decomposed into two phase rotations and single-mode squeezing¹⁹. We will later exploit this fact in order to demonstrate the existence of local operations needed for our protocol.

Decoupling a single bosonic mode

The sequence S⁽⁴⁾L⁽³⁾S⁽³⁾L⁽²⁾S⁽²⁾L⁽¹⁾S⁽¹⁾ used for decoupling relies on several properties that can be derived directly from the definition of symplectic matrices. We start with four multi-mode symplectic matrices S⁽¹⁾, S⁽²⁾, S⁽³⁾, and S⁽⁴⁾ that are generic — i.e., as mentioned above in our discussion of the edge cases, fulfilling the constraints on the feasible form of the given symplectic matrices imposed by the ensuing discussion in this section. One can assume randomly sampled symplectic matrices are generic, since the probability of failure is statistically trivial.

Our claim is that carefully engineered local Gaussian operations L⁽¹⁾, L⁽²⁾, and L⁽³⁾ can be used to construct the following symplectic matrices of the form

$$\begin{array}{ll}{{{{\boldsymbol{T}}}}}^{(k)}&={{{{\boldsymbol{S}}}}}^{(2k)}{{{{\boldsymbol{L}}}}}^{(2k-1)}{{{{\boldsymbol{S}}}}}^{(2k-1)}\\ &=\left(\begin{array}{lllll}0&1&0&\ldots &0\\ -1&{T}_{22}^{(k)}&{T}_{23}^{(k)}&\ldots &{T}_{2,2N}^{(k)}\\ 0&{T}_{32}^{(k)}&{T}_{33}^{(k)}&\ldots &{T}_{3,2N}^{(k)}\\ \vdots &\vdots &\vdots &\ddots &\vdots \\ 0&{T}_{2N,2}^{(k)}&{T}_{2N,3}^{(k)}&\ldots &{T}_{2N,2N}^{(k)}\end{array}\right),\end{array}$$

(6)

with k ∈ {1, 2}. The matrices T^(k) correlate the output Q-quadrature of the first mode to its input P-quadrature only. By constructing T⁽¹⁾ and T⁽²⁾ via the above “sandwiching” of GUCs and local operations, a resultant symplectic matrix R can be obtained from the whole sequence

(7)

which is also depicted graphically in Fig. 1. The resulting R consists of a single-mode Gaussian operation on mode 1 (upper diagonal block), and a multi-mode Gaussian operation on the remaining N − 1 modes (lower diagonal block). Effectively, R decouples mode 1: i.e., it induces no interactions between this mode and the rest of the system.

To demonstrate the mechanism behind Eqs. (6)–(7), we can introduce a helpful geometric interpretation. The following simple fact reflects the definition of symplectic matrices: the rows (or columns) of a symplectic matrix form an orthonormal symplectic basis¹⁴. Specifically, for an arbitrary 2N × 2N symplectic matrix S, we can denote its rows by ${{{\boldsymbol{S}}}}={({{{{\bf{u}}}}}_{1},{{{{\bf{v}}}}}_{1},\ldots ,{{{{\bf{u}}}}}_{N},{{{{\bf{v}}}}}_{N})}^{T}$ and its columns by S = ( x₁, y₁, …, x_N, y_N), where u_k, v_k, x_k and y_k are 2N-dimensional column vectors, with 1 ≤ k ≤ N. Since the matrix S describes a physical unitary process that preserves the canonical commutation relations, it must satisfy the matrix equation SΩS^T = Ω. This results in an explicit set of orthogonality relations between the rows of S: ${{{{\bf{u}}}}}_{i}^{T}{{{\boldsymbol{\Omega }}}}{{{{\bf{u}}}}}_{j}={{{{\bf{v}}}}}_{i}^{T}{{{\boldsymbol{\Omega }}}}{{{{\bf{v}}}}}_{j}=0$ and ${{{{\bf{u}}}}}_{i}^{T}{{{\boldsymbol{\Omega }}}}{{{{\bf{v}}}}}_{j}={\delta }_{ij}$, where i, j ∈ {1, 2, …, N}. Additionally, the columns satisfy $\,{{{{\bf{x}}}}}_{i}^{T}{{{\boldsymbol{\Omega }}}}\,{{{{\bf{x}}}}}_{j}={{{{\bf{y}}}}}_{i}^{T}{{{\boldsymbol{\Omega }}}}{{{{\bf{y}}}}}_{j}=0$ and $\,{{{{\bf{x}}}}}_{i}^{T}{{{\boldsymbol{\Omega }}}}{{{{\bf{y}}}}}_{j}={\delta }_{ij}$. Comparing these two sets of relations to the similar properties of orthogonal matrices, one can then think of symplectic matrices as geometric transformations on the spaces spanned by the row (or column) vectors.

With this in mind, we can interpret the general idea of decoupling an individual mode from the others as building up a certain destructive interference between the quadratures (via the geometric orthogonality relations above). The first step, where we construct T⁽¹⁾, is understood as finding the suitable local operation L⁽¹⁾ such that

$$\begin{array}{ll}{{{{\boldsymbol{T}}}}}^{(1)}&={{{{\boldsymbol{S}}}}}^{(2)}{{{{\boldsymbol{L}}}}}^{(1)}{{{{\boldsymbol{S}}}}}^{(1)}\\&=\left(\begin{array}{lll}-&{{{{\bf{u}}}}}_{1}^{T}&-\\ -&{{{{\bf{v}}}}}_{1}^{T}&-\\ \vdots &\vdots &\vdots \\ -&{{{{\bf{u}}}}}_{N}^{T}&-\\ -&{{{{\bf{v}}}}}_{N}^{T}&-\end{array}\right){{{{\boldsymbol{L}}}}}^{(1)}\left(\begin{array}{lllll}\vert &\vert &\cdots \,&\vert &\vert \\ {{{{\bf{x}}}}}_{1}&{{{{\bf{y}}}}}_{1}&\cdots \,&{{{{\bf{x}}}}}_{N}&{{{{\bf{y}}}}}_{N}\\ \vert &\vert &\cdots \,&\vert &\vert \end{array}\right)\end{array}$$

(8)

is of the form shown in Eq. (6). Note: we have expressed S⁽¹⁾ = ( x₁, y₁, …, x_N, y_N) in terms of its column vectors, and ${{{{\boldsymbol{S}}}}}^{(2)}={({{{{\bf{u}}}}}_{1},{{{{\bf{v}}}}}_{1},\ldots ,{{{{\bf{u}}}}}_{N},{{{{\bf{v}}}}}_{N})}^{T}$ in terms of its row vectors.

Now, suppose there exists an L⁽¹⁾ that transforms the first column of S⁽¹⁾ such that L⁽¹⁾ x₁ = Ωu₁. Then, we claim that Eq. (8) indeed takes the form of Eq. (6) as desired (the existence of such an operation will be discussed later). The reason for this claim is as follows: by the geometric properties above, x₁ is naturally orthogonal to each of the other columns except y₁, and thus L⁽¹⁾ x₁ will be orthogonal to each of the modified columns except for L⁽¹⁾y₁. Furthermore, since L⁽¹⁾ x₁ = Ωu₁ by construction, it will also be orthogonal to each of the rows of S⁽²⁾ except for v₁. Thus T⁽¹⁾ = S⁽²⁾L⁽¹⁾S⁽¹⁾ will be of the expected form. By an almost identical calculation, we can show that a suitable choice of L⁽³⁾ will result in T⁽²⁾ = S⁽⁴⁾L⁽³⁾S⁽³⁾ having the desired form of Eq. (6); we simply use the row and column symplectic bases of S⁽⁴⁾ and S⁽³⁾ respectively.

With the matrices T⁽¹⁾ and T⁽²⁾ constructed, we now proceed to the second step of our protocol to fully isolate the first mode from the others. Simply speaking, we repeat the construction above, only replacing S^(1)/(2) with T^(1)/(2) and slightly modifying the form of the local operation L⁽²⁾. As before, let us start by expressing ${{{{\boldsymbol{T}}}}}^{(2)}={({{{{\boldsymbol{\alpha }}}}}_{1},{{{{\boldsymbol{\beta }}}}}_{1},\ldots ,{{{{\boldsymbol{\alpha }}}}}_{N},{{{{\boldsymbol{\beta }}}}}_{N})}^{T}$ in terms of its row vectors, and T⁽¹⁾ = (χ₁, γ₁, …, χ_N, γ_N) in terms of its column vectors.

We now construct a “sandwich” R = T⁽²⁾L⁽²⁾T⁽¹⁾ by choosing L⁽²⁾ in such a way that it transforms the first two columns χ₁, γ₁ of T⁽¹⁾ respectively to L⁽²⁾χ₁ = Ωα₁ and ${{{{\boldsymbol{L}}}}}^{(2)}{{{{\boldsymbol{\gamma }}}}}_{1}=\left({{{{\boldsymbol{\alpha }}}}}_{1}\cdot {{{{\boldsymbol{\beta }}}}}_{1}-{{{{\boldsymbol{\gamma }}}}}_{1}\cdot {{{{\boldsymbol{\chi }}}}}_{1}\right){{{\boldsymbol{\Omega }}}}{{{{\boldsymbol{\alpha }}}}}_{1}-{{{\boldsymbol{\Omega }}}}{{{{\boldsymbol{\beta }}}}}_{1}$, which can be satisfied simultaneously simply by letting ${{{{\boldsymbol{L}}}}}_{1}^{(2)}$ be the symplectic form ω. Since L⁽²⁾χ₁ and L⁽²⁾γ₁ are linearly independent, the two-dimensional plane spanned by the pair of vectors Ωα₁, Ωβ₁ is identical to that spanned by the vectors L⁽²⁾χ₁, L⁽²⁾γ₁. Consequently, this plane is orthogonal to every other row vector Ωα_j, Ωβ_j, and column vector L⁽²⁾χ_j, L⁽²⁾γ_j for j ≥ 2, as guaranteed by the geometrical orthogonality relations. Putting these together, we find the resulting R indeed takes the form shown in Eq. (7)– thus decoupling the first mode as desired.

It remains to be shown that the appropriate local operations L⁽¹⁾, L⁽²⁾, and L⁽³⁾ can in fact be constructed. By definition, each of these operations is the direct sum of N individual single-mode Gaussian operations: i.e., ${{{{\boldsymbol{L}}}}}^{(i)}={{{\rm{diag}}}}({{{{\boldsymbol{L}}}}}_{1}^{(i)},\ldots ,{{{{\boldsymbol{L}}}}}_{N}^{(i)})$, where the ${{{{\boldsymbol{L}}}}}_{k}^{(i)}$ are 2 × 2 symplectic matrices. Thus, in order to transform one 2N-dimensional vector (e.g., x₁) to another (e.g., Ωu₁) using N single-mode operations (e.g., L⁽¹⁾), it suffices to show that we can transform any generic 2-dimensional vector to another using one single-mode operation (e.g., ${{{{\boldsymbol{L}}}}}_{1}^{(1)}$). As demonstrated in Fig. 4, this can be satisfied generically – that is, for any pair of vectors $(q_i,p_i)\ne (0,0)$. The required single-mode operation is realized using a sequence of three elementary operations: (i) rotation to the Q–axis, (ii) dilation, and (iii) rotation to the final direction. In the language of quantum optics, rotation and dilation correspond to phase-shifting and finite squeezing, respectively. Thus, the existence of L⁽¹⁾, L⁽²⁾, and L⁽³⁾ is always guaranteed, unless for the initial quadrature u and the final quadrature v, there exists a mode i such that $\left({u}_{2i-1}^{2}+{u}_{2i}^{2}\right)\left({v}_{2i-1}^{2}+{v}_{2i}^{2}\right)=0$ and ${u}_{2i-1}^{2}+{u}_{2i}^{2}+{v}_{2i-1}^{2}+{v}_{2i}^{2}\,\ne \, 0$.

Fig. 4: Geometric argument for the existence of local operations.

Given any non-zero single-mode quadrature vector (q₁, p₁), it is possible to transform to another non-zero quadrature vector (q₂, p₂) using only phase-space rotations and finite squeezing. This local transformation is constructed (in the Heisenberg picture) via L_i = R(−θ)Z(r)R(φ) with θ, φ ∈ [0, 2π) and r a non-negative real number.

The entire decoupling procedure above can be easily generalized. For example, we could apply the same protocol to 16 randomly-sampled symplectic matrices in order to isolate the first two modes from the rest of the system. In particular, if we have four R-type matrices of the form given in Eq. (7), we can apply the decoupling protocol on the (N − 1)-mode sub-matrices, while leaving the first mode intact up to local operations. This will result in a new symplectic matrix with the first and second modes decoupled. Since each of the R-type matrices are themselves constructed using 4 randomly-sampled symplectic matrices, we are effectively performing a sequence S⁽¹⁶⁾L⁽¹⁵⁾ ⋯ S⁽²⁾L⁽¹⁾S⁽¹⁾ to isolate the 2 modes. At this point, we can proceed inductively in order to sequentially decouple any ℓ modes from the N-mode system. We do so by using 4^ℓ multi-mode symplectic matrices interspersed with 4^ℓ − 1 carefully engineered local Gaussian operations.

Note: although we demonstrated here how to decouple the first mode, this choice was for convenience only. Different sets of local operations will allow us to decouple any mode from the system (by transforming the appropriate column vectors into the appropriate row vectors); equivalently, we are free to re-label the modes arbitrarily. Finally, we highlight this decoupling scheme can be applied inductively because our results are independent of the total number of modes N.

Structure of the general protocol

Let us now discuss how we arrive at our main result. Suppose we wish to construct a specific ℓ-mode target Gaussian operation S^⊙. Without loss of generality, we may assume this acts on the first ℓ modes of the N mode system, so that the desired Gaussian operation is of the block diagonal form S^⊙ ⊕ S^*. Here S^⊙ is the 2ℓ × 2ℓ target symplectic matrix, and S^* is a 2(N − ℓ) × 2(N − ℓ) symplectic matrix, representing some arbitrary operation on the remaining N − ℓ modes (which we are not concerned with).

We can engineer this interaction using an interference-based sequence consisting of 4^ℓ copies of a given GUC S. We start by fictitiously decomposing the first copy of S in the sequence into the product of two matrices:

$${{{\boldsymbol{S}}}}={{{\boldsymbol{S}}}}^{\prime} \left(\begin{array}{ll}{{{{\boldsymbol{S}}}}}^{\odot }&0\\ 0&{{{{\boldsymbol{I}}}}}_{2(N-\ell )}\end{array}\right),$$

(9)

where I_2(N−ℓ) represents a 2(N − ℓ) × 2(N − ℓ) identity matrix. We then apply our decoupling protocol on 4^ℓ − 1 copies of S and one copy of ${{{\boldsymbol{S}}}}^{\prime}$ (as before, interspersed with 4^ℓ − 1 local operations ${{{{\boldsymbol{L}}}}}^{(1)},{{{{\boldsymbol{L}}}}}^{(2)},\ldots ,{{{{\boldsymbol{L}}}}}^{({4}^{\ell }-1)}$). This results in a Gaussian operation that isolates the first ℓ modes, i.e., a symplectic matrix of the form

$$\left(\begin{array}{lllll}{{{{\boldsymbol{L}}}}}_{1}&0&\cdots \,&0&0\\ 0&{{{{\boldsymbol{L}}}}}_{2}&0&\cdots \,&0\\ \vdots &0&\ddots &\ddots &\vdots \\ 0&\vdots &\ddots &{{{{\boldsymbol{L}}}}}_{\ell }&0\\ 0&0&\cdots \,&0&{{{{\boldsymbol{S}}}}}^{* }\end{array}\right)$$

(10)

with L_k, for 1 ≤ k ≤ ℓ, the 2 × 2 symplectic matrices corresponding to the ℓ single-mode local operations on the isolated ℓ bosonic modes. At last, we just need to apply one final local Gaussian “recovery” operation L^(r) of the form

$${{{{\boldsymbol{L}}}}}^{(r)}=\left(\begin{array}{lllll}{({{{{\boldsymbol{L}}}}}_{1})}^{-1}&0&\cdots \,&0&0\\ 0&{({{{{\boldsymbol{L}}}}}_{2})}^{-1}&0&\cdots \,&0\\ \vdots &0&\ddots &\ddots &\vdots \\ 0&\vdots &\ddots &{({{{{\boldsymbol{L}}}}}_{\ell })}^{-1}&0\\ 0&0&\cdots \,&0&{{{{\boldsymbol{I}}}}}_{2(N-\ell )}\end{array}\right)$$

(11)

to finish the construction of the whole sequence. Putting the above steps all together, the process for constructing a desired ℓ-mode Gaussian operation (isolated from the remaining N − ℓ modes of the N mode system) can be summarized using the following equations:

$$\begin{array}{ll}{{{{\boldsymbol{S}}}}}^{{{{\rm{eff}}}}}&={{{{\boldsymbol{L}}}}}^{(r)}{{{\boldsymbol{S}}}}{{{{\boldsymbol{L}}}}}^{({4}^{\ell }-1)}{{{\boldsymbol{S}}}}\cdots {{{\boldsymbol{S}}}}{{{{\boldsymbol{L}}}}}^{(1)}{{{\boldsymbol{S}}}}\\ &={{{{\boldsymbol{L}}}}}^{(r)}\left[{{{\boldsymbol{S}}}}{{{{\boldsymbol{L}}}}}^{({4}^{\ell }-1)}{{{\boldsymbol{S}}}}\cdots {{{\boldsymbol{S}}}}{{{{\boldsymbol{L}}}}}^{(1)}{{{\boldsymbol{S}}}}^{\prime} \right]\left(\begin{array}{ll}{{{{\boldsymbol{S}}}}}^{\odot }&0\\ 0&{{{{\boldsymbol{I}}}}}_{2(N-\ell )}\end{array}\right)\\ &=\left(\begin{array}{lllll}{({{{{\boldsymbol{L}}}}}_{1})}^{-1}&0&\cdots \,&0&0\\ 0&{({{{{\boldsymbol{L}}}}}_{2})}^{-1}&0&\cdots \,&0\\ \vdots &0&\ddots &\ddots &\vdots \\ 0&\vdots &\ddots &{({{{{\boldsymbol{L}}}}}_{\ell })}^{-1}&0\\ 0&0&\cdots \,&0&{{{{\boldsymbol{I}}}}}_{2(N-\ell )}\end{array}\right)\left[\begin{array}{ccccc}{{{{\boldsymbol{L}}}}}_{1}&0&\cdots \,&0&0\\ 0&{{{{\boldsymbol{L}}}}}_{2}&0&\cdots \,&0\\ \vdots &0&\ddots &\ddots &\vdots \\ 0&\vdots &\ddots &{{{{\boldsymbol{L}}}}}_{\ell }&0\\ 0&0&\cdots \,&0&{{{{\boldsymbol{S}}}}}^{* }\end{array}\right]\left(\begin{array}{ll}{{{{\boldsymbol{S}}}}}^{\odot }&0\\ 0&{{{{\boldsymbol{I}}}}}_{2(N-\ell )}\end{array}\right) =\left(\begin{array}{ll}{{{{\boldsymbol{S}}}}}^{\odot }&0\\ 0&{{{{\boldsymbol{S}}}}}^{* }\end{array}\right)\end{array}$$

(12)

Clearly any arbitrary Gaussian operation on the first ℓ modes can be constructed, since there is no constraint on the form of the target symplectic matrix S^⊙ in the above calculation.

Explicit formulae for the local Gaussian operations

As we have seen, the local Gaussian operations L^(k) for k ≥ 2 are determined by the mode-decoupling protocol. It is notable that the values of each element of these local symplectic matrices can be calculated explicitly using the given multi-mode symplectic matrices S^(k). As a matter of fact, there exist closed-form expressions for the single-mode-decoupling local operations, given four randomly-sampled symplectic matrices. We provide an example of these formulae below, but stress that this is not the unique solution.

Let ${{{{\boldsymbol{S}}}}}^{(k)}=\left({S}_{ij}^{(k)}\right)$ for k ∈ {1, 2, 3, 4} be four 2N × 2N randomly-sampled symplectic matrices. To isolate the first mode from the system, we need to construct the sequence S⁽⁴⁾L⁽³⁾S⁽³⁾L⁽²⁾S⁽²⁾L⁽¹⁾S⁽¹⁾ with L^(k) the local Gaussian operations. We first decompose each of the local operations into N single-mode Gaussian operations matrices:

$${{{{\boldsymbol{L}}}}}^{(k)}=\left(\begin{array}{llll}{{{{\boldsymbol{L}}}}}_{1}^{(k)}&0&\cdots \,&0\\ 0&{{{{\boldsymbol{L}}}}}_{2}^{(k)}&\ddots &\vdots \\ \vdots &\ddots &\ddots &0\\ 0&\cdots \,&0&{{{{\boldsymbol{L}}}}}_{N}^{(k)}\end{array}\right).$$

(13)

where ${{{{\boldsymbol{L}}}}}_{m}^{(k)}$ are 2 × 2 symplectic matrices. Then according to the aforementioned geometric interpretation, to transform an arbitrary random vector u to Ωv with v another arbitrary random vector, we can simply let

$${\left({{{{\boldsymbol{L}}}}}_{m}^{(k)}\right)}_{ij}=\frac{{(-1)}^{j+1}{v}_{\bar{i}}{u}_{\underline{j}}}{{({v}_{2m})}^{2}+{({v}_{2m-1})}^{2}}+\frac{{(-1)}^{i}{v}_{\underline{i}}{u}_{\bar{j}}}{{({u}_{2m})}^{2}+{({u}_{2m-1})}^{2}},$$

(14)

with 1 ≤ m ≤ N and i, j ∈ {1, 2}. Here, if i = 1, then $\bar{i}=2m-1$, $\underline{i}=2m$. Meanwhile if i = 2, then $\bar{i}=2m$, $\underline{i}=2m-1$. Note that when both denominators are zero (and thus both numerators also zero), the formula above should be calculated by taking the limit. Since L⁽¹⁾ is meant to transform the first column of S⁽¹⁾ into the first row of S⁽²⁾ multiplied by the symplectic form Ω, we let

$$\begin{array}{l}{\left({{{{\boldsymbol{L}}}}}_{m}^{(1)}\right)}_{ij} =\frac{{(-1)}^{j}{S}_{1,\bar{i}}^{(2)}{S}_{\underline{j},1}^{(1)}}{{({S}_{1,2m}^{(2)})}^{2}+{({S}_{1,2m-1}^{(2)})}^{2}}+\frac{{(-1)}^{i+1}{S}_{1,\underline{i}}^{(2)}{S}_{\bar{j},1}^{(1)}}{{({S}_{2m,1}^{(1)})}^{2}+{({S}_{2m-1,1}^{(1)})}^{2}}\end{array}$$

(15)

For the same reason, the elements of L⁽³⁾ are thus given by:

$$\begin{array}{l}{\left({{{{\boldsymbol{L}}}}}_{m}^{(3)}\right)}_{ij} =\frac{{(-1)}^{j}{S}_{1,\bar{i}}^{(4)}{S}_{\underline{j},1}^{(3)}}{{({S}_{1,2m}^{(4)})}^{2}+{({S}_{1,2m-1}^{(4)})}^{2}}+\frac{{(-1)}^{i+1}{S}_{1,\underline{i}}^{(4)}{S}_{\bar{j},1}^{(3)}}{{({S}_{2m,1}^{(3)})}^{2}+{({S}_{2m-1,1}^{(3)})}^{2}}\end{array}$$

(16)

With these formulae, we can carry out the matrix multiplications in order to calculate T⁽¹⁾ = S⁽²⁾L⁽¹⁾S⁽¹⁾ and T⁽²⁾ = S⁽⁴⁾L⁽³⁾S⁽³⁾. Therefore, according to the protocol, we only need to set

$${u}_{k}={T}_{k,2}^{(1)},\quad {v}_{k}={T}_{2,k}^{(2)}+\mathop{\sum }\limits_{k=1}^{2N}\left({T}_{k,1}^{(1)}{T}_{k,2}^{(1)}-{T}_{1,k}^{(2)}{T}_{2,k}^{(2)}\right){T}_{1,k}^{(2)},$$

(17)

for 1 ≤ k ≤ 2N, and use Eq. (14) to obtain the remaining local operation L⁽²⁾.

General mechanism of identifying the color sets

The graph-theory inspired language makes it possible for us to come up with the following mechanism to properly color the bosonic modes involved in an arbitrary multi-mode Gaussian interaction:

• 1. Set up the graph (G_S) corresponding to the given Gaussian interaction (S)

• 2. Pick a vertex (e.g., ν) that is the starting point of at least two distinctive arrows.

• 3. Find all the immediate successors of this vertex ν, i.e., those vertices {μ_i} that are linked with ν by arrows (from ν → μ_i). Then, contract all of these successors {μ_i} into a single vertex, while removing any redundant arrows from the graph. (That is to say, if we have two arrows starting and ending with the same pair of vertices, only one of the arrows will be kept.)

• 4. Repeat the above two steps, if possible, until no further contractions can be made (i.e., there is no vertex in the resulting graph with at least two distinctive outgoing arrows).

• 5. The total number of colors $\gamma \left({G}_{{{{\boldsymbol{S}}}}}\right)$ is equal to the number of the vertices in the final resultant graph, after all possible contractions have been performed. Each vertex k of this resultant graph represents a color set (consisting of all the vertices in the original graph G_S contracted to form k). Universal interference-based Gaussian operations can then be constructed between any subset of vertices within the same resultant color set; but not between vertices that end up in different color sets.