1. Introduction

Multiport interferometers are a crucial technology for optical communication and information processing, both in classical and in quantum optics. Classical applications include mode (de)multiplexers for few-mode fibers [1,2], self-aligning coupling into fiber [3], and spatial-mode and polarisation converters [4]. On-chip multiport interferometers, consisting of an array of reconfigurable beam splitters (BSs) and phase shifters (PSs), are well suited for manipulation of photonic quantum states in e.g. quantum teleportation [5], quantum key distribution [6] or photonic qubit gates [7], due to their inherent phase stability, reconfigurability and ease of fabrication.

One particular quantum-optical task which multiport interferometers are well suited for is boson sampling [8]. The boson sampling task consists of sampling from the output photon number distribution of a large interferometer, which is fed with single-photon inputs. Since the first demonstrations [7,9–12], many advances have been made, by devising alternative sampling schemes that are easier to implement [13,14] and by improving the efficiency of single-photon sources [15]. A direct implementation of this task in quantum hardware outperforms simulations on a classical computer for a not unreasonable number of photons, making it a promising technique for an unambiguous demonstration of a quantum advantage.

The hardness proof of boson sampling requires that the unitary matrix that governs the mode transformation $a_{\text{out}}^{†}=U{a}_{\text{in}}^{†}$ is randomly chosen according to the Haar measure, and that the number of modes is much larger than the number of input photons. This has created interest in implementing large random unitary matrices in multiport interferometers [16], and the accuracy with which this can be done.

In this work, we study the implementation of unitary (i.e. lossless, coherent) N-to-N optical mode transformations, in planar multiport interferometers with realistic fabrication tolerances. We use a recently developed decomposition algorithm [17], that implements a unitary transformation in a square array of beam splitter-phase shifter pairs. It can be shown that this decomposition has superior loss tolerance to an older decomposition [18], which uses a triangular arrangement.

First, we generate many Haar-random matrices and study the resulting distribution of beam splitter and phase shifter settings. We find that the nonlinearity of the decomposition means that while matrices which are drawn according to the Haar measure are the natural generalization of a uniform distribution, the resulting distribution of beam splitter settings is highly nonunniform. In particular, at the center of the network, the distribution of beam splitter setitngs is stronlgy skewed towards low values. For moderate interferometer sizes (20 modes) and realistic errors in fabrication this has the effect that neither decomposition can implement any unitary transformation faithfully. Moreover, our results show that the allowable fabrication tolerances decrease with the size of the interferometer, meaning that any level of fabrication tolerance sets a limit on the size of a reconfigurable interferometer.

We also study techniques to mitigate this effect, by adding reduncancy to the system in the form of an additional layer of beam splitters and by numerically optimizing the beam splitter settings. We find that this for small networks, this technique restores fidelity of the required tranformations, for realistic imperfections. This result points the way to the study of the robust creation of random matrices in photonic networks.

2. Results

Figure 1 shows the problem under study. We start with a unitary transformation we need to implement. The decomposition algorithms translate the unitary matrix into a set of beam splitter (BS) reflectivities (R_k) and phase shifts (φ_k). These can then be implemented in a multiport interferometer, of which each node is a beam splitter and phase shifter pair (see inset).

 

Fig. 1 A unitary matrix can be implemented into a multiport interferometer via a mathematical decomposition. The interferometer consists of pairs of beam splitters and phase shifters (see inset). The square decomposition results in the structure of the interferometer shown.

Download Full Size | PDF

Our first goal is to understand what the implementation of random unitary matrices looks like in terms of reflectivities and phase shifts. To do this, we performed the square decomposition on a large ensemble, consisting of 5000 random unitary matrices of size 20-by-20 [19].

Figure 2(a) shows the highly nonuniform spatial distribution of average reflectivity. Each grayscale square in the figure represents a beam splitter at the same location in the underlying interferometer, through which light travels from left to right. The modes are labeled along the y-axes and the depth along the x-axes. The colour indicates the average reflectivity, which ranges from 0 to 0.5. It is surprising that the centre of the interferometer has low values of reflectivity. In fact, the majority of beam splitters have low reflectivity and the overall average is 0.18. Note that low reflectivity means most light is transmitted, and thus travels along diagonal lines across the interferometer. Similar results can be found for the Reck decomposition by using the expressions for reflectivity distributions presented in [16].

 

