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

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:

All suggestions or ideas are greatly appreciated, thanks.

• Sure it is possible, just work it out by substituting the action on the creation/annihilation oeprators as you discuss above into the Fock states. Does this answer your question? Or are you after a simple closed expression for such a transformation? (Educated guess: It can't be that hard to guess it and prove it by induction or the like.) – Norbert Schuch Jun 6 at 19:21
• Thanks for the suggestion. Glad to hear it's possible, but are you able to show an example please as I can't quite see how I'd do that. I'm looking for a matrix expression / approach, such that I could code it up in a common programming language and then obtain the output photon number distribution for any input photon number distribution. e.g. if I put input as 2 photons in MODE1 and 1 photon in MODE2 like my example, I'd want the output to let me see the coefficients of each possible output photon number state (so I can square them to find the probability of occuring). Thanks – Qconfused9102 Jun 6 at 19:28
• Well, for the case of at most one photon in each incident beam, you have worked it out yourself above! (Sorry, checked again: Not exactly what you did, but you can immediately change what you did for that case.) – Norbert Schuch Jun 6 at 19:36
• I've tried working it out but can't get a single expression for the output. Specifically, the beamsplitter is a 2 x 2 matrix, and acts on a (2 x 1) vector, thus resulting in a (2 x 1) vector. This means that the b1 and b2 operators are separated. However, some terms of the output are known to require BOTH b1 and b2 together, e.g. |1>|2>, so how can a (2 x 1) output vector contain this information? As there are 4 possible output states (|3>|0>, |2>|1>, |1>|2> and |0>|3>, do I need to somehow extend the dimension of my problem? Sorry if I'm missing something obvious? Thanks – Qconfused9102 Jun 6 at 19:51
• Indeed, you want to extend at least you output space to the full space you can span with 3 photons, i.e. the 4 states you mention. (You can also extend the input space to that dimension, but you don't have to.) Each photon can go either way, so a priori, you can get any number of photons in any output (unless there is destructive interference such as in HOM.) – Norbert Schuch Jun 6 at 19:58

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.

• Thanks for the very helpful response. First approach seems to work nicely (2nd one using Dmatrix is nice too but I need to read more to understand angular momentum states mappings). My only remaining question is: how can we use this approach with inputs in superposition states. Let's consider an example where an input photon is in a superposition state across both input modes: |INPUT> = 1/sqrt(2) * (|0>|1> +|1>|0>). I ultimately want to scale this analysis to consider more modes (2 time bins + 2 input paths) so I guess I can just use more indices over more modes as you suggest. Thanks – Qconfused9102 Jun 7 at 17:37
• @Qconfused9102 Best place to read on the Schwinger rep. is Gordon Baym's Lectures on Quantum Mechanics. Yes handing linear combos of arbitrary states becomes a burden if you do the multiplication method, and you have act on each part of the combo separately, then keep track of the coefficient of the combos. From that perspective the $D$-function trick helps as the example in your comment would be like rotating $\frac{1}{\sqrt{2}} \left(\vert 1/2,1/2> + \vert 1/2,-1/2>\right)$. – ZeroTheHero Jun 7 at 18:23
• thanks again for this. Out of interest, do you think this approach could work with coherent states too? i.e. replacing creation operators with displacement operators... in either the initial or Jordan-Schwinger approach? – Qconfused9102 Jun 10 at 19:42
• This Schwinger trick is basically for $SU(n)$ so if you are generating an angular momentum coherent state then yes. In the case of the standard coherent state it is possible to expand the 2-mode coherent state into $SU(2)$ subspace basically using the Jordan map, i.e. take a state $\vert n_1n_2\rangle$ in your 2-mode coherent state and think of this as an angular momentum state, but otherwise I'm not 100% of what you man by "replacing creation operators with displacement operators". You could presumably replace rotations by displacements... – ZeroTheHero Jun 10 at 20:41
• ... and then you'd get $D$-functions for the group of translations (they are Bessel functions I think). – ZeroTheHero Jun 10 at 20:41