© The Author(s) 2022. Published by China Science Publishing & Media Ltd. and EDP Sciences.

INTRODUCTION

Unitary operators are powerful tools in quantum information processing for both investigating non-classical phenomena and exploring quantum computational resources [1]. Great efforts have been made to apply universal linear operation protocols on high-dimensional optical bases, for extending information capacity and enhancing computing versatility [2–4]. Among them, the orbital angular momentum (OAM) [5] modes provide abundant high-dimensional quantum entanglement resources [6–8], while the spiral phase distributions of OAM states have been employed for quantum imaging [9] and remote sensing [10] techniques.

It has been known that any lossless optical setup can be described by a unitary operator [11]. Thus, generally speaking, the generating [12,13], multiplexing [14–16], sorting [17–19] and manipulating [20,21] OAM modes can be considered certain unitary transformations. Here, we are focused on the arbitrary unitary operator within the OAM-domain. It means that both the input and output states are encoded on a series of optical OAM modes, namely the basis. In our scenario, the weight coefficients of input state would be transformed to those of the output with a unitary operator while the basis keeps. Such a setup or device to perform the unitary operation is quite crucial, since it is the fundamental tool to process the high-dimensional information encoded on OAM modes. However, as analyzed in [21], although there have been great achievements about generating and measuring the OAM modes, how to design an on-demand lossless implementation for an arbitrary unitary operator within OAM-domain has not been well-solved yet. The requirement of closed operation has made it tough to demonstrate unitary operations within OAM-domain, especially with the demand for near-unity success probability at the same time. With the technique of OAM parity sorter [17,18], the four-dimensional Pauli $\widehat{X}$ gate and $\widehat{Z}$ gate have been reported [20]. However, it would be complicated to construct arbitrary unitary operators based on $\widehat{X}$ gates and $\widehat{Z}$ gates. Recently, with multi-plane light conversion (MPLC) [21,22] method, cyclic and Fourier transformations are reported with dimensionality up to five. Though a broad variety of unitary operations can be performed with MPLC approach, the transformation efficiency would decrease significantly when the number of OAM modes exceeds four according to the reports [21], thus near-unity efficiency has not been achieved for higher dimensionalities. For techniques with spin-orbital coupling [23,24], only two OAM modes can be utilized due to the two-dimensional encoding space of spin angular momentum so that such scheme is limited and not scalable. Thus, to achieve an OAM-encoded unitary operation scheme with universality, scalability, and near-unity success probability simultaneously is still an open question.

In this work, inspired by the Fourier-diagonal decomposition theorem [25,26] and the demonstration of frequency-domain unitary operation protocols [27,28], we extend this technique to optical encoding basis of OAM domain. A scalable and intrinsic near-lossless scheme for programmable unitary gates based on matrix decomposition [25] is proposed and demonstrated with dimensionality of $3\times 3$. Moreover, the parallelism of our proposed scheme is also presented with two 3×3 matrices. Limited by the experimental conditions, only numerical simulations are presented for arbitrary OAM-encoded unitary matrices. It should be mentioned that our proposed scheme could also be applied on other two photonic encoding bases that belong to a Fourier transform pair. Thus, as an alternative to verify our proposal, proof of principle experiments have been performed on path domain, in which average fidelity value of $0.97\pm 0.02$ and transformation efficiency of $0.96\pm 0.03$ are evaluated through 80 experimental results with matrix dimensionality of$3\times 3$.

THEORY

Without loss of generality, our target is to implement an arbitrary unitary gate $U$ acting on optical states, which could be expressed as

$\begin{array}{c}|\beta \rangle =U|\alpha \rangle ,\end{array}$(1)

where $|\alpha \rangle$ and $|\beta \rangle$ denote the input and output states. Usually $|\alpha \rangle$ and $|\beta \rangle$ are encoded on a series of optical modes, namely the basis, which could be any degree of freedom, including optical path, frequency, temporal-mode, and transverse-spatial-mode. It should be mentioned that here, the bases of both input and output states are considered the same. Thus, with a unitary operator of $U$, the weight coefficients of $|\alpha \rangle$ would be transformed to those of $|\beta \rangle$ while the basis is not changed. For most practical applications, $U$ is a finite-dimensional operator and can be mathematically expressed as a unitary matrix with respect to the utilized encoding basis. It has been presented [25] that an arbitrary matrix of dimensionality $N$ can be mathematically decomposed into the production of diagonal matrices and cyclic matrices alternately, where the total number of matrices is M, or equivalently saying, $M$ diagonal matrices $D$ spaced by discrete Fourier transformation (DFT) matrix $F$ and its inverse matrix, which is also its Hermitian conjugate $F^{†}$:

$\begin{array}{c}{U}_{N\times N}=F{D}_M{F}^{†}\cdots F{D}_5{F}^{†}{D}_{4}F{D}_3{F}^{†}{D}_{2}F{D}_1{F}^{†}\\ =(F{D}_M{F}^{†})\cdots (F{D}_5{F}^{†})D_4(F{D}_3{F}^{†})D_2(F{D}_1{F}^{†}),\end{array}$(2)

where $M\le 2N-1$ [26] denotes the number of diagonal matrices required for the decomposition and the actual value would vary according to the target matrices. Specifically, to establish a unitary operator $U$, the diagonal matrices $D$ in eq. (2) would also be unitary. To clearly present our idea, the right-hand term in the first line of eq. (2) is modified as the form of the second one. Particularly, the terms in brackets are diagonal matrices sandwiched by DFT and inverse DFT matrices. Actually, they are also diagonal matrices but acting on a transformed domain, which is determined by the DFT. Thus, eq. (2) could be understood as that a unitary operator can be decomposed into a series of diagonal matrices in two domains, which are linked with Fourier transformation. A diagonal matrix, which is also referred to as a spectral shaper, would not lead to cross-mode effects. Actually, the required cross-components for unitary operations are achieved by DFT matrices in eq. (2). Thus, if two photonic degree-of-freedoms (DOFs) are Fourier transform pair, any unitary matrix can be implemented by performing diagonal matrix operations on these two DOFs in succession for the investigated photonic qudits.