Fig. 2 Interferometers implementing Haar-random unitary matrices show a specific distribution of beam splitter reflectivities. a) shows the spatial distribution of the average reflectivity in a size 20 interferometer. b) Shows the underlying histograms for three regions in the interferometer, the first column, top row and centre. c) shows how the centre-of-interferometer histogram scales with the size of the interferometer. Note the change of scale in c).

Download Full Size | PDF

For Fig. 2(b), we have selected the three regions from the interferometer which are marked in subfigure a: the first column, top row and the interferometer centre, a square with sides of 20% the interferometer size. For each of these we show the distribution of reflectivities that underlies the average of Fig. 2(a), plotted in their corresponding colours. Most interesting is the distribution for the centre, which is peaked at low values and for which we found no values of the reflectivity larger than 0.4 in our sample.

Figure 2(c) shows that this effect becomes more pronounced as the size of the interferometer increases. The figure shows how the distribution of the centre of the interferometer changes with interferometer size. We have plotted the corresponding distribution for sizes 20, 50 and 100. The distribution becomes more sharply peaked at low values when increasing the size and the average reflectivity becomes lower. The distributions for the first column and top row do not change with size, thus the overall average reflectivity becomes lower as the interferometer size increases. From a similar analysis we found that the Reck decomposition also has this scaling property.

These results can be understood intuitively through the properties of a Haar-random unitary matrix. Given a matrix U that describes an interferometer, the amount of light that travels from input j to output i in a classical experiment is |U_i,_j|². For Haar-random unitaries, the mean magnitude of every element is the same, 〈|U_i,_j|²〉 = 1/N. There is only one path, however, that light can take from the top input to the bottom output, or vice versa. On these paths, which correspond to the diagonals of our square mesh, light is transmitted at each BS (transmission T = 1 − R): thus transmission has to be high and, correspondingly, reflectivity has to be low.

Summarizing our results thus far, we have found that when taking randomly chosen unitary matrices and implementing them in a particular arrangement of beam splitters and phase shifters, the resulting distribution is highly nonuniform. This result is at first sight surprising, given that the distribution of Haar-random matrices is the generalization of the uniform distribution, in the sense that it is translation invariant. The crucial point is that the mapping between the unitary matrix and the BS and PS causes the latter not to have this invariance.

We now introduce the problem of interferometer imperfections. In particular, we investigate one type of imperfection that stands out when implementing Haar-random unitary matrices in on-chip photonic networks: limited beam splitter tunability. Most reconfigurable realisations of multiport interferometers use Mach-Zehnder interferometers (MZIs) to implement variable beam splitters [7,10]. These interferometers contain two static 50:50 beam splitters. In practice, these beam splitters are not exactly 50:50, which means the MZI can generally not reflect or transmit all light. As shown above, low reflectivities are needed for the majority of MZIs in a large interferometer implementing random unitaries, thus this is problematic.

We quantify the error resulting from this limitation. Our procedure is to decompose a unitary matrix into a set of BS settings, implement the constraints, and tranform back to see how close the resulting unitary is to our initial unitary. First, we generate a Haar-random unitary and decompose it assuming a perfect interferometer. Next, we modeled the BS error by drawing the reflectivities of the static BSs, which are drawn from a normal distribution with standard deviation σ and mean 0.5, where we refer to σ as the fabrication error. We calculate calculated the minimum and maximum reflectivity of the corresponding MZI. We then compute the resulting unitary matrix, setting the BS values to the achievable value closest to their ideal values as given by the decomposition.

To quantify the degree of similarity between the unitary U_exp which is achieved by this process and the target unitary U_t, we evaluate the fidelity, which is defined as $F\left(U,U_{exp}\right)=\left(1/N\right)\text{Tr}\left(|U_t^{†}{U}_{\exp }|\right)$, where Tr(A) is the trace of A. This fidelity measure computes the single photon quantum state fidelities at the output of the interferometer, averaged over the photon inputs. It provides an upper bound for the amount of crosstalk via an inequality that bounds the Kolmogorov (trace) distance between the ideal and achieved distribution [20]. Since the fidelity runs from 0 to 1, at values close to 1 we typically plot the infidelity 1 − F (as in Fig. 3).

 

Fig. 3 The effect of unbalanced MZIs on the fidelity of the decompositions as a function of the size of the fabrication error. The triangular Reck (open symbols) decomposition and the square decomposition by Clements et al. filled symbols are used for various interferometer dimensions. a) us the fraction of the random unitaries that are affected by imperfections. b) the fidelity between the target and the effective unitary for the affected matrices when using our adapted version of the decompositions. The error bars show the standard deviation of all data points used in the average.

Download Full Size | PDF

