1. Introduction

The ability to arbitrarily transform an optical mode has applications spanning communications, imaging, and information processing. General variable control of modal unitary transformations has numerous applications across optics, ranging from fiber optic communications [1–3] to information proccessing with optical networks [4] and imaging [5–7]. In the area of quantum information, arbitrary unitary transformations will be used in quantum cryptography [8], simulation of quantum systems with random walks [9,10], quantum information processing with qudits [11,12], quantum neural networks [13], quantum computing with photonic waveguide modes [14] or boson sampling [15].

An optical mode transformation that is lossless and linear is described by a unitary matrix $U$, mapping a basis of $N$ input modes onto a basis of $N$ output modes. Since any such matrix has $N^{2}$ free parameters, a method for its implementation must have at least $N^{2}$ controllable parameters, which is an experimentally challenging scaling. One implementation method is based on optical Fourier transforms (FT) [16–18]. In this paper, we show that only $N^{2}$ controllable parameters are needed to implement an arbitrary unitary transformation on $N$ modes using FTs. What is more, unlike previous methods, which were optimization based, we introduce a deterministic algorithm to design an arbitrary unitary transformation based on this method.

In a seminal paper [19], Reck et al. gave a prescription to implement any unitary on an array of modes by using a triangle-shaped lattice of beamsplitters interleaved with phase shifters. A more compact square lattice was proposed and implemented by Clements et al. [20]. Since these first works, implementations have been shown in a range of integrated optical platforms [21–23]. Photonic circuits with up to $64^{2}$ controllable parameters have been demonstrated [24], and a fully programmable arbitrary unitary acting on $N=12$ waveguides has been realized [25,26]. However, the fabrication and control complexity associated with this method makes it challenging to implement larger unitaries.

Before these integrated optical implementations, a different approach was based on a lens or curved mirror. These enact an approximate FT of the spatial field distribution [27]. Using these, a unitary is decomposed into a series of FTs interleaved with phase shifters [16–18], as shown in Fig. 1. The phase shifters are varied to implement a given unitary, whereas the FTs do not change. A series of experiments using multiple reflections from a curved mirror and a phase-shifting spatial light modulator array (SLM) successfully demonstrated a variety of unitary transformations [16,18]. Fourier transforms can also be realized for other types of modes, for instance waveguides modes and spectral temporal modes, in an efficient manner [28–30].

 

Fig. 1. In this work we show that any linear unitary transformation between $N$ channels can be implemented by means of a succession of $6N+1$ phase masks (diagonal operators) and $6N$ Fourier transforms.

Download Full Size | PDF

The FT method is the focus of this paper. In particular, we give a deterministic algorithm to find the requisite phase shifts in the FT method. Rather than using the full continuous FT, we use the discrete Fourier transform (DFT). While there is an existence proof showing that a unitary could be decomposed into a sequence of FTs alternating with phase shifters [31], there is no prescription for doing so with a sequence of realizable length. Previous work [32,33] has provided methods to decompose an arbitrary complex matrix as a sequence of Fourier transforms and non-unitary diagonal matrices. However, since their prescriptions require non-unitary diagonal matrices, they cannot be implemented using phase masks. That said, an optimization algorithm to determine these phase shifts, wavefront matching, was recently introduced and experimentally validated [34]. While practical, iterative optimization has a number of drawbacks for the FT method: 1. The output of the algorithm is just an approximate solution. 2. It does not prescribe the design parameters. Thus, it is unknown what is required to achieve a unitary of a given dimension, level of optical loss, or amount of error. 3. Relative to the Reck et al. deterministic algorithm, it is computationally slow. Consequently, there is a need for the deterministic algorithm we introduce here.

2. Decomposition method

Any lossless linear transformation on a closed system of $N$ optical modes is described by a unitary matrix $U\in \mathcal {U}_{N}(\mathbb {C})$. Reck et al. showed that any unitary transformation between optical modes can be implemented as a lattice of beam splitters [19], a multiport interrferometer. A beam splitter is an optical element that mixes two modes $i$ and $j$ according to unitary matrix $T(\theta,\phi )\in \mathcal {U}_{N}(\mathbb {C})$ parametrized by two angles $\phi,\theta \in [0,2\pi )$