It should be mentioned that for the Fourier terms in eq. (2), it is not required to actually construct the corresponding matrix operation. These Fourier and inverse Fourier operators are naturally inserted in eq. (2) when the spectral shaping is alternatively performed on two DOFs linked by Fourier transform. Actually, there are some Fourier transform pairs, which have been well investigated in photonic domain, including linear basis (x-basis) and transverse wave vector basis (k-basis) [29,30], OAM and azimuthal angle [31], and frequency and time-bin. Specifically, arbitrary unitary operators with OAM encoded states can be implemented with diagonal matrices in the OAM and azimuthal angle domain.

With well-defined OAM basis, an arbitrary state vector $|\alpha \rangle$ can be expanded as

$\begin{array}{c}|\alpha \rangle ={\displaystyle \sum _l{c}_l|l\rangle },\end{array}$(3)

where $\left\{c_l\right\}$ represent OAM spectrum coefficients and $|l\rangle$ denotes the $l\text{th}$ eigenstate with OAM of $l\hslash$ carried per photon. As mentioned above, in our scenario, the spectrum coefficients of $\left\{c_l\right\}$ would be transformed with a unitary operator of U while the basis of $|l\rangle$ keeps. To perform the decomposition of eq. (2), both the OAM domain and corresponding Fourier transformation should be employed. It has been investigated that transverse OAM basis $|l\rangle \langle l|$ and angle basis $|\phi \rangle \langle \phi |$ are connected via Fourier transformation [31]:

$\begin{array}{c}\begin{array}{l}|l\rangle ={\displaystyle \frac{1}{\sqrt{2\text{π}}}}{\displaystyle {\int }_{-\text{π}}^{\text{π}}\exp (-\text{i}l\phi )}|\phi \rangle \text{d}\phi ,\\ |\phi \rangle ={\displaystyle \frac{1}{\sqrt{2\text{π}}}}{\displaystyle \sum _{l=-\infty }^{\text{+}\infty }\exp (\text{i}l\phi )}|l\rangle .\end{array}\end{array}$(4)

Thus, it is feasible to implement well-designed diagonal matrices within the OAM-domain and angle-domain alternately to achieve any target OAM-domain matrix operator with decomposition form shown in eq. (2). The question is how to implement the diagonal matrices for both OAM and angle domains. Figure 1 shows an example of OAM-domain unitary operator with the three-layer scheme corresponding to the case of $M=3$ in eq. (2). There are two layers to implement diagonal matrices in angle domain and one layer for OAM domain. As shown later, this scheme could perform most unitary operations with dimensionality of $3\times 3$, and some particular operations of $4\times 4$ and $6\times 6$.

Figure 1 The scheme of OAM-domain unitary operation with three-layer.

In the angle-domain, a diagonal matrix can be readily achieved by a single wave front modulation element such as a spatial light modulator (SLM) [32] or metasurfaces [33]. Here, the SLM is considered and noted as angle modulator in Figure 1. A typical angle modulation function $f(\phi )$ can be achieved by the pattern ${\text{Φ}}_{\text{ANG}}(r,\phi )$ settled on SLM:

$\begin{array}{c}{\Phi }_{\text{ANG}}(r,\phi )=f(\phi ),\end{array}$(5)

where $r$ and $\phi$ indicate the radius and azimuthal angle under polar coordinates in the transverse plane according to beam propagation direction, respectively. Thus, the only question is how to realize the required diagonal matrices in OAM-domain, which is also referred to as an OAM spectral shaper [34].

Here, a straight-forward method is considered to achieve a programmable and arbitrary OAM spectral shaper. With the help of OAM mode-sorter based on Cartesian to log-polar coordinate transformation [35], each OAM eigenstate can be mapped to different position in the focal plane. As shown in Figure 1, two static optical elements ${\text{Φ}}_1(x,y)$ and ${\text{Φ}}_2(u,v)$ are employed to perform the mode sorter [35,36]:

$\begin{array}{c}\begin{array}{l}{\Phi }_1(x,y)=\frac{2\text{π}a}{\lambda f}\left[y\text{arctan}\left(\frac{y}{x}\right)-x\text{ln}\left(\frac{\sqrt{x^2+y^2}}{b}\right)+x\right],\\ {\Phi }_2(u,v)=-\frac{2\text{π}ab}{\lambda f}\text{exp}\left(-\frac{u}{a}\right)\text{cos}\left(\frac{v}{a}\right),\end{array}\end{array}$(6)

where $a$ and $b$ are adjustable parameters, which can be designed according to the beam width of the exact transverse optical field. After passing through these two static optical elements together with an auxiliary lens of focal length $f$, the $l\text{th}$ OAM eigenstate is sorted to a single spot, in which the center coordinate $x$ in the focal plane reads,

$\begin{array}{c}{C}_x(l)=\frac{f\lambda }{2\text{π}a}l,\end{array}$(7)

while the center coordinate $y$ is only determined by parameter $b$ in eq. (6) and irrelevant to $l$. According to eq. (7), different OAM modes are spatially separated in the focal plane of the auxiliary lens after the mode sorter consists of ${\text{Φ}}_1$ and ${\text{Φ}}_2$ shown in eq. (6). The coordinate $x$ of the focus center of the $l\text{th}$ OAM mode is proportional to the value of $l$. As is shown in Figure 1, an SLM is employed and spatially divided into several blocks along the $x$-axis, which are aligned with the coordinates of the focus centers of each OAM mode. Different phase modulations are settled on these blocks to achieve spectral shaping on OAM modes. In practical implementations, the OAM spectral shaper is determined by the target unitary matrix of U. Suppose that the desired spectral shaping function is noted as $g(l)$, the OAM spectral shaper can be implemented with the following wave front modulation function:

$\begin{array}{c}{{\displaystyle \Phi }}_{\text{OAM}}(x,y)=g(l),l=\left[\frac{2\text{π}ax}{f\lambda }\right],\end{array}$(8)

