Universal multiport interferometer

Universal unitary (mode transformations)

In quantum mechanics, an operator ${\textstyle {\hat {U}}}$ is unitary if it and its Hermitian conjugate, ${\textstyle {\hat {U}}^{\dagger }}$, satisfy the condition ${\textstyle {\hat {U}}{\hat {U}}^{\dagger }={\hat {I}}={\hat {U}}^{\dagger }{\hat {U}}}$, where ${\textstyle {\hat {I}}}$ is the identity operator [1]. By extension, a unitary transformation describes the change of a state, ${\textstyle |\psi \rangle \rightarrow |\psi '\rangle ={\hat {U}}|\psi \rangle }$, or the change an operator, ${\textstyle {\hat {A}}\rightarrow {\hat {A}}'={\hat {U}}{\hat {A}}{\hat {U}}^{\dagger }}$, under the action of a unitary operator. Hence, a platform capable of imposing completely general unitary transformations is referred to, itself, as a universal unitary. Importantly, any lossless experimental setup can be described by a unitary operator [2]. However, designing an experimental setup that implements any given unitary operator is, unfortunately, a more challenging task.

In general, fully determining any ${\textstyle N}$-dimensional unitary requires specifying ${\textstyle N^{2}}$ independent real parameters [3]. For the simple case of transforming a two-beam array, a universal unitary can be implemented using a variable beam splitter and three phase-shifters [3,4]. In 1994, Michael Reck and Anton Zeilinger generalized this well-known approach by proving that variable beam splitters and phase-shifters, when arranged in an interferometric mesh with ${\textstyle N^{2}}$ arms, can be used to impose any ${\textstyle N{\text{x}}N}$ unitary mode transformation [2]. Using their deterministic algorithm to decompose a given unitary into a triangular network of these two optical elements, it is possible to experimentally realize a discrete universal unitary, specifically for ${\textstyle N{\text{x}}N}$ mode transformations. The resulting device is commonly referred to as a universal multiport interferometer.

In 2016, Clements et al. introduced a variation of Reck and Zeilinger's decomposition, again using beam splitters and phase shifters, but arranged in a symmetrically crossing network as opposed to a triangular network. Importantly, this variation has a smaller optical depth - the longest path through the interferometric mesh - and thus experiences lower propagation losses.

Transformations beyond the subset of linear ${\textstyle N{\text{x}}N}$ mode transformations cannot be realized, in general, using such decompositions. Hence, they are not universal for the complete set of unitary transformations. This distinguishes it from other universal unitary platforms which implement transformations existing within another subset. For example, the universal unitary commonly discussed in quantum computing is the universal gate whereby any ${\textstyle N}$-qubit gate can be realized by a circuit of single qubit gates and CNOT gates. The classical analog of such universality is the idea that an arbitrary Boolean function can be realized using a combination of NOT gates and any one of the two-bit gates (e.g. AND, OR).

Definitions

Unitary operator and transformations

In quantum mechanics, the inverse of an operator ${\textstyle {\hat {O}}}$ is represented by ${\textstyle {\hat {O}}^{-1}}$ , and is defined such that ${\textstyle {\hat {O}}{\hat {O}}^{-1}={\hat {I}}}$ where ${\textstyle {\hat {I}}}$ is the identity operator. Among the complete set of operators, those that satisfy the additional condition that ${\textstyle {\hat {O}}^{-1}={\hat {O}}^{\dagger }}$ are called unitary and are denoted ${\textstyle {\hat {U}}}$. This is where the superscript "$\dagger$" signifies the Hermitian conjugate (or Hermitian adjoint, Hermitian transpose, conjugate transpose) which is the conjugate transpose of the given operator, i.e., ${\textstyle {\hat {O}}^{\dagger }=\left({\hat {O}}^{*}\right)^{T}}$. Therefore, any unitary operator ${\textstyle {\hat {U}}}$ must, by definition, satisfy the following condition:$${\hat {U}}{\hat {U}}^{\dagger }={\hat {I}}={\hat {U}}^{\dagger }{\hat {U}}.$$Unitary operators play an important role in quantum mechanics as most physical processes evolve unitarily. An exception to this is measurement, which is an inherently non-unitary process.