Figure 3 shows us the effect of the imperfection when using this imperfect decomposition. We have performed the adapted decomposition while varying the fabrication error of static BSs, which is displayed on the x-axis. This we have done for various interferometer sizes up to size 50 and for both decompositions.

In Fig. 3(a), we show what fraction of the random unitaries is affected by the error. We have repeated the procedure described above 1000 times, and we record in which fraction of the cases we are hampered by the constraints caused by our imperfect physical beam splitters. We see that, for larger interferometer sizes, unbalanced MZIs affect the fidelity for even small fabrication error. This means that as MZI multiport interferometers grow in size, they are inevitably affected by the error at some point. The ratio is the same for both decompositions.

Figure 3(b) shows the average fidelity for those random unitary matrices that cannot be implemented perfectly. The y-axis shows one minus the fidelity, which means that a value of 0 implies the effective unitary is equal to the target. With current state-of-the-art fabrication tolerance (0.025, [21,22]), we are limited to 0.9995 fidelity when building a 50-mode interferometer. To relate this value to experiment, we compare the results of single photon experiments of the effective matrix to the target unitary. We define P^exp as the set of single-photon transition probabilities of this implementation and P as the same set for the target unitary. Then, for 0.9995 fidelity, $〈|P_i^{\exp }-P_i|〉/〈P_i〉=0.02$: probabilities are off by 2% on average, with maximum averaging 25%. The triangular interferometer is slightly less robust to these imperfections than the square interferometer.

To compare the reflectivity distribution of random unitary matrices to those of other interesting interferometer applications, we have also performed the new decomposition on the Fourier transform. The Fourier transform is used in several quantum algorithms, such as Shor’s algorithm [23], and is, like Boson Sampling, well-defined for any number of input modes. Moreover, it has the property that the transmitted intensity is equal to the average transmitted intensity for a Haar-random matrix.

Figure 4 shows the decomposition of a 20-by-20 Fourier matrix using the square decomposition. Similar to the decomposition of a Haar-random matrix, the resulting reflectivity distribution has low reflectivity on the diagonals and high values at the edges. However, low values of the reflectivity occur only exactly on the diagonal and not in the general center of the system. This means that the number of interferometers which is at risk of being affected by an imperfect beam splitter goes roughly as N² for a Haar-random unitary, but as N for a Fourier matrix. Fig. 4(b) shows that this is reflected in the fraction of times that the matrix fails to implement, given a particular error: Fourier matrices are significantly more robust to fabriction errors than random matrices. to obtain this result for Fourier matrices, we have generated 1000 random realizations of the physical BS reflectivities, and computed the constraints on the tuning range of each MZI arising from those.

 

Fig. 4 Decomposition of Fourier matrices. The top panel shows the decomposition of a 20-by-20 fourier matrix. The bottom panel shows the fraction of times a random matrox (black points) and a fourier matrix (red points) fails to implement, where we have averaged over many realizations of our model of the imperfections.

Download Full Size | PDF

This result has two implications. Firsts, it implies that Fourier matrices are in some sense ’easier’ to implement than Haar-random matrices: they are affected by beamsplitter imperfections, but since there are fewer instances of low reflectivity, the probability of succesfully implementing a Fourier transformation on a system with given fabrication errors is higher than for a Haar-random unitary. This means that while Fourier matrices have been used as benchmarks for boson sampling [24], they cannot be used as a benchmark the tunability of the interferometer. For Hadamard matrices, the opposite is true: since these require many of the beam splitters to be tuned to 0 reflectivity [25], implementing a Hadamard matrix fails for any level of beam splitter imperfection.

Second, it implies that our heuristic argument based on intensity transmission along the diagonals is incomplete: it holds for beam splitters which are exactly on the diagonals, but for those beam splitters which are slightly off-diagonal the phase relations must be considered as well.

3. Numerical optimization and redundancy

While both the triangular and square decompositions are unique, our strategy of dealing with imperfections in the interferometer locally is non-optimal. For cases where unity fidelity cannot be achieved, it is still possible to increase the fidelity by numerical optimization [26, 27]. Moreover, a natural strategy to mitigate the effect of imperfect components is to add redundancy to the network, in the form of additional beam splitters. This latter strategy is motivated by the observation that both the square and triangular decomposition contain the optimal number of optical elements, in the sense that for both, the number of free parameters in the BS-PS setting is equal to the number of free parameters in a unitary of the corresponding size. If we consider the BS-PS settings as a parameter space, the subspace of optimal BS-PS settings is therefore a point. When we add constraints, parts of that parameter space become inaccessible, and that point may not be in the accesible region. By adding redundancy, we make the system underdetermined, meaning that there are now many solutions which are equally satisfactory. Our goal is to find under which conditions, with high probabily, one of these solutions will lie in the accessible region of the parameter space.

