Programmable holographic technique for implementing unitary and nonunitary transformations

Beyond the possibilities of linear transformations in polarization space, whose dimensionality is constrained by limited orthogonal states, we propose a technique for implementing both unitary and nonunitary transformations with higher dimensionality. Any high-dimensional matrix can be decomposed into a product of two processes realizable by utilizing spatial phase modulation and free-space propagation, in a simple, fixed, and scalable setup. Given that perfect power transmission for an arbitrary matrix may not be possible, the method is optimized to reach the theoretical best. Projected applications of the method described here include a means of restricting the infinite-dimensional Hilbert space to a finite-dimensional basis for information processing purposes, simultaneous multichannel optical routing, and a method of optical orbital angular momentum sorting and generation.


I. INTRODUCTION
Linear optics represents an outstanding platform for information processing and communication [1]. At the singlephoton level, it gives a direct access to many quantummechanical concepts such as single-qubit gates and projective measurements, and these transform smoothly to classical wave phenomena as the intensity increases. Twodimensional operations can be conveniently implemented in the polarization space, but to access higher dimensionality or additional quantum bits, one has to resort to other methods. Two main current approaches are the spatially encoded optical information [2] and the orbital angular momentum (OAM) degrees of freedom [3][4][5]. Due to their substantial differences, each has its own benefits and drawbacks.
Of the two methods, spatial encoding is more widely adopted. Two limitations preventing a broader adoption of OAM for high-dimensional operations are decreasing efficiency for manipulating larger OAM values and a nontrivial mapping of the potentially infinite range of values to finite-dimensional basis. An existing proposal for OAM implementation requires too many optical elements and is not fully scalable [6]. Another method has achieved a certain matrix transformation within a finite range of OAM values, but is hard to expand for other transformations [7].
The principal current method for processing spatially encoded optical information is the beam splitter or, more generally, an optical multiport [8,9]. These building blocks suffer, however, from low fixed dimensionality and a narrow parametric space. Typically, parametric behavior is introduced by using polarization-sensitive elements or phase delays, regulating, in either case, one of the internal continuous parameters of the setup per added element. To achieve a full control over the overall transformations applied to N inputs and outputs, one has to insert O(N 2 ) such elements, resulting in alignment challenges and complicated mapping between their control values and matrix elements of the desired transformation [2].
In this work, we focus on an N-dimensional optical system, in which an arbitrary linear mapping can be decomposed into a product of two processes realizable by utilizing spatial phase modulation and free-space propagation, making a full range of both unitary and nonunitary matrices, in a simple, fixed, and scalable setup. Using this technique, we reach, as a special case, a natural conversion between spatial encoding and OAM degrees of freedom, the latter constrained to N modular equivalence classes. Our method does not provide perfect power transmission for an arbitrary matrix, but we prove that within this realization, the efficiency is optimal. It reaches unit efficiency for some extremal cases that are of practical interest, and does not (in theory) fall below a universal lower bound depending only on N.