where $[\cdot ]$ rounds the value to the nearest integer. As indicated in Figure 1, such shape function of ${\text{Φ}}_{\text{OAM}}(x,y)$ is achieved with another SLM, which is placed at the focal plane of the auxiliary lens after the second optical element of ${\text{Φ}}_2$. After OAM spectral shaping, another two static optical elements are required to recover the OAM modes, i.e., to transform the Cartesian coordinate back to log-polar coordinate. The mode-sorter has been proved to be reversible [37] and the static optical elements to perform the reverse mode-sorting are exactly the same as those shown in eq. (6). Since the static elements for OAM mode-sorter can be replaced by customized spherical lens and cylindrical lens [36], only one programmable SLM performing ${\text{Φ}}_{\text{OAM}}(x,y)$ shown in eq. (8) would be enough for any OAM spectral shaping task.

SIMULATION

We have mathematically modeled the scheme shown in Figure 1 and simulated the field evolution according to Huygens-Fresnel principle under paraxial approximation. According to the parameters of SLMs (Holoeye Pluto), the simulation area is set to $1080\times 1080$ pixels with a pixel size of $8\text{ μm}\times 8\text{ μm}$. The operation wavelength is considered 1550 nm.

Firstly, the performance of our proposed OAM spectral shaper is evaluated. As shown in Figure 2A, L1 and L2 are lenses with focal length of $f$. A superposed state with 11 OAM modes is shown in Figure 2C as the input state. The mode interval between adjacent OAM encoding channels is chosen as 4 to compress the mode crossing induced by mode-sorter [38]. From left to right, Figure 2B shows field evolution of an initial superposed OAM state passing through the spectral shaper step-by-step (from index 1 to index 6), while the corresponding positions are marked in Figure 2A. Here, brightness and color indicate the amplitude and phase distributions, respectively. After the OAM mode-sorter, the initial angular momentum state is mapped to linear state. With the help of lens L1, linear momentum is transformed to linear coordinate in the focal plane. Then the OAM spectral shaping is performed by the SLM settled on the focal plane. After that, the optical vortex is rolled back again by another reversely operating mode-sorter.

Figure 2 Simulations of OAM spectral shaper with Huygens-Fresnel principle. (A) The scheme of OAM spectral shaper. (B) Step-by-step field evolution of an initial superposed OAM state passing through the spectral shaper. (C) Initial OAM spectrum. (D) Output OAM spectrum. The brightness and color indicate the amplitude and phase distributions, respectively.

The OAM spectral shaper under test is expressed as $g\left(l\right)=\text{π}l/4$. Since only lossless scheme and unitary matrices are concerned, the required diagonal operators in matrix decomposition in eq. (2) are phase-only and thus referred to as $D\left(g\left(l\right)\right)=\text{Diag}(\text{exp}(\text{i}g(l)))$ for simplicity. In the following sections, $f\left(\phi \right)$ and $g\left(l\right)$ are adopted instead of$\text{ exp}(\text{i}f(\phi ))$ and $\text{exp}(\text{i}g(l))$ for angular phase modulation and OAM phase spectral shaping, respectively.

In the case shown in Figure 2, $g\left(l\right)=\text{π}l/4$ means $0$ and $\text{π}$ phase modulation for adjacent OAM channels. Actually, it is the most challenging case since modulation function is “sharpest” as well as the phase crosstalk is maximum. Such extreme example is considered to examine the crosstalk tolerance of our proposed OAM spectral shaper. Besides, $g\left(l\right)=\text{π}l/4$ is a typical clock matrix and could also be achieved by field rotation of $\text{π}/4$ rad with a Dove prism [39]. Thus, the OAM spectral shaping effect can be observed intuitively through the rotation between the output and initial field in Figure 2B. We also provided the OAM spectrum after spectral shaper in Figure 2D. The complex OAM spectrum is calculated with mode-matching method [40]. With this method, the eigenmode expansions of the input and output OAM states during simulation can be evaluated and are expressed as $|{\alpha }_{\text{in}}\rangle ={{\displaystyle \sum }}^{\text{}}{c}_{l,\text{in}}|l\rangle$ and $|{\alpha }_{\text{out}}\rangle ={{\displaystyle \sum }}^{\text{}}{c}_{l,\text{out}}|l\rangle$, where $\left\{c_{l,\text{in}}\right\}$ and $\left\{c_{l\text{,out}}\right\}$ are input and output OAM spectrum coefficients, respectively. Thus, the spectral shaping function $g\left(l\right)$ can be deduced with $g\left(l\right)=c_{l,\text{out}}/c_{l,\text{in}}$. From these simulations, the fidelity value of output OAM spectrum compared with the target one is 0.988, which is estimated according to eq. (9). It can be seen that the OAM spectral shaper is valid for high-dimensional input states and the unitary operations could be achieved in succession.

Secondly, the OAM spectral shaping function $g\left(l\right)$ and two required angular modulation functions of $f_1\left(\phi \right)$ and $f_2\left(\phi \right)$ for a given target matrix are calculated through an optimization algorithm [20]. Our target is to achieve the maximum success probability. The success probability is proportional to the trace value of $V{V}^{†}$ subject to $\text{fidelity}\left(V,U\right)\ge 0.999$, where $U$ and $V$ denote target unitary matrix and designed matrix, respectively. The success probability of implemented matrix after optimization is calculated by $\text{Tr}(V{V}^{†})/N$. The fidelity is defined as [41]

$\begin{array}{c}\text{fidelity}(U,V)={\left|\frac{\text{Tr}(U^{†}V)}{\sqrt{\text{Tr}(U^{†}U)\cdot \text{Tr}(V^{†}V)}}\right|}^2.\end{array}$(9)

The designed matrix $V$ in eq. (9) is calculated by OAM spectral shaping function $g\left(l\right)$ and two angular modulation functions $f_1\left(\phi \right)$ and $f_2\left(\phi \right)$. Here, we give an example with dimensionality of $2\times 2$:

$\begin{array}{c}{\left(\begin{array}{ccc}\ddots & & \\ & \begin{array}{cc}{V}_{2\times 2}{}^{11}& V_{2\times 2}{}^{12}\\ V_{2\times 2}{}^{21}& V_{2\times 2}{}^{22}\end{array}& \\ & & \ddots \end{array}\right)}_{K\times K}=F^{†}D(f_1(\phi ))F\times D(g(l))\times F^{†}D(f_2(\phi ))F.\end{array}$(10)

The number $K$ in eq. (10) is the dimensionality of numerical optimization. The central $2\times 2$ region in the left side of eq. (10) is the optimized matrix $V$. We note that $K=N+2d$, where $N$ is the dimensionality of target matrix ($N=2$ for the case in eq. (10)) and $d$ is the number of guard bands. The guard bands are required to eliminate cut-off effects, which arises from the differences between the finite-dimension discrete Fourier transform required in eq. (2) and the infinite-dimension nature of the Fourier relation of OAM eigenstates and angle eigenstates expressed in eq. (4). In this work, both for the OAM-encoding and path-encoding, the number $K$ is chosen to be $K=64$ during the numerical calculation of spectral shaping functions. The DFT matrices in eq. (10) are normalized with $1/\sqrt{K}$ in amplitude so that $\text{Tr}(V{V}^{†})/N\le 1.$ Besides, we have $\text{Tr}(V{V}^{†})/N=1$ if and only if all other elements that share the same rows or columns with the central $N\times N$ block are zero.

The angular modulation functions $f_1\left(\phi \right)$ and $f_2\left(\phi \right)$ are quasi-continuous since the $[0,2\text{π})$ interval is divided into $K$ bins. The angular modulation functions are represented by sine expansion:

$\begin{array}{c}f(\phi )={\displaystyle \sum _{n=1}^p{A}_n\text{sin}(n\phi +{\theta }_n)}.\end{array}$(11)

The parameters $\left\{A_n\right\}\in [0,m\text{π})$ and$\left\{{\theta }_n\right\}\in [0,2\text{π})$ are to be determined. The scalar $m$ may grow as the size $N$ of target matrices increases, to ensure enough interactions between farther OAM channels. Here $m=1$ and $m=2$ are chosen for 2-dimensional and 3-dimensional unitary matrices, respectively. By choosing finite series expansion of sine expression, a relatively smooth modulation function can be obtained. Thus, the corresponding phase masks could have lower spatial resolution according to the Nyquist criterion and can be realized with more facility. The number $p$ in eq. (11) is determined by the dimension $N$ of target matrix since higher orders of sine components in angular modulation functions allow interactions between farther OAM eigenstates. $p=3$ is chosen in our simulations. The entire optimization is done with the MATLAB “fmincon” function. To further explain the optimization procedure and provide an executable example, the functions $g\left(l\right)$, $f_1\left(\phi \right)$ and $f_2\left(\phi \right)$ required for $2\times 2$ and $3\times 3$ Hadamard matrices $H_2$ and $H_3$ as well as modulation functions settled on SLMs are shown in Figure 3.

Figure 3 Optimized approach for OAM-domain Hadamard matrices with dimensionality of $2\times 2$ and $3\times 3$. (A) Modelled setup for achieving OAM-domain matrix transformation. (B) Optimized angular modulation functions $f_1\left(\phi \right)$ and $f_2\left(\phi \right)$ for Hadamard matrix $H_2$. (C) OAM spectral shaper $g\left(l\right)$ for $H_2$. (D) and (E) Optimized modulation functions for Hadamard matrix $H_3$.

$\begin{array}{c}{H}_2=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right),H_3=\frac{1}{\sqrt{3}}\left(\begin{array}{ccc}1& 1& 1\\ 1& {\text{e}}^{\text{i}2\text{π}/3}& {\text{e}}^{\text{i}4\text{π}/3}\\ 1& {\text{e}}^{\text{i}4\text{π}/3}& {\text{e}}^{\text{i}2\text{π}/3}\end{array}\right).\end{array}$(12)

The schematic setup for OAM-domain matrix transformation is shown in Figure 3A. Three SLMs labelled as SLM1, SLM2, and SLM3 perform the first angular modulator, OAM spectral shaper, and the second angular modulator, respectively. The colored phase modulation patterns corresponding to each SLM are also shown. For Hadamard matrix $H_2$, the optimized angular modulation functions are shown under polar coordinates in Figure 3B, where the red and blue solid lines indicate $f_1\left(\phi \right)$ and $f_2\left(\phi \right)$, respectively. The unit of polar axis in Figure 3B is rad. The optimized OAM spectral shaper noted as $g\left(l\right)$ is shown in Figure 3C, where the OAM eigenstates of $l=0$ and $l=1$ are chosen as encoding channels for $H_2$. In Figure 3C, 10 guard bands are plotted. Since the optimized spectral shaping functions tend to converge for more guard bands, much less OAM channels need to be modulated even when number $K$ in eq. (10) is chosen to be as large as 64.

Numerical estimation is made to quantize the influence of the number of guard bands on transformation fidelity and success probability. For those encoding channels beyond guard bands, there are two situations that have to be considered. The signals that enter those channels are left unchanged for OAM-encoding since those OAM channels always exist objectively. For path-encoding employed in experiments in section 4, the channels beyond guard bands are filtered out since the number of accessible encoding channels is determined by the size of SLMs (Experiment). For both cases, the fidelity and probability of unitary transformation as functions of numbers of guard bands are plotted in Figure 4. The target matrices are chosen as $H_2$ and $H_3$ in Figures 4A and 4B and Figures 4C and 4D, respectively.

Figure 4 Fidelity and probability of unitary transformation as functions of numbers of guard bands. The target matrices are chosen as $H_2$ and $H_3$ in (A, B) and (C, D), respectively. The channels beyond guard bands are unchanged in (A, C) and filtered out in (B, D).