It acts as the identity matrix on all the other channels. An arbitrary beam splitter $T_{ij}(\theta,\phi )$ can be factorized in the following way

Instead of a beamsplitter-based method, here we investigate a factorization method based on Fourier transforms. As a starting point, we consider the Discrete Fourier Transform (DFT), whose action is described by a unitary matrix whose elements are given by $F_{jk}=\frac {1}{\sqrt {N}}e^{i2\pi jk/N}$. Our design is built as a succession of Fourier transforms and phase masks, as depicted in Fig. 1:

The phase masks $\left \{ D^{(i)}\right \} _{i\in \left \{ 0,..,L\right \} }$ are the only element in this setup that we need to control. Each phase mask on $N$ modes is described by a diagonal matrix parametrized by $N$ angles, $D_{jk}^{(i)}=e^{i\alpha _{j}^{(i)}}\delta _{jk}$. Thus, one needs at least $N$ of them to construct an arbitray unitary with its $N^{2}$ parameters such as the Reck multiport interferometer.

We present a way to find a decomposition of an arbitrary unitary matrix in the form displayed in Eq. (3), consisting of $L+1=6N+1$ unitary diagonal matrices and $6N$ DFT matrices. In our factorization method, we start from the decomposition into beam splitters given in [20], i.e. a mesh of beam splitters arranged in $N$ consecutive layers:

 

Fig. 2. Any unitary matrix can be realized by a mesh of beam splitters, as described in [20]. Our decomposition method is based on replacing each layer of beam splitters by a succession of discrete Fourier transforms and phase masks, as schematically depicted here at the top. Each layer of beam splitters requires six Fourier transforms (grey rounded rectangles) and 6 phase-mask diagonal matrices (coloured rectangles). Only two diagonal matrices per layer (red and yellow rectangles) depend on the unitary matrix that is being implemented, while the rest (blue rectangles) are fixed. The general expression for the phase masks is given in Eqs. (7), (8), and (11).

Download Full Size | PDF

In the following, we will show how to decompose a given unitary $U$ by prescribing the seven distinct phase masks in Fig. 2 ($E,\Gamma (\phi ),H,\Gamma (\theta ),G,\Gamma (\eta$) and $\Gamma (\chi )$) that are required by each bilayer. In the decomposition displayed in Eq. (4), each term of the form $\prod _{k=1}^{N/2-1}T_{2k}(\chi _{k}^{(i)},\eta _{k}^{(i)})$ represents a layer of beam splitters connecting each even channel $2j$ with the odd channel $2j+1\ (\mathrm {mod}\ {N})$, whereas each term $\prod _{j=1}^{N/2}T_{2j-1}(\theta _{j}^{(i)},\phi _{j}^{(i)})$ represents a layer of beam splitters connecting each even channel $2j$ with the odd channel $2j-1\ (\mathrm {mod}\ {N})$. It is convenient to relabel the indices as $\left \{ 0,1,2,3,\ldots \right \} \rightarrow \left \{ 0,\frac {N}{2},1,\frac {N}{2}+1,\ldots \right \}$. We see that $U$ can also be decomposed as a succession of layers of beam splitters such that in the odd layers each beam splitter connects each channel $j\in \left \{ N/2,\ldots,N-1\right \}$ with the channel $j-N/2$, whereas in the even layers each beam splitter connects the channel $j\in \left \{ N/2,\ldots,N-1\right \}$ with the channel $j-N/2+1\ (\mathrm {mod}\ {\frac {N}{2}})$.

