Composing beam splitters

Question

Let $a, b$ and $c$ be independent modes in a system $S$ and in environments $E_1$, $E_2$ respectively. Suppose $a$ goes through a beam-splitter characterized by a parameter $\theta$ coupling it to mode $b$, so that first this first interaction we may write the unitary $$U_\theta = \exp(i\theta(a^\dagger b + b^\dagger a)) $$ (I'm forgetting about relative phases, global signs and what-not; this should work for any kind of beam splitter, polarizing or otherwise).

Now suppose that immediately after the action $U_\theta$, the system modes (which evolved as $a'=U_\theta^\dagger aU_\theta$) go into another beam splitter with parameter $\eta$, this time coupling $S$ to the mode $c$ in $E_2$: this time, the interaction may be described by another unitary of the form $$U_\eta= \exp(i\eta(a'^\dagger c + c^\dagger a')). $$ Intuitively, if we compose the two transformation we should get a single unitary $U_{f(\theta, \eta)}=U_\eta$ parametrized by some function $f$ of the original parameters, acting on the mode $a$ and the tensor product $b\otimes c$. However, since the exponents of $U_\theta$ and $U_\eta$ do not commute, I'm not finding a straightforward way to derive the explicit expression for $U_{f(\theta, \eta)}$. Can someone provide an answer?

---

I've looked into applying the BCH formula to this problem, and I think I might be onto something. The first few terms of the expansion (I'm following the boxed identity at Wikipedia's page) other than the one proportional to $1/24$ are nonzero and give for the exponent $Z$ the sum \begin{equation}\tag{*} Z:=i\left(f(\theta, \eta)(a^\dagger b + b^\dagger a) + f(\eta, \theta)(a^\dagger c + c^\dagger a)\right) -i\frac{\theta\eta}{2}(b^\dagger c - c^\dagger b) \end{equation} where $f(x,y):=x-\frac{1}{12}xy^2-\frac{1}{720}xy^4-\frac{1}{360}x^3y^2$. I am not sure about my calculations (those pesky signs!), but it seems to me that most higher terms involve at some point the commutator between $a^\dagger b + b^\dagger a$ (or $a^\dagger c + c^\dagger a$) and itself, and thus do not contribute to $(*)$, so there may be hope to simplify the whole thing.

Answer 1

The easiest way to deal with concatenated beam splitters is to evaluate their action on creation or annihilation operators as $$U\begin{pmatrix}a\\b\\c\end{pmatrix}U^\dagger=\mathsf{U}\begin{pmatrix}a\\b\\c\end{pmatrix}$$ for some unitary matrix $\mathsf{U}$ (belonging to SU(3)). That is, each of $a$, $b$, and $c$ transform into a linear combination of $a$, $b$, and $c$ following a beam splitter.

For example, taking $U_\theta=\exp(i\theta(a^\dagger b+ a b^\dagger))$, we have $$\mathsf{U}_\theta=\begin{pmatrix}\cos\theta&-i\sin\theta&0\\-i\sin\theta&\cos\theta &0\\0&0&1\end{pmatrix}$$ such that $a\to a\cos\theta-b\sin\theta$, etc. ($c\to c$ because the $\theta$ beam splitter doesn't affect mode $c$; the minus sign might be in the wrong spot but I think it's correct). This is known as a beam splitter matrix.

Now the nice thing is that matrix multiplication is straightforward. So the composition of two unitaries (that effect linear optical transformations according to the group SU(3)) is $$U_\eta U_\theta\begin{pmatrix}a\\b\\c\end{pmatrix}U^\dagger_\theta U^\dagger_\eta= U_\eta\left[\mathsf{U}_\theta\begin{pmatrix}a\\b\\c\end{pmatrix}\right]U^\dagger_\eta=\mathsf{U}_\eta\mathsf{U}_\theta\begin{pmatrix}a\\b\\c\end{pmatrix}.$$ All you need is matrix multiplication $\mathsf{U}_\eta\mathsf{U}_\theta$ to discover how every operator evolves. In our case, $$\mathsf{U}_\eta=\begin{pmatrix}\cos\eta&0&-i\sin\eta\\0&1&0\\-i\sin\eta&0&\cos\eta \end{pmatrix},$$ so you just need to multiply together the two matrices.

Why are these input-output relations useful? First, if you are evaluate expectation values of operators for the evolved state, you can use the Heisenberg picture to evolve the operators like $a$ and $b$ instead of evolving your state. Second, if you write your state as some analytic function of creation operators acting on the [three-mode] vacuum, the unitaries move around nicely: $$U [f(a^\dagger,b^\dagger,c^\dagger)|\rm{vac}\rangle]=f(Ua^\dagger U^\dagger,Ub^\dagger U^\dagger,Uc^\dagger U^\dagger)U|\rm{vac}\rangle,$$ and $U|\rm{vac}\rangle=|\rm{vac}\rangle$ for linear optical unitaries, so knowing how the operators update can tell you how any state updates.

As another tack, you are allowed to just say $U=U_\eta U_\theta$ and that is an explicit expression for a unitary as a function of the two parameters. You can use the Baker Campbell Hausdorff formula to rewrite this as a beam splitter first acting on modes $b$ and $c$ followed by $U_\eta$, but that doesn't help much in general other than for specific initial states.

Answer 2

Just to make sure that I understand the question, I represent the process in terms of unitaries applied to an initial (pure) state: $$ |\psi_{fin}\rangle = U_{ac} U_{ab} |S\rangle|E_1\rangle|E_2\rangle . $$ If the initial states are not pure, we can replace them by density operators and apply the adjoint unitaries at the back. But, we are interested in the channel which is given here as the product of the two unitaries. So the question is $$ U_{ac} U_{ab} = \exp[i\theta(a^{\dagger}b+b^{\dagger}a)] \exp[i\eta(a^{\dagger}c+c^{\dagger}a)] = ? $$ Note that $a$ remains $a$ after the first beam splitter.

One can try to use Baker-Campbell-Haussdorff, but the problem is that every time we compute a commutator, we end up with another operator that doesn't commute with original two. So it is necessary to a perform sequence of such commutations and hope that they form a pattern that can be summed.

Another approach is to produce a differential equation for the result. It is done by inserting a variable $t$ into the exponents of the operators. We make an educated guess for the form of the final result. So the equation is $$ \exp[it\theta X] \exp[it\eta Y] = \exp[if(t)X +ig(t)Y +ih(t)Z] , $$ where $X=a^{\dagger}b+b^{\dagger}a$, $Y=a^{\dagger}c+c^{\dagger}a$, $Z=b^{\dagger}c-c^{\dagger}b$, and $f(t)$, $g(t)$ and $h(t)$ are unknown functions to be determined. They all become zero for $t=0$. The form of the result is based on the commutations $[X,Y]=Z$, $[Z,X]=-Y$, and $[Y,Z]=-X$, which means that these combinations form a closed algebra.

Next, we take the derivative with respect to $t$, and then remove as many of the unitary operators by applying the adjoints of the respective operators on the right hand side of both sides of the equation. The resulting equation may need some remaining commutations of the form $\exp(X)Y\exp(-X)$ to be evaluated, and can then be separated into different differential equations for the three functions. By solving the differential equations with the boundary condition of being zero at $t=0$, we end up with the functions. The final values for the operator is obtained by setting $t=1$. I am not going to work through this calculation, unless there is a real need to show how it goes.