According to the numerical results shown in Figure 4, the condition of $d\ge N-1$ would be enough for most practical cases. Thus, only two or three guard bands are required for practical $H_2$ gate. Due to the value of $K=64$, the optimized angular modulation function shown in Figure 3B seems quasi-continuous. Those optimization results for Hadamard matrix $H_3$ are shown in Figures 3D and 3E, where OAM eigenstates $l=-1$, $l=0$ and $l=1$ are chosen as encoding channels. According to the optimization, the fidelity above 0.999 and success probability above 0.970 are achieved simultaneously for $H_2$ and $H_3$. After the phase patterns settled on SLMs are obtained, the field evolution of Hadamard matrix $H_2$ has been calculated according to Huygens-Fresnel principle. The OAM eigenstates of $l=-3$ and $l=3$ are chosen to be the encoding channels. The transformation fidelity is calculated as 0.998.

In Figure 5, each row represents the field evolution of different initial OAM states and six columns are corresponding to the initial field, field after SLM1, field after mode-sorter, field before SLM2, field before L2, and output field (from left to right). All rows are achieved with the same angular phase patterns so that the same OAM phase spectral modulation functions are obtained. Moreover, in each map, the brightness and color are corresponding to the amplitude and phase, respectively. It should be mentioned that there is no need for post-selection to perform OAM-domain closed transformation with our proposal and thus the near-unity success probability is preserved. Thus, our results for implementing OAM-domain Hadamard matrix (or OAM beam splitter) may inspire some new applications such as observation of OAM-domain Hong-Ou-Mandel (HOM) effect [42] and OAM-encoded measurement-device-independent quantum key distribution (MDI-QKD) [43].

Figure 5 Step-by-step field evolution of OAM-domain Hadamard matrix.

Additionally, another important feature of our proposal is parallel computing, which is considered an advantage of optical information processing. With our scheme, the same matrix could be applied on different OAM channels simultaneously without any extra hardware cost. To further explain this, the optimized OAM spectral shaping function $g\left(l\right)$ is shown in Figure 6A for two $H_2$ gates acting on two groups of OAM channels with $l=-3\&-2$ and $l=4\&5$. For SLM modulation functions, the only change is to replace the OAM spectral shaping function $g\left(l\right)$ by $\begin{array}{c}\text{conv}(g(l),\delta (l-l_1)+\delta (l-l_2)+\cdots ),\end{array}$(13) where $\text{conv}(\cdot )$ denotes convolution and $\delta \left(l-l_1\right)$ denotes the Dirac delta function. The OAM charge $l_1,l_2,\dots$ is the central OAM channels utilized by each identical matrix. It is unnecessary to change the angular modulation functions $f_1\left(\phi \right)$ and $f_2\left(\phi \right)$ for such parallel processing. Two $H_2$ gates acting on OAM channels with $l=-14\&-10$ and $l=10\&14$ are simulated and the entire $4\times 4$ matrix operator is evaluated and shown in Figure 6B. Furthermore, the simulated results of dual $H_3$ gates are shown in Figure 6C. The utilized OAM channels are $l=-18\&-14\&-10$ and $l=10\&14\&18$. The $H_3$ gate is also noted as $3\times 3$ DFT. In Figure 6C, the $H_3$ gate is tested with input superposed OAM states of $|{\omega }_1\rangle ,|{\omega }_2\rangle$ and $|{\omega }_3\rangle$ one-by-one, while $|{\omega }_n\rangle$ denotes the nth column of $H_3{}^{†}$. With this test, the phase accuracy of the unitary matrix can be evaluated. The colored bars and empty bars in Figure 6C indicate input complex OAM spectrum and output OAM spectrum, respectively. It should be mentioned that the simulation in Figure 6 is made under the same three-layer model as shown in Figure 3A. It has been observed that the OAM-domain parallel gates work well without any extra hardware cost.

Figure 6 The OAM-domain parallel $H_2$ gates and $H_3$ gates achieved by the three-layer structure. (A) OAM spectral shaping function for dual $H_2$ gates. (B) The entire $4\times 4$ matrix operator made up of two $H_2$ gates. (C) Dual $H_3$ gates tested with input superposed OAM states of $|{\omega }_1\rangle ,|{\omega }_2\rangle$ and $|{\omega }_3\rangle$ one-by-one. (D) and (E) Performances of parallel $H_2$ and $H_3$ gates as functions of OAM channels separation, respectively.

In Figures 6D and 6E, the performances of parallel $H_2$ and $H_3$ gates are plotted as functions of OAM modes separation, respectively. It could be found that the parallel gates operate independently with negligible crosstalk if the interval of utilized OAM channels is larger than two times of the number of necessary guard bands.

Actually, the number of parallel matrices is only limited by the number of available OAM channels. Here, the mode interval of adjacent OAM channels is settled as 4 to avoid mode crosstalk induced by OAM mode-sorters. The OAM channel resources would be utilized more efficiently if some efforts are made to reduce the crosstalk [44]. According to the mathematic decomposition form in eq. (2) as well as the numerical optimization form in eq. (10), it is quite straight forward to extend this three-layer structure in Figure 1 to larger scales. With alternately cascaded angular spectral shaper layers and OAM spectral shaper layers, the system arrangement would grow with complexity of $(2N-1)$ [26] according to the dimensionality $N$ of the target unitary matrix.

In our design, the free space focus systems shown in the conceptual setup in Figure 3A only serve for the OAM mode-sorters. Besides this, there is no need for free-space propagation or extra Fourier transformation blocks between two adjacent layers of spectral shapers. For OAM domain unitary operations, this is one of the main differences to the MPLC method [21,22], in which some free-space propagation between adjacent layers is required. It should be mentioned that our proposal is not constrained by the actual method of OAM spectral shaping. In our simulation, the OAM spectral shaping is achieved with OAM mode-sorters, which would induce extra experimental difficulties. Actually, the system arrangement shown in Figure 1 could be largely simplified with compact OAM spectral shaper. For some specific target matrices, the recently reported single-step OAM shaper [34,45] may be employed.