II. PRINCIPLE
To construct the desired high-dimensional Hilbert space, we take inspiration from the eigenstates of the conjugate pair of OAM and angle operators [10]. It is simplest to consider the OAM state |l and the angle state |ϕ , which are connected via Fourier transform [11]: A more careful analysis, however, which includes an angle operator, usually requires a finite basis of angular momentum states together with an appropriate limiting procedure [12,13]. For all practical purposes, including those considered here, there will be a natural cutoff at large positive and negative values of l. So while the OAM states would, in principle, give an infinite-dimensional basis, the conjugate "pure" angle states are impossible to implement as unbounded OAM states are necessarily required with nondecreasing amplitudes. Angle states model infinitely thin angular sections of the beam, which can only be approximated in practice.
To overcome these limitations, we work with the experimentally realizable "quasi"-angle states, in exchange for lowering the dimensionality from infinite to finite. A quasiangle state is considered as a Gaussian beam spot radially dislocated from the center with certain azimuthal angle. Here we will consider |ϕ n = u(r − r n ), n = 1,2, . . . ,N, where r = (x,y), r n = (r 0 cosϕ n ,r 0 sinϕ n ), with r 0 being the constant dislocation from the center, ϕ n = 2π (n − 1)/N the azimuthal angle, and u the complex amplitude for Gaussian beam at its waist, where w 0 r 0 is the size of the beam waist. The corresponding conjugate basis of quasi-OAM states is found via the discrete Fourier transform (DFT), where l m = m − 1.
Here, N is arbitrary but limited by the condition on separation of the Gaussian spots, and gives the dimension to both quasi-OAM and quasiangle bases. We provide an intuitive comparison of angle and OAM states as well as the quasiangle and quasi-OAM states with N = 6 in Fig. 1. For a given N , we can then represent an arbitrary N -dimensional state vector α = (α n ) using the quasiangle states for the basis as a ket, Let us consider a generic state of this form that is to be transformed by a fixed matrix T to another state, |β = N n=1 β n |ϕ n , (6) such that β = T α. In an operational formalism, it can be denoted by |β = T |α . In particular when T † T = T T † = I , the matrix is unitary, preserving the total intensity of the pattern.
It is important to emphasize that the set of Gaussian spots depicted in Fig. 1 are each a superposition of a large number of true OAM states [14], but these are close to being orthogonal and can be thought of as a set of quasiangle states. Considering, in this way, |ϕ n to be the dislocated Gaussian beam spot, the transformation could be simply implemented by splitting all the original Gaussian beams, which form the pattern of the state |α , into several unequal beams and recombining the latter to form the new pattern as the state |β . Thanks to the modulation convenience provided by the spatial light modulator (SLM), the aforementioned process can be implemented with a programmable hologram technique. We will discuss the details of such implementation in the following. Figure 2 shows the implementation principle with the simplest case of N = 2, allowing a planar projection for clarity. As the first step, by applying an appropriate diffraction hologram on SLM1, the Gaussian beam spots forming the state |α can be split into designated directions with distinct complex weights. These beams are then recombined as a physical implementation of a linear transformation. A diffraction pattern of exp(−ik x x) on a SLM mimics the behavior of a blazed grating. Specifically, the diffracted beam would be tilted with respect to the incident by a small angle of tan −1 (k x /k) (provided k x k in paraxial approximation) in the x direction. In order to split any Gaussian beam spot of state |α into several directions with certain complex weights, a superposition of gratinglike diffraction patterns as the hologram on SLM1 is required,

A. Splitting and recombining
where a mn stands for the transition amplitudes from the nth spot in state |α to the mth spot in state |β , and k mn = (k mn x ,k mn y ) represents tilting the beams with angles of tan −1 (k mn x /k) in x and tan −1 (k mn y /k) in y. With paraxial approximation, the components of k mn are calculated as where R m = [r 0 cos(ϕ m + π ),r 0 sin(ϕ m + π )]. It is important to notice that due to the 2f system adopted in the subsequent process, as shown in Fig. 2, the Gaussian beam spots on SLM2 are positioned in a centrally reversed orientation, which the extra rotation angle of π represents. The more the beam is tilted, the longer optical path it has during propagation from SLM1 to SLM2. Therefore, different phase compensation is necessary in each split direction, To prevent the crosstalk between different Gaussian beam spots on SLM2 as a result of natural divergence during propagation, a phase pattern encoding one thin lens per spot is superimposed on SLM1. In the paraxial approximation scheme, neglecting the constant phase shift, the transmission function of a focusing lens with a focal length of f is The binary function χ (r), valued one when |r| is smaller than a certain threshold and zero elsewhere, is present, in Eq. (7), to avoid the spatial coincidences of the grating patterns for different n. The threshold is determined by the size and separation of Gaussian beam spots on the SLM. Due to the superposition of different gratings being taken, the hologram F diff1 (r) requires in general modulations in both amplitude and phase. Fortunately, there have been many experimentally verified methods to achieve both amplitude and phase modulation on a phase-only SLM [15]. Some are implemented with two cascading SLMs, while others only depend on one SLM. In this work, only one phase-only SLM is adopted to mimic the amplitude and phase modulation, as discussed in Sec. III.
An amplitude-modulating hologram will inevitably be lossy, as the inequality |F diff1 (r)| 1 (11) must be satisfied for all (x,y). As soon as more than one term is present in the rightmost sum in Eq. (7), amplitude modulation is required, which inevitably comes with energy losses. The summation is upper bounded by for ∀n ∈ {1,2, . . . ,N}, and we consider the slightly stronger inequality as the condition for a physically realizable hologram.

