Compact Gaussian quantum computation by multi-pixel homodyne detection

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 measurement-based quantum computing [2–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–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 multi-pixel 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.

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.

Zoom In Zoom Out Reset image size

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 the $\hat {p}$ quadrature, corresponding hence to the $\hat {p}$-eigenstate |s〉_p, being $\hat {p}|s\rangle _p = s|s \rangle _p$, with zero eigenvalue, i.e. $|0\rangle _p = 1/\sqrt {2 \pi } \int \mathrm {d}s|s\rangle _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 $\hat a_{u_i}$ and $\hat {a}^{\dagger }_{u_i}$ with $\hat {a}_{u_i} = \int \mathrm {d} \rho \hat {a}(\rho ) u^{*}_{i} (\rho )$ (where $\int \mathrm {d} \rho$ stands for $\int \hspace {-0.2cm} \int _{\rho \in \mathbb {R}^2} {\mathrm { d}}^2 \rho$ if ρ is a 2D variable), and $\hat {a}(\rho ) = \sum _{k=1}^{\infty } \hat {a}_{u_k} u_k (\rho )$. These input modes can be grouped in a vector $\boldsymbol{\hat{\rm a}_{\rm u}} = (\vec {a}_{u}, \vec {a}_{u}^{\dagger })^{\rm{T}} = (\hat {a}_{u_1},\hat {a}_{u_2}, \ldots ,\hat {a}_{u_1}^{\dagger }, \hat {a}_{u_2}^{\dagger }\ldots )^{\rm{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 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

$$\begin{equation} v_{i} (\rho)= \left\{ \begin{array}{@{}ll@{}} \kappa_i \, \, u_{\mathrm{LO}} (\rho), & \forall \rho \in S_{i}, \\ 0 & {\rm elsewhere} \end{array} \right. \end{equation} \tag{ 1 }$$

and correspondingly $\hat {a}_{v_i} = \int \mathrm {d} \rho \hat {a}(\rho ) v_{i}^* (\rho )$, where κ_i is a normalizing constant such that $\int _{\rho \in \mathbb {R}} \mathrm {d} \rho v_{i}^2 (\rho ) = \kappa _i^2 \int \hspace {-0.2cm} \int _{\rho \in S_{i}} \mathrm {d}^2 \rho \, | u_{\mathrm {LO}} (\rho ) |^2 = 1$. Equation (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: $\vec {a}_{v} = U_{\rm{T}} \vec {a}_{u}$, which brings the input modes $\hat {a}_{u_i}$ onto the pixel modes $\hat {a}_{v_i}$, with ${U_{\rm{T}}}_{ij} = \kappa _i \int _{\rho \in S_i} \mathrm {d} \rho \,u^{*}_{\mathrm {LO}}(\rho ) u_j(\rho )$. 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₁(ρ). 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

$$\begin{equation} U_{\rm{T}} = \frac{1}{2} \left( \begin{array}{@{}cccccccc@{}} 1 & 1 & -1 & 1 \\ 1 & 1 & 1 & -1 \\ 1 & -1 & 1 & 1 \\ 1 & -1 & -1 & -1 \\ \end{array} \right). \end{equation} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

Above we have assumed that the input state is $\hat {p}$-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 Δ_OPOi,j = e^iϕ_iδ_ij. The input modes entering the detection system are u_i'(ρ) = e^−iϕ_iu_i(ρ) and $\vec {a}_{u'} = \Delta_{\rm OPO}\vec {a}_{u}$. The transformation of these input modes to the pixel set as defined above is now $\vec {a}_{v} = U_{\rm{T}} \Delta _{\mathrm {OPO}}^{*} \vec {a}_{u'}$, as Δ*_OPO moves the phase shifted modes to the modes $\vec {a}_{u}$ and U_T is the transformation from $\vec {a}_{u}$ to $\vec {a}_{v}$.

With the same approach, the ability to shape the phase of the local oscillator by phase-shifting the pixel modes can be modeled by the action of a diagonal matrix $\vec {a}_{v'} = \Delta _{\mathrm {LO}} \vec {a}_{v}$ with elements Δ_LOi,j = e^iφ_iδ_i,j, leading to the modes v'_i(ρ) = e^−iφ_iv_i(ρ).

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. $\vec {a}_{v^{''}} = O \vec {a}_{v'} \equiv \vec {a}_{\mathrm {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

$$\begin{equation} \vec{a}_{\mathrm{out}} = O \Delta_{\mathrm{LO}} G \, \, \vec{a}_{u} \equiv U_{\mathrm{MPHD}} \, \, \vec{a}_{u}, \end{equation} \tag{ 3 }$$

where Δ_LO = diag(e^iφ₁, e^iφ₂, ..., e^iφ_N), O is real orthogonal and G = U_TΔ*_OPO is a fixed matrix. Note that since $\Delta _{\mathrm {LO}} \Delta _{\mathrm {LO}}^\dagger = \mathcal I$ and $G G^{\dagger } = \mathcal I$ the product OΔ_LOG = U_MPHD is unitary. 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].

We now address the characterization of the subset of unitary operations, acting on the annihilation operator $\vec {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

$$\begin{equation} U_{\mathrm{th}} = U_{\mathrm{MPHD}}. \end{equation} \tag{ 4 }$$

It is easily seen that a necessary and sufficient condition for equation (4) is

$$\begin{equation} {U_{\mathrm{th}}'}^{\rm{T}} U_{\mathrm{th}}' = D, \end{equation} \tag{ 5 }$$

with U_th' = U_thG^† 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' ^TU'_th upon substitution of equation (3), which yields U_th' ^TU'_th = G*G^TΔ^T_LOO^TOΔ_LOGG^† = Δ²_LO, where we have used Δ^T_LO = Δ_LO since this matrix is diagonal, and where Δ²_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

$$\begin{equation} \Delta_{\mathrm{LO}} = D^{\frac{1}{2}} = ({U_{\mathrm{th}}'}^{\rm{T}} U_{\mathrm{th}}')^{\frac{1}{2}}. \end{equation} \tag{ 6 }$$

Then, it is automatically ensured that, for each solution Δ_LO, the matrix

$$\begin{equation} O = U_{\mathrm{th}}' \Delta_{\mathrm{LO}}^{-1} \end{equation} \tag{ 7 }$$

is orthogonal and fulfils equation (4). Indeed by direct calculation we obtain U_MPHD = OΔ_LOG = U_th'Δ⁻¹_LOΔ_LOG = 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.

3.1. 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]:

$$\begin{equation} U_{\mathrm{th}} = U_{\mathrm{lin}} = \left( \begin{array}{@{}cccccccc@{}}\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{10}} & \frac{2 {\rm{i}}}{\sqrt{10}} & 0 \\ \frac{\rm{i}}{\sqrt{2}} & -\frac{\rm{i}}{\sqrt{10}} & \frac{2}{\sqrt{10}} & 0 \\ 0 & -\frac{2}{\sqrt{10}} & \frac{\rm{i}}{\sqrt{10}} & \frac{\rm{i}}{\sqrt{2}} \\ 0 & -\frac{2{\rm{i}}}{\sqrt{10}} & -\frac{1}{\sqrt{10}} & \frac{1}{\sqrt{2}} \\ \end{array}\!\!\right). \end{equation} \tag{ 8 }$$

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

$$\begin{equation} {U_{\mathrm{lin}}'}^{\rm{T}} U_{\mathrm{lin}}' = \left( \begin{array}{@{}cccccccc@{}} \frac{-2 - {\rm{i}}}{\sqrt{5}} & 0 & 0 & 0 \\ 0 & \frac{2 - {\rm{i}}}{\sqrt{5}} & 0 & 0 \\ 0 & 0 & \frac{2 + {\rm{i}}}{\sqrt{5}} & 0 \\ 0 & 0 & 0 & \frac{-2 + {\rm{i}}}{\sqrt{5}} \end{array}\!\!\right). \end{equation} \tag{ 9 }$$

Since the matrix in equation (9) is diagonal, the condition in equation (5) is satisfied, which implies that the corresponding cluster state matrix (8) can be implemented experimentally. Among the 2^N = 16 possible solutions for Δ_LOlin and the corresponding orthogonal matrix O_lin satisfying equations (6) and (7) for U_th = U_lin we can chose for instance the set

$$\begin{eqnarray} \fl { \Delta_{{\rm LOlin}1} = \left(\begin{array}{@{}cccccccc@{}}{\rm e}^{{\rm{i}}\frac{\gamma + \pi}{2}} & 0 & 0 & 0 \\ 0 & {\rm e}^{-{\rm{i}}\frac{\gamma}{2}} & 0 & 0 \\ 0 & 0 & {\rm e}^{{\rm{i}}\frac{\gamma}{2}} & 0 \\ 0 & 0 & 0 & {\rm e}^{-{\rm{i}}\frac{\gamma + \pi}{2}} \end{array}\!\!\right) = \left(\begin{array}{@{}cccccccc@{}} -0.23+0.97 {\rm{i}} & 0 & 0 & 0 \\ 0 & 0.97-0.23 {\rm{i}} & 0 & 0 \\ 0 & 0 & 0.97+0.23 {\rm{i}} & 0 \\ 0 & 0 & 0 & -0.23-0.97 {\rm{i}} \end{array} \right)}, \nonumber \\ \fl { O_{{\rm lin}1} = \frac{1}{\sqrt{2}} \left( \begin{array}{cccc} -\cos \alpha & \sin \alpha & \sin \alpha & -\cos \alpha \\ \sin \alpha & -\cos \alpha & \cos \alpha & -\sin \alpha \\ \sin \alpha & \cos \alpha & \cos \alpha & \sin \alpha \\ -\cos \alpha & -\sin \alpha & \sin \alpha & \cos \alpha \end{array}\!\!\right) = \left(\begin{array}{@{}cccccccc@{}} -0.16 & 0.69 & 0.69 & -0.16 \\ 0.69 & -0.16 & 0.16 & -0.69 \\ 0.69 & 0.16 & 0.16 & 0.69 \\ -0.16 & -0.69 & 0.69 & 0.16 \end{array}\!\!\right)} \end{eqnarray} \tag{ 10 }$$

with $\gamma \equiv {\rm arctan} \frac{1}{2}$ and $\alpha = {\rm cos}^{-1} \sqrt{\frac{1}{2}(1-\frac{2}{\sqrt{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].

Zoom In Zoom Out Reset image size

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. $({\hat{a}_{\rm{in}}' \atop \hat{a}_{1}'})= U^{{\rm{in}},1}_{\mathrm {BS}}({\hat{a}_{\rm{in}} \atop \hat{a}_{1}})$ with $U^{{\rm{in}},1}_{\mathrm {BS}} = \frac{1}{\sqrt{2}}({1 \ \ \mathrm {i} \atop \mathrm {i} \ \ 1})$. Measuring the quadratures $\hat{q}'_{\rm{in}} \sin \theta_{\rm{in}} + \hat{p}'_{\rm{in}} \cos \theta _{\rm{in}}$ and $\hat{q}'_{1} \sin \theta_{1} + \hat{p}'_{1} \cos \theta_{1}$ allows to transfer the input state on the second mode of the cluster state modulo a transformation M_tele(θ_in,θ₁) controllable via the choice of the measurement angles θ _in,θ₁ [21].⁵ This corresponds to the effective transformation of the quadratures

$$\begin{equation} \fl \left(\begin{array}{@{}cc@{}} \hat{q}^{'} \\ \hat{p}^{'} \\ \end{array}\!\!\right) = M_{\rm{tele}}(\theta_{\rm{in}},\theta_{1}) \left(\begin{array}{@{}cc@{}} \hat{q} \\ \hat{p} \\ \end{array}\!\!\right)\quad \mbox{with}\; M_{\rm{tele}}(\theta_{\rm{in}},\theta_{1}) = - \left(\begin{array}{@{}cc@{}}\displaystyle\frac{\cos \theta_+}{\cos \theta_-} & \displaystyle\frac{\sin \theta_- - \sin \theta_+}{\cos \theta_-} \\ \displaystyle\frac{\sin \theta_- + \sin \theta_+}{\cos \theta_-} & \displaystyle\frac{\cos \theta_+}{\cos \theta_-} \\ \end{array} \right) \end{equation} \tag{ 11 }$$

and θ_± = θ_in ± θ₁ (the difference from equation (8) of [21] being due to a different definition of the angles θ_in,θ₁). In particular, the choice θ_in = θ₁ = π/2 teleports exactly the input state onto the second mode of the cluster, i.e.

$$\begin{equation} M_{\rm{tele}} \left( \theta_{\rm{in}} = \frac{\pi}{2},\theta_1 = \frac{\pi}{2} \right)= \left(\begin{array}{@{}ll@{}} 1 & 0 \\ 0 & 1 \\ \end{array}\!\right). \end{equation} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

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) = {\mathrm { e}}^{{\rm{i}} s \hat {q}}$ and the shear $D_{2,q}(s) = {\mathrm { e}}^{{\rm{i}} s \hat {q}^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 $({\hat{q} \atop \hat{p}})$ by rotating them, i.e. $F = \left(\begin{array}{@{}c@{\;\;}c@{}}0& -1\\ 1 & 0\end{array}\right)$. For a general cluster state it has been shown [6, 21] that measuring the quadrature

$$\begin{equation} D^{\dagger}_{q,2}(s) \hat{p} D_{q,2}(s) = \hat{p} + s \hat{q} = g (\hat x \sin \theta + \hat{p} \cos \theta) \end{equation} \tag{ 13 }$$

on one of its modes, with $g = \sqrt {1 + s^2}$ and $\theta = \arctan s$, effectively transfers to the adjacent (say, right hand) mode a state transformed according to the matrix

$$\begin{equation} M(s) = \left(\begin{array}{@{}cc@{}} -s & -1 \\ 1 & 0 \\ \end{array}\!\right), \end{equation} \tag{ 14 }$$

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 quadrature $\hat {p}$ after the transformation $\hat {a} \rightarrow {\mathrm { e}}^{{\rm{i}} \theta } \hat {a}$ (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₂ = 0 we obtain from equation (14)

$$\begin{equation} M(s_2 = 0) = \left(\begin{array}{@{}cc@{}} 0 & -1 \\ 1 & 0 \\ \end{array}\!\right), \end{equation} \tag{ 15 }$$

leading to $M(s_2 = 0) M_{\rm{tele}}(\theta _{\rm{in}} = \frac {\pi }{2},\theta _1 = \frac {\pi }{2}) = F$. As expressed by equation (13), this corresponds to the measurement of the quadrature $\hat {p}$ 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 $(\theta _{\rm{in}},\theta _1,\theta _2,\theta _3) = (\frac {\pi }{2},\frac {\pi }{2},0,\theta _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 quadrature $\hat {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θ₁ ,e^iθ₂,e^iθ₃) = (i,i,1,e^iθ₃). Hence the unitary transformation which transforms the three $\hat {p}$-squeezed modes plus the input state $(\hat a_{\rm{in}},\hat a_{1},\hat a_{2},\hat a_{3})$ into the three-mode cluster state coupled with the initial state by a beam-splitter interaction reads

$$\begin{equation} U_{\mathrm{ft}} = D_{\mathrm{meas}} (U^{{\rm{in}},1}_{\mathrm{BS}} \otimes I^{2,3}) (I^{\rm in} \otimes U^{1,2,3}_{\mathrm{lin}}) \equiv D_{\mathrm{meas}} U_{\mathrm{BS}} U^{1,2,3}_{\mathrm{lin}}, \end{equation} \tag{ 16 }$$

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ⁱⁿ⊗U^1,2,3_lin. The matrix U^1,2,3_lin can be chosen as (see appendix B)

$$\begin{equation} U^{1,2,3}_{\mathrm{lin}}= \left( \begin{array}{@{}ccc@{}} 0 & -\sqrt{\frac{2}{3}} & -\frac{{\rm{i}}}{\sqrt{3}} \\ -\frac{{\rm{i}}}{\sqrt{2}} & -\frac{\rm{i}}{\sqrt{6}} & -\frac{1}{\sqrt{3}} \\ -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{6}} & -\frac{\rm{i}}{\sqrt{3}} \end{array}\!\right), \end{equation} \tag{ 17 }$$

which, setting θ₃ = 0, results in the 'desired' matrix

$$\begin{equation} U_{\mathrm{ft}} = \left( \begin{array}{@{}cccc@{}} \frac{\rm{i}}{\sqrt{2}} & 0 & \frac{1}{\sqrt{3}} & \frac{\rm{i}}{\sqrt{6}} \\ -\frac{1}{\sqrt{2}} & 0 & -\frac{\rm{i}}{\sqrt{3}} & \frac{1}{\sqrt{6}} \\ 0 & -\frac{\rm{i}}{\sqrt{2}} & -\frac{\rm{i}}{\sqrt{6}} & -\frac{1}{\sqrt{3}} \\ 0 & -\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{6}} & -\frac{\rm{i}}{\sqrt{3}} \end{array}\!\right). \end{equation} \tag{ 18 }$$

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.⁶ 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

$$\begin{equation} {U_{\mathrm{ft}}'}^{\rm{T}} U_{\mathrm{ft}}' =\left( \begin{array}{@{}cccc@{}} -\frac{{\rm{i}}+\sqrt{2}}{\sqrt{3}} & 0 & 0 & 0 \\ 0 & \frac{{\rm{i}}+\sqrt{2}}{\sqrt{3}} & 0 & 0 \\ 0 & 0 & \frac{-{\rm{i}}+\sqrt{2}}{\sqrt{3}} & 0 \\ 0 & 0 & 0 & -\frac{-{\rm{i}}+\sqrt{2}}{\sqrt{3}} \end{array}\!\right). \end{equation} \tag{ 19 }$$

Among the 2^N = 16 possible solutions for Δ_LOft and the corresponding orthogonal matrix O_ft satisfying equations (6) and (7) for U_th = U_ft we can choose for instance the set

$$\begin{eqnarray} \fl { {{\Delta_{{\rm LO}_{\rm ft}1}}}\!=\!\left(\begin{array}{@{}cccc@{}} \mathrm{e}^{\frac{1}{2} \mathrm{i} (\zeta +\pi )} & 0 & 0 & 0 \\ 0 & \mathrm{e}^{\frac{\mathrm{i} \zeta }{2}} & 0 & 0 \\ 0 & 0 & \mathrm{e}^{-\frac{\mathrm{i} \zeta }{2}} & 0 \\ 0 & 0 & 0 & \mathrm{e}^{-\frac{1}{2} \mathrm{i} (\zeta +\pi )} \end{array}\right)\!=\!\left(\begin{array}{@{}cccc@{}} -0.3+0.95 \mathrm{i} & 0 & 0 & 0 \\ 0 & 0.95+0.3 \mathrm{i} & 0 & 0 \\ 0 & 0 & 0.95-0.3 \mathrm{i} & 0 \\ 0 & 0 & 0 & -0.3-0.95 \mathrm{i} \end{array}\!\right),}\nonumber\\ \end{eqnarray} \tag{ 20 }$$

where we have set $\zeta = \arctan (1/\sqrt {2})$. The corresponding orthogonal matrix is

$$\begin{equation} \fl {O_{\mathrm{ft}}}_1\!=\!\frac{1}{\sqrt{2}}\left(\begin{array}{@{}cccc@{}}\cos \beta & \sin \beta & -\sin \beta & -\cos \beta \\ \sin \beta & -\cos \beta & -\cos \beta & \sin \beta \\ -\cos \beta & -\sin \beta & -\sin \beta & -\cos \beta \\ \sin \beta & -\cos \beta & \cos \beta & -\sin \beta \end{array}\right)\!=\!\left(\begin{array}{@{}cccc@{}} 0.67 & 0.21 & -0.21 & -0.67 \\ 0.21 & -0.67 & -0.67 & 0.21 \\ -0.67 & -0.21 & -0.21 & -0.67 \\ 0.21 & -0.67 & 0.67 & -0.21 \end{array}\right) \end{equation} \tag{ 21 }$$

with $\beta = \cos ^{-1}(\sqrt {\frac {1}{2}(1+\sqrt {\frac {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Δ_LOft1G = 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 a $\hat {p}$ 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 ${\mathrm { e}}^{{\rm{i}} \hat {q} s}$. 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 ${\mathrm { e}}^{-{\rm{i}} h(\hat {q})} \hat {p}\, {\mathrm { e}}^{{\rm{i}} h(\hat {q})}$ on this mode allows implemention of the gate ${\mathrm { e}}^{-{\rm{i}} h(\hat {q})}$, where $h(\hat {q})$ is an arbitrary function of the quadrature $\hat {q}$ [6]. Then, to implement the gate ${\mathrm { e}}^{{\rm{i}} s \hat {q}}$ the quadrature $\hat {p}_{s \hat {q}} = {\mathrm { e}}^{-{\rm{i}} s \hat {q}} \hat {p} \,{\mathrm { e}}^{{\rm{i}} s \hat {q}} = \hat {p} + s$ has to be measured. This can be achieved by measuring $\hat {p}$ 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 quadrature $\hat {p}$ on mode 2 of the cluster⁷.

We have also checked that the other elementary operations of the Gaussian universal set {${\mathrm { e}}^{{\rm{i}} \hat {q} s}, {\mathrm { e}}^{{\rm{i}} \hat {q}^2s}, F, C_Z$}, namely ${\mathrm { e}}^{{\rm{i}} \hat {q}^2s}$, C_Z may be implemented approximately, i.e. a solution exists for matrices U_th close (in the sense of the matrix distance d(M₁,M₂) = ∥M₁ − M₂∥_F , where $\parallel{\kern-3pt}M{\kern-3pt}\parallel_{F} = \sqrt {\sum _{i,j} |M_{i,j}|^2 }$ 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].

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–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.

We acknowledge useful discussions with J-F Morizur and D Markham. The research is supported by the ERC starting grant FREQUAM, and by the ANR project NAMOCS. CF is a member of the Institut Universitaire de France.

A.1. Basis and notations

We introduce here the notations and conventions that we are going to use. The local bosonic operator $\hat {a}(\rho )$ satisfies

$$\begin{equation} [\hat{a}(\rho),\hat{a}^{\dagger}(\rho') ] = \delta(\rho - \rho') \end{equation} \tag{ A.1 }$$

and is related to the local quadrature operators by

$$\begin{eqnarray} \hat{q} (\rho) &=& \hat{a}^{\dagger}(\rho) + \hat{a}(\rho), \nonumber \\ \hat{p} (\rho) &=& {\rm{i}} \left( \hat{a}^{\dagger}(\rho) - \hat{a}(\rho) \right), \end{eqnarray} \tag{ A.2 }$$

so that $[\hat {q}(\rho ),\hat {p}(\rho ') ] = 2\mathrm { i} \delta (\rho - \rho ')$. Given an orthonormal basis of wave functions u_k(ρ) satisfying the orthogonality relation

$$\begin{equation} \int \mathrm{d} \rho \, u^{*}_k (\rho) u_l (\rho) = \delta_{k,l}, \end{equation} \tag{ A.3 }$$

the local bosonic operator can be decomposed as

$$\begin{equation} \hat{a}(\rho) = \sum_{k=0}^{\infty} \hat{a}_{u_k} u_k (\rho) \end{equation} \tag{ A.4 }$$

with

$$\begin{equation} \hat{a}_{u_k} = \int \mathrm{d} \rho \hat{a}(\rho) u^{*}_k (\rho). \end{equation} \tag{ A.5 }$$

A.2. Multimode homodyne detection

We proceed now with the modelization of the change of basis which the optical modes undergo as an effect of the MPHD, as reported in [18]. The beam splitter mixes the input modes with the local oscillator, yielding the modes

$$\begin{eqnarray} \hat{a}_{A,u_k} &=& \frac{1}{\sqrt{2}} \left( \hat{a}_{\mathrm{LO},f_k} + \hat{a}_{u_k} \right), \nonumber \\ \hat{a}_{B,u_k} &=& \frac{1}{\sqrt{2}} \left( \hat{a}_{\mathrm{LO},f_k} - \hat{a}_{u_k} \right) \end{eqnarray} \tag{ A.6 }$$

with

$$\begin{eqnarray} \hat{a}_A(\rho) &=& \sum_k \hat{a}_{A,u_k} u_k (\rho), \nonumber \\ \hat{a}_{A,u_k} &=& \int_{\rho \in \mathbb{R}} \hspace{-0.4cm} \mathrm{d} \rho \, u^{*}_k (\rho) \hat{a}_A(\rho). \end{eqnarray} \tag{ A.7 }$$

Then, light is propagated to the pixels. Let us introduce the output intensities on each pixel

$$\begin{eqnarray} \hat{i}_{A, i} &=& \int_{\rho \in S_{i}} \mathrm{d} \rho \, \hat{a}^\dagger_A(\rho) \hat{a}_A(\rho), \nonumber \\ \hat{i}_{B, i} &=& \int_{\rho \in S_{i}} \mathrm{d} \rho \, \hat{a}^\dagger_B(\rho) \hat{a}_B(\rho). \end{eqnarray} \tag{ A.8 }$$

Taking the differences of the intensity signals between pairs of corresponding pixels gives

$$\begin{eqnarray} \fl \hat{s}_{i } &=& \hat{i}_{A,i} - \hat{i}_{B,i} \nonumber \\\fl &=& \frac{1}{2} \sum_{k,l} \left[ (\hat{a}^\dagger_{\mathrm{LO},u_k} + \hat{a}^\dagger_{u_k}) (\hat{a}_{\mathrm{LO},u_l} + \hat{a}_{u_l}) - (\hat{a}^\dagger_{\mathrm{LO},u_k} - \hat{a}^\dagger_{u_k}) (\hat{a}_{\mathrm{LO},u_l} - \hat{a}_{u_l}) \right] \int_{\rho\in S_{i}} \hspace{-0.4cm} \mathrm{d} \rho \, u^{*}_k (\rho) u_l (\rho) \nonumber \\\fl &=& \sum_{k,l} \left( \hat{a}^\dagger_{\mathrm{LO},u_k}\hat{a}_{u_l} + \hat{a}^\dagger_{u_k}\hat{a}_{\mathrm{LO},u_l} \right) \int_{\rho\in S_{i}} \hspace{-0.4cm} \mathrm{d} \rho \, u^{*}_k (\rho) u_l (\rho). \end{eqnarray} \tag{ A.9 }$$

Now we assume that the local oscillator is in an intense coherent state in a certain mode that we will indicate by $\hat {a}_{\mathrm {LO},u_{\mathrm {LO}}}$ (u_LO is indeed one of the modes of the basis {u_i}), and we keep only the terms implying such a mode. Furthermore, we use $\hat {a}^\dagger _{\mathrm {LO},u_{\mathrm {LO}}} \simeq \langle \hat {a}^\dagger _{\mathrm {LO},u_{\mathrm {LO}}} \rangle = \alpha _0^{*}$, with $\alpha _0 \gg \langle \hat {a}^\dagger _{u_k} \rangle \, \forall k$. With this we obtain

$$\begin{eqnarray} \fl \hat{s}_{i } &\simeq& \sum_{l} \hat{a}^\dagger_{\mathrm{LO},u_{\mathrm{LO}}}\hat{a}_{u_l} \int_{\rho\in S_{i}} \hspace{-0.4cm} \mathrm{d} \rho \, u_{\mathrm{LO}}^{*} (\rho) u_l (\rho) + \sum_{k} \hat{a}^\dagger_{u_k}\hat{a}_{\mathrm{LO},u_{\mathrm{LO}}} \int_{\rho\in S_{i}} \hspace{-0.4cm} \mathrm{d} \rho \, u^{*}_k (\rho) u_{\mathrm{LO}} (\rho) \nonumber \\ \fl &\simeq& \sum_{k} \left( \alpha^{*}_0 \hat{a}_{u_k} \int_{\rho\in S_{i}} \hspace{-0.4cm} \mathrm{d} \rho \, u_{\mathrm{LO}}^{*} (\rho) u_k (\rho) + \alpha_0 \hat{a}^\dagger_{u_k} \int_{\rho\in S_{i}} \hspace{-0.4cm} \mathrm{d} \rho \, u^{*}_k (\rho) u_{\mathrm{LO}} (\rho) \right). \end{eqnarray} \tag{ A.10 }$$

According to the definition in equation (A.5) and equation (1) we can write

$$\begin{equation} \hat{a}_{v_{i}} = \int_{\rho \in \mathbb{R}} \mathrm{d} \rho \hat{a}(\rho) v_{i}^{*} (\rho) = \kappa_i \int_{\rho\in S_{i}} \hspace{-0.4cm} \mathrm{d} \rho \hat{a}(\rho) u_{\mathrm{LO}}^{*} (\rho). \end{equation} \tag{ A.11 }$$

Substituting equation (A.4) in equation (A.11) we obtain

$$\begin{equation} \hat{a}_{v_{i}} = \kappa_i \sum_j \int_{\rho\in S_{i}} \hspace{-0.4cm} \mathrm{d} \rho u_{\mathrm{LO}}^{*}(\rho) u_j(\rho) \hat{a}_{u_j} = \sum_j {U_{\rm{T}}}_{i,j} \hat{a}_{u_j}, \end{equation} \tag{ A.12 }$$

where we have defined the matrix realizing the change of bases

$$\begin{equation} U_{{\rm{T}} \, i, j} = \kappa_i \int_{\rho\in S_i} \hspace{-0.4cm} \mathrm{d} \rho u_{\mathrm{LO}}^{*}(\rho) u_j(\rho). \end{equation} \tag{ A.13 }$$

Remark that the matrix U_T is neither necessarily square nor has finite dimension. In terms of the pixel modes defined in equation (1), (A.10) becomes

$$\begin{eqnarray} \hat{s}_{i } &=& \left( \alpha'^{*}_0 \hat{a}_{v_{i}} + \alpha'_0 \hat{a}^\dagger_{v_{i}} \right) \nonumber \\ &=& |\alpha_0'| \frac{1}{2} \left( {\rm e}^{{\rm{i}} \varphi_0} (\hat{q}_{v_{i}} + \mathrm{i} \hat{p}_{v_{i}}) + {\rm e}^{-{\rm{i}} \varphi_0} (\hat{q}_{v_{i}} - \mathrm{i} \hat{p}_{v_{i}}) \right) \nonumber \\ &=& |\alpha_0'| (\hat{q}_{v_{i}} \cos \varphi_0 - \hat{p}_{v_{i}} \sin \varphi_0 ) = |\alpha_0'| \hat{q}^{\varphi_0}_{v_{i}}, \end{eqnarray} \tag{ A.14 }$$

where we have set α₀' = |α₀|e^−iφ₀/κ. We see that the measurement of the quadrature $\hat {p}$ on each output mode requires setting the phase of the local oscillator to φ₀ = 3π/2.

A.3. Signal recombination

The pixel modes can be digitally multiplied by real gains and recombined, leading to the modes $v''_{i} (\rho ) = \sum _{j} g^i_{j} v'_{j} (\rho )$. 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ⁱ_j. Then the corresponding annihilation operators are

$$\begin{equation} \hat{a}_{v''_{i}} = \sum_{j} O_{i,j} \hat{a}_{v'_{j}}, \end{equation} \tag{ A.15 }$$

and we can obtain information on a quadrature as

$$\begin{eqnarray} \hat{s}''_{i} &=& \sum_{j} O_{i,j} \hat{s}'_{j } = \left( \alpha'^{*}_0 \sum_{j} O_{i,j} \hat{a}_{v'_{j}} + \alpha_0' \sum_{j} O_{i,j} \hat{a}^\dagger_{v'_{j}} \right) \nonumber \\ &=& \left( \alpha'^{*}_0 \hat{a}_{v''_i} + \alpha_0' \hat{a}^\dagger_{v''_i} \right) \nonumber \\ &=& |\alpha_0'| (\hat{q}_{v''_{i}} \cos \varphi_0 - \hat{p}_{v''_{i}} \sin \varphi_0 ) \nonumber \\ &=& |\alpha_0'| (\hat{q}_{v'_{i}} \cos (\varphi_0 +\varphi_i) - \hat{p}_{v'_{i}} \sin (\varphi_0 + \varphi_i) ). \end{eqnarray} \tag{ A.16 }$$

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]

$$\begin{eqnarray} & U = X + \mathrm{i} Y, \nonumber \\ & Y - V X = 0, \nonumber \\ & X X^{\rm{T}} + Y Y^{\rm{T}} = \mathcal{I}, \nonumber \\ & X^{\rm{T}}Y = Y^{\rm{T}} X, \quad X Y^{\rm{T}} = Y X^{\rm{T}}. \end{eqnarray} \tag{ B.1 }$$

Setting A = XX^T the previous equations lead to the linear system

$$\begin{equation} V A V = \mathcal{I} - A. \end{equation} \tag{ B.2 }$$

For a three-mode linear cluster state, the adjacency matrix is given by

$$\begin{equation} V = \left( \begin{array}{@{}lll@{}} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \\ \end{array} \right). \end{equation} \tag{ B.3 }$$

The linear system in equation (B.2) yields in this case the solution

$$\begin{equation} A = \left( \begin{array}{@{}ccc@{}} \frac{2}{3} & 0 & -\frac{1}{3} \\ 0 & \frac{1}{3} & 0 \\ -\frac{1}{3} & 0 & \frac{2}{3} \end{array} \right). \end{equation} \tag{ B.4 }$$

From equation (B.4) we easily obtain the symmetric solution X_s such that X²_s = A, and the corresponding U_s and V_s satisfying equation (B.1)

$$\begin{eqnarray*} \fl X_{\rm{s}} = \left( \begin{array}{@{}ccc@{}} \frac{\left(3+\sqrt{3}\right)}{6} & 0 & \frac{\left(-3+\sqrt{3}\right)}{6} \\ 0 & \frac{1}{\sqrt{3}} & 0 \\ \frac{\left(-3+\sqrt{3}\right)}{6} & 0 & \frac{\left(3+\sqrt{3}\right)}{6} \end{array} \right); \quad Y_{\rm{s}} = \left(\begin{array}{@{}ccc@{}} 0 & \frac{1}{\sqrt{3}} & 0 \\ \frac{1}{\sqrt{3}} & 0 & \frac{1}{\sqrt{3}} \\ 0 & \frac{1}{\sqrt{3}} & 0 \end{array} \right); \quad U_{\rm{s}} = \left(\begin{array}{@{}ccc@{}} \frac{\left(3+\sqrt{3}\right)}{6} & \frac{\rm{i}}{\sqrt{3}} & \frac{\left(-3+\sqrt{3}\right)}{6} \\ \frac{\rm{i}}{\sqrt{3}} & \frac{1}{\sqrt{3}} & \frac{\rm{i}}{\sqrt{3}} \\ \frac{\left(-3+\sqrt{3}\right)}{6} & \frac{\rm{i}}{\sqrt{3}} & \frac{\left(3+\sqrt{3}\right)}{6} \end{array} \nonumber \right). \end{eqnarray*}$$

From this symmetric solution, all the other solutions can be obtained. Indeed, it is easy to see from equation (B.1) that if X_s is a solution, also $X = X_{\rm{s}} \mathcal {O}$ is a solution for any real orthogonal matrix $\mathcal {O}$. Hence, any matrix $U^{1,2,3}_{\mathrm {lin}} = X + \mathrm {i} Y = X_{\rm{s}} \mathcal {O} + \mathrm {i} V X_{\rm{s}} \mathcal {O} = (\mathcal {I} + \mathrm {i} V) X_{\rm{s}} \mathcal {O} = U_{\rm{s}} \mathcal {O}$ can be used to build a cluster state. With this we obtain $U^{1,2,3}_{\mathrm {lin}} = U_{\rm{s}} \mathcal {O}$. An example is provided in equation (17) of the main text.