EXPERIMENT

Due to the lack of OAM mode-sorter elements, which requires advanced 3D printing technology, proof-of-principle experiments in path domain are performed to verify our proposed scheme. To utilize the matrix decomposition in eq. (2), we notice that there is a one-to-one mapping from transverse wave vector in the object plane to transverse coordinate in the focal plane. To encode qudits, the path basis can be defined by a series of transverse coordinates lying on y-axis in the focal plane. Thus, the optical field $\widetilde{E_o}$ in the object plane is linked to the aforementioned path-encoded state vector through the Fourier transformation and expressed as

$\begin{array}{c}\widetilde{E_o}(y)={\displaystyle \sum _{n=-\infty }^{\text{+}\infty }{c}_n\text{exp}(-\text{i}n{k}_{y}y)},\end{array}$(14)

where $c_n$ denotes the complex amplitude of the nth path-encoding channel in the focal plane. It should be mentioned that $\widetilde{E_o}$ and the path-encoded state vector of ${{\displaystyle \sum }}^{\text{}}{c}_n|n\rangle$ are equivalent expressions of the same state under momentum space and position space, respectively. Therefore, unitary operations can be achieved by alternately performing spectral shaping in momentum space and position space. The distance between adjacent paths is noted as $L$, while the distance between adjacent SLMs is set to $2f$. This leads to a corresponding transverse wave vector of

$\begin{array}{c}{k}_y=\frac{L}{2f}\left(\frac{2\text{π}}{\lambda }\right)=\frac{\text{π}L}{\lambda f}.\end{array}$(15)

Figure 7 (A) Conceptual set-up of path-encoded matrix transformation, and the modulation functions settled on SLMs. (B) Optimized phase modulation function settled on SLM2 for implementing $H_3$ gate. (C) One period of optimized phase gratings settled on SLM1 and SLM3 for $H_3$ gate.

The conceptual setup is shown in Figure 7A. The modulation functions on SLM1 and SLM3 are settled as the following to introduce the required transverse wave vector components: $\begin{array}{c}f(y)={\displaystyle \sum _{n=1}^p{A}_n\text{sin}(n{k}_{y}y+{\theta }_n)}.\end{array}$(16) Eq. (16) is very similar to the azimuthal phase modulation function in eq. (11) for OAM-domain matrix transformation. The parameters $\left\{A_n\right\}\in [0,2\text{π})$ and $\left\{{\theta }_n\right\}\in [0,2\text{π})$ are to be determined by optimizations and possess the same meanings as shown in eq. (11). The only difference is that the azimuthal coordinate $\phi$ is replaced by linear coordinate $y$. The modulation function on SLM2 is settled as $\begin{array}{c}\Phi (x,y)=g\left(\left[\frac{y}{L}\right]\right),\end{array}$(17) which is also a direct mapping of eq. (8). The function of $g(n)$ represents the phase modulation applied to the nth path channel. Practically, $g(n)$, $f_1(y)$, $f_2(y)$ are optimized with the same method as $g(l)$, $f_1(\phi )$, $f_2(\phi )$ in OAM-domain unitary transformation. The $H_3$ gate is shown as an example to clarify it. The optimized $g(n)$ is shown in Figure 7B, where path channels $n=-1\&0\&1$ are chosen for encoding. The $f_1(y)$ and $f_2(y)$ functions are shown in Figure 7C as phase-only gratings programmed on SLM1 and SLM3, respectively. Gratings can be considered spectral shaping in transverse momentum space. Here, only one period of grating is plotted for simplicity and clarity. A lens with focal length of $f$ is added on SLM2 to maintain the Fourier relation between optical fields on adjacent SLM planes. It should be mentioned that the one-dimensional Fourier relation is employed here, which can be expressed as a Fourier transform matrix and is different from the two-dimensional Fourier transform that is achieved by free space propagation. Two extra lenses noted as L1 and L2 in Figure 7A are employed to compensate for transverse wave vector so that the setup can be scalable. As indicated in eq. (14), we are only interested in the field variations along y-axis in the transverse plane. The field along x-axis is simply set as Gaussian shape and some ancillary lenses are programmed on SLMs to compensate the natural expansion of Gaussian beams.

Eighty different unitary transformations have been performed with the experimental setup as shown in Figure 8, while 75 of them are randomly generated. For each different target matrix, a particular set of modulation functions has to be optimized. The number of guard bands is set as $d=2$. As shown in Figure 8, a laser beam of zero-order Gaussian profile operating at 1550 nm (RIO Orion) is incident to the optical system through a collimator. Three SLMs (Holoeye Pluto) operating with reflective mode are utilized and the focal length $f$ mentioned before is set as 40 cm. The distance $L$ between adjacent paths is designed as 2 mm considering the reflective area of SLMs. Specially, the input state is both generated and switched by SLM1. The transverse k-space expression $\widetilde{E_o}$ in eq. (14) of the input state, or equivalently saying, the interference pattern of the input beams at SLM1 plane in conceptual setup in Figure 7A is pre-calculated and directly settled on SLM1 to modulate the incident Gaussian beam from the collimator in experimental setup shown in Figure 8.

Figure 8 Sketch of the experimental set-up for path-encoded matrix transformation. HWP: half wave plate, SLM: spatial light modulator, CCD: charge-coupled device camera.

The experimental results are summarized in Figures 9A and 9B. The distributions of fidelity and success probability of implemented matrices are shown as solid bars, respectively. Here, the fidelity is evaluated by a charge coupled device (CCD) camera placed at the output port in the experimental setup shown in Figure 8. Since the target matrix $U$ is unitary, the value of $\text{fidelity}\left(V,U\right)$ between $U$ and implemented matrix $V$ equals $\text{fidelity}\left(V{U}^{†},U{U}^{†}\right)$ according to eq. (9). $\begin{array}{c}\text{fidelity}(U,V)=\text{fidelity}(U{U}^{†},V{U}^{†})={\left|\frac{\text{Tr}(V{U}^{†})}{\sqrt{N\text{Tr}[{(V{U}^{†})}^{†}(V{U}^{†})]}}\right|}^2.\end{array}$(18)