B. Realigning and spatial filtering
After splitting and recombining, any reformed Gaussian beam spot is actually a superposition of tilted beams from different directions. However, when considering cascaded matrix transformations, it is essential that the incoming beams from different directions are realigned to the original longitudinal direction z. Once again, this realignment can be ensured by applying a corresponding diffraction hologram on SLM2. But any such realignment will inevitably result in several unwanted diffractions as a by-product. These could be cut off in the near field, requiring N carefully positioned apertures, or in the far field using a 2f focusing lens system with a central pinhole in the Fourier plane, which is employed in this work and shown in Fig. 2. Note that the latter reverses the image but this is compensated by the reversed design of SLM1 and the corresponding SLM2 appropriately.
In the simplistic example of Fig. 2, the lower Gaussian spot formed right before SLM2 is composed of two beams with a relative phase. In order to realign the upper of these two beams (green dashed line) to the z direction, a gratinglike diffraction pattern that tilts the beam in the +x direction is required. However, this pattern would also unselectively tilt the lower incoming beam (red solid line) to the +x direction as an unwanted side effect. Meanwhile, a flat 033827-3 pattern that simply transmits without changing direction is also required, as the lower incoming beam could be kept in its propagation direction. Once again, this flat pattern would keep the upper incoming beam in its direction as well, and would introduce another unwanted diffraction. Thus any incoming beam would be partially realigned, but also partially diffracted in some unwanted directions. The pinhole filtering in the 2f system is designed to block any components of this resulting field propagating in directions except for the z axis by transferring to the Fourier plane and back.
In a real application, it makes sense to integrate the transmission function of the lens L1 into SLM2 for convenience. Generally, the hologram on SLM2 is then described by where b mn stands for the complex weights for the partial phase ramp realigning a beam from the nth spot in state |α to the mth spot in state |β . Again the hologram will involve both amplitude and phase modulation and will be constrained by the condition Consequently, both the hologram itself and the spatial filtering can introduce additional loss at this stage. Finally, to ensure the Gaussian beam spots in state |β at their waists a focal length behind lens L2, an additional phase pattern encoding several thin lenses is also necessary, which is similar with those on SLM1. Therefore, we implement the pattern together with the transmission function of lens L2 on SLM3, As the lenses shown in Fig. 2 are encoded in the hologram patterns on SLM2 and SLM3, the whole setup to implement unitary and nonunitary transformations is finally shown in Fig. 3 (we take N = 5 for this illustration).

C. Efficiency
Taking the two above steps together, in order to obtain a generic transformation between N -dimensional states |α and |β described by the matrix T = (t mn ), we need to find matrices (a mn ) and (b mn ) such that a mn b mn = ηt mn , ∀m,n ∈ {1,2, . . . ,N}, where η is the efficiency, with conditions of Let us first assume for now that T only has nonzero elements. We note that starting with any solution to the problem, the inequalities holding for a mn and b mn can always be leveraged so that one set of the inequalities becomes equalities. Without loss of generality, we assume that Similarly, we are free to assume that a mn is real and positive for all indices as b mn can always absorb the necessary phase information. Under these assumptions, we can express b mn = η t mn a mn , ∀m,n ∈ {1,2, . . . ,N}.
Finding the most efficient realization of the matrix T corresponds to determining the maximum value of the parameter η without violating Eq. (19) or, equivalently, finding the lowest value of the right-hand side of Interestingly, these inequalities can also be converted to equalities while further decreasing the upper bound. Using pairwise modification of amplitudes of a mn such that Eq. (20) stays satisfied in an iterative algorithm, the optimum is reached exactly when all the left-hand sides of Eq. (22) This can easily be done using Lagrange's method, resulting in where u m and v n satisfy Or, donating S = ( √ |t mn |), which is actually an eigenvalue equation for the vector v: The particular choice of the largest eigenvalue can guarantee the existence of an eigenvector with all positive elements, and in turn the positivity of values of a mn , due to the Perron-Frobenius theorem [16]. Choosing the largest eigenvalue goes against the goal of maximizing η, but no other choice is guaranteed to comply with the prescribed constraints, and in fact cannot, due to the orthogonality of eigenvectors corresponding to different eigenspaces of the symmetric matrix S T S.
It follows that by identifying κ as the largest eigenvalue of S T S with a corresponding eigenvector v of positive elements, we can decompose the matrix T into an elementwise product of a mn and b mn as given in Eq. (17), with an efficiency factor of 1/κ. (We note that κ is equivalently the squared largest singular value of S and v is a corresponding right-singular vector.) By arguments of continuity, this result can be extended to matrices containing vanishing elements as well, and can be easily generalized even to rectangular matrices when mapping between unequally dimensional states.
As a special case that S T S is the identity, the matrices T can, in principle, be realized with unit efficiency. One might expect, based on arguments of energy conservation, that this should cover all unitary matrices. However, this condition is independent of unitarity (indeed, it states that a matrix formed of square roots of magnitudes of elements of T is unitary, rather than T itself) and for most unitary matrices the statement is not true. This is due to the limitations of what amplitude and phase modulation alone can achieve in splitting and recombining a beam, in particular due to the inevitable presence of unwanted diffractions, as illustrated in Fig. 2. Indeed, the only unitary matrices with η = 1 are generalized permutation matrices with arbitrary complex units, which correspond to mapping each input to exactly one output spot. The worst case is for the matrices in which all elements have the same magnitude, e.g., the DFT. For matrices of order N satisfying this criterion, the efficiency according to the above theory reaches a universal lower bound of η = N −3/2 .
Finally, in this discussion, we should note that while our scheme allows sequential applications of matrices, decomposing a matrix into a product can never make the total efficiency better. This is a direct consequence of an inequality valid for the largest singular values of two factors and of their product shown in [17].