The odd layers can be written as $X\Xi {}_{1}^{(i)}X\Xi _{2}^{(i)}$, where $\Xi _{1}^{(i)},\Xi _{2}^{(i)}$ are diagonal matrices. The even layers have the same structure as the odd layers after a cyclic shift of the first half of the channels. Thus, each even layer can be expressed as $P^{T}X\Omega _{1}^{(i)}X\Omega _{2}^{(i)}P$, where $\Omega _{1}^{(i)},\Omega _{2}^{(i)}$ are also diagonal matrices, and $P$ is the permutation matrix given by

It follows that any unitary matrix admits the following decomposition

At this point, we still need to decompose $X$ and $P$ as a product of phase masks and Fourier transforms. We start by factorizing them into products of circulant and diagonal matrices. A circulant matrix is a matrix such that each row is obtained by applying a cyclic shift by one slot to the right to the previous row. Since any circulant matrix is diagonalized by the DFT matrix $F$, a product of circulant and diagonal matrices can always be re-expressed as a product involving only $F$, $F^{\dagger }$ and diagonal matrices.

Define the diagonal matrix $G:=\begin {pmatrix}I & 0\\ 0 & iI \end {pmatrix}$ and the circulant matrix $Y:=\frac {1}{\sqrt {2}}\begin {pmatrix}I & -iI\\ -iI & 1 \end {pmatrix}$. First, we note that $X=GYG$. Second, the permutation matrix $P$ can be factorized as a product of three circulant matrices and four diagonal matrices:

Since $\Pi D\Pi$ is diagonal whenever $D$ is diagonal, we can decompose $U$ using only $F$ and diagonal matrices.

In the end, diagonalizing all the circulant matrices we obtain the following expression

Finally, $p:\mathcal {U}_{N}\rightarrow \mathcal {U}_{N}$ is just the map $p(U):=\Pi U\Pi$. Note that when applied on a diagonal matrix, it just inverts the order of the diagonal entries after the first one: $p\left (diag(a_{0},a_{1},\ldots,a_{N-1})\right )=diag(a_{0},a_{N-1},\ldots,a_{1}).$

In summary, our procedure to create any unitary $U$ using phase masks and DFTs is based on using Eqs. (7) and (8) to express a unitary according to Eq. (3). First, we permute the channels of the unitary as described above, which corresponds to computing the matrix $U_{P}=P^{T}UP$, with $P$ being the permutation matrix of Eq. (5). Then, we find the decomposition of $U_{P}$ as a lattice of beam splitters by the procedure described in Clements et al. [20]. That is, we find the parameters $\left \{ (\chi ^{(i)},\eta ^{(i)},\theta ^{(i)},\phi ^{(i)})\right \} _{i=\left \{ 0,\ldots,N/2-1\right \} }$ for each lattice layer $i$ such that $U_{P}$ is factorized in the form of Eq. (4). The procedure for finding these parameters is explained in [20], but the general idea is to null, one by one, all the off-diagonal elements of $U_{P}$ by means of an appropriate succession of beam splitters. We then apply these parameters $\left \{ (\chi ^{(i)},\eta ^{(i)},\theta ^{(i)},\phi ^{(i)})\right \} _{i}$ as phase masks along with other fixed phase masks, all interleaved with DFTs, to replace layer $i$. In Fig. 2, we indicate all seven different diagonal matrices (e.g., phase masks), $E,$ $H$, $G$, $p(G)$, and $\Gamma (\mathbf {v})$ (labelled by the value of $\mathrm {\mathbf {v}}=\chi,\eta,\theta,$ and $\phi$), at the location of their application within one layer $i$ of our method. All the control parameters are contained in the diagonal matrices $\Gamma (\mathbf {v})$, whereas the rest of the diagonal matrices are fixed. In summary, applying the structure in Fig. 2 in place of each the $N/2$ beamsplitter lattice layers results in an implementation of an arbitrary unitary using only phase masks and Fourier transforms.

We should observe that fabrication imperfections can pose a serious challenge in the case of large photonic circuits. In any implementation of our scheme, some form of fine-tuning may be required in order to correct for such fabrication imperfections. Several proposals for efficient error correction already exist in the case of meshes of MZIs [35,36], but we leave their generalization to the FT method for future work. Nevertheless, we remark that a deterministic solution is always a good starting point for fine-tuning. Even a strategy based on parameter shifting may be enough if the deviation from the intended behaviour is small.

