Quantum light with finite efficiency using a “beam splitter model”

Detectors (in my case I'm interested in homodyne detectors) with imperfect efficiency and losses are said to be able to be modeled using a beamsplitter and are to be in a statistical mixture with vacuum.

I am struggling to reproduce this "statistical mixture" result by modeling the system as a beamsplitter.

Just to overview how the homodyne works without losses: A homodyne detector measures the subtracted intensity output of two detectors. Since a quantum beamsplitter returns $c^\dagger = a^\dagger + b^\dagger$ and $d^\dagger = a^\dagger - b^\dagger$.

The quantum operator becomes $c^\dagger c$ - $d^\dagger d$ = $a^\dagger b + b^\dagger a$. With a coherent state as our mode-b: this reduces to being $|\beta| X(\theta)$.

Where $X(\theta)$ is the "generalized quadrature", a quantum observable of the form:

$$X(\theta) = a^\dagger e^{i\theta} + a e^{-i \theta}$$

Theta($\theta$) is the relative phase between the two fields (or can be thought as the complex phase associated with $\beta$.

First I want to understand how losses in the quantum channel ($a^\dagger$) affect our output observable.

So to do this I am going to have my ($a^\dagger$) input mode to interfere with an additional mode (which will remain as vacuum) through a beamsplitter. Writing the beamsplitter's transmission as $\eta$ (my losses) and "reflection" as (1-$\eta$) we have my output state (after losses) as: \begin{bmatrix} \eta & \sqrt{\eta(1-\eta)} \\ \sqrt{\eta(1-\eta)} & (1-\eta)\\ \end{bmatrix}

My output state is:

$|\psi \rangle = \sqrt{\eta}|1, 0\rangle + \sqrt{1-\eta} |0, 0\rangle$

Now I think the way to move forward is to partial trace out the second vacuum mode, so that I can obtain a "single-mode" that I can use as the new quantum input-mode of the homodyne. So I'm going to get the density matrix, and trace it over this second vacuum mode.

The tensored density-matrix for this two-mode state is:

$|\psi \rangle \langle \psi| = (\sqrt{\eta}|1, 0\rangle + \sqrt{1-\eta} |0, 0\rangle) (\sqrt{\eta}\langle 1, 0| + \sqrt{1-\eta} \langle 0, 0|)$

$= \eta|1, 0\rangle \otimes \langle 1, 0| + \sqrt{(1-\eta) \eta} (|0, 1\rangle \otimes \langle 0, 0| + |0, 0\rangle \otimes \langle 1, 0|) + (1-\eta) |0, 0\rangle \otimes \langle 0, 0|$

You can write this out and sum over one of the two dimensions. Tracing over the second mode, it looks like I just end up back at the normal density matrix of a superposition state (I don't get a mixture):

\begin{bmatrix} \eta & \sqrt{\eta(1-\eta)} \\ \sqrt{\eta(1-\eta)} & (1-\eta)\\ \end{bmatrix}

This suggest to me that I can just perform a substitution and swap $a^\dagger$ with $\sqrt{\eta} a^\dagger|0\rangle + \sqrt{1-\eta} |0\rangle$

I'm not sure if this is result is consistent with the literature. If I'm understanding any losses are "inefficiencies" which can all be multiplied together as an effective mixture with vacuum. I believe in this resource they include one-mode losses in their measurement of efficiency, which they model as a statistical mixture of the state with vacuum.

Two-mode losses, on the other hand, in which both outputs of the two ends of the detector have losses, I think works out differently. Doing the beam splitter model twice on modes c and d and extracting terms proportional to b (all others are small relative to the local oscillator), I get something looking like:

$$\eta \beta X(\theta) + \sqrt{\eta(1-\eta)}((L1*c+c*L1)-(L2*d-d*L2)$$

Where I've made L1 and L2 to be loss channels of c and d respectively, using the same beam-splitter model as above (and assuming both have the same losses). I'm a bit stuck in how to work out that second term. Are the two losses correlated? Can I combine terms? Does this end up working out to be a statistical mixture?

To quote:

Remarkably, all imperfections of the experiment losses in transmission of the signal photon, quantum efficiency of the HD, trigger dark counts, mode matching of the signal photon and the local oscillator, and spatiotemporal coherence of the signal photon had a similar effect on the reconstructed state: admixture of the vacuum 0 to the ideal Fock state 1.

Finally, I'm also trying to see how this can be used to find out how this changes the probabilities of measuring certain quadrature values (given what the quantum state is. I understand the normal projection operator provides the probabilites in the fock-state basis (which can be found here).

With this statistical mixture model, it's said that this projection operator can be modified to look like:

$$\sum_{n,m} B_{m+k, m}(\eta) B_{n+k, n}(\eta) \langle n| \theta, x \rangle \langle \theta, x| m \rangle |n+k\rangle \langle m+k|$$

I'm struggling to see how to derive this (particularly the sums).