III. SIMULATION OF EXAMPLE MATRIX TRANSFORMATIONS
The implementation of both amplitude and phase modulation is crucial for this work. In this simulation, we have adopted the checkerboard approach investigated by [18], which only relies on one SLM for simplicity, in exchange for lowering spatial modulation resolution. In this approach, two pixels on a SLM are paired to act as a superpixel. This principle is illustrated in Fig. 4. For any complex variable within the unit circle on the complex plane, one could always find two variables on the circle edge to represent it via its arithmetic average, Then, for any amplitude and phase modulation hologram, whose maximum magnitude would never exceed one, we can find a phase-only hologram to mimic it approximately. Without loss of generality, adjacent pixels on the SLM are paired in the x direction (horizontal) to be the superpixels, marked by red outlines in Fig. 4(b). Intuitively, the resolution of modulation in the x direction would be reduced by half, while that in the y direction (vertical) is unchanged.
All the electromagnetic fields in this work are considered in a paraxial approximation regime. The field propagation is described by utilizing the Fresnel-Huygens integral, with free-space wavelength of 1.  pixel size is 8 μm (referring to the PLUTO phase-only SLM by Holoeye Photonics, though transmissive SLM is assumed in this work for simplicity, instead of the actual reflective one). Note that an additional constant phase delay is applied on SLM3 to illustrate the simulation results more clearly. The dimension in the following examples is N = 5, which could be extended without difficulty to other integers. In what follows, we discuss the results of our simulation for several important cases.

A. Pauli shift and clock matrices
Ubiquitous in quantum physics, Pauli matrices are three 2 × 2 matrices that are complex, unitary, and Hermitian. Pauli matrices are extremely useful in quantum information and computation, doubling as a base for decomposing qubit states and a set of generators for describing operations on them. In discrete phase-space formulation, the matrix σ 1 acts as a displacement in position and σ 3 in momentum. Inspired by these properties, a possible extension of these two Pauli matrices to higher dimensions (also known as Weyl and Schwinger matrices) [19][20][21] is the shift matrix, and the clock matrix where ω = exp(i2π/N). Only some properties of Pauli matrices σ 1 and σ 3 are retained, and we note, in particular, that 1 and 3 are not Hermitian for N > 2 and there is no direct analog of σ 2 to form a representation of su(N ).
Thus, the shift matrix 1 corresponds to a translation operation, while the clock matrix corresponds to a delay operation, both in a quasiangle state basis. A concrete example is shown in Fig. 5(a) In this basis, 1 effectively becomes a clock matrix, while 3 is a shift matrix. We illustrate the two matrix transformations with l n = 1 in Fig. 5(b).
Generalized Pauli matrices are implemented perfectly with our proposal as each spot is just deflected to a new location and/or phase shifted. The theoretical efficiency as discussed above is 100%. In our numerical simulation, the efficiency is calculated as where |β id is the idealized output state and |β cal is the calculated state. In the above four cases, the calculated efficiency is 91.78 ± 2.21%, with the small drop due to the rasterization effects in the phase pattern on the SLM and the field itself.