Property 1: Preserving inner product

Suppose there exists a unitary operator ${\textstyle {\hat {U}}}$ which transforms the state ${\textstyle |v\rangle }$ according to ${\textstyle |v\rangle \rightarrow {\hat {U}}|v\rangle =|v'\rangle }$ and which transforms the state ${\textstyle |w\rangle }$ according to ${\textstyle |w\rangle \rightarrow {\hat {U}}|w\rangle =|w'\rangle }$. Then, the inner product of the transformed states, $|v'\rangle$ and $|w'\rangle$, in bra-ket notation is given by,$${\begin{array}{lcl}\langle v'|w'\rangle &=&\left({\hat {U}}|v\rangle \right)^{\dagger }{\hat {U}}|w\rangle \\&=&\left(\langle v|{\hat {U}}^{\dagger }\right){\hat {U}}|w\rangle \\&=&\langle v|{\hat {U}}^{\dagger }{\hat {U}}|w\rangle \\&=&\langle v|{\hat {I}}|w\rangle \\\langle v'|w'\rangle &=&\langle v|w\rangle .\end{array}}$$Hence, the inner product is preserved under the unitary transformation of states. The consequence of this property is information is conserved under unitary transformation.

Property 2: Diagonalization of Hermitian operators

If an operator ${\textstyle {\hat {O}}}$ is Hermitian such that ${\textstyle {\hat {O}}={\hat {O}}^{\dagger }}$, then according to finite-dimensional spectral theorem, ${\textstyle {\hat {O}}}$ can be diagonalized by some unitary operator ${\textstyle {\hat {U}}}$. Mathematically, this property is represented by,$${\hat {O}}={\hat {U}}{\hat {\lambda }}{\hat {U}}^{\dagger },$$where ${\textstyle {\hat {\lambda }}}$ is a diagonal matrix with ${\textstyle \lambda _{i}}$ as elements, i.e.,$${\hat {\lambda }}\Rightarrow \sum _{i}^{d}\lambda _{i}|\lambda _{i}\rangle \langle \lambda _{i}|.$$It follows that,

${\begin{array}{lcl}{\hat {U}}^{\dagger }{\hat {O}}{\hat {U}}&=&{\hat {U}}^{\dagger }\left({\hat {U}}{\hat {\lambda }}{\hat {U}}^{\dagger }\right){\hat {U}}\\&=&{\hat {I}}{\hat {\lambda }}{\hat {I}}\\{\hat {U}}^{\dagger }{\hat {O}}{\hat {U}}&=&{\hat {\lambda }},\end{array}}$

or, in other words, the action of ${\textstyle {\hat {U}}}$ can be described as a mere change-of-basis of ${\textstyle {\hat {O}}}$ to one in which it is diagonal, namely, the set$\{|\lambda _{i}\rangle \langle \lambda _{i}|\}_{d}$. In the physical world, unitary operators are the generalizations of rotations in 3-dimensional Euclidean space.

Property 3: Dimensionality

Suppose there exists a unitary operator ${\textstyle {\hat {U}}}$ of dimension ${\textstyle N}$, then the number of elements contained in its matrix representation is ${\textstyle N^{2}}$. The elements are, in general, complex so this corresponds to a total of ${\textstyle 2N^{2}}$ real independent parameters. However, the conditions of unitarity reduce the number of real independent parameters.

According to the definition of a unitary operator, ${\textstyle {\hat {U}}}$ is required to satisfy ${\textstyle {\hat {U}}^{\dagger }{\hat {U}}={\hat {I}}}$, but ${\textstyle {\hat {I}}}$ is Hermitian. Therefore, the diagonal elements of ${\textstyle {\hat {U}}^{\dagger }{\hat {U}}}$ must be real valued. This corresponds to ${\textstyle N}$ real independent parameters.