Figure 9 Experimental results of path-encoded matrix transformation. The distributions of (A) fidelity and (B) success probability for 80 different tests. The typical experimental results for (C) Hadamard matrix H3, (D) H2, and (E–H) some other unitary operators. For each group, there are the amplitude of target matrix U, the measured amplitude and the measured phase of implemented matrix V from left to the right. The color bar is also shown.

The amplitude of $V{U}^{†}$ can be directly measured by the CCD camera. Actually, to measure $V{U}^{†}$ is a phase test. During the test, the input states of $|{\psi }_1\rangle$, $|{\psi }_2\rangle$ and $|{\psi }_3\rangle$ are settled as the columns of target matrix $U^{†}$. The results would also form an identical matrix if the implemented matrix $V$ is sufficiently accurate, as shown in Figures 9C–9H. The phase differences among diagonal elements of $V{U}^{†}$ are assumed to be zero during this estimation because subsequent phase shifters could be used for error correction [46]. Besides, according to eq. (2), the implementing errors of any diagonal matrix $D$ would influence both relative phase within matrix rows and relative phase between matrix rows simultaneously. Thus, if the relative phase within matrix rows is accurate evidently, it is unlikely to remain a large relative phase between matrix rows.

The amplitude of $V$ is scanned with input states $|1\rangle$, $|2\rangle$ and $|3\rangle$, providing an intuitive vision of absolute value of matrix elements. The success probability is calculated by the output optical power from the three path-encoding channels normalized by the total output power received by the CCD camera.

The average value of fidelity is $0.97\pm 0.02$ among 80 different tests. The success probability is shown in Figure 9B with an average value of $0.96\pm 0.03$. The main reason of deteriorated probability is that part of optical energy entered guard bands and failed to vanish through interference. The experimental results of $H_3$ gate and $H_2$ gate (acting on the first and third channel) are shown in Figures 9C and 9D, respectively. The original data captured by CCD camera are also shown. As concrete examples, the data of another four typical unitary matrices are shown in Figures 9E–9H, while the fidelity value ranges from 0.96 to 0.99.

A brief error analysis for experimental results is as follows. First, the conditioned maximum algorithm is utilized to determine the modulation functions of $g(n)$, $f_1(y)$, $f_2(y)$. The fidelity values are conditioned by $\text{fidelity}\ge 0.999$. As a result, the average probability value tends to be lower than the average fidelity value. Second, a five-layer structure is sufficient for arbitrary $3\times 3$ matrices according to the decomposition theory [26]. Though the three-layer structure employed in this work is enough for most cases, the fidelity and probability could not both achieve 1 at the same time for some particular target matrices, hence there is a tradeoff. Third, experimental errors in terms of optical misalignment and nonlinear effects induced by CCD camera would also deteriorate both the fidelity and probability. The saturated effect of some CCD images will result in lower fidelity and efficiency estimation. The fidelity is calculated through phase test shown in Figures 9C and 9D. The saturated effect induces reduction of intensity of main lobs but has little influence on side lobs, thus lowers the weight of main lobs as well as the value of $\text{Tr}(V{U}^{†})$. The detailed evaluation considering a $2\times 2$ case is below.

After normalization, suppose that an averaged reduction of $E_{\text{reduct}}$ is induced for all diagonal elements of $V{U}^{†}$ due to saturation. At the same time, the difference among diagonal elements of $V{U}^{†}$ would decrease by $E_{\text{even}}$. The $E_{\text{even}}$ value is comparable with $E_{\text{reduct}}$ value and we have $E_{\text{even}}\le E_{\text{reduct}}$ if the derivative value of the saturation function is less than 1. According to the definition of fidelity in eq. (18), the numerator term of $\text{Tr}(V{U}^{†})$ would decrease by $2E_{\text{reduct}}$, while the denominator term of $\text{Tr}[{(V{U}^{†})}^{†}V{U}^{†}]$ would decrease from about ${(1+E_{\text{even}})}^2+{(1-E_{\text{even}})}^2$ to $2$, assuming that the values of diagonal elements of $\text{Tr}(V{U}^{†})$ are about 1. The denominator term would decrease by $2E_{\text{even}}{}^2$, which is one order smaller than the reduction of numerator. The square root calculation in denominator further shrinks the reduction. Hence, the saturated effect will result in lower fidelity estimation.

The amplitude measurement shown in Figures 9C and 9D is used to evaluate transformation efficiency (success probability). The saturated effect lowers the measured optical power from the three path-encoding channels and decreases the efficiency value.

It could be concluded that the three-layer structure is valid to achieve both high fidelity and high success probability for most 3-dimensional unitary matrices. Thus, the experimental results of path-encoded unitary transformation are important evidence to prove the versatility of our OAM domain proposal.

DISCUSSION

First, the relations and differences between our work and the frequency-domain protocols are discussed. Our method shares the same basic mathematical theorem [25] of matrix decomposition with Lougovski’s work [27]. In their work, the frequency-time Fourier relation is employed. We notice that frequency encoding is different from other photonic encoding since energy conservation is not satisfied for unitary linear operations such as the Pauli $\widehat{X}$ gate. Thus, it cannot be achieved by only passive setup. In this work, we focus on OAM and path domains employing the OAM-angle Fourier relation and the position-momentum Fourier relationship. Correspondingly, the unitary operations within both OAM and path domains could be implemented with passive optical elements.

