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

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 various physical platforms to accomplish different tasks relevant to quantum information science [1,2].The major challenge is to develop a coherent quantum interface across different physical platforms.Most prior investigations have focused on quantum state transfer (i.e., SWAP operations) between different physical 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 in hybrid systems, it is thus desirable to have the capability of processing quantum information directly over multiple physical platforms.
As essential to various information processing tasks, in this paper, we aim at the general construction of multimode Gaussian operations over relevant modes in 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 [16].By contrast, however, hybrid bosonic systems lack ondemand beam-splitter operations, and thus the existing protocols are often inapplicable.Indeed, bosonic modes hosted on disparate physical platforms (e.g.microwaveoptical 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 singlemode squeezing [19]; 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 [31].
A similar setup to the one described above was studied in Ref. [32], where the authors introduced the notion of using interference for hybrid bosonic 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 quan-tum information between the two involved bosonic modes without any additional pre-or post-processing.However, while powerful, their results focus only on the SWAP operation (i.e.transduction) rather than more general Gaussian operations.More importantly, the methods presented in [32] are specific to two-mode systems and cannot be directly generalized.This hinders the applicability of their scheme since, in practice, hybrid devices involve many interacting modes.
In this paper, we resolve the aforementioned challenges and develop a novel 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 quadratureby-quadrature.With slight modification, this 'modedecoupling' 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 can be identified using an efficient graph algorithm.

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 {q k , pk }.Consider an N -mode system described by a vector of quadrature operators x := (q 1 , p1 , . . ., qN , pN ) T .In the Heisenberg picture, any Gaussian unitary transfor-mation ÛS mapping xk → ÛS xk Û † S can be equivalently characterized by a 2N × 2N real symplectic matrix S mapping x → Sx [19] (see Methods for details).Without loss of generality, we can work entirely in terms of these symplectic scattering matrices: given two Gaussian unitaries ÛR and ÛS , we have ÛRS = ÛR Û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 whole 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 [33].Now, for GUC's 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 (i) [1 ≤ i ≤ 4 ] be the symplectic matrices associated with 4 local Gaussian operations (here "local" means each of the L (i) consists of individual single-mode operations).If we carefully engineer the L (i) 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 We consider a sequence of 4 multi-mode symplectic matrices S (j) (blue) interspersed with local operations k (yellow) are singlemode 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 âin k and âout k 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 (1) and T (2) .(b) As an example, we demonstrate the decoupling of mode 1.In the first layer, T (1) 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 (1) of the specified form, where * denotes an arbitrary matrix element or sub-block.Note T (2) has the same structure.(c) The second recursive layer of the sequence involves sandwiching another set of local operations between T (1) and T (2) (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.
interference-based sequence: The specific choice of local operations L (i) 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 (1) , S âin (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 S = S S ⊕ I N − , where I N − is the (N − )-mode identity matrix and S (red) is another N -mode operation.(b) Recursive layers of the universal sequence.We organize the 16 multi-mode operations (either S or S) into groups of four, and apply the decoupling protocol to each group to yield four matrices R (i) with the first mode decoupled.We then apply the decoupling protocol again using R (i) in order to additionally decouple the second mode.(c) The recursive step results in a matrix of the form L1 ⊕ L2 ⊕ S * , where L1, L2 (brown) are single-mode operations, and S * (grey) is an arbitrary (N − )-mode operation on the remaining modes.Finally, we apply local "recovery" operations We are then left with the desired effective operation of the form S ⊕ S * .S (2) , S (3) , and S (4) , we can construct an interferencebased sequence by interspersing these matrices with carefully chosen local Gaussian operations L (1) , L (2) , and L (3) such that the resulting Gaussian operation S (4) L (3) S 3 L (2) S (2) L (1) S (1) 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 (1) = S (2) L (1) S (1) and T (2) = S (4) L (3) S (3) 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 (2) L (2) T (1) .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 (i) .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.
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 S = S S ⊕ I N −2 , where I N −2 is the (N − 2)-mode identity matrix and S is another N -mode operation.This decomposition is purely mathematical.As shown in Fig. 2(b), we then apply the mode-decoupling protocol to the interference-type sequence of 15 copies S and one copy of S (interspersed with 15 local operations L (i) ) in order to yield a two-mode decoupled intermediate operation as shown in Fig. 2(c).After that, one additional set of local Gaussian "recovery" operations L (r) k is applied in order to cancel the resulting single-mode operations L 1 and L 2 from the mode-decoupling.This is done by choosing L The resulting sequence in Fig. 2(d) 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.

Dealing with the edge cases
The Eqs. (14)(15)(16)(17), as we will discuss in more detail later, for the construction of the local Gaussian operations may not be applicable when for a GUC S and a certain mode m, both S k,2m−1 and S k,2m (or S 2m−1,k and S 2m,k ) are zero, where k represents the quadrature we hope to engineer.Although this is a rare situation in the practical settings, it calls for a more careful look into the applicability of the general scheme.First of all, as 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, we usually can resolve the issue by randomization and saturation, i.e. by simply replacing the given S with L (2) S L (1) where L (1) and L (2) are randomly sampled local Gaussian operations, so that the vanishing elements will disappear.
However, there exist certain exceptional situations, which we refer to as the edge cases, where the vanishing elements cannot be removed in a similar way as above using any randomized interference-based sequence of the form (1) consisting of multiple copies of S interspersed by randomly-sampled local Gaussian operations L (i) , considering that the only resources we are granted are the free access to local Gaussian operations and multiple uses of the given Gaussian interaction.The simplest example of an edge case is the permutation of modes up to additional single-mode Gaussian operations, which can be realized physically as circulators.Obviously, no Gaussian operation but permutation of modes can be obtained via any interferencebased sequences, since the local Gaussian operations, which themselves can be interpreted as trivial permutations, cannot introduce more complicated couplings between modes.Therefore, general-purpose Gaussian operations involving more than one-modes are not obtainable 3 , 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 the swap the two color sets (right panel).This demonstrates the invariant grouping behavior of the modes of same colors.Here, the L (i) matrices are randomly-sampled local Gaussian operations.
with such a special Gaussian interaction as the given input.It turns out this example also extends to a general categorization of symplectic matrices, only that a general multi-mode Gaussian operation permutes aggregation of modes, which we call the color sets, instead of individual modes.That is to say, randomized interference-based sequences as introduced in the beginning of this paragraph can only yield permutations of color sets, all of which belong to the same permutation group generated by the one determined by a single S. As a result, with a sufficiently long randomized interference-based sequence, we can replace a given edge case S with the new symplectic matrix a collection of mutually non-interfering {S c } acting on the color sets labelled by the color c with γ(G S ) being the total number of colors [34], where each S c is a fully-randomized symplectic matrix with no vanishing elements.Thus, one can immediately see that any Gaussian operations can be constructed on any set of bosonic modes inside the same color set and correspondingly cannot be constructed on modes belonging to different color sets.
We now introduce a systematic way of identifying the color sets as well as the permutation specified by the given S by representing S as a graph.Here, a graph means a collection of vertices and vertex-connecting arrows determined by the following rules: We assign a vertex (e.g., i) to each bosonic mode (e.g., âi ) involved in the operation; and the vertex i is linked to a vertex j through an arrow with the arrowhead pointing to j, if and only if the sub-block the 2 × 2 zero matrix [35].Once the graph has been set up, the color sets can be figured out by recursively contracting the vertices according to a definite set of rules listed in the Methods.A non-trivial example is shown in Fig. 3, which yields a two-vertex and two-arrow graph as shown in the right panel of Fig. 3(c).Each of the two colored vertices in this simplified graph, as a result of vertex contraction, is a color containing two bosonic modes, with the arrows connecting them indicating how they are permuted by the given S. Fig. 3(d) shows interferencebased sequences generate a permutation group as alluded to in the previous paragraph.The mechanism we introduced for identifying the color sets can be efficiently calculated on a classical computer, since the vertex contraction steps can be efficiently executed and therefore, the overall complexity depends polynomially on the number of modes involved.More notably, the physical overhead of randomizing and saturating the given Gaussian interaction S to suit it to the general scheme, i.e. the number of copies of S needed for the the randomized interference-based sequence needed, will not exceed 2N 2 , where N is the total number of the involved bosonic modes [34].

DISCUSSION
In comparison with other existing ideas aimed at addressing similar Gaussian control problems in hybrid systems, our scheme has several manifest upsides.In addition to benefits we alluded to earlier, Ref. [32] will not lead to a systematic categorization of symplectic matrices (see [34] for more detail), which is an indispensable piece for the puzzle, as we developed here, due to its rigid focus on the two-bosonic-mode situation.On the other hand, the method in Ref. [31] is also compatible with multi-mode situations and can yield a variety of Gaussian operations by tuning locally accessible parameters without changing the given GUC.However, this method requires some overly demanding resources such as infinite squeezing and perfect homodyne measurement, which are absent from the list of the requirements of the scheme presented here; in the meantime, whether and how this method can yield arbitrary desired Gaussian operation are still unknown.
Beyond Gaussian operations, our scheme bears similarity to the quantum approximate optimization algorithm (QAOA) [36] -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,37], 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 novel 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 more detail in the Supplementary Material [34].

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 [19].Nevertheless, for the sake of completeness, we review the salient details below.
In this work, we study Gaussian unitary operations of the form Û = exp(−i Ĥt) where Ĥ is bilinear in the field operators.In the Heisenberg picture, such operations will realize the transformation â → Û † â Û .This is equivalently characterized by the scattering matrix transforming the quadrature operators x → Sx.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: (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 1 , . . ., 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 R(θ) = exp −iθâ † â is represented in the quadrature basis by the symplectic matrix For single-mode squeezing Ẑ(r) = exp r(â 2 − â †2 )/2 , the associated symplectic matrix representation is given in the quadrature basis by Any 2 × 2 (local) symplectic matrix can be decomposed into two phase rotations and single-mode squeezing [19].
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 (4) L (3) S (3) L (2) S (2) L (1) S (1) 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 (1) , S (2) , S (3) , and S (4) 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 (1) , L (2) , and L (3) can be used to construct the following symplectic matrices of the form 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 (1) and T (2) via the above "sandwiching" of GUC's and local operations, a resultant symplectic matrix R can be obtained from the whole sequence 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 [14].Specifically, for an arbitrary 2N × 2N symplectic matrix S, we can denote its rows by S = (u 1 , v 1 , . . ., u N , v N ) T and its columns by S = (x 1 , y 1 , . . ., 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: u T i Ωu j = v T i Ωv j = 0 and u T i Ωv j = δ ij , where i, j ∈ {1, 2, . . ., N }.Additionally, the columns satisfy x T i Ωx j = y T i Ωy j = 0 and x T i Ωy j = δ 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 (1) , is understood as finding the suitable local operation L (1)  such that is of the form shown in Eq. ( 6).Note: we have expressed S (1) = (x 1 , y 1 , . . ., x N , y N ) in terms of its column vectors, and S (2) = (u 1 , v 1 , . . ., u N , v N ) T in terms of its row vectors.Now, suppose there exists an L (1) that transforms the first column of S (1) such that L (1) x 1 = Ωu 1 .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 1 is naturally orthogonal to each of the other columns except y 1 , and thus L (1) x 1 will be orthogonal to each of the modified columns except for L (1) y 1 .Furthermore, since L (1) x 1 = Ωu 1 by construction, it will also be orthogonal to each of the rows of S (2) except for v 1 .Thus T (1) = S (2) L (1) S (1) will be of the expected form.By an almost identical calculation, we can show that a suitable choice of L (3) will result in T (2) = S (4) L (3) S (3) having the desired form of Eq. ( 6); we simply use the row and column symplectic bases of S (4) and S (3) respectively.
With the matrices T (1) and T (2) 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 (2) .As before, let us start by expressing T (2) = (α 1 , β 1 , . . ., α N , β N ) T in terms of its row vectors, and T (1) = (χ 1 , γ 1 , . . ., χ N , γ N ) in terms of its column vectors.
We now construct a "sandwich" R = T (2) L (2) T (1) by choosing L (2) in such a way that it transforms the first two columns χ 1 , γ 1 of T (1) respectively to L (2) which can be satisfied simultaneously simply by letting L (2) 1 be the symplectic form ω. Since L (2) χ 1 and L (2) γ 1 are linearly independent, the two-dimensional plane spanned by the pair of vectors Ωα 1 , Ωβ 1 is identical to that spanned by the vectors L (2) χ 1 , L (2) γ 1 .Consequently, this plane is orthogonal to every other row vector Ωα j , Ωβ j , and column vector L (2) χ j , L (2) γ 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 (1) , L (2) , and L (3) can in fact be constructed.By definition, each of these operations is the direct sum of N individual single-mode Gaussian operations: i.e.
k are 2 × 2 symplectic matrices.Thus, in order to transform one 2Ndimensional vector (e.g.x 1 ) to another (e.g.Ωu 1 ) using N single-mode operations (e.g.L (1) ), it suffices to show that we can transform any generic 2-dimensional vector to another using one single-mode operation (e.g.L 1 ).As demonstrated in Fig. 4, this can be satisfied generically -that is, for any pair of vectors (q i , p i ) = (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 (1) , L (2) , and L (3) is always guaranteed unless unless for the initial quadrature u and the final quadrature v, the there exists a mode i, such u 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 now 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 randomlysampled symplectic matrices, we are effectively performing a sequence S (16) L (15) • • • S (2) L (1) S (1) 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 (we which are not concerned with).
We can engineer this interaction using an interferencebased 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: 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 S (as before, interspersed with 4 − 1 local operations L (1) , L (2) , . . ., L (4 −1) ).This results in a Gaussian operation that isolates the first modes, i.e. a symplectic matrix of the form 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 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: 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

FIG. 2 .
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 L (j) k and L (r) k

FIG. 3 .
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 grey 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 âin k (blue arrows) are first injected into a four-port circulator.Then from the outputs âout k (red arrows), those of mode 2 and 3 are routed into a 50:50 beam-splitter.The resulting outputs, also labelled âout 2 and âout3 , 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 the swap the two color sets (right panel).This demonstrates the invariant grouping behavior of the modes of same colors.Here, the L (i) matrices are randomly-sampled local Gaussian operations.

FIG. 4 .
FIG.4.Geometric argument for the existence of local operations.Given any non-zero single-mode quadrature vector (q1, p1), it is possible to to another non-zero quadrature vector (q2, p2) using only phase-space rotations and finite squeezing.This local transformation is constructed (in the Heisenberg picture) via Li = R(−θ)Z(r)R(ϕ) with θ, ϕ ∈ [0, 2π) and r a non-negative real number.