3. Implementation of the DFT

Our approach requires the capability to optically perform the DFT. While the standard continuous FT is routinely approximated in optical setups [28,37,38], implementing the DFT is nontrivial. There are some proposals to realize the DFT based on AWGs [39], star couplers [40,41] and other devices for orthogonal frequency division multiplexing [42,43]. These devices can couple each input mode to each output mode with equal strength. In order to realize a true DFT, one also has to match the complex phase of each entry of the transmission matrix, which is a non-trivial problem. Moreover, many of these architectures are complicated compound photonic devices themselves.

Here we discuss one possible procedure to optically compute the DFT that is based on the phenomenon of self-imaging inside a multimode waveguide [44]. The idea of using multimode intereference (MMI) couplers to realize the DFT is not new, and it was first proposed in [45]. However, their prescription uses an $2N \times 2N$ MMI coupler to output two copies of the $N$ dimensional DFT on half of the input modes, which is not amenable to our goal. Here we describe a method to implement the DFT on N modes with an $N\times N$ MMI coupler, using all $N$ modes. As it only needs a rectangular planar waveguide and phase shifts, we believe that our proposal could be easily scaled to a large number of modes. In fact, the idea of programming arbitrary unitaries in a rectangular dielectric has attracted recent interest [46]. Another previous work has also considered the possibility of implementing programmable unitaries as a cascade of MMIs and phase layers, but did not include an explicit prescription for doing so. [47].

Consider a planar waveguide of width $w$ and index of refraction $n$. We parametrize the transversal coordinate as $x$ and the longitudinal coordinate as $z$. Let us assume hard wall boundary conditions, so that it supports guided modes of the form $\psi _{n}(x)=\sin (k_{xn}x)$, where $k_{xn}=\frac {\pi (n+1)}{w}$. Furthermore, let us assume that the length of the waveguide is much larger than its width. Then, in the paraxial limit we can approximate $k_{zj}=nk_{0}-k_{xj}^{2}/2nk_{0}$, where $k_{0}:=\frac {2\pi }{\lambda }$.

Consider now that at $z=0$ we input a wavepacket $f(x-x_{j}^{in})$ centered at $x_{j}^{in}$. For simplicity, let us assume that $x_{j}:=\left (j+\frac {1}{2}\right )\frac {w}{N}$, for some integer $j\in \left \{ 0,\ldots,N-1\right \}$. This defines a vector basis for our target DFT matrix in terms of $N$ wavepacket modes. Under the assumptions listed above, it has been shown in [44] that when the propagation length is set to be equal to $z_{N}=\frac {2nk_{0}}{\pi N}w^{2}$, the output field is given by

In other words, the output field is a superposition of $N$ repetitions of the input wavepacket at $N$ distinct positions and weighted by complex phases. The wavepackets $f\left (x-x_{k}^{out}\right )$ define the output mode basis, where $x_{j}^{out}=\left (N-j-\frac {1}{2}\right )\frac {w}{N}$, and the complex phase weights compose a unitary matrix, $S_{jk}:=\frac {1}{\sqrt {N}}e^{i\chi _{jk}}$. In [44] these weights were shown to be

Equation (14) implies a possible optical implementation of the DFT. The setup would consist of a planar multimode waveguide with length $z_{N}$ coupled to $N$ input channels and $N$ output channels as in Fig. 3. For an input field $E_{in}=\sum _{j=0}^{N-1}a_{j}f\left (x-x_{j}^{in}\right )$, the output field is $E_{out}=\sum _{j=0}^{N-1}\tilde {a}_{j}f\left (x-x_{j}^{out}\right )$, where the coefficients of the output are related to the coefficients of the input by $\tilde {\mathbf {a}}=S\mathbf {a}$.

 