The Hermiticity of ${\textstyle {\hat {I}}}$ implies ${\textstyle {\hat {U}}^{\dagger }{\hat {U}}={\hat {U}}{\hat {U}}^{\dagger }}$, which further implies the lower triangular elements of ${\textstyle {\hat {U}}^{\dagger }{\hat {U}}}$ can be expressed in terms of the upper triangular elements such that the off-diagonal elements correspond to only ${\textstyle 2\left(N(N-1)/2\right)}$ real independent parameters instead of ${\textstyle 2\left(N(N-1)\right)}$.

The total number of real independent parameters is thus reduced from ${\textstyle 2N^{2}}$ to:$$N_{\text{ind}}=2N^{2}+N+2\left(N(N-1)/2\right),$$ and hence the number of free parameters is:$$N_{\text{Free}}=2N^{2}-\left(N+2\left(N(N-1)/2\right)\right)=N^{2}.$$Since ${\textstyle {\hat {U}}^{\dagger }{\hat {U}}}$ has the same dimension as ${\textstyle {\hat {U}}}$, the number of free parameters describing ${\textstyle {\hat {U}}^{\dagger }{\hat {U}}}$ is the same number of free parameters describing ${\textstyle {\hat {U}}}$.

For example, suppose ${\textstyle {\hat {U}}}$ is two-dimensional.$${\hat {U}}^{\dagger }{\hat {U}}\Rightarrow {\begin{bmatrix}u_{11}^{*}&u_{21}^{*}\\u_{12}^{*}&u_{22}^{*}\end{bmatrix}}\cdot {\begin{bmatrix}u_{11}&u_{12}\\u_{21}&u_{22}\end{bmatrix}}={\begin{bmatrix}|u_{11}|^{2}+|u_{21}|^{2}&u_{11}^{*}u_{12}+u_{21}^{*}u_{22}\\u_{12}^{*}u_{11}+u_{22}^{*}u_{21}&|u_{12}|^{2}+|u_{22}|^{2}\end{bmatrix}}={\begin{bmatrix}1&0\\0&1\end{bmatrix}}$$$${\hat {U}}{\hat {U}}^{\dagger }\Rightarrow {\begin{bmatrix}u_{11}&u_{12}\\u_{21}&u_{22}\end{bmatrix}}\cdot {\begin{bmatrix}u_{11}^{*}&u_{21}^{*}\\u_{12}^{*}&u_{22}^{*}\end{bmatrix}}={\begin{bmatrix}|u_{11}|^{2}+|u_{12}|^{2}&u_{11}u_{21}^{*}+u_{12}u_{22}^{*}\\u_{21}u_{11}^{*}+u_{22}u_{12}^{*}&|u_{21}|^{2}+|u_{22}|^{2}\end{bmatrix}}={\begin{bmatrix}1&0\\0&1\end{bmatrix}}$$The Hermiticity condition gives, for the upper and lower triangular elements:$${\begin{cases}u_{11}^{*}u_{12}+u_{21}^{*}u_{22}=u_{11}u_{21}^{*}+u_{12}u_{22}^{*}\\u_{12}^{*}u_{11}+u_{22}^{*}u_{21}=u_{21}u_{11}^{*}+u_{22}u_{12}^{*},\end{cases}}$$which are equivalent equations up to complex conjugation. By setting the right-hand side to the matrix representation of ${\textstyle {\hat {I}}}$, one arrives at the following system of three equations:$${\begin{cases}|u_{11}|^{2}+|u_{21}|^{2}=1\\|u_{12}|^{2}+|u_{22}|^{2}=1\\u_{11}^{*}u_{12}+u_{21}^{*}u_{22}=0.\end{cases}}$$As there are only three equations, one of the parameters is free. Using the formula derived above for ${\textstyle N_{\text{free}}}$, one finds the expected number of free parameters for a ${\textstyle 2{\text{x}}2}$ unitary is one, which is in agreement with the example.

Photon modes

