Direct dialling of Haar random unitary matrices

Random unitary matrices find a number of applications in quantum information science, and are central to the recently defined boson sampling algorithm for photons in linear optics. We describe an operationally simple method to directly implement Haar random unitary matrices in optical circuits, with no requirement for prior or explicit matrix calculations. Our physically-motivated and compact representation directly maps independent probability density functions for parameters in Haar random unitary matrices, to optical circuit components. We go on to extend the results to the case of random unitaries for qubits.

Here we present a simple procedure for choosing a HRU on an optical circuit, implemented in terms of recursive decompositions of a unitary operator [21,22], by choosing values of the physical parameters independently from simple distributions. This procedure is useful for applications where the exact unitary description of the implemented circuit is less important than a guarantee that it is drawn from the correct distribution. While similar parameterisations exist in the mathematical literature [23], an operational application within linear optics is not widely appreciated. We extend the result to systems of qubits, by deriving a mapping between a linear-optical circuit on m = 2 n modes and a circuit operating on n qubits. Note that constructions for pseudo-HRUs on qudit and qubit systems are also available, serving as a general framework to investigate randomising operations in complex quantum many-body systems [24,25].
Choosing a HRU is analogous to choosing a random number from a uniform distribution, in that it should be unbiased. The probability of selecting a particular unitary matrix from some region in the space of all unitary matrices should be in direct proportion to the volume of the region as defined by the Haar measure, which is the unique translation-invariant measure on the space of uni- (c) A variable reflectivity beamsplitter can be effectively implemented with an MZI composed of a phase shift, θ between two 1/2 reflectivity beamsplitters. The MZI phases may be chosen directly from the distributions shown, with increasing bias towards θ = π. The integrated optics implementation with directional couplers results in a bias towards θ = 0. Line colours in (b) and (c) correspond to shading in (a).
As we will show, this approach is particularly relevant to recursive circuit decompositions that allow any unitary matrix to be implemented over m optical modes, by choosing appropriate values for beamsplitter reflectivities and phase shifters. We first consider the triangular scheme [16] shown in figures 1, 2a, which is a variant of that proposed by Reck et al. [21], and which represents an m × m unitary matrix U as a product of unitary operators labeled R n , U = m−1 i=0 R m−i [31]. Each block R n is chosen to transform the mode j = m − n, or the corresponding basis state, denoted by |Ψ m−n , into the n-dimensional unit vector |v n , over modes j = m − n to k = m − 1, i.e., This vector undergoes further transformations under subsequent blocks R i (n < i ≤ m) to finally produce |f n that occupies all m modes: Orthogonality between each of the m |f i vectors is guaranteed. Further, if the vector |v n is chosen from the unbiased distribution of unit vectors in n dimensions, the property of left invariance ensures that |f n does not become biased by the operation of the subsequent R i . The next and main task is therefore determining how an unbiased vector in n dimensions may be implemented with R n by choosing values for the linear optical components from which it is constructed, according to the expansion shown in figure 1(b). To achieve this, consider the complex Gaussian vector in n dimensions: where |Ψ i denotes the ith basis state and the z i are independent and identically distributed normal random variables with the probability density function (pdf), P zi (z) = 1/π exp − |z| 2 . This independence means that the pdf for v n is the product of the pdfs for these elements and depends only on the magnitude of the vector: We now show how this basis x, which we call the Cartesian basis, can be mapped to a new basis, r. We call the latter the physical basis, since, as we demonstrate below, it contains the variables corresponding directly to components in a physical realisation of the vector in linear optics. Namely, we denote by r 0 the power of the input to the given block R n , while the other r i stand for the reflectivities of beamsplitters (see also figure 1(b)). Next, combining the definition (1) and equation (3), we find where the matrix (in the Pauli basis) B(r) = √ rσ z + √ 1 − rσ x has been used to describe a beamsplitter as a function of its reflectivity. Finally, taking into account that x i = |z i | 2 , we find, We must show that the pdfs for the vector v n are separable in the physical basis so that the experimental parameters can be chosen independently. We also need to derive the form of the marginal distributions for the r i and φ i , from which experimental parameters must be chosen to obtain a Haar unitary. Since there is no functional dependence on the α i parameters in equation (4) and there is a one-to-one mapping α i → φ i , these phases can be chosen uniformly and independently from the interval [0, 2π).
Finding the pdfs for the beamsplitter reflectivities requires a more careful change in bases, using the Jacobian, The pre-factor from (4) is expressed in the r basis simply as exp(−r 0 ), so is trivially separable. We therefore consider the Jacobian matrix with For the four cases where the variable r n = 1 has been introduced for convenience.
We show that this form of matrix (lower Hessenberg) can always be transformed into a lower triangular matrixfor which the determinant is simply the product of the diagonal elements-by elementary operations, which do not change the absolute value of the determinant.
The first step is to perform a set of operations on the j = 0 column, c 0 , that set the upper n − 1 terms to zero, as follows: where k runs from 1 to m − 1, c k,0 are those quantities after k operations, and c k is the kth column. We can then place the column c 0 as the rightmost column, at which point the matrix is lower triangular. After the procedure is complete the element J 0,n−1 = 1 (see appendix for detailed proof).
The Jacobian determinant is given by multiplying the diagonal elements of the shifted matrix which are given by equation (16). The explicit form of the pdf in the r basis is which is manifestly separable.
It can be verified by explicit integration that this expression is appropriately normalised. Since the pdf is separable in this basis, the variables r i are independent, and can be chosen according to their marginal distributions, where, for clarity, r n,i denotes the reflectivity of the ith beamsplitter in the nth rotation, R n . We now integrate over r 0 to obtain a compact form for the pdf of n-dimensional unit vectors, and express the pdf for the full circuit of beamsplitters, P C (r) as the product of the pdfs for the diagonal arrays of beamsplitters: Recalling the beamsplitter transformation B(r) = √ rσ z + √ 1 − rσ x , we note that a variable reflectivity beamsplitter can be constructed as a Mach-Zehnder interferometer (MZI), from a variable phase shifter θ between two 1/2 reflectivity beamsplitters, H = B(1/2), to give B v (θ) = cos θ 2 I + i sin θ 2 σ x (up to a global phase). It is then useful to re-express the pdfs in terms of MZI phase shifts. The further change of variables, r = cos 2 θ 2 , gives In the setting of integrated optics, where beamsplitters are implemented with directional couplers on waveguides according to D(1/2) = 1 √ 2 (I + iσ x ) for reflectivity of 1/2, the pdfs are given by (24) but with sin and cos functions interchanged, i.e., In practical terms, an optical circuit composed of beamsplitters and variable phase shifters can directly dial up a configuration corresponding to a HRU, by choosing phase shifter values from the derived pdfs. A six mode example is given in figure 2.
We note that the version of the triangular scheme used here, in which each beamsplitter in every block R n couples two adjacent modes, differs from the original scheme [21], in which the first mode is consecutively coupled with modes 2, 3, ..., n. It is easy to check, however, that the mapping (7)-(9) can be applied to the original scheme as well, by replacing r with 1 − r and relabelling the output modes: {x 0 , x 1 , ..., x m−1 } → {x m−1 , x 0 , ..., x m−2 }. Such a change of variables does not affect the Jacobian determinant and the final expression for the reflectivity pdfs for the original scheme is obtained by replacing r with 1 − r in equation (21) (the phases are again chosen uniformly and independently from the interval [0, 2π)).
Next, we analyse the alternative decomposition of unitary matrices, proposed by Clements et al. [22], which corresponds to a rectangular mesh of beamsplitters and phase shifters, as shown in figure 3 for six modes. While the triangular scheme might be more resilient to loss and other errors in experiments in which only a small proportion of its (upper) input ports are accessed, the rectangular scheme is likely to be beneficial for experiments that involve accessing most of its inputs. The more compact rectangular scheme may also fit a greater number of modes on standard wafers used in the fabrication of integrated photonic circuits.

