Matrix representation of beamsplitter, for numerical computation of output based on given photon number (Fock state) input

Question

Question: Is it possible to express the effect of a simple 50% beamsplitter on photon number states using matrices, such that the output can be computed by matrix calculations rather than manual substitution of equations?

To explain the problem, consider a 50% beamsplitter and define: $a_{1,2}^{+}$ = creation operators for input and $b_{1,2}^{+}$ = creation operators for the output. For a balanced 50% beamsplitter, we can write how the beamsplitter maps the input to output creation operators:

$$a_1^{+} \rightarrow \frac{1}{\sqrt{2}} \left( b_1^{+} + b_2^{+} \right)$$ $$a_2^{+} \rightarrow \frac{1}{\sqrt{2}} \left( b_1^{+} - b_2^{+} \right)$$

We can then compute the output photon number distribution for given input photon number states by writing the Fock states in terms of creation operators, followed by substituting in the above equations.

For example, consider 2 photons incident on the beamsplitter in MODE1 and 1 photon incident in MODE2 at the input. We write them first as photon number states, then express in terms of creation operators acting on the vacuum.

$$ |\psi_{in}\rangle = |2\rangle_1|1\rangle_2 = \frac{1}{\sqrt{2}} a_1^{+2} a_2^{+}|0\rangle_1|0\rangle_2 $$

We can find the output state by rewriting $a_1^{+}$ & $a_2^{+}$ in terms of $b_1^{+}$ & $b_2^{+}$ (using the above equations):

$$ |\psi_{out}\rangle = \frac{1}{2\sqrt{2}} \left( b_1^{+} + b_2^{+} \right)^2 \left( b_1^{+} - b_2^{+} \right) |0\rangle_1|0\rangle_2 $$

which can then be evaluated by expanding the brackets and evaluating the effect of creation operators on the vacuum (i.e. $b_1^{+} |0\rangle_1 = |1\rangle_1$ and $b_1^{+2} |0\rangle_1 = \sqrt{2}|2\rangle_1$ etc.). The result of this example is thus:

$$ |\psi_{out}\rangle = \frac{1}{4} \left( \sqrt{6} |3\rangle_1 |0\rangle_2 + \sqrt{2} |2\rangle_1 |1\rangle_2 - \sqrt{2} |1\rangle_1 |2\rangle_2 - \sqrt{6} |0\rangle_1 |3\rangle_2 \right) $$

The above method works, but requires algebraic equation substitution, which is not a mathematical operation. My question is, how can we frame the problem so that the input could be entered as a vector/matrix and the output (e.g. as a vector to show the coefficients of each possible photon number output state) calculated by matrix operations, which are well suited for performing on a computer. I've commonly seen the beamsplitter effect expressed as a matrix: e.g.

$$B = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix}$$ [ref: https://physics.stackexchange.com/questions/330888/how-to-write-the-output-state-of-a-beam-splitter]

so one could potentially write:

$$\begin{pmatrix} b_1^{+} \\ b_2^{+} \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix} \begin{pmatrix} a_1^{+} \\ a_2^{+} \end{pmatrix} $$

but this only expresses the relationship in terms of single creation operators and it's unclear how to formulate a matrix problem to support arbitrary inputs (i.e. where each input may contain multiple photons).

Possible Solutions: From my research I note:

• a similar question has been asked before (Output of a beamsplitter with photon number (Fock) state inputs), but the answer talks about substituting expressions that relate input to output creation operators - i.e. this is not a simple matrix operation.
• a paper that looks promising is: https://journals.aps.org/pra/abstract/10.1103/PhysRevA.40.1371 (unpaywalled: http://people.bu.edu/teich/pdfs/PRA-40-1371-1989.pdf) but this is very complex and fully general, and I can't work out how to apply this / if this is a valid approach for me.
• A Python package QuTip is available which can handle Fock states and creation operators as matrices/vectors, so perhaps this would be a good framework for implementing this?

All suggestions or ideas are greatly appreciated, thanks.

Answer 1

I'm not completely sure what you mean but I think this is easily done.

What you need to remember that is that your beam splitter acts on products of creation operators as a product transformation. Thus, if your beamsplitter is represented by the matrix \begin{align} U=\left(\begin{array}{cc} U_{11}&U_{12}\\ U_{21}&U_{22} \end{array}\right) \tag{1} \end{align} then: \begin{align} a_i^\dagger &\to \sum_{j}a_j^\dagger U_{ji}\, , \\ a_1^\dagger &\to a_1^\dagger U_{11} +a_2^\dagger U_{21} \, , \tag{2} \\ a_2^\dagger &\to a_2^\dagger U_{12} +a_2^\dagger U_{22} \end{align} and thus on a product state \begin{align} \vert \hbox{in}\rangle =\frac{1}{\sqrt{2}} a_1^\dagger a_1^\dagger a_2^\dagger \vert 0\rangle \end{align} you get the output \begin{align} \vert\hbox{out}\rangle =\frac{1}{\sqrt{2}}\sum_{ijk}a_i^\dagger a_j^\dagger a_k^\dagger U_{i1}U_{j1}U_{k2}\, . \tag{3} \end{align} Note that nothing limits (1) to $2\times 2$ matrices. In general the indices $ijk$ in (3) will run over out modes so in a $3$-channel interferometer they would run from $1$ to $3$ for instance.

