Extending pairwise beam splitters operators to an overall transformation

Question

Let $\bar x_s, \sigma_s$ be the vector of first moments and the covariance matrix of an $N$-mode Gaussian state $\rho_s$. If $\rho_s$ is paired with an empty environment of similar dimension, the total state is described by the moments \begin{equation} \bar x = \begin{pmatrix} \bar x_s \\ 0_{\scriptstyle{N}} \end{pmatrix}, \quad \sigma = \begin{pmatrix} \sigma_s & 0 \\ 0 & \sigma_e \end{pmatrix} \end{equation} where $\sigma_e = \frac12 I_{2N}$. Now consider the action of a family of beam splitters $\{S_j\}$ coupling each system mode $j$ with the corresponding $j$-th environmental mode with transmittance $\sqrt{t_j}$ and reflectivity $r_j$. Since we are only interested in the system, we take the partial trace right after the transformation, in order to obtain $\bar x_s'$ and $\sigma_s'$.

My question is: what is the correct expression for the evolved subsystem? My main doubt arises from the way I should be writing the total transformation $S$. Each beam splitter mixing two modes can be individually written as a matrix $$S_j = \begin{pmatrix}\sqrt{t_j}I_2 & \sqrt{r_j}I_2 \\ -\sqrt{r_j}I_2 & \sqrt{t_j}I_2\end{pmatrix}$$ and I originally wrote $S=\bigoplus_j S_j$, but I think I'm making a mistake since that would correspond to the interleaved ordering of the quadratures $$(q_1^s, p_1^s, q_1^e, p_1^e,..., q_N^s, p_N^s, q_N^e, p_N^e)$$ instead of the ordering $$(q_1^s, p_1^s, ..., q_N^s, p_N^s, q_1^e, p_1^e,...,q_N^e, p_N^e)$$ implied by the original form of $\sigma$. I tried writing down the matrices explicitly for $N=2$ but the corrected version doesn't seem to admit a nice closed form for general $N$, making the computation of the partial trace somewhat harder.

---

Now consider the general case $\sqrt t_j\in \mathbb C$. Write $\sqrt t_j = |\sqrt t_j|e^{i\phi_j}$; then the total transformation should be equivalent to a phase shift $F = \bigoplus_{j=1}^N R(\phi_j)$ with $$R(\phi_j) = \begin{pmatrix}\cos\phi_j & \sin\phi_j \\ -\sin\phi_j & \cos\phi_j \end{pmatrix} $$ followed by a pure loss with real transmittances $|\sqrt t_j|$, yielding $$\sigma_s' = (DF)\sigma(DF)^\top + \frac12(I - D^2)$$ where now we define $D := \text{diag}(|\sqrt t_1|, |\sqrt t_1|,...,|\sqrt t_N|, |\sqrt t_N|)$.

Answer 1

Using formulae from the quick reference, a beam splitter acts with matrix (their $F$ is OP's $S$) $$S=\begin{pmatrix}\sqrt{t}{I}_2&\sqrt{r}{I}_2\\-\sqrt{r}{I}_2&\sqrt{t}{I}_2\end{pmatrix}$$ in the basis of $\{q^s,p^s,q^e,p^e\}$. Given that the covariance matrix for those modes begins as a block-diagonal matrix with $\sigma^s$ and $\sigma^e$ on the diagonal, we have the updated covariance matrix $S\begin{pmatrix}\sigma^s &0\\0&\sigma^e\end{pmatrix}S^\dagger=\begin{pmatrix}t \sigma^s+r \sigma^e&\sqrt{tr}(\sigma^e-\sigma^s)\\\sqrt{tr}(\sigma^e-\sigma^s) &r\sigma^s+t\sigma^e\end{pmatrix}$, from which we trace out the environment to find $\sigma^s\to t\sigma^s+r\sigma^e$ (the 0 parts are 0 matrices of appropriate size). What happens to the multimode version?

Acting with a beam splitter on one mode and then tracing out the corresponding environment is an operation that commutes with the same being done on another mode with another environment. As such, we can inspect the individual action of each of these operations and put them together to achieve the cumulative action.

When this individual operation acts on a multimode system, it is convenient to include only the one environment and perform this operation on the final mode; the rest will be taken care of accordingly. Since it leaves the other modes unchanged, we write $$S_N=\begin{pmatrix}I_{2N-2}&0&0\\ 0&\sqrt{t_N}{I}_2&\sqrt{r_N}{I}_2\\0&-\sqrt{r_N}{I}_2&\sqrt{t_N}{I}_2\end{pmatrix}$$ in the basis of $\{q_1^s,p_1^s,q_2^s,p_2^s,q_3^s,p_3^s\cdots,q_N^s,p_N^s,q_N^e,p_N^e\}$. Act with this one the covariance matrix of the system and its one environment $\Sigma=\begin{pmatrix}\sigma^s&0\\0&\sigma^e\end{pmatrix},$ where we break $\sigma^s=\begin{pmatrix}A&B\\B^\dagger&C\end{pmatrix}$ into blocks of size $2N-2\times 2N-2$, $2N-2\times 2$, and $2\times 2$ for $A,B,C$ respectively. Then \begin{equation} S_N \Sigma S_N^\dagger=\left( \begin{array}{ccc} A & B \sqrt{t_N} & -B \sqrt{r_N} \\ B^\dagger \sqrt{t_N} & C t_N+r_N \sigma^e & \sqrt{r_N} \sqrt{t_N} \left(\sigma^e-C\right) \\ -B^\dagger \sqrt{r_N} & \sqrt{r_N} \sqrt{t_N} \left(\sigma^e-C\right) & C r_N+t_N \sigma^e \\ \end{array} \right) \end{equation} and ignoring the environment yields \begin{equation} \sigma^s=\begin{pmatrix}A&B\\B^\dagger &C\end{pmatrix}\to\left( \begin{array}{ccc} A & B \sqrt{t_N} \\ B^\dagger \sqrt{t_N} & C t_N+r_N \sigma^e \\ \end{array} \right). \end{equation} What is the upshot? After mode $N$ has this operation, the $2\times 2$ part of the covariance matrix corresponding to that mode gets multiplied by $t_N$ and then $r_N I_2/2$ is added to it. As well, all of the off-diagonal components of the covariance matrix with one of the indices corresponding to either $q^s_N$ or $p^s_N$ get multiplied by $\sqrt{t_N}$. The same will occur when any of the other modes has this operation.

Indeed, collecting the results, the covariance matrix of the system transforms as $$\sigma^s\to D\sigma^s D^\dagger+(I_{2N}-D^2)/2$$ for $D=\text{diag}(\sqrt{t_1},\sqrt{t_1},\sqrt{t_2},\sqrt{t_2},\cdots,\sqrt{t_N},\sqrt{t_N})$. Each diagonal $2\times 2$ block gets updated is it would in the absence of other modes, while each off-diagonal block gets multiplied by a corresponding $\sqrt{t_i t_j}$.