In the second-quantization of quantum mechanics, the electromagnetic (EM) field is quantized and this allows one to define a photon. In this sense, a photon is an excitation in a specific mode of the quantized EM field. For example, the mode of a plane wave is characterized by its wavevector ${\textstyle {\vec {k}}}$ and a photon associated with a given plane wave mode must be labelled by such wavevector in order to differentiate it. In addition, two frequently studied spatial mode sets are the set of Hermite-Gaussian (HG) modes and the set of Laguerre-Gaussian (LG) modes. These are two transverse modal bases representing the eigenfunctions of the paraxial wave equation in Cartesian and cylindrical coordinates, respectively.

To analyze a given set of modes, the electromagnetic field is treated as a set of quantum harmonic oscillators. The states are then built-up using a creation operator, ${\textstyle {\hat {a}}^{\dagger }}$, and an annihilation operator, ${\textstyle {\hat {a}}}$, both of which are defined for the mode basis being studied. These operators obey bosonic commutation relations:$$\left[{\hat {a}},{\hat {a}}\right]=0=\left[{\hat {a}}^{\dagger },{\hat {a}}^{\dagger }\right],\left[{\hat {a}},{\hat {a}}^{\dagger }\right]=1$$and they transform the photon number state ${\textstyle |n\rangle }$, which contains ${\textstyle n}$ photons (again, existing in the mode being studied), according to the following relations:$${\begin{array}{lcl}{\hat {a}}|n\rangle &=&{\sqrt {n}}|n-1\rangle \\{\hat {a}}^{\dagger }|n\rangle &=&{\sqrt {n+1}}|n+1\rangle .\end{array}}$$The photon number states ${\textstyle |n\rangle }$ are orthonormal, i.e., ${\textstyle \langle n|m\rangle =\delta _{n,m}}$, and form a basis for the Hilbert space of such ${\textstyle n}$ photons. In general, an ${\textstyle n}$ photon state ${\textstyle |\psi _{n}\rangle }$ can be expressed as:$$|\psi _{n}\rangle =\sum _{n}c_{n}|n\rangle$$where ${\textstyle \{c_{n}\}}$​ are complex coefficients satisfying the normalization condition ${\textstyle \sum _{n}|c_{n}|^{2}=1}$. In other words, the general photon state ${\textstyle |\psi _{n}\rangle }$ can be decomposed into a superposition of photon number states with weights given by ${\textstyle c_{n}=\langle \psi |n\rangle }$.

Unitary mode transformations

In the previous section, the creation and annihilation operators were defined for a given basis mode basis described by ${\textstyle {\hat {a}}}$. When ${\textstyle {\hat {a}}}$ acts on a specific mode ${\textstyle j}$, this labels both the creation and annihilation operators as a subscript. The set of all modes ${\textstyle \{{\hat {a}}_{j}\}}$ where${\textstyle {\ce {j=\{0,1,...,J\}}}}$ forms the basis, with ${\textstyle J}$ being the maximum number of modes in the system. Using a unitary transformation, the creation operators of can be expressed in a rotated mode basis related to the original basis as:$${\hat {b}}_{k}^{\dagger }=\sum _{j=0}^{J}U_{k,j}{\hat {a}}_{j}^{\dagger },$$where ${\textstyle U_{k,j}}$ are elements of the unitary matrix ${\textstyle U}$. This transformation allows photon states to be expressed in a new basis, facilitating tasks such as simplifying descriptions or enabling experimental realizations in specific bases.

For example, states with ${\textstyle n}$-photons in a mode ${\textstyle k}$ are generated by applying the operator ${\textstyle \left({\hat {b}}_{k}^{\dagger }\right)^{n}}$ to the vacuum state:$$|n\rangle _{k}={\frac {1}{\sqrt {n!}}}\left({\hat {b}}_{k}^{\dagger }\right)^{n}|0\rangle .$$Substituting the expression for ${\textstyle \left({\hat {b}}_{k}^{\dagger }\right)^{n}}$​, this can be written as:$$|n\rangle _{k}={\frac {1}{\sqrt {n!}}}\left(\sum _{j=0}^{J}U_{k,j}{\hat {a}}_{j}^{\dagger }\right)^{n}|0\rangle .$$The vacuum state ${\textstyle |0\rangle }$, which contains no photons, remains invariant under unitary transformation.