B. Some sparse matrices
A permutation matrix is a further generalization of Pauli matrices that is unitary and can also be achieved with a theoretical efficiency of 100%. We illustrate this with a specific permutation matrix as The calculated results are shown in Fig. 6(a). Comparing state |β with state |α , the first and second spots (counting from the positive x semiaxis counterclockwise) are exchanged. The same is with the fourth and fifth spots, while the third spot is retained. The whole process fulfills the matrix transformation implied by T 1 and the calculated efficiency is 91.99%.
The efficiency would be worse when more than one element per row or column of the transformation matrix is nonzero. It can be expected to stay moderately high as long as the matrix is sufficiently sparse, or close to a sparse matrix. An example inspired by classical random walk on a line could be the tridiagonal Toeplitz matrix The calculated results are shown in Fig. 6(b). Though this requires multiplexing in most inputs and outputs and is not unitary, it was still achieved in our simulation with calculated efficiency of 52.78%, compared to the theoretical efficiency of 1/ √ 3 ≈ 57.74%. Furthermore, we sampled fields at two key locations between states |β and |α to give more information about the whole process. One is right before SLM1 (second column in Fig. 6), where the original Gaussian beam spots are diverging due to free-space propagation. The other is right before SLM2 (third column in Fig. 6), where one could see the amplitude pattern of beams from different directions.

C. DFT matrix
To investigate a more involved case of a matrix without any zero elements, we illustrate the potential for practical use of the theory with a matrix describing DFT, As indicated by Eq. (4), the Fourier relationship between the quasiangle and quasi-OAM states could be equivalently given as |l n = W −1 |ϕ n , ∀n ∈ {1,2, . . . ,N} |ϕ n = W |l n , ∀n ∈ {1,2, . . . ,N}.
In particular, the second equality indicates that each quasi-OAM state would be projected to one certain Gaussian beam spot after the DFT matrix transformation. This represents As there are N quasi-OAM states whose total azimuthal phase variation is given by 2πl n , each of them approximates a vortex beam with topological charge of l ≡ l n modN . In fact, this correspondence is manifested if the quasi-OAM state is simply spatially filtered from the respective vortex beam. It is important to mention that vortex beams with topological charges of l n + mN(m ∈ Z) become indistinguishable in the projection with the condition of w 0 r 0 . In our simulation shown in Fig. 7, Laguerre-Gaussian beams with topological charge of l [22] are spatially filtered by a fixed mask to reach the quasi-OAM states as state |α . The radius of the pinholes on this fixed mask is 1.5 times the size of beam waist w 0 . Then the latter are projected to the corresponding Gaussian beam spots applying the DFT matrix, which clearly sorts them according to their topological charges. The field before spatial filtering is also illustrated, where one could see the unwanted diffractions after focusing by SLM2. The pinhole is schematically illustrated by a white dashed circle in the figure, enabling the field inside it to pass through.
The efficiency is calculated to be 8.41 ± 1.10% (over the five studied cases), compared to the theoretical efficiency of 5 −3/2 ≈ 8.94%. Some calculated efficiencies even surpass the theoretical expectation due to the fact that state |α is made up of binary cutoff beam spots with slightly larger size, rather than ideal Gaussian beams as assumed in the theoretical model.
Conversely, the inverse DFT matrix could be viewed as a generator for quasi-OAM states and even OAM states themselves [23], if one starts from the quasiangle states.  By introducing the DFT matrix, we have illustrated a topological charge measurement method of a vortex beam. Furthermore, the vortex beam could also be generated with the inverse DFT matrix.

IV. DISCUSSION
By applying the experimentally realizable and highdimensional quasiangle state, the programmable holographic technique presented here successfully demonstrates the implementation of unitary and nonunitary transformations. The efficiency of the transformation depends on the matrix applied and has a universal lower bound of η = N −3/2 . Though only square transformation matrices are investigated in this work, it is important to state that our proposed scheme could be readily extended for any rectangular matrix. Furthermore, though an angular pattern is considered here as a choice inspired by optics, the distribution of Gaussian beam spots in the plane is largely irrelevant and could also be tailored for a particular purpose.
It is worthwhile to mention that all the optical elements need to be interferometrically stabilized with a high precision in experiment. The SLMs may be further replaced by custom refractive elements for a fixed transformation to obtain a higher efficiency [24]. A possible application of the technique would be an integrated implementation with nanophotonics technologies to further reduce the size of the setup and improve the dynamic characteristics.
Inspired by the angle and OAM states, the quasiangle state expands the dimensionality remarkably when compared to two polarization states and paves the way for high-dimensional Hilbert space manipulations. The measurement of OAM topological charge is demonstrated in detail, indicating a metrology method for optical vortices. Our proposed scheme also has potential for optical free-space communication involving OAM multiplexing and optical computing.