In your specific example the coefficient of the properly normalized state output state $\frac{1}{\sqrt{3!}}\vert 3\rangle_1\vert 0\rangle$ is found by restricting the sum to all the combinations of $(ijk)$ that multiply the factor $(a_1^\dagger)^3$. Here these are obviously $(ijk)=(1,1,1)$ so the coefficient is $\frac{1}{\sqrt{2}}U_{11}U_{11}U_{12}=\frac{1}{4}$. Thus, the coefficient of $\vert 3\rangle_1\vert 0\rangle_2=\sqrt{6}/4$.

The coefficient of the properly normalized state $\frac{1}{\sqrt{2}}\vert 2\rangle_1\vert 1\rangle_2$ will be a sum of terms with $(ijk)=(1,1,2), (1,2,1)$ or $(2,1,1)$: \begin{align} \frac{1}{\sqrt{2}}\left(U_{11}U_{11}U_{22}+ U_{11}U_{21}U_{12}+U_{21}U_{11}U_{12}\right) =\frac{1}{\sqrt{2}}\frac{1}{2\sqrt{2}}\left(1+ (-1)+(-1) \right)= -\frac{1}{4} \end{align} so the coefficient of $\vert 2\rangle_1\vert 1\rangle_2$ is $-\frac{\sqrt{2}}{4}$. I have a sign difference with you probably because of the way I place my indices in (2): your definition of the transformation might have the indices $(ij)$ flipped w/r to my way of writing things. If you flip indices you'd get $\frac{\sqrt{2}}{4}$.

Note that, in my notation, the creation operators are mapped to vectors \begin{align} a^\dagger_1\vert 0\rangle \to \vert 1\rangle &=\left(\begin{array}{c} 1 \\ 0\end{array}\right)\, ,\\ a^\dagger_2\vert 0\rangle\to \vert 2\rangle &=\left(\begin{array}{c} 0 \\ 1\end{array}\right)\, , \end{align} so that \begin{align} U\vert 1\rangle &= \vert 1\rangle \langle 1\vert U\vert 1\rangle + \vert 2\rangle \langle 2\vert U\vert 1\rangle\, ,\\ &= \vert 1\rangle U_{11} +\vert 2 \rangle U_{21}\, ,\\ &= a^\dagger_1 \vert 0\rangle U_{11} + a^\dagger_2\vert 0\rangle U_{21} \end{align} so that indeed (2) looks correct.

There is more...

You can simplify the problem further by using the Schwinger representation (or Jordan map)

so that the angular momentum state \begin{align} \vert jm\rangle \Leftrightarrow \frac{(a_1^\dagger)^{j+m}(a_2^\dagger)^{j-m}}{\sqrt{(j+m)!(j-m)!}}\vert 0\rangle\, . \end{align} Thus your input state is mapped to the angular momentum state \begin{align} \vert\psi_{in}\rangle \to \vert \textstyle\frac{3}{2},\frac{1}{2}\rangle\, . \end{align} Your beamsplitter transformation then corresponds (see Yurke, B., McCall, S.L. and Klauder, J.R., 1986. SU (2) and SU (1, 1) interferometers. Physical Review A, 33(6), p.4033.) to a rotation $U=R_z(0)R_y(\pi/2)R_z(0)$ so that you can then, in a single shot, write \begin{align} U\vert\psi_{in}\rangle &\to R_y(\pi/2)\vert \textstyle\frac{3}{2},\frac{1}{2}\rangle\, ,\\ &=\sum_m \vert \textstyle\frac{3}{2},m\rangle D^{3/2}_{m,1/2}(0,\pi/2,0) \end{align} where $D^j_{mm'}(\alpha,\beta,\gamma)$ is a Wigner D-matrix .

It's then easy to get the output in terms of angular momentum states: \begin{align} \vert \textstyle\frac{3}{2},\frac{3}{2}\rangle \left(-\frac{\sqrt{3}}{2\sqrt{2}}\right) +\vert \textstyle\frac{3}{2},\frac{3}{2}\rangle \left(-\frac{1}{2 \sqrt{2}}\right)\vert \textstyle\frac{3}{2},\frac{1}{2}\rangle +\frac{1}{2 \sqrt{2}}\vert \textstyle\frac{3}{2},-\frac{1}{2}\rangle + \vert \textstyle\frac{3}{2},-\frac{3}{2}\rangle \left(\frac{\sqrt{3}}{2\sqrt{2}}\right) \end{align} and then convert back to occupation number states.