FIG. 3:
A 6 × 6 unitary operator implemented with a six-mode linear optical circuit according to the rectangular scheme. Here,rn,i stands for the reflectivity of the ith beamsplitter of the blockRn. Within eachRn (n = 2, ..., 6), we enumerate the beamsplitters according to the sequence s, which consists of n − 1 indices, with odd (even) numbers arranged in descending order and followed by even (odd) numbers, arranged in ascending order (see also main text). In this figure, the beamsplitters in the ith row mix the modes i − 1 and i (e.g., the beamsplitters of the third row,r4,3,r6,1 and r5,2, couple the modes 2 and 3).
The rectangular scheme obeys the blocked structure, analogous to the triangular scheme described above. That is, an m × m unitary matrix U can be written down as a product of blocksR n (hereafter the tilde refers to the decomposition of reference [22]). Each of these blocks, as previously, transforms the mode m − n into a vector over modes m−n up to m−1 (see also figure 1(b)). More precisely, for odd (even) m, U = m/2 j=1R 2j−1 i=0R m−2i ). Moreover, the mapping of equations (7)-(9) for the operator R n for the triangular scheme can be used forR n as well, by a simple substitution. Namely, for even (odd) m we replacer n,i by r n,s(i) , ∀i, where s is a sequence of n − 1 indices, with odd (even) numbers arranged in descending order and followed by even (odd) numbers, arranged in ascending order (e.g., for m = n = 6, we have s = {5, 3, 1, 2, 4}).
This substitution leaves the corresponding Jacobian determinant unaffected. Therefore, the pdfs given in equation (21) for reflectivities r n,i for the triangular scheme correspond to that of the rectangular scheme, but for r n,s(i) . In other words, P rn,i (r) =Pr n,s(i) (r). Subsequently, we find Alternatively, one can reorder the reflectivitiesr n,i according to the sequence s, as is done in figure 3, yielding P rn,i (r) =Pr n,i (r). Finally, the phases of the rectangular scheme, analogous to the triangular scheme, are chosen uniformly and independently from the interval [0, 2π).
Given the above parameterisation of Haar-random optical circuits, we now address the effects of errors, caused by imperfections in integrated photonics manufacturing. Before going into detail, we emphasize the important feature of our approach: due to the separability of the derived probability distributions, errors on a given component of the circuit do not propagate to other independently chosen parameters. In turn, a major source of individual errors is the imperfection of directional couplers. Used to implement the balanced beamsplitters of MZIs, directional couplers should ideally couple 1/2 of the light between waveguides so that each MZI can achieve the full reflectivity range. Fabrication tolerances, however, introduce errors and limitations on this range. Furthermore, we note that upper MZIs in the triangular scheme and central MZIs in the rectangular scheme are those most sensitive to errors, according to their polynomially growing pdfs.
Although schemes exist to minimise the effect of such errors and produce near perfect MZIs [27,28] it is worthwhile considering the influence of many small errors over a large circuit (this simple model is also useful to the qubit picture that we develop below). As an estimate to this effect we address the range of unitary operations covered by the proposed parameterisation, which we evaluate in terms of the coverage of the unitary space (see references [23,29] for more details), That is, cov m is the ratio between the reachable and full unitary spaces, assuming that the phase shifters cover their full range [0, 2π). The range of MZI reflectivities, in turn, is [|ε|, 1 − |ε|], where ε is a small random error.
In figure 4 we plot the coverage versus the circuit size m, which shows that for such moderate errors our parameterisation achieves high coverage rates. Since the pdfs for the triangular and rectangular schemes have been shown to be equivalent and independent, the coverage plotted in figure 4 is valid for both.
We now briefly show how these results may be extended to the scenario of quantum information processing with qubits, independently of any particular physical implementation. We suggest a mapping between a unitary operation on m = 2 p optical modes and the same unitary operation on p qubits, such that the pdfs derived above can be directly applied to systems of qubits. Labelling the optical modes as qubit basis states {|0...00 , |0...01 , |0...10 , ..., |1...11 } we map the optical beamsplitters and phase shifters to single qubit Hadamard gates, and n-qubit logic gates where the state of a single target qubit is transformed depending on the states of n − 1 control qubits. The target qubit operations are the NOT gate or qubit-flip operator, σ x , and the qubit-phase gates, Φ = e iφσz and Φ = σ x Φσ x . Each optical phase is mapped to a n-qubit Φ or Φ logic gate, and each optical beamsplitter is mapped as an MZI to a n-qubit Φ or Φ logic gate between two single qubit Hadamard gates on the respective target qubit.
The mapping can be understood with reference to figure 5, which explicitly details the case for 3 qubits and 8 optical modes and present a full circuit example for 2 qubits and 4 optical modes. The target for the n-qubit phase operations is always the final qubit; the conditioning configuration of the control qubits determines which element in the qubit space receives the phase. The addition of Hadamard operations on the final qubit allows the mapping of 1/2 reflectivity beamsplitters, and therefore MZIs, that operate between pairs of optical modes that differ in labelling only by the final bit. The further addition of n-qubit NOT gates allows MZIs to be mapped from pairs of optical modes that may differ in labelling by more than one bit. Any subset of the MZI operations may be implemented on qubits by simply omitting controlled phases where appropriate.
While not designed to be optimal, this one-to-one mapping between n-qubit phase gates and optical MZIs illustrates one way in which the distributions expressed in figure 2(c) may be used to directly implement a HRU on qubits.
We have presented a recipe to directly generate HRUs in linear optics with a proof that is straightforward in comparison to previous works [6,23]. Experimental conformation of these results can make use of tomography that does not require further optical circuitry [30]. The for- mula in its general form is applicable to boson sampling where Haar unitaries are required, and the extension to systems of qubits invites wider applications.
This work was completed shortly after the tragic death of one of the authors, Nick Russell. Those of us who knew him are grateful for his contributions to our work and our lives, which we continue to miss. We acknowledge support from the Engineering and Physical Sciences Research Council (EPSRC), the European Research Council (ERC), including QUCHIP (H2020-FETPROACT-3-2014: Quantum simulation), the U.S. Army Research Office (ARO) grant W911NF-14-1-0133. A.L. acknowledges support from an EPSRC early career fellowship.