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 $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'$ from the evolved state $(S\bar x, S\sigma S^\top)$.
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}t_jI_2 & r_jI_2 \\ -r_jI_2 &t_jI_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.