Using this formalism, one may mathematically test the effect of a given unitary on a set of modes (e.g., spatial modes of light).

Motivation

Givens rotations

A Givens rotation, also known as a Jacobi rotation, is a well-known operation in linear algebra that performs a rotation in a two-dimensional subspace of a higher-dimensional space. Mathematically, it the Givens rotation has the following matrix representation:$$G(i,j,\phi )={\begin{bmatrix}1&\cdots &0&\cdots &0&\cdots &0\\\vdots &\ddots &\vdots &&\vdots &&\vdots \\0&\cdots &\cos(\phi )&\cdots &\sin(\phi )&\cdots &0\\\vdots &&\vdots &\ddots &\vdots &&\vdots \\0&\cdots &-\sin(\phi )&\cdots &\cos(\phi )&\cdots &0\\\vdots &&\vdots &&\vdots &\ddots &\vdots \\0&\cdots &0&\cdots &0&\cdots &1\end{bmatrix}},$$where ${\textstyle \{i,j\}}$ denote the rows in which the rotation terms appear. The left multiplication of ${\textstyle G(i,j,\phi )}$ on another matrix ${\textstyle M}$ results in only rows ${\textstyle i}$ and ${\textstyle j}$ of ${\textstyle M}$ being affected. The effect of the Givens operation thus reduces to the transformation of two input amplitudes, ${\textstyle a}$ and ${\textstyle b}$ (where ${\textstyle a}$ and ${\textstyle b}$ are elements of the ${\textstyle i}$- and ${\textstyle j}$-th rows of ${\textstyle M}$, respectively), into the new amplitudes, ${\textstyle c}$ and ${\textstyle d}$, as follows:$${\begin{bmatrix}c\\d\end{bmatrix}}={\begin{bmatrix}\cos \phi &\sin \phi \\-\sin \phi &\cos \phi \end{bmatrix}}{\begin{bmatrix}a\\b\end{bmatrix}}.$$The Givens rotation can be used to zero out a specific element of a vector (e.g., making ${\textstyle d=0}$) or systematically triangularize a matrix, making it essential for linear algebra algorithms like matrix factorization and solving systems of equations.

This is the same matrix that defines the Jacobi rotation, but the choice of angle ${\textstyle \phi }$ differs by a factor of approximately 2.

Experimental Givens rotation device

In 1986, Mirsalehi et al. proposed a lossless integrated-optical implementation of a Givens rotation device using electro-optic grating diffraction and phase modulators to perform the necessary operations for efficient and high-speed data processing.

The proposed device operates with two coherent, monochromatic input waves representing amplitudes ${\textstyle a}$ and ${\textstyle b}$. The phase modulators adjust the relative phase of these inputs, while the diffraction grating computes the sine and cosine components. The outputs ${\textstyle c}$ and ${\textstyle d}$ are coherently combined to produce the desired rotation.

1. The input and output light signals are guided in waveguides. The use of z-cut lithium niobate waveguides ensures low-loss and high-speed operation.
2. A thick diffraction grating modulated by a voltage generates the sine and cosine multiplications naturally. The input wave amplitudes ${\textstyle a}$ and ${\textstyle b}$ are processed through the grating to produce transmitted, ${\textstyle \cos(\phi )}$, and diffracted, ${\textstyle \sin(\phi )}$, components.
3. Electro-optic phase shifters adjust the phases of the optical waves to ensure coherent addition and subtraction, corresponding to the operations required for the rotation matrix.

The final implementation achieves the desired outputs:$$c=a\cos \phi +b\sin \phi ,\quad d=-a\sin \phi +b\cos \phi .$$Mirsalehi et al. proposed using such a Givens device as a building block in lattice filters and wavefront processors. With this in mind, it was already known that such meshes could perform useful operations, but it was not until nearly a decade later, when Reck et al. published their work that these meshes gained significant attention.

