Arbitrary optical wave evolution with Fourier transforms and phase masks

A large number of applications in classical and quantum photonics require the capability of implementing arbitrary linear unitary transformations on a set of optical modes. In a seminal work by Reck et al. it was shown how to build such multiport universal interferometers with a mesh of beam splitters and phase shifters, and this design became the basis for most experimental implementations in the last decades. However, the design of Reck et al. is difficult to scale up to a large number of modes, which would be required for many applications. Here we present a constructive proof that it is possible to realize a multiport universal interferometer on N modes with a succession of 6N Fourier transforms and 6N+1 phase masks, for any even integer N. Furthermore, we provide an algorithm to find the correct succesion of Fourier transforms and phase masks to realize a given arbitrary unitary transformation. Since Fourier transforms and phase masks are routinely implemented in several optical setups and they do not suffer from the scalability issues associated with building extensive meshes of beam splitters, we believe that our design can be useful for many applications in photonics.

Introduction.-Theability to arbitrarily transform an optical mode has applications spanning communications, imaging, and information processing.Such a 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) [2][3][4].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, we introduce a deterministic algorithm to design an arbitrary unitary transformation based on this method.
General variable control of modal unitary transformations will have applications across optics.For example, in fiber optic communications, spatial multiplexing will require transforming between an array of Gaussian profile modes from, say, a ribbon of single-mode fibers to the non-Gaussian spatial modes of one multimode fiber [5].Routing of optical channels requires a reconfigurable network, described by a unitary [6,7].Information processing with optical networks takes advantage of the ultralow latency and ultra-high clock speed of photonic waveguides.Capitalizing on this allows for one to, for instance, concatenate two unitary transformations to quickly multiply two matrices, a key ingredient in a neural network [8,9].Another area of application, imaging, is, at its heart, a unitary spatial transformation.General unitary transformations would enable novel imaging functionalities, such as cancelling the optical scattering that inhibits imaging through human tissue [10,11].Image process-ing, such as noise reduction, sharpening, or compression, could be done on the field itself, rather than the intensity recorded by the sensor [12].Turning to the area of quantum information, a generalization of the N = 2 qubit is a photon in a superposition of N modes, a qudit [13,14].In quantum cryptography, a protocol that uses qudits (and unitaries on them) rather than qubits improves the robustness to noise [15].Quantum computing logic gates, such as the controlled-NOT gate, can be implemented using unitary transformations on photonic waveguide modes [16].Moreover, random walks in waveguide-network transformations simulate a variety of quantum systems [17], such as molecules [18].The problem of sampling the output probability distribution when multiple photons traverse such networks is hard to simulate in a classical computer and hence it may be a viable path to achieve quantum supremacy with photonic devices, as proved in [19].The underlying reason is that sampling the bosonic statistics of the output photons is linked to the problemn of estimating the permanent of a large unitary matrix [20][21][22][23], which is #P-hard [24].These applications motivate why the field is spending considerable effort to develop controllable unitary transformations.
We now briefly outline these efforts and methods.While the first methods to create arbitrary transformations were developed during early radio and microwave engineering, Reck et al. introduced them to optics in a seminal paper in 1994 [1].One of the simplest tools available in optics is a phase shifter.However, by themselves modal phase-shifters are insufficient to build a general unitary transformation.In addition, one must use mode-mixing elements such as beamsplitters.Reck et al.