To study these ideas, we took the square interferometer design and added an additional layer if beam splitters and phase shifters, respecting the staggered layout. Similar to the techniques used to obtain Fig. 3, we then produce many random realistic networks and Haar-random matrices.

We then numerically optimize our fidelity using the quasi-Newtonian optimization function fminunc in MATLAB, which implements a BFGS algorithm [28]. We apply the standard technique of using a function with unbounded domain but limited range to convert our constrained problem into an unconstrained one. As our initial guess for the numerical optimization, we implement the square decomposition in the first N by N beamsplitters, and set the rest to be as reflective as possible given our constraints. We also performed direct numerical optimization on a square set of beam splitters, using the same techniques.

Figure 5 shows the enhancement in fidelity which is obtained by this procedure, given by the reduction of the infidelity, expressed as a ratio to the infidelity obtained with the direct approach. The black points show the enhancement given by numerical optimization, the red points show the enhancement given by adding an additional layer and then optimizing. The inset shows the raw fidelities. We sampled a minimum of 100 combinations of matrices and sets of imperfections in each point.

 

Fig. 5 Numerical optimization and redundancy in a photonic network with realistic errors. We plot the fidelity enhancement, defined as (1 − F_before)/(1 − F_after)). The black points show the fidelity enhancement when only numerical optimization is applied, and the red points show enhancement when an additional layer of beam splitters is added, followed by numerical optimization. The inset shows the raw fidelities, where the black squares are the unoptimized values, the red circles are the optimized values for a square network, and the blue triangles are the optimized value for a network with one redunant row of beam splitters.

Download Full Size | PDF

For the case where we optimize the square network, we observe that the fidelity enhancement increases with network size. As the network gets bigger, there is more scope for numerical optimization. In the extreme case of one beam splitter, there is no scope for optimization at all.

Conversely, for the case where we add another layer, the enhancement decreases with network size. Again, this can be understood by considering the case of two modes. In that case, our scheme is identical to a ’double Mach-Zehnder interferometer’ (see [25,29] for recent discussion and [30] for the demonstration of such a concept), where the effect of unbalanced MZIs is largely circumvented by using two imperfect MZIs to implement one perfect variable BS. In this case, one expects complete fidelity. As we go to larger network sizes, the effect of a single additional layer becomes less noticable. We note that it is an interesting open problem how many additional layers are required as a function of network size to achieve constant fidelity, i.e. whether our redundancy approach is efficient in the number of BS.

Finally, we note that our optimization method is not efficient. The run-time of our numerical algorithm scales strongly with increasing network size. In the case of networks with inbuilt redundancy, the solver performs worse due to the fact that the unconstrained problem has more than one solution. We leave the question of an efficient algorithm to mitigate network imperfections as an open problem. One interesting approach would be to apply self-aligning techniques [1–4,25,30,31] to this problem

4. Discussion & conclusion

Arkhipov showed [32] a bound on the overall distance between the target unitary and the unitary achieved by a network of optical components, given accuracy with which each component is set to its target value. However, this approach does not take into account the fact that not all points in the parameter space of component settings may be achieved equally easily in experiment, nor does it take into account the fact that the required settings may themselves be a function of network size, inducing further scaling behaviour. We have demonstrated that as the network grows, the required settings tend to move towards parts of the parameter space which - with the widely used implementation of tunable beam splitters by MZI’s - are not easily accessible.

In conclusion, we showed that the reflectivities in a multiport interferometer implementing Haar-random unitary matrices are such that fidelities are severely limited by unbalanced Mach-Zehnder Interferometers. We showed that, using optimisation of the parameters, some fidelity can be regained. More importantly, we found that slightly increasing the depth of the interferometer can improve the fidelity even in the presence of considerable error. This approach may also prove useful to mitigate the effects of other types of fabrication imperfections, such as unbalanced loss. The next step is to find a closed-form or low overhead method of finding these settings for a realistic interferometer with added depth. With such a solution in hand, one can greatly increase the fidelity of future large multiport interferometers.

Note

During the preparation of this manuscript, closely related results were published in [33]. [33] gives analytical expressions for the distribution of beam splitter reflectivity at each site in both the triangular and square mesh, but does not consider redundancy.