Reck and Zeilinger Scheme

Reck et al. showed that a triangular arrangement of ${\textstyle 2{\text{x}}2}$ beam splitters and phase shifters could be systematically programmed, using a straightforward analytical approach, to implement any unitary transformation across a set of optical channels.

Beam splitters:

A ${\textstyle 2{\text{x}}2}$ beam splitter mixes two input modes ${\hat {a}}^{\dagger }$and ${\hat {b}}^{\dagger }$, producing two output modes ${\hat {c}}^{\dagger }$ and ${\hat {d}}^{\dagger }$. The transformation is given by:$${\hat {c}}^{\dagger }=\cos \theta \,{\hat {a}}^{\dagger }+i\sin \theta \,e^{-i\phi }{\hat {b}}^{\dagger },$$and$${\hat {d}}^{\dagger }=i\sin \theta \,e^{i\phi }{\hat {a}}^{\dagger }+\cos \theta \,{\hat {b}}^{\dagger }.$$It is important to note that the above transformations, when written in a matrix form, can be represented as:$${\begin{bmatrix}{\hat {c}}\\{\hat {d}}\end{bmatrix}}={\begin{bmatrix}\cos \left(\theta \right)&ie^{-i\phi }\sin \left(\theta \right)\\ie^{i\phi }\sin \left(\theta \right)&\cos \left(\theta \right)\end{bmatrix}}\cdot {\begin{bmatrix}{\hat {a}}\\{\hat {b}}\end{bmatrix}},$$which is a modified Givens rotation. The transmittance of the beam splitter appears in the above matrix as the term ${\textstyle T=\cos \left(\theta \right)}$.

The universal unitary for ${\textstyle 2{\text{x}}2}$ beam transformations is more commonly written in the following form:$$U(\theta ,\phi )=\exp {\left(i\phi _{c}\right)}\cdot {\begin{bmatrix}\cos \left(\theta \right)\exp {\left(i\phi _{a}+i\phi _{b}\right)}&-\sin \left(\theta \right)\exp {\left(i\phi _{b}\right)}\\\sin \left(\theta \right)\exp {\left(i\phi _{a}\right)}&\cos \left(\theta \right)\end{bmatrix}},$$which is a combination of the modified Givens rotation matrix seen above and three phase-shifters, namely ${\textstyle \exp {\left(i\phi _{a}\right)}}$, ${\textstyle \exp {\left(i\phi _{b}\right)}}$, and ${\textstyle \exp {\left(i\phi _{c}\right)}}$. These are the four free parameters which must be set to fully characterize the ${\textstyle 2{\text{x}}2}$ unitary matrix (as expected, ${\textstyle N^{2}=4}$). The third phase-shifter, ${\textstyle \exp {\left(i\phi _{c}\right)}}$, represents a global offset which can usually be neglected in most practical applications, though it does play an important role when considering geometric phase.

Phase Shifters:

A phase shifter adds a phase ${\textstyle \phi }$ to the state of a photon passing through it. In terms of creation operators, it performs the transformation:

${\hat {a}}_{\text{out}}^{\dagger }=e^{i\phi }{\hat {a}}_{\text{in}}^{\dagger }.$

The same phase ${\textstyle \phi }$ can be achieved by propagating through a material with linear refractive index ${\textstyle n}$ and thickness ${\textstyle d}$, where:

$\phi =n{\frac {2\pi }{\lambda }}d.$

Algorithm:

The objective is to determine the set of matrices ${\textstyle \{T_{N,q}\}}$ such that:$$U(N)\cdot T_{(N,N-1)}\cdot T_{(N,N-2)}\cdots T_{(2,1)}\cdot D=I(N),$$where $N$ and ${\textstyle q=N-1,...,1}$ are the port numbers in the triangular mesh. The matrix ${\textstyle T_{N,q}}$ is a modified Givens rotation matrix. That is, it is an identity matrix with the elements ${\textstyle T_{N,q}}$

Step 1: Initial multiplication