Next, a brief comparison between our proposal and previous approaches is made for both OAM-domain and path-domain. Different from other OAM-domain approaches, the cut-off effects from infinite OAM-encoding basis to finite subset are largely alleviated since ancillary guard bands are employed. During the unitary operation, part of the input energy would enter the guard bands, but such bands could be managed to vanish eventually through interference according to meticulous design. The number of required spectral shaper layers grows in O(N) scale according to dimensionality N of unitary operation. The reason is that each spectral shaper could provide O(N) modulation on all the channels in parallel, thus the whole setup can fulfill $N^2$ independent real numbers required by a unitary operation. Compared with the MPLC approach [21,22,47], our design has similar system arrangement of O(N) layers as dimensionality $N$ grows. However, near-unity transformation efficiency and near-perfect fidelity can be achieved theoretically at the same time with our proposal according to the results of numerical simulations. Besides, our proposed scheme can perform parallel unitary operations on different OAM encoding channels simultaneously without increasing system arrangement. It should be mentioned that although our method is based on the mathematic decomposition of the Fourier/diagonal matrix, there is no DFT matrix required. Specifically, our proposed unitary operation scheme for OAM-domain is performed within OAM-domain and angle-domain alternately, while there is no actual DFT transformation between them. Actually, the OAM shaper is employed to perform the diagonal operation, while the DFT matrix is required for the MPLC method. However, since there is no requirement of DFT transformation, our scheme has to be performed within two optical DOF linked by Fourier transformation, e.g., the OAM-domain vs. angle-domain and path-domain vs. wavevector-domain.

For path-domain operation, both of our protocol and Reck’s scheme can be adopted and keep similar features of programmability, scalability, and near-unity success probability in theory, although they are based on different matrix decomposition. In this work, our experiment is based on three SLMs to perform the transformation. In our setup, the insertion loss of utilized three SLMs is about 11.9 dB. Thus, the scale of our experimental scheme would suffer from the insertion loss of the SLM since more SLMs are required for transformation with higher dimensionalities. It should be mentioned that our proposal is also manageable in both free space optics and integrated circuits. The SLMs could be replaced by three-dimensional set of metasurfaces [48] to achieve a compact and low-loss module for unitary operations. Due to the similar alignment of path channels in our proposal and Reck’s scheme, an on-chip version of our design could also be implemented. The key function of Fourier transform could be achieved by the structures of Rowland circle according to the recently reported work [49]. Note that in this design no MMI couplers [29] nor fully connective $N\times N$ blocks [30] are needed. Our proposal would provide an alternative choice in addition to the reported integrated architectures [29,30].

In conclusion, we have proposed and demonstrated a scheme for unitary operation in OAM and optical path domains. A three-layer setup is implemented to evaluate the performance. Field simulation of OAM-domain Hadamard matrices and experimental implementation of randomly generated path-domain unitary matrices have been carried out. For path domain, the fidelity value of $0.97\pm 0.02$ and transformation efficiency of $0.96\pm 0.03$ are experimentally evaluated simultaneously.

Acknowledgments

The authors would like to thank Dr. Peng Zhao in Appotronics Corporation Ltd. for valuable discussions and helpful comments.

Funding

This work was supported by the National Key Research and Development Program of China (Grant Nos. 2018YFB2200402 and 2017YFA0303700), the National Natural Science Foundation of China (Grant No. 61875101), the Beijing Academy of Quantum Information Science, the Beijing National Research Center for Information Science and Technology (BNRist), the Beijing Innovation Center for Future Chip, and the Tsinghua University Initiative Scientific Research Program.

Conflict of interest

The authors declare that they have no conflict of interest.

References

All Figures

	Figure 1 The scheme of OAM-domain unitary operation with three-layer.
In the text

	Figure 2 Simulations of OAM spectral shaper with Huygens-Fresnel principle. (A) The scheme of OAM spectral shaper. (B) Step-by-step field evolution of an initial superposed OAM state passing through the spectral shaper. (C) Initial OAM spectrum. (D) Output OAM spectrum. The brightness and color indicate the amplitude and phase distributions, respectively.
In the text

Figure 3 Optimized approach for OAM-domain Hadamard matrices with dimensionality of $2\times 2$ and $3\times 3$. (A) Modelled setup for achieving OAM-domain matrix transformation. (B) Optimized angular modulation functions $f_1\left(\phi \right)$ and $f_2\left(\phi \right)$ for Hadamard matrix $H_2$. (C) OAM spectral shaper $g\left(l\right)$ for $H_2$. (D) and (E) Optimized modulation functions for Hadamard matrix $H_3$.

In the text

	Figure 4 Fidelity and probability of unitary transformation as functions of numbers of guard bands. The target matrices are chosen as $H_2$ and $H_3$ in (A, B) and (C, D), respectively. The channels beyond guard bands are unchanged in (A, C) and filtered out in (B, D).
In the text

	Figure 5 Step-by-step field evolution of OAM-domain Hadamard matrix.
In the text

Figure 6 The OAM-domain parallel $H_2$ gates and $H_3$ gates achieved by the three-layer structure. (A) OAM spectral shaping function for dual $H_2$ gates. (B) The entire $4\times 4$ matrix operator made up of two $H_2$ gates. (C) Dual $H_3$ gates tested with input superposed OAM states of $|{\omega }_1\rangle ,|{\omega }_2\rangle$ and $|{\omega }_3\rangle$ one-by-one. (D) and (E) Performances of parallel $H_2$ and $H_3$ gates as functions of OAM channels separation, respectively.

In the text

	Figure 7 (A) Conceptual set-up of path-encoded matrix transformation, and the modulation functions settled on SLMs. (B) Optimized phase modulation function settled on SLM2 for implementing $H_3$ gate. (C) One period of optimized phase gratings settled on SLM1 and SLM3 for $H_3$ gate.
In the text

	Figure 8 Sketch of the experimental set-up for path-encoded matrix transformation. HWP: half wave plate, SLM: spatial light modulator, CCD: charge-coupled device camera.
In the text

Figure 9 Experimental results of path-encoded matrix transformation. The distributions of (A) fidelity and (B) success probability for 80 different tests. The typical experimental results for (C) Hadamard matrix H3, (D) H2, and (E–H) some other unitary operators. For each group, there are the amplitude of target matrix U, the measured amplitude and the measured phase of implemented matrix V from left to the right. The color bar is also shown.

In the text