N
Fourier Transform e iϕ (2)   1 e iϕ (2)   2 e iϕ (2)   3 e iϕ (2) 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.
gave a prescription to implement any chosen unitary on an array of beam modes by using a triangle-shaped lattice of variable-reflectivity beamsplitters interleaved with phase shifters.However, the complexity of this method meant it was not demonstrated until an integrated optical implementation over twenty years later [25].Soon after, a more compact square lattice of beamsplitters and phase shifters was proposed and implemented by Clements et al. [26].Since the first implementation, a range of integrated optical platforms have hosted proofof-principle applications of these lattices [27][28][29].However, the fabrication and control complexity associated with this method has, so far, limited demonstrations of arbitrary unitaries to N = 6 waveguides [25].Before these integrated optical implementations, a different type of mixer was proposed and implemented, a lens or curved mirror.The latter elements enact an approximate Fourier transform of the spatial fielddistribution [30].Using these, a unitary is decomposed into a series of FTs interleaved with phase-shifters [2][3][4], as shown in Fig. 1.The phase shifters are varied to implement a given unitary, whereas the FTs do not change.An SLM of N pixels per side will inherently contain the N 2 control parameters required to implement a unitary for a spatial mode that varies along a row of pixels.Moreover, since they are based on television technology (e.g., 4K resolution) they currently have up to N 2 = 8 million pixels.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 [2,4].Fourier transforms can also be realized for other types of modes, for instance waveguides modes and spectral temporal modes, in an efficient manner [31][32][33].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) since it is more amenable to matrix algebra.While there is an existence proof showing that a unitary could be decomposed into a sequence of FTs alternating with phase-shifters [34], there is no prescription for doing so with a sequence of realizable length.In [35], a method was found to decompose an arbitrary matrix as a sequence of Fourier transforms and non-unitary diagonal matrices.However, their method is not adequate for linear optical setups, since their prescription makes use of non-unitary diagonal masks.Furthermore, the length scales as N 3 , which is very far from the optimal scaling L ∼ N .That said, an optimization algorithm to determine these phase shifts, 'wavefront matching', was recently introduced and experimentally validated [36].While practical, iterative optimization has a number of drawbacks for the FT method: 1.There is no guarantee of a solution nor its global optimality 2. It does not prescribe the design parameters, e.g., the required number of required FT-phase shift iterations, resolution (i.e., SLM pixel size), and range (i.e., number of pixels).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.'s deterministic algorithm, it is computationally slow.4. It does not provide physical insight into how to develop improved methods.Consequently, there is a need for the deterministic algorithm we introduce here.
In the first section, we map the Clements et al. lattice onto the FT method.That is, we decompose a layer of beamsplitters in the lattice into a short sequence of FTs and fixed phase shifts.We use this to adapt their deterministic algorithm to find the requisite variable phase shifts in the FT method.Consequently, we give an explicit prescription for how to design a unitary transformation with the FT method.While the FT is often numerically computed using the DFT, in the second section, we show that the DFT can also directly occur in optical systems.We show that, for example, modal propagation in a box waveguide is described by a DFT.In the appendix, we give a detailed derivation of our FT method decomposition.More broadly, the FT method is not restricted to optics, being directly applicable to many other setups such as neutral atoms in optical traps or phonon modes in ion chains.
Decomposition method.-Anylossless, noiseless, linear transformation on a closed system of N optical modes is described by a unitary matrix U ∈ U N (C).Reck et al. showed that any unitary transformation between optical modes can be implemented as a lattice of beam splitters [1].Such a lattice is also known as a multiport interferometer.A beam splitter is an optical element that mixes two modes i and j according to unitary matrix T (θ, φ) ∈ U N (C) parametrized by two angles := e iφ cos(θ) −sin(θ) e iφ sin(θ) cos(θ) .
(1) It acts as the identity matrix on all the other channels.An arbitrary beam splitter T ij (θ, φ) can be factorized in the following way where X represents a 50-50 beam splitter, i.e., X := Hence, one only needs controllable phase shifters and fixed 50-50 beam splitters to build the lattice of beam splitters designed by Reck et al.
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 = 1 √ N e i2πjk/N .Our design is built as a succession of Fourier transforms and phase masks: The phase masks D (i) i∈{0,..,L} are the only element in this setup that we can control.Each phase mask on N modes is described by a diagonal matrix parametrized by N angles, D (i) jk = e iα (i) j δ jk .Thus, it is clear that one needs at least N of them to simulate an arbitrary multiport interferometer such as the Reck scheme.
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, rather than the Reck et al. method we start from the decomposition in beam splitters given by Clements et al. in [26].Their design consists of a mesh of beam splitters arranged in N consecutive layers.The composition of the action of all the beam splitters in the mesh has the following form where T j (θ, φ) is a beam splitter mixing channels j and j + 1.As pointed out above in Eq.( 2), any beam splitter can be implemented with two 50-50 beam splitters and two phase shifters.Therefore, the design proposed in [26] can be understood as a succession of layers of 50-50 beam splitters and phase masks.The procedure to translate the lattice of beam splitters into a composition of phase masks and DFT's is schematically depicted in Figure (2).In a nutshell, our decomposition builds on ϕ (1)  1 , θ (1)   1 ϕ (1)  2 , θ (1)   2 ϕ (1)  3 , θ (1)   3 Any unitary matrix can be realized by a mesh of beam splitters, as described in [26].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 (6,7,8).
this by factoring each layer of 50-50 beam splitters in the mesh as a product of Fourier transforms and phase masks.Since the proof is constructive, it automatically gives a method to find the parameters in terms of U .Our decomposition needs only N 2 free controllable parameters, which is optimal.However, an improvement by a constant factor in the optical length (i.e. the total number of phase masks required) may still be possible.
Consider that we want to decompose a given unitary matrix U .In the decomposition displayed in Eq.( 4), each term of the form k ) represents a layer of beam splitters connecting each even channel 2j with the odd channel 2j + 1 (mod N ), whereas each term ) represents a layer of beam splitters connecting each even channel 2j with the odd channel 2j − 1 (mod N ).It will turn out to be more convenient for us if we relabel the indices as {0, 1, 2, 3, ...} → 0, N 2 , 1, N 2 + 1, ... .Then, we can 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 ∈ {N/2, ..., N − 1} with the channel j − N/2, whereas in the even layers each beam splitter connects the channel j ∈ {N/2, ..., N − 1} with the channel j − N/2 + 1 (mod N 2 ).On the one hand, the odd layers can be written as 2 are diagonal matrices.On the other hand, the even layers have the same structure as the odd layers after a cyclic shift of the first half of the channels.In other words, each even layer of beam splitters can be expressed as P T XΩ (i) 1 XΩ (i) 2 P , where 2 are also diagonal matrices, and P is just the permutation matrix given by Then, it follows that any unitary matrix admits the following decomposition At this point, we only have to find how to decompose the matrices X and P as a product of phase masks and Fourier transforms.In order to do this, we will find how to factorize them in 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 † and diagonal matrices.
Let us define the diagonal matrix G := I 0 0 iI and the circulant matrix . First, we note that X = GY G. Second, we observe that the permutation matrix P can be factorized as a product of three circulant matrices and four diagonal matrices in the following way where C is the cyclic shift matrix C = δ j,j+1 (mod N/2) of size N 2 × N 2 .Therefore, we have shown how to decompose any unitary matrix U as a product of diagonal and circulant matrices.If then we diagonalize the circulant matrices, we immediately obtain a factorization of U involving only F , F † and diagonal matrices.But the inverse of the DFT matrix is just F † = ΠF = F Π, where Π is the following permutation matrix Since ΠDΠ 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 where the terms B (i) , A (i) are given by where we made use of the notation The diagonal matrices E, H are defined as , and the diagonal matrix Γ(v) is defined as a function of a real vector v ∈ R N/2 : Finally, p : U N → U N is just the map p(U ) := ΠU Π.
Note that when applied on a diagonal matrix, it just inverts the order of the diagonal entries after the first one: p (diag(a 0 , a 1 , ..., a N −1 )) = diag(a 0 , a N −1 , ..., a 1 ).
We now summarize the procedure to create any unitary transformation U using phase-masks and DFTs.This procedure is based on using Eqs.(6,7) to express a unitary matrix according to the factorization in Eq.( 3).First, we permute the channels of the unitary transformation as described in the proof, which corresponds to computing the matrix U P = P T U P , with P being the permutation matrix defined in Eq.( 5).Then, we find the decomposition of U P as a lattice of beam splitters by the procedure described in Clements et al. [26].That is, we find the parameters (χ (i) , η (i) , θ (i) , φ (i) ) i={0,...,N/2−1} for each lattice layer i such that U P is factorized in the form given by Eq.( 4).The procedure for finding these parameters is explained in [26], 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 (χ (i) , η (i) , θ (i) , φ (i) ) 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 Γ(v) (labelled by the value of v = χ, η, θ, and φ), at the location of their application within one layer i of our method.All the control parameters are contained in the diagonal matrices Γ(v), whereas the rest of the diagonal matrices are fixed.In the end, the computation of all the diagonal matrices is quite efficient, as it only requires O(N 2 ) operations.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.
Physical Implementation of the DFT.-Our decomposition method requires the capability to optically perform the DFT.Although the standard continuous Fourier transform is routinely approximately performed in optical experiments with several setups, such as lenses [37], curved mirrors [38], or arrayed waveguide gratings [31], implementing the DFT by optical means is not a trivial problem.In waveguide systems a 'star coupler', sometimes called a N × N symmetric multiport or splitter, is sometimes said to perform a DFT-like operation in that any given input mode is transformed to a flat distribution of output modes.That is, the unitary matrix that describes the star coupler matches the magnitudes of the DFT matrix, |F jk |, but the phases will likely be incorrect.Since setting the magnitude only removes half of the free parameters in a unitary matrix, N 2 /2 parameters still must be adjusted to match a DFT.This cannot be accomplished by sandwiching the star coupler between two phase-masks since they only have 2N parameters together.Consequently, to implement our method there is a need for an optical DFT in waveguide systems.
Here we discuss one possible procedure to optically compute the DFT that is based on the phenomenon of self-imaging inside a multimode waveguide [39].The idea of using multimode intereference (MMI) couplers to realize the DFT is not new, and it was first proposed in [40].However, their prescription uses an 2N × 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 × N MMI coupler, using all N modes.As it only needs a planar waveguide and phase shifts, we believe that our proposal could be easily scaled to a large number of modes.Furthermore, it can be generalized to many other setups, such as neutral atoms in optical traps.
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 ψ n (x) = sin(k xn x), where k xn = π(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 2 xj /2nk 0 , where k 0 := 2π λ .
Consider now that at z = 0 we input a wavepacket f (x − x in j ) centered at x in j .For simplicity, let us assume that x j := j + 1 2 w N , for some integer j ∈ {0, ..., N − 1}.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 [39] that when the propagation length is set to be equal to z N = 2nk0 π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 (x − x out k ) define the output mode basis, where x out j = N − j − 1 2 w N , and the complex phase weights compose a unitary matrix, S jk := 1 √ N e iχ jk .In [39] these weights were shown to be It is straightforward to check that in fact the unitary matrix S is nothing else than the DFT matrix left and right multiplied by a diagonal matrix and a permutation matrix where the permutation matrix R and the diagonal matrix Θ are given by Eq.( 11) 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 figure (3).For an input field , where the coefficients of the output are related to the coefficients of the input by ã = Sa.
Consider that we want to implement an arbitrary unitary matrix U directly using such an MMI.Let us define the unitary matrix U R = RU R T .We can use the previous results to find its factorization, U R = D (0) L i=1 F D (i) .Then, writing F in terms of S using Eq.( 11), the unitary matrix U can be factorized as U = D(0) L i=1 S D(i) , where we have defined the new diagonal matrices D(0 Note that all matrices D(i) are indeed diagonal matrices, since R is just a permutation matrix and Θ 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.
This idea is not limited to optical modes in multimode waveguides.In fact, any physical system with confined modes of the form ψ n (x) = sin(k xn x) and a parabolic dispersion relation E j = 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 ϕ(x, t = 0) = N j=1 a j φ x − x in N −j+1 .Now, the state at any time is given by ϕ(x, t) = n c n e iEj t ψ n (x).For free propagation, the dispersion relation is parabolic.Consequentially, all the mathematical expressions are formally equivalent to the ones that describe multimode interference in a waveguide.In particular, one could apply this protocol to neutral atoms confined in an optical trap.
Conclusion.-We have given an explicit, analytical, and deterministic procedure to design an implementation of an arbitrary unitary transformation of dimension N using only discrete Fourier transforms and controllable phase-masks.The number of control parameters is the minimum possible, N 2 .For mode sorting or multiplexing, where the global phase of each output state is irrelevant, the required number of DFT-phase-mask layers needed is 6N .One additional phase-mask is required for a completely arbitrary unitary.Thus, the scaling of the number of layers with the dimension is linear and is optimal, up to an overall factor.We have also prescribed the first practical method to implement a DFT in integrated optics and even in systems outside optics, such as ion traps.The unitary matrix factorization we give could also be useful in quantum computation theory for decomposing quantum algorithms in terms of just two types of operations, the quantum Fourier transform and diagonal operators.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.

Figure 3 .
Figure3.Free propagation inside a multimode waveguide can be used to realize the DFT (in this case, a 6 × 6 DFT).When a wavepacket is injected in the waveguide, the output at a particular value of the propagation length is a superposition of N copies of the wavepacket weighted with complex phases.It can be shown that he output field is nothing else than the DFT of the input field modulo phase masks and permutation of the indices of the input and output modes.Therefore, a combination of free propagation inside multimode waveguides and controllable phase masks is enough to realize arbitrary unitary transformations.