Compact Gaussian quantum computation by multi-pixel homodyne detection

We study the possibility of producing and detecting continuous variable cluster states in an optical set-up in an extremely compact fashion. This method is based on a multi-pixel homodyne detection system recently demonstrated experimentally, which includes classical data post-processing. It allows to incorporate the linear optics network, usually employed in standard experiments for the production of cluster states, in the stage of the measurement. After giving an example of cluster state generation by this method, we further study how this procedure can be generalized to perform gaussian quantum computation.


Introduction
Quantum computing is a very promising subject despite the formidable challenges it must overcome [1]. One of those challenges lies in the difficulty of finding a system allowing for a large number of entangled states to be manipulated and for many quantum gates to be implemented successively. In this context, an interesting avenue consists in measurementbased quantum computing [2][3][4] (MBQC) protocols which rely on the availability of a large, multipartite entangled state on which a series of measurements is performed. For each operation that we wish to implement, a specific sequence of measurements on its modes has to be chosen. Successive measurement outcomes are in general used to determine which measurement has to be performed afterwards, and to perform a final correction stage on the output state. Optical systems are promising candidates for the experimental realization of MBQC [5,6] and convincing demonstrations have already been made, both in the discrete [7] and continuous variable (CV) regime [8]. CV systems are especially promising since they allow for the creation and processing of cluster states in a deterministic fashion, contrarily to the discrete case [9].
Most often, the creation of multimode entangled states, such as cluster states in a CV quantum optical setup, is implemented using a series of squeezers followed by a network of beam-splitters and dephasers, which transform the squeezed input modes into entangled output modes [8,[10][11][12]. The configuration of this network varies considerably with the state to be generated, and its complexity grows rapidly with the number of modes, which renders this method poorly scalable.
In a recent experiment at the Australian National University [13], some of us have experimentally demonstrated the possibility of performing exceptionally compact operations on optical modes using spatial (transverse) modes of light. In this experiment, multiple transverse modes were generated and individually squeezed in optical cavities, and then mixed together to produce a multimode squeezed beam. This beam was then measured by means of a multipixel homodyne detection (MPHD) system, followed by digital post-processing of the acquired signals, e.g. multiplication of each signal by a suitable gain [14]. This allowed them to measure simultaneously an arbitrary quadrature on all the modes in a chosen mode set. The use of this method to emulate cluster states statistics was also sketched in that work.
It is also possible to generate multimode quantum states of light in the frequency domain. First experiments in this context have been performed with continuous pumping of an optical parametric oscillator (OPO) to generate cluster states between symmetric resonant frequencies of the cavity [15]. Alternatively, an intrinsically entangled multimode source of frequency modes is obtained by pumping an optical oscillator with a femtosecond pump, which possesses a frequency-comb spectrum, and which realizes a synchronously pumped optical parametric oscillator (SPOPO). In this context, the formation of three-mode squeezed non-classical states has already been demonstrated [16]. A multimode detection system analogous to the one employed in [13] can be implemented for the detection of these frequency modes, where the spatial pixels are replaced by frequency bins. An appealing advantage of this setup is that the source is now scalable, as all the modes are produced within only one optical device, and potentially hundreds of modes can be entangled. The remaining question is thus to assess to what extent such a method can be applied to generate cluster states and to MBQC.
This motivates our study of the ensemble of equivalent unitary operations which can be readily implemented by the multi-pixel detection followed by digital signal recombination. In particular, we give a procedure to assess whether and how a given cluster state can be realized with this method. Then, we apply this idea to perform Gaussian quantum computation in the measurement-based model and give some specific examples.

Modelization of the operations which can be performed on the modes
Let us first model an experiment that would consist in generating a multimode squeezed state and in performing a MPHD from the point of view of canonical transformations on the optical input modes. Such an experiment can be split into four steps, as can be seen in figure 1: (a) generating a composite beam made of independent orthogonal squeezed modes, which can be either spatial or frequency modes; (b) mixing it on a beam-splitter with a shaped local oscillator; (c) propagating it to multi-pixel detectors; and (d) recombining digitally the acquired signals.
As discussed above, there are several ways of generating multimode squeezed beams. Nevertheless, independently of the choice adopted, the result can be modeled as an optical field composed of N orthogonal and normalized modes of the electromagnetic field u i (ρ), where ρ could either be spatial coordinates (x, y) or frequency components ω depending on the realization. We assume that in each mode the state of the field is infinitely squeezed, say along thep quadrature, corresponding hence to thep-eigenstate |s p , beingp|s p = s|s p , with zero eigenvalue, i.e. |0 p = 1/ √ 2π ds|s q [17]. A brief discussion of the effect of a finite squeezing degree is carried out at the end of section 4.
These N modes are accompanied by vacuum modes in all the other modes of the complete modal basis {u i (ρ)}. With each mode are associated a pair of bosonic operatorsâ u i andâ † u i withâ u i = dρâ(ρ)u * i (ρ) (where dρ stands for ρ∈R 2 d 2 ρ if ρ is a 2D variable), andâ(ρ) = ∞ k=1â u k u k (ρ). These input modes can be grouped in a vectorâ u = ( a u , a † u ) T = (â u 1 ,â u 2 , . . . ,â † u 1 ,â † u 2 . . .) T . Once the multimode squeezed beam is generated, it is mixed on the beam-splitter with a local oscillator. Then, it is propagated to two arrays of detectors composed of pixels, each of surface S i , that collect the intensity of the incoming light from the two respective beams, and the differences of the two signals are taken pixel by pixel. This propagation is governed by diffraction in the spatial case, and by frequency dispersion realized with the help of a (a) Sketch of the experiment performed with frequency modes. The (infinitely)p-squeezed modes u i (ρ) generated by SPOPO are mixed with a local oscillator on a 50/50 beam-splitter. Both outputs are then dispersed and each frequency band is sent to a different pixel of the two multi-pixel detectors. The differences of the acquired signals are taken pixel by pixel. These are later processed by an ordinary computer which multiplies every trace by a suitable gain, for instance to obtain a cluster state. The additional correction step indicated in the dashed box is only required when performing a quantum computation; it reinterprets the measurement result on the mode of interest depending on the outcomes of the measurements on the other modes. Note that an analogue device can similarly be implemented for the detection of spatial input modes. (b) Functional modelization of the same experiment. OPO models squeezing quadratures of source modes, U T , light propagation from source to multi-pixel detector, LO , local oscillator phases in the pixel basis, and O, digital data processing. Refer to main text for details. prism in the frequency case (see figure 1). To describe the effect of the detection of the beam it is convenient to introduce another set of modes, that we shall call pixel modes; these are defined as and (1) implies that each pixel mode is a 'slice' of the local oscillator, i.e. it shares its shape but only on a window with width given by the pixel dimension (and zero outside). Note that the pixel modes form a complete basis, allowing to reconstruct all the possible modes only in the limit of infinite pixel number. The detection process can be modeled as an effective transformation corresponding to light propagation from the squeezed modes set to the pixel modes one, and to a homodyning in that basis with an intense local oscillator u LO,1 (ρ). Here for consistency we assume that ρ is, from the beginning, the variable that spatially describes the field in the multi-pixel detector plane. For instance, in figure 1 it does indeed correspond to frequency ω, but in the spatial case it depends on the actual imaging system between squeezing generation and the detectors. This transformation is described by a matrix U T acting on the bosonic operators: a v = U T a u , which brings the input More details can be found in appendix A.2 and in [13,18]. Note that in principle the detection matrix U T does not need to be square, since the number of pixels P may differ from the number of input modes-the latter may be potentially infinite.
To set the ideas, consider as an example that we will use later, the simple case in which four input squeezed modes are detected with the help of four pixel modes. Let us assume that the intensity distribution of the local oscillator mode is the first mode of this basis, i.e. u LO (ρ) = u 1 (ρ). We take a toy model for the input modes, describing the mode n by a square profile in which we introduce n − 1 π -phase shifts (or 'flips') at regular intervals, as illustrated in figure 2 (top panels). These modes are however qualitatively similar to the modes generated by a SPOPO, which in turn are close to Hermite-Gauss modes [16]. In this simple case, the matrix U T is unitary and takes the simple form Above we have assumed that the input state isp-squeezed on each mode. Experimentally, it is in principle possible to modify the squeezing quadrature. For simplicity, we model that by changing the input mode phase instead of its state, which is completely equivalent in an homodyne detection scheme. This is done by introducing a diagonal dephasing matrix OPO with OPO i, j = e iφ i δ i j . The input modes entering the detection system are u i (ρ) = e −iφ i u i (ρ) and a u = OPO a u . The transformation of these input modes to the pixel set as defined above is now a v = U T * OPO a u , as * OPO moves the phase shifted modes to the modes a u and U T is the transformation from a u to a v .
With the same approach, the ability to shape the phase of the local oscillator by phaseshifting the pixel modes can be modeled by the action of a diagonal matrix Finally, after the detection, it is possible to digitally recombine the electronic signals coming from each pixel by multiplying them with real gains [14]. This allows us to detect the field according to desired modes (e.g. as we shall see later, according to the modes which correspond to a cluster state), and amounts to applying a P-dimensional orthogonal matrix, i.e. a v = O a v ≡ a out (see also the appendix).
In summary, we find that the series of the possible transformations that we can perform on the input modes are As implicit in the previous discussion, the matrices O and LO are easily tunable in the experiments. The dephasing matrix OPO can also be adjusted to some extent by changing the spectral shape of the pump [19].

Characterization of the feasible transformations
We now address the characterization of the subset of unitary operations, acting on the annihilation operator a u , which can be emulated by means of equation (3). The question: 'which class of unitary transformations can we emulate in the laboratory by means of the MPHD plus data processing'? can be recast into the following question: 'given a unitary matrix U th , can it be written in the form of equation (3)'? Answering this questions requires finding a solution for the parameters of O and LO (given a certain fixed detection matrix G) such that It is easily seen that a necessary and sufficient condition for equation (4) is with U th = U th G † and where D is a diagonal matrix with unit modulus complex elements. The necessity is simply proved by assuming that equation (4) holds, and by direct calculation of the product U th T U th upon substitution of equation (3), which yields U th where we have used T LO = LO since this matrix is diagonal, and where 2 LO is obviously a diagonal matrix with unit modulus complex elements. For the sufficiency, we can show that if equation (5) holds, we can always find at least a set of 'experimental' matrices LO and O such that their combination of the form of U MPHD in equation (3) satisfies equation (4). This is achieved by taking as a realization of the diagonal matrix LO any of the solutions of the equation Then, it is automatically ensured that, for each solution LO , the matrix Four-mode linear cluster state. Each circle represent a cluster mode and the edges represent entanglement between modes, namely that the two modes have been acted upon with a C Z -type interaction e iq i ⊗q j [6].
is orthogonal and fulfils equation (4). Indeed by direct calculation we obtain LO LO G = U th , which completes our proof. As a major result of our work this shows that, given a certain unitary matrix U th , if a solution of equation (5) exists, with the use of equations (6) and (7) we can determine suitable experimental parameters in terms of the gain coefficients O and the local oscillator phase shaping LO . These allow the statistics of the quadrature outcomes to be obtained with the MPHD measurement as if the matrix U th -usually corresponding to a complex optical network-had been applied to the input squeezed modes, and the same quadratures were measured one by one. Let us also note that our proposal is versatile, in the sense that different matrices U th can be produced with little or no reconfiguration of the experimental setup.

Cluster states in the multimode beam
To provide an illustration of our method, we consider the simple case in which four input squeezed modes are detected with the help of four pixel modes, which entails the mode transformation given in equation (2). Let us consider as a concrete example the case in which the desired unitary transformation to be emulated is the one which transforms N squeezed modes in a linear cluster state.
In the case N = 4 (figure 3), such a transformation is represented by the following matrix [10]: In order to find a simple algebraic solution for the problem we consider the dephasing matrix * OPO = diag(1, i, i, 1). With the use of equations (5), (6), (8) and (2) we find Since the matrix in equation (9) is diagonal, the condition in equation (5) is satisfied, which implies that the corresponding cluster state matrix (8)   The input mode is coupled via a beam-splitter (dashed line) to the first mode of a three-mode cluster state, and then these two modes are simultaneously measured, which teleport the input state onto the second mode. Then, a suitable quadrature measurement is performed on the second mode of the cluster, projecting the last mode onto the desired output. With the notation used in this picture [20], the angles appearing on each mode specify which quadraturê q i sin θ i +p i cos θ i shall be measured on it in order to implement the desired operation.
satisfying equations (6) and (7) for U th = U lin we can chose for instance the set with γ ≡ arctan 1 2 and α = cos −1 1 2 (1 − 2 √ 5 ). Hence, we conclude that shaping the phase of the local oscillator on each pixel mode as prescribed by the matrix LOlin1 in equation (2) and choosing the digital recombination gains according to O lin1 allows us to obtain statistics for the quadrature measurement as if a network of beam-splitters generating a linear cluster state had been applied to the input squeezed modes. We have checked that a numerical solution exists for other choices of the dephasing matrix, namely for the one corresponding to the physical dephasings of the cavity modes in the experiment of [16].

Simple examples of quantum computations
We now address the emulation of a simple example of quantum computation in the measurement-based model (MBQC). Since we are performing homodyne detection on all the modes, we are restricted to Gaussian operations [6]. This is a particularly simple framework to start with, because the implementation of Gaussian operations in the MBQC model does not require effective adaptivity-usually necessary to implement operations deterministically-and the measurements can be performed on all the modes simultaneously. The outcomes are recorded, and the 'adaptation' of the measurement basis depending on the outcomes of previous measurements can be recovered in a further classical processing step (e.g. by electronically correcting the photocurrents, see dashed box in figure 1) [6].
We consider here single-mode operations. A quantum computation implies first coupling an initial quantum state to be processed to a cluster state, and then suitable measurements of the quadratures on each mode of the cluster (figure 4). The measurement of the last mode is not part of the computation but provides a read-out of the output state, the latter being the 'result' of the computation. The coupling between the initial state and the first mode of the cluster can be realized by a generalized teleportation [21]: these two modes are coupled by a beam-splitter interaction, e.g. (ˆa in . Measuring the quadratureŝ q in sin θ in +p in cos θ in andq 1 sin θ 1 +p 1 cos θ 1 allows to transfer the input state on the second mode of the cluster state modulo a transformation M tele (θ in , θ 1 ) controllable via the choice of the measurement angles θ in , θ 1 [21]. 5 This corresponds to the effective transformation of the quadratures and θ ± = θ in ± θ 1 (the difference from equation (8) of [21] being due to a different definition of the angles θ in , θ 1 ). In particular, the choice θ in = θ 1 = π/2 teleports exactly the input state onto the second mode of the cluster, i.e.
Let us consider as an example of quantum computation the Fourier transform of the initial state. Together with the quadrature displacement D 1,q (s) = e isq and the shear D 2,q (s) = e isq 2 this is one of the elementary operations for universal one-mode Gaussian quantum computation [6], and it has been recently implemented experimentally with the use of a beam-splitter network generating a four-mode cluster state [8]. The Fourier transform that we wish to implement should act on the quadratures of the input state (qˆp ) by rotating them, i.e. F = 0 −1 1 0 . For a general cluster state it has been shown [6,21] that measuring the quadrature D † q,2 (s)p D q,2 (s) =p + sq = g(x sin θ +p cos θ) on one of its modes, with g = √ 1 + s 2 and θ = arctan s, effectively transfers to the adjacent (say, right hand) mode a state transformed according to the matrix apart from 'by-product' displacements which are not accounted for by equation (14). These may be corrected in the end of the computation, although this is not strictly necessary [6]. Note that the measurement of the quadrature in equation (13) corresponds to the measurement of the quadraturep after the transformationâ → e iθâ (modulo the scale factor g). We see that in the simple case under consideration a single computational measurement, to be performed on the mode 2 of the cluster state in figure 4, is necessary after the teleportation step. Indeed, with s 2 = 0 we obtain from equation (14) M(s 2 leading to M(s 2 = 0)M tele (θ in = π 2 , θ 1 = π 2 ) = F. As expressed by equation (13), this corresponds to the measurement of the quadraturep on the mode 2. As well as for the measurements implementing the teleportation, the corresponding outcome should be recorded to implement the final correction step above mentioned [6,21].
Hence, to summarize, in order to implement a Fourier transform on an initial state, we have to couple it via a beam-splitter interaction on a three-mode cluster state, and then to measure on each mode the quadrature specified by the angles (θ in , θ 1 , θ 2 , θ 3 ) = ( π 2 , π 2 , 0, θ 3 ), respectively (see figure 4). The measurement on the third mode of the cluster yields an outcome probing the implementation of the desired transformation, and hence it can be performed along any quadrature.
As anticipated, we take the perspective of measuring (simultaneously) the quadraturê p on all the output modes after suitable rotations have been applied, here expressed by the matrix D meas = diag(e iθ in , e iθ 1 , e iθ 2 , e iθ 3 ) = (i, i, 1, e iθ 3 ). Hence the unitary transformation which transforms the threep-squeezed modes plus the input state (â in ,â 1 ,â 2 ,â 3 ) into the three-mode cluster state coupled with the initial state by a beam-splitter interaction reads where for simplicity we have indicated with the same notation the matrix U 1,2,3 lin , which builds the three-mode linear cluster state from three (infinitely) p-squeezed modes, and the same matrix times the identity on the fourth mode I in ⊗ U 1,2,3 lin . The matrix U 1,2,3 lin can be chosen as (see appendix B) which, setting θ 3 = 0, results in the 'desired' matrix In order to be able to implement the matrix U ft by means of the MPHD method, we have to find a combination of LO and O which satisfies equation (4) with U th = U ft . 6 We may assume here appropriate dephasing between the modes, e.g. * OPO = diag(1, 1, −i, i). With the use of equation (6) we obtain Among the 2 N = 16 possible solutions for LO ft and the corresponding orthogonal matrix O ft satisfying equations (6) and (7) for U th = U ft we can choose for instance the set where we have set ζ = arctan(1/ √ 2). The corresponding orthogonal matrix is with β = cos −1 ( 1 2 (1 + 2 3 )). Note that all the modes are simultaneously measured according to the procedure described in the previous section. It is readily verified that O ft1 LO ft 1 G = U ft . Hence, with the MPHD we can emulate the formation of a three-mode cluster state, the interaction of its first mode with the input mode via a beam-splitter, and the measurement of the suitable quadratures to propagate the quantum computation along the cluster. Indeed, shaping the phase of the local oscillator on the pixel modes as prescribed by the matrix LOft1 in equation (2) and choosing the digital recombination gains according to O ft1 allows us to obtain a statistics for the quadrature measurement of mode 4 as if this would contain the Fourier transform of the state in the first mode (here in ap squeezed state), modulo reinterpreting these results taking into account the by-product operators. Experimentally, addressing the variety of possible input states could be realized via proper seeding of the multimode cavity with spatially and temporally shaped light.
The same procedure also allows to implement the quadrature displacement e iqs . This is not a symplectic operation, therefore the formalism of [21] cannot be used to determine the angles to measure. In reference to the scheme of figure 4, after teleporting the input state onto the second mode of the cluster state, the measurement of the quadrature e −ih(q)p e ih(q) on this mode allows implemention of the gate e −ih(q) , where h(q) is an arbitrary function of the quadratureq [6]. Then, to implement the gate e isq the quadraturep sq = e −isqp e isq =p + s has to be measured. This can be achieved by measuringp and adding s to the result. Hence, everything goes exactly as for the Fourier transform in terms of the quadratures to measure, except that we have to add 's' to the measurement of the quadraturep on mode 2 of the cluster 7 .
We have also checked that the other elementary operations of the Gaussian universal set {e iqs , e iq 2 s , F, C Z }, namely e iq 2 s , C Z may be implemented approximately, i.e. a solution exists for matrices U th close (in the sense of the matrix distance d( is the Frobenius norm) to the ones associated with these operations. The physical meaning of these approximate operations deserves however a detailed study and will be addressed elsewhere.
Experimentally, the degree of squeezing is necessarily finite. This does not change the operations to perform in order to obtain a desired cluster state, with respect to the infinite squeezing case [28]. Though this entails errors in a quantum computation, which manifest as extra noise that modifies the result of the quadrature measurement in the last mode coding the result of the computation [6]. These extra noise terms affect a quantum computation in our proposal in the same way as in other proposals for Gaussian quantum computation, and can be precisely addressed with the technique sketched in the supplementary information section of [8].
Issues related to the possibility of providing a full fault-tolerant quantum computation are given in [24,25].

Conclusion and prospective views
In this work we have demonstrated the possibility of measuring cluster states and implementing Gaussian quantum computation in an extremely compact fashion. The method is based on the simultaneous measurement of all the (highly) squeezed optical modes produced in a single cavity by MPHD, and on the classical post-processing of the acquired signals, which are multiplied by adequate electronic gains. This procedure requires the determination of the suitable phase shape of the local oscillator to be employed in the MPHD, as well as the gains to apply to the traces recorded in each mode. In particular, as a first example we have provided the explicit solution in terms of these experimental parameters which allows the formation of a four-mode linear cluster state to be mimicked. As a simple example of quantum computation, we have considered the implementation of a Fourier transform on an input mode.
We remark that the state preparation and the computation in the setup that we have presented are performed in a single step, yielding hence a quantum depth equal to 1-the quantum depth being indeed defined as the number of measurement steps that a computation task requires, when trying to do as many measurements as possible at the same time. This implies that our method is not suitable for achieving universal quantum computation [27]. However, it is interesting to explore sub-universal quantum operations, following a recent trend in the field [29][30][31]. Furthermore, depth 1 computations in the measurement-based setting may provide an advantage with respect to the circuit model for performing the same operation [26].
To go in this direction, it is important to stress that by the Gottesmann-Knill theorem [6,22] all the results of manipulations with Gaussian elements only, such as the squeezed states and homodyne detection discussed in this work, can be efficiently simulated with a classical computer. In order to outperform classical processing capability at least one non-Gaussian operation is necessary. This can be readily done by extracting a controlled fraction of a given mode and placing a photon resolving detector on this extracted mode. Mode extraction can be performed either by a slightly reflecting beam-splitter in front of one of the detection pixels, or by using mode dependent nonlinear effects [23]. In these cases, no other effective adaptation of the measurement basis is required than the one performed at the stage of the post-processing, due to the fact that the homodyne detection to be performed on the remaining modes is a Gaussian measurement [6]. Though, as already stated, not sufficient to achieve universality, this would allow a specific class of non-Gaussian operations to be performed, obtained as the combination of a non-Gaussian gate followed by a Gaussian one, the latter belonging to the set that we have characterized. The precise characterization of the resulting non-Gaussian operation stems as an interesting perspective of this work.

A.3. Signal recombination
The pixel modes can be digitally multiplied by real gains and recombined, leading to the modes v i (ρ) = j g i j v j (ρ). Since the output functions v i (ρ) should form an orthonormal basis, we are restricted to gains which are elements of an orthogonal matrix, i.e. O i, j = g i j . Then the corresponding annihilation operators arê and we can obtain information on a quadrature aŝ

Appendix B. Three-mode linear cluster state matrix
In this appendix we derive the matrix U 1,2,3 lin which builds-up the three-component linear cluster state. For a general cluster state identified by an adjacency matrix V the corresponding unitary matrix which brings from the squeezed modes to the given cluster state can be chosen as any possible matrix satisfying [28] U = X + iY,