A symmetric 50:50 beamsplitter can be represented by a unitary 2x2 matrix acting on a vector of creation operators for the input states:
$$ \begin{pmatrix} \hat{a}_c^\dagger\\ \hat{a}_d^\dagger \end{pmatrix} = \begin{pmatrix} \frac{i}{\sqrt{2}} & \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}} & \frac{i}{\sqrt{2}} \end{pmatrix}\begin{pmatrix} \hat{a}_a^\dagger\\ \hat{a}_b^\dagger \end{pmatrix}.$$
This gives rise to interesting effects such as the Hong-Ou-Mandel effect where a (1,1)-output is impossible for a (1,1)-input.
Is it possible to generalize this two-mode beamsplitter to a three-mode beamsplitter where each input photon has the same probability to end up in any of the three output modes?
My own attempt: I tried to construct a 3x3 matrix representation myself following the principles of the two-mode beamsplitter but was unsure how to make it unitary. The way I interpret the two-mode matrix is each element $U_{i,j}$ is the amplitude for a photon in mode $i$ to end up in mode $j$, meaning for example that the only difference between $U_{1,1}$ and $U_{1,2}$ is a $\frac{\pi}{2}$ phase $\exp(i\frac{\pi}{2})=i$, since we are dealing with a symmetric beamsplitter. I wanted to use roots of unity for the 3x3 beamsplitter matrix, to represent an equal probability for a photon to end up in any of the output modes, meaning only a phase difference between their amplitude:
$$\begin{pmatrix}\frac{1}{\sqrt{3}} & \frac{e^{i\frac{2\pi}{3}}}{\sqrt{3}} & \frac{e^{i\frac{4\pi}{3}}}{\sqrt{3}}\\ a & b & c\\ d & e & f\end{pmatrix}.$$
However, from this point on I was unable to fill out the rest of the matrix as I realised I had no idea if this was going to be a unitary matrix or not.