Fig. 3. We show that free propagation inside a multimode slab waveguide can be used to realize the DFT (here, a $6\times 6$ DFT). When a wavepacket is injected, the output at a suitable distance is a superposition of $N$ copies of the wavepacket weighted with complex phases, representing the DFT of the input field (modulo phase shifts and permutation of input and output modes). Thus, a combination of free propagation inside such waveguides and controllable phase masks is enough to realize arbitrary unitary transformations. In the plot, we show the phase and the amplitude of the electric field inside the MMI, $E(x,z)=A(x,z)e^{i\theta (x,z)}$ (the amplitude has been normalized so that the maximum is $A_\textrm {max}=1$. The phase is represented by the color, while the amplitude is represented by the pixel intensity.)

Download Full Size | PDF

Consider that we want to implement an arbitrary unitary matrix $U$ directly using such an MMI. Let us define the unitary matrix $U_{R}=RUR^{T}$. We can use the previous results to find its factorization, $U_{R}=D^{(0)}\prod _{i=1}^{L}FD^{(i)}$. Then, writing $F$ in terms of $S$ using Eq. (14), the unitary matrix $U$ can be factorized as $U=\tilde {D}^{(0)}\prod _{i=1}^{L}S\tilde {D}^{(i)}$, where we have defined the new diagonal matrices $\tilde {D}^{(0)}:=R^{T}D^{(0)}\Theta ^{*}R$, $\tilde {D}^{(i)}:=\Theta ^{*}R^{T}D^{(i)}R\Theta ^{*},\,i=1,\ldots,L-1$, and $\tilde {D}^{(L)}:=\Theta ^{*}R^{T}D^{(L)}R$. Note that all matrices $\tilde {D}^{(i)}$ are indeed diagonal matrices, since $R$ is just a permutation matrix and $\Theta$ is also diagonal. In summary, to use such an MMI in place of an exact DFT one simply needs to modify the $L+1=6N+1$ phase-masks in our method.

The simplicity of our scheme for implementing the DFT comes at the cost of a relatively large footprint. In our scheme, the length of the slab waveguide scales linearly with $N$. Hence the overall area of the fully programmable unitary with $O(N)$ DFTs scales like $O(N^{3})$. As a side remark, we observe that there are plenty of applications for an optical implementation of the DFT or other operations involving the composition of a few DFTs (such as a convolution) [41]. For those transformations that require only $O(1)$ DFTs, our scheme could still result in a very compact photonic circuit.

This idea is not limited to optical modes in multi- mode waveguides. In fact, any physical system with confined modes of the form $\psi _n(x) = sin(k_{xn} x)$ and a parabolic dispersion relation $\omega _j = \hbar \frac {k^{2}_j}{2m}$ can be used to realize the DFT. In this case, instead of propagating modes in a waveguide we consider a wavefunction that evolves inside a rectangular well according to the Schrödinger equation. We start with an input field of the form $\phi (x,t = 0) = \sum ^{N}_{j=1} a_j \phi (x - x^{\textrm{in}}_{N-j+1} )$. Now, the state at any time is given by $\phi (x,t) = \sum _n c_n e^{i \omega _j t} \psi _n(x)$. For free propagation, the dispersion relation is parabolic. Consequentially, all the mathematical expressions are formally equivalent to the ones that describe multimode interfer- ence in a waveguide. In particular, one could apply this protocol to neutral atoms confined in an optical trap. Indeed, it is now possible to realize arbitrary 1d potentials in optical traps for neutral atoms, and in particular one can realize a rectangular well potential [48].

4. Conclusion

With this result there is now an analytical, and deterministic procedure to design an implementation of an arbitrary unitary transformation. An even dimension unitary requires 6N DFTs and 6N+1 controllable phase masks. Thus, the scaling of the number of layers with the dimension is optimal, up to an overall constant factor. We have also described the first practical method to implement a DFT in integrated optics and even in systems outside optics, such as ion traps. We expect these results to be useful in a variety of classical and quantum information applications in photonics using various optical degrees of freedom including frequency-time, orbital angular momentum, and position-momentum.