Multiply ${\textstyle U(N)}$ from the right by a succession of matrices ${\textstyle T_{(N,q)}(\omega _{N,q},\phi _{N,q})}$ for ${\textstyle q=N-1,...,1}$. This is where the matrix ${\textstyle T_{(N,q)}}$ is an ${\textstyle N}$-dimensional identity matrix with the elements $I_{NN},I_{Nq},I_{qN},$ and $I_{qq}$ replaced by the corresponding ${\textstyle 2{\text{x}}2}$ beam transformation matrix elements. Hence, it represents a modified Givens rotation matrix.

By the properties of the Givens rotation matrix, ${\textstyle \omega _{Nq}}$ and ${\textstyle \phi _{Nq}}$ in ${\textstyle T_{(N,q)}(\omega _{N,q},\phi _{N,q})}$ can be chosen such that, upon multiplication with ${\textstyle U(N)}$, the resulting matrix element at ${\textstyle (N,q)}$ vanishes. Changing the index ${\textstyle q\rightarrow q'<q}$ and performing another multiplication with specially chosen values of ${\textstyle \omega _{N,q'}}$ and ${\textstyle \phi _{N,q'}}$, the resulting matrix element at ${\textstyle (N,q')}$ vanishes. Repeating successive multiplications until the index ${\textstyle q=1}$ is reached will result in the last row vanishing (expect the on-diagonal element which remains 1). Due to the unitarity of each transformation, the rightmost column will also vanish (again, expect the on-diagonal element which remains 1). This step reduces the effective dimension of ${\textstyle U}$ to ${\textstyle N-1}$.$$U(N)\cdot T_{(N,N-1)}\cdot T_{(N,N-2)}\cdots T_{(N,1)}={\begin{bmatrix}{\begin{bmatrix}U(N-1)\end{bmatrix}}&0\\0&e^{i\alpha }\end{bmatrix}}.$$Step 2: Recursive multiplication

Multiply the reduced matrix from the right by a succession of matrices ${\textstyle T_{(N-1,q)}(\omega _{N-1,q},\phi _{N-1,q})}$ for ${\textstyle q=N-2,...,1}$. Following the same thought-process as in step 1, this will result in the second-to-last row vanishing and by unitarity, the second-to-rightmost column vanishing (except for the on-diagonal element). The resulting reduced matrix is of the following form:$$U(N)\cdot T_{(N,N-1)}\cdot T_{(N,N-2)}\cdots T_{(N,1)}\cdot T_{(N-1,N-2)}\cdot T_{(N-1,N-3)}\cdots T_{(N-1,1)}={\begin{bmatrix}{\begin{bmatrix}U(N-2)\end{bmatrix}}&0&0\\0&e^{i\beta }&0\\0&0&e^{i\alpha }\end{bmatrix}}.$$Repeating this step in a recursive fashion until the matrix multiplication involves ${\textstyle T_{(2,1)}(\omega _{2,1},\phi _{2,1})}$ will result in a transformed diagonal matrix. Notice that the elements along the diagonal have modulus 1.

Step 3: Recovering the unitary

The final step is to separate the unitary ${\textstyle U(N)}$ from the successive ${\textstyle U(2)}$ transformations. This is accomplished by multiplying the transformed diagonal matrix by another diagonal matrix ${\textstyle D}$ whose elements are also modulus 1 such that the outcome is the identity matrix:$$U(N)\cdot T_{(N,N-1)}\cdot T_{(N,N-2)}\cdots T_{(2,1)}\cdot D=I(N).$$In practice, ${\textstyle D}$ represents a set of phase shifters that compensate for the phases appearing along the diagonal of the transformed matrix.

By the properties of the identity matrix, the product of the final transformed matrix and ${\textstyle D}$ represents the inverse of ${\textstyle U(N)}$,$$U(N)=\left(T_{(N,N-1)}\cdot T_{(N,N-2)}\cdots T_{(2,1)}\cdot D\right)^{-1}.$$

Variants and similar decompositions

Maybe add this section. I haven't decided yet.