t}+\delta(t)\right)\ast\hat{x}_{\text{in}}(t)\right)+\sqrt{2\kappa}\hat{X}_{\text{\tiny{M}}}(0)\Theta_{\kappa_{\text{\tiny{M}}}}(t)+\text{Noise}.$
(97)

We assume a rectangular pulse. Thus

$\displaystyle\hat{\mathbf{X}}^{\prime}{}^{\text{out}}=\frac{1}{\sqrt{\tau}}\int_{0}^{\tau}\hat{x}^{\prime}_{\text{out}}(t)dt.$
(98)

To obtain a symmetric gate we need to shift one of the quadratures. Thus,
after detection we shift of the atomic quadrature by the feedforward procedure
with $\mathbf{K}_{\text{f}}$ coefficient:

$\displaystyle\hat{X}_{\text{\tiny{A}}}(\tau)\rightarrow\hat{X}_{\text{\tiny{A}}}(\tau)+\mathbf{K}_{\text{f}}\hat{\mathbf{X}}^{\prime}{}^{\text{out}},\;\mathbf{K}_{\text{f}}=\frac{\sqrt{2\eta\tau}g_{\text{\tiny{M}}}\left(e^{-\kappa\tau}(\kappa\tau+2)+\kappa\tau-2\right)}{\sqrt{\kappa}(\kappa\tau-1+e^{-\kappa\tau})}.$
(99)

The symmetric QND Gate is characterized by the equations that read:

$\displaystyle\begin{cases}&\hat{X}_{\text{\tiny{A}}}(\tau)=\hat{X}_{\text{\tiny{A}}}(0)+\mathfrak{G}\hat{X}_{\text{\tiny{M}}}(0)+\mathfrak{N}_{\text{X}_{\text{\tiny{A}}}},\\\
&\hat{P}_{\text{\tiny{A}}}(\tau)=\hat{P}_{\text{\tiny{A}}}(0),\\\
&\hat{X}_{\text{\tiny{M}}}(\tau)=\hat{X}_{\text{\tiny{M}}}(0)+\mathfrak{N}_{\text{X}_{\text{\tiny{M}}}},\\\
&\hat{P}_{\text{\tiny{M}}}(\tau)=\hat{P}_{\text{\tiny{M}}}(0)-\mathfrak{G}\hat{P}_{\text{\tiny{A}}}(0)+\mathfrak{N}_{\text{P}_{\text{\tiny{M}}}},\end{cases},\text{where}\qquad\mathfrak{G}=g_{\text{\tiny{A}}}g_{\text{\tiny{M}}}\sqrt{\eta}\frac{2}{\kappa^{2}}\left(e^{-\kappa\tau}(\kappa\tau+2)+\kappa\tau-2\right)$
(100)

Here, the noises are defined as:

$\displaystyle\mathfrak{N}_{\text{P}_{\text{\tiny{A}}}}=0,$ (101)
$\displaystyle\mathfrak{N}_{\text{X}_{\text{\tiny{M}}}}=\sqrt{\tau
2\Gamma_{\text{\tiny{M}}}}\;\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}}},$
$\displaystyle\mathfrak{N}_{\text{X}_{\text{\tiny{A}}}}=-\frac{g_{\text{\tiny{A}}}\sqrt{2}}{\mathbf{K}_{2}\sqrt{\kappa}}\hat{\mathbf{X}}^{\text{in}}+\frac{\mathbf{K}_{\text{f}}}{\mathbf{K}_{5}}\sqrt{\frac{\eta}{\tau}}\hat{\mathbf{X}}^{\text{in}}_{\text{f}}+\frac{\mathbf{K}_{\text{f}}}{\mathbf{K}_{1}}\sqrt{\frac{1-\eta}{\tau}}\hat{\mathbf{X}}^{\text{vac}}+\frac{g_{\text{\tiny{M}}}\mathbf{K}_{\text{f}}}{\kappa\mathbf{K}_{3}}\sqrt{\frac{4\Gamma_{\text{\tiny{M}}}}{\kappa\tau}}\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}}}_{\text{f}}+$
$\displaystyle+\left(\mathbf{K}_{\text{f}}\sqrt{\frac{2\eta}{\kappa\tau}}(1-e^{-\kappa\tau}(2\kappa\tau+1))-g_{\text{\tiny{A}}}(1-e^{-\kappa\tau})\right)\hat{x}_{\text{c}}(0)+\mathbf{K}_{\text{f}}\sqrt{\frac{2}{\kappa\tau}}(1-e^{-\kappa\tau})\hat{x}^{\prime}_{\text{c}}(0)$
$\displaystyle\mathfrak{N}_{\text{P}_{\text{\tiny{M}}}}=-g_{\text{\tiny{M}}}(1-e^{-\kappa\tau})\hat{p}^{\prime}_{\text{c}}(0)+\sqrt{\tau
2\Gamma_{\text{\tiny{M}}}}\;\bm{\hat{\zeta}}^{\text{P}_{\text{\tiny{M}}}}-\sqrt{2\kappa\eta}\;\frac{g_{\text{\tiny{M}}}}{\kappa\mathbf{K}_{6}}\hat{\mathbf{P}}^{\text{in}}-$
$\displaystyle-2\sqrt{\eta}g_{\text{\tiny{M}}}\frac{(1-e^{-\kappa\tau}(1+\kappa\tau))}{\kappa}\hat{p}_{\text{c}}(0)-\frac{g_{\text{\tiny{M}}}}{\mathbf{K}_{2}}\sqrt{2\kappa(1-\eta)}\hat{\mathbf{P}}^{\text{vac}}.$

$\displaystyle\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}},\text{P}_{\text{\tiny{M}}}}=\frac{1}{\sqrt{\tau
2\Gamma_{\text{\tiny{M}}}}}\int_{0}^{\tau}dt\hat{\zeta}_{\text{X}_{\text{\tiny{M}}},\text{P}_{\text{\tiny{M}}}}(t),$
(102)
$\displaystyle\langle\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}}}\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}}}\rangle=\bm{\hat{\zeta}}^{\text{P}_{\text{\tiny{M}}}}\bm{\hat{\zeta}}^{\text{P}_{\text{\tiny{M}}}}\rangle=1,\;\langle\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}}}\bm{\hat{\zeta}}^{\text{P}_{\text{\tiny{M}}}}\rangle=\bm{\hat{\zeta}}^{\text{P}_{\text{\tiny{M}}}}\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}}}\rangle=0$
(103)
$\displaystyle\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}}}_{\text{f}}=\frac{\mathbf{K}_{3}}{\sqrt{2\Gamma_{\text{\tiny{M}}}}}\left((\kappa\tau-1+e^{-\kappa\tau})\ast\hat{\zeta}_{\text{X}_{\text{\tiny{M}}}}(\tau)\right),\qquad\langle\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}}}_{\text{f}}\bm{\hat{\zeta}}^{\text{X}_{\text{\tiny{M}}}}\rangle=\frac{\mathbf{K}_{3}\mathbf{K}_{4}}{\sqrt{\tau}},$
(104)
$\displaystyle\mathbf{K}_{3}=\sqrt{\frac{6\kappa}{3+2\kappa\tau(3+\kappa\tau(\kappa\tau-3)))-3e^{-2\kappa\tau}-12e^{-\kappa\tau}\kappa\tau}}$
(105)
$\displaystyle\mathbf{K}_{4}=\frac{2(1-\kappa\tau)-2e^{-\kappa\tau}+\kappa^{2}\tau^{2}}{2\kappa}$
(106)
$\displaystyle\hat{\mathbf{X}}^{\text{vac}}=\mathbf{K}_{1}\left(\left(1-2e^{-\kappa\tau}\right)\ast\hat{x}_{\text{vac}}(\tau)\right),\qquad\mathbf{K}_{1}=\sqrt{\frac{\kappa}{\kappa\tau-2+4e^{-\kappa\tau}-2e^{-2\kappa\tau}}}$
(107)
$\displaystyle\hat{\mathbf{P}}^{\text{vac}}=\mathbf{K}_{2}\left(\left(1-e^{-\kappa\tau}\right)\ast\hat{p}_{\text{vac}}(\tau)\right),\qquad\mathbf{K}_{2}=\sqrt{\frac{2\kappa}{4e^{-\kappa\tau}+2\kappa\tau-3-e^{-2\kappa\tau}}}$
(108)
$\displaystyle\hat{\mathbf{P}}^{\text{in}}=\mathbf{K}_{6}\left((1-e^{-\kappa\tau}(2\kappa\tau+1))\ast\hat{p}_{\text{in}}(\tau)\right),$
(109)
$\displaystyle\mathbf{K}_{6}=\sqrt{\frac{2\kappa}{(2\kappa\tau-7)+4e^{-\kappa\tau}(2\kappa\tau+3)-e^{-2\kappa\tau}(5+4\kappa\tau(2+\kappa\tau))}}$
(110)
$\displaystyle\hat{\mathbf{X}}^{\text{in}}=\mathbf{K}_{2}\left(\left(1-e^{-\kappa\tau}\right)\ast\hat{x}_{\text{in}}(\tau)\right)$
(111)
$\displaystyle\hat{\mathbf{X}}^{\text{in}}_{\text{f}}=\mathbf{K}_{5}\left(\left(1-4\kappa\tau
e^{-\kappa\tau}\right)\ast\hat{x}_{\text{in}}(\tau)\right),$ (112)
$\displaystyle\mathbf{K}_{5}=\sqrt{\frac{\kappa}{(\kappa\tau-4)+8e^{-\kappa\tau}(1+\kappa\tau)-4e^{-2\kappa\tau}(1+2\kappa\tau(1+\kappa\tau))}}$
(113)
$\displaystyle\mathbf{K}_{7}=\frac{2e^{-\kappa\tau}(3+2\kappa\tau)+(\kappa\tau-4)-2e^{-2\kappa\tau}(1+2\kappa\tau)}{\kappa},\qquad\langle\hat{\mathbf{X}}^{\text{in}}\hat{\mathbf{X}}^{\text{in}}_{\text{f}}\rangle=\mathbf{K}_{2}\mathbf{K}_{5}\mathbf{K}_{7}$
(114)

The physical parameters of the systems are the pulse duration $\tau$, the
optical decay rates of the cavity $\kappa_{\text{\tiny{A,M}}}$, the coupling
strengths $g_{\text{\tiny{A,M}}}$, the gate efficiency $\eta$, the mechanical
damping coefficient $\gamma_{\text{\tiny{M}}}$ that shows how good the
mechanics is isolated from the thermal bath with average phonon number
$n_{\text{th}}$ (the two latter parameters are combined in the reheating rate
$\Gamma_{\text{\tiny{M}}}=\gamma_{\text{\tiny{M}}}(2n_{\text{th}}+1)$).

Figures show the dependencies of the HOM element on different parameters. We
assume that Hong-Ou-Mandel interference takes place if the HOM-element of the
output state of the gate lies above its corresponding input threshold.

a

b

Figure 6:  Matrix element
$\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$ of the output state of
the gate as a function of the coupling strength for the pulse length
$\tau\kappa_{\text{\tiny{A}}}=100$: a) Light-atom QND gate with the
independent single-photon (light) and single-polariton (atom) states at the
input. Dependence on coupling strength $g_{\text{\tiny{A}}}$ with efficiencies
$\eta=1,\;0.9,\;0.1$. Well pronounced maximum shifts to the left with
decreasing efficiency. The inset demonstrates HOM-element as the function of
the pulse duration $\tau\kappa_{\text{\tiny{A}}}$. b) Light-mechanical QND
gate with the independent single-photon (light) and single-phonon (mechanics)
states at the input. Dependence on coupling strength $g_{\text{\tiny{M}}}$ for
the rethermalization rates
$\Gamma_{\text{\tiny{M}}}=0.01\kappa_{\text{\tiny{M}}},\;0.001\kappa_{\text{\tiny{M}}}$
for the different efficiencies ($\eta=1,\;0.1$). Note,
$\Gamma_{\text{\tiny{M}}}=0.01\kappa_{\text{\tiny{M}}}$ is already too high
for $\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$ to surpass the
output classical threshold, but even for $\eta=0.1$ the HOM-element lies quite
close to the input classical threshold. For both (a) and (b), the dashed gray
line is the output threshold and the blue curves of the corresponding
thickness and dashing are the input thresholds (phase randomized).

a

b

Figure 7:  Matrix element
$\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$ of the output state
for the non-adiabatic atom-mechanical QND gate with the independent single-
boson states at the input: a) dependence on coupling strength for the
different efficiencies $\eta=1,\;0.8,\;0.1$. Blue curves of the same dashing
indicate the corresponding input thresholds. Efficiency $\eta=0.1$ is too
small to allow HOM-element to surpass the input thresholds. b) dependence on
the squeezing $S$ of the light pulse for the different efficiencies
$\eta=1,\;0.8,\;0.1$. For both (a) and (b),
$\kappa_{\text{\tiny{A}}}=\kappa_{\text{\tiny{M}}}=\kappa,\;g_{\text{\tiny{A}}}=g_{\text{\tiny{M}}}=g$.

##  Appendix B Ideal QND gate. Nonclassicality thresholds

Here, we first describe the ideal case of the Quantum NonDemolition (QND)
transformation and compare it with the beam splitter (BS) transformation. We
then proceed to introduce the matrix elements of the output quantum state,
that correspond to the bunching of excitations, and describe how to calculate
them in different cases.

###  B.1 Comparison of a beam splitter with a QND transformation

A beam splitter (BS) transformation, characterized by the Hamiltonian
$H_{\text{BS}}=\Theta\hbar i\left(a^{\dagger}b-b^{\dagger}a\right)$, describes
an evolution of two quantum oscillators a and b. The only parameter of this
transformation is the transmittance coefficient $\mathsf{T}=\cos^{2}\Theta$.

A quantum Non-Demolition (QND) gate, characterized by the Hamiltonian
$H_{\mathsf{G}}=\mathsf{G}\hbar i(a+a^{\dagger})(b^{\dagger}-b)/2$, describes
another type of evolution of the two oscillators. Gain $\mathsf{G}$ is the
only parameter characterizing the ideal QND gate transformation.

Hamiltonians and unitary transformations corresponding to the BS and QND
transformations are the following:

$\displaystyle H_{\text{BS}}=\Theta\hbar
i\left(a^{\dagger}b-b^{\dagger}a\right),$ $\displaystyle
U_{\text{BS}}=\sum\limits_{n=0}^{\infty}\frac{\Theta^{n}(a^{\dagger}b-b^{\dagger}a)^{n}}{n!}),$
(115) $\displaystyle H_{\mathsf{G}}=\frac{\mathsf{G}}{2}\hbar
i(a+a^{\dagger})(b^{\dagger}-b),$ $\displaystyle
U_{\mathsf{G}}=\sum\limits_{n=0}^{\infty}\frac{g^{n}(a^{\dagger}+a)^{n}(b^{\dagger}-b)^{n}}{n!}.$
(116)

For $U_{\mathsf{G}}$ we can use the Zassenhaus formula and derive it in the
normal order:

$\displaystyle
U_{\mathsf{G}}=\sum\limits_{n=0}^{\infty}g^{n}n!\sum\limits_{i,j,l,m=0}^{\infty}\frac{(a^{\dagger})^{i}a^{j}}{i!j!\left(\frac{n-i-j}{2}\right)!2^{\frac{n-i-j}{2}}}\frac{(-1)^{m}(b^{\dagger})^{l}b^{m}}{l!m!\left(\frac{n-l-m}{2}\right)!(-2)^{\frac{n-l-m}{2}}},\;\frac{n-l-m}{2},\;\frac{n-i-j}{2}\in\mathbb{N}_{0},$

where $\mathbb{N}_{0}$ is the set of natural numbers including $0$.

$\displaystyle
e^{t(a^{\dagger}+a)}=e^{ta^{\dagger}}e^{ta}e^{-\frac{t^{2}}{2}[a^{\dagger},a]}=e^{ta^{\dagger}}e^{ta}e^{\frac{t^{2}}{2}},\qquad
e^{t(b^{\dagger}-b)}=e^{tb^{\dagger}}e^{-tb}e^{-\frac{t^{2}}{2}[b^{\dagger},-b]}=e^{tb^{\dagger}}e^{-tb}e^{-\frac{t^{2}}{2}}$
$\displaystyle\sum\limits_{n=0}^{\infty}\frac{(a^{\dagger}+a)^{n}}{n!}t^{n}=\sum\limits_{i,j,k=0}^{\infty}\frac{(a^{\dagger})^{i}a^{j}}{i!j!k!2^{k}}t^{i+j+2k},\qquad\sum\limits_{n=0}^{\infty}\frac{(b^{\dagger}-b)^{n}}{n!}t^{n}=\sum\limits_{l,m,k=0}^{\infty}\frac{(-1)^{m}(b^{\dagger})^{l}b^{m}}{l!m!k!(-2)^{k}}t^{l+m+2k},$
(117)
$\displaystyle\frac{(a^{\dagger}+a)^{n}}{n!}=\sum\limits_{i,j,k=0}^{\infty}\frac{\delta_{n,i+j+2k}}{i!j!k!2^{k}}(a^{\dagger})^{i}a^{j}=\sum\limits_{i,j=0}^{\infty}\frac{(a^{\dagger})^{i}a^{j}}{i!j!\left(\frac{n-i-j}{2}\right)!2^{\frac{n-i-j}{2}}}$
(118)
$\displaystyle\frac{(b^{\dagger}-b)^{n}}{n!}=\sum\limits_{l,m,k=0}^{\infty}\frac{(-1)^{m}\cdot\delta_{n,l+m+2k}}{l!m!k!(-2)^{k}}(b^{\dagger})^{l}b^{m}=\sum\limits_{l,m=0}^{\infty}\frac{(-1)^{m}(b^{\dagger})^{l}b^{m}}{l!m!\left(\frac{n-l-m}{2}\right)!(-2)^{\frac{n-l-m}{2}}}$
(119)

In the article we are interested in
$\langle\text{HOM}|U_{\text{G,BS}}|\varphi\rangle_{\text{in}}$, where
$\langle\text{HOM}|=(1/\sqrt{2})\left(\langle 0|_{\text{a}}\langle
2|_{\text{b}}-\langle 2|_{\text{a}}\langle 0|_{\text{b}}\right)$. Then we need
$\langle\text{HOM}|U_{\text{G,BS}}$. Note that even $\langle
1|_{\text{a}}(a^{\dagger})^{2}=0$ (and $\langle
2|_{\text{a}}(a^{\dagger})^{3}=0$). Then, using $\langle
0|_{\text{a}}(a)^{k}=\sqrt{k!}\langle k|_{\text{a}},\;\langle
1|_{\text{a}}(a)^{k}=\sqrt{(k+1)!}\langle k+1|_{\text{a}},\;\langle
2|_{\text{a}}(a)^{k}=\sqrt{(k+2)!/2}\langle k+2|_{\text{a}}$, we can easily
obtain:

I. For the QND transformation:

$\displaystyle\langle\text{HOM}|\;U_{\mathsf{G}}=\frac{1}{\sqrt{2}}\left(\langle
0|_{\text{a}}\langle 2|_{\text{b}}-\langle 2|_{\text{a}}\langle
0|_{\text{b}}\right)\;U_{\mathsf{G}}=$ (120)
$\displaystyle=\sum\limits_{n,m,j=0}^{\infty}\frac{(-1)^{m}g^{n}n!}{\sqrt{j!m!}}(\frac{\sqrt{(m+2)(m+1)}\langle
m+2|_{\text{b}}\langle j|_{\text{a}}+\sqrt{(j+2)(j+1)}\langle
m|_{\text{b}}\langle j+2|_{\text{a}}}{2\cdot\left(\frac{n-j}{2}\right)!\cdot
2^{\frac{n-j}{2}}\cdot\left(\frac{n-m}{2}\right)!\cdot(-2)^{\frac{n-m}{2}}}-$
$\displaystyle-\frac{\sqrt{j+1}\langle m|_{\text{b}}\langle
j+1|_{\text{a}}}{\left(\frac{n-1-j}{2}\right)!\cdot
2^{\frac{n-1-j}{2}}\cdot\left(\frac{n-m}{2}\right)!\cdot(-2)^{\frac{n-m}{2}}}-\frac{\langle
m|_{\text{b}}\langle j|_{\text{a}}}{2\cdot\left(\frac{n-2-j}{2}\right)!\cdot
2^{\frac{n-2-j}{2}}\cdot\left(\frac{n-m}{2}\right)!\cdot(-2)^{\frac{n-m}{2}}}+$
$\displaystyle+\frac{\sqrt{m+1}\langle m+1|_{\text{b}}\langle
j|_{\text{a}}}{\left(\frac{n-j}{2}\right)!\cdot
2^{\frac{n-j}{2}}\cdot\left(\frac{n-1-m}{2}\right)!\cdot(-2)^{\frac{n-1-m}{2}}}+\frac{\langle
m|_{\text{b}}\langle j|_{\text{a}}}{2\cdot\left(\frac{n-j}{2}\right)!\cdot
2^{\frac{n-j}{2}}\cdot\left(\frac{n-2-m}{2}\right)!\cdot(-2)^{\frac{n-2-m}{2}}}).$

Let us use
$|\varphi\rangle_{\text{in}}=|1\rangle_{\text{b}}|1\rangle_{\text{a}}$ and
obtain matrix element
$|\langle\text{HOM}|\;U_{\mathsf{G}}\;|1\rangle_{\text{b}}|1\rangle_{\text{a}}|^{2}$:

$\displaystyle|\langle\text{HOM}|\;U_{\mathsf{G}}\;|1\rangle_{\text{b}}|1\rangle_{\text{a}}|^{2}=|\sum\limits_{k=0}^{\infty}\frac{(-1)^{\text{k}}g^{2k+1}(2k+1)!}{(k!)^{2}2^{2k-1}}+\sum\limits_{k=1}^{\infty}\frac{(-1)^{\text{k}}g^{2k+1}(2k+1)!}{k!(k-1)!2^{2k-1}}|^{2}=\frac{16\mathsf{G}^{2}(-8+\mathsf{G}^{2})^{2}}{(4+\mathsf{G}^{2})^{5}},$
$\displaystyle(g=\mathsf{G}/2).$ (121)

II. For the BS transformation:

$\displaystyle\langle\text{HOM}|\;U_{\text{BS}}=\frac{1}{\sqrt{2}}\left(\langle
2|_{\text{b}}\langle 0|_{\text{a}}-\langle 0|_{\text{b}}\langle
2|_{\text{a}}\right)U=$
$\displaystyle=\frac{1}{\sqrt{2}}\sum\limits_{n=0}^{\infty}\left(\frac{(-1)^{n+1}\Theta^{2n+1}(2\sqrt{2})^{n+1}(\sqrt{2})^{\text{n}}\langle
1|_{\text{a}}\langle
1|_{\text{b}}}{(2n+1)!}+\frac{(-1)^{n}\Theta^{2n}(2\sqrt{2})^{\text{n}}(\sqrt{2})^{\text{n}}\left(\langle
2|_{\text{b}}\langle 0|_{\text{a}}-\langle 0|_{\text{b}}\langle
2|_{\text{a}}\right)}{(2n)!}\right).$ (122)

Let us check the result by obtaining
$|\langle\text{HOM}|\;U_{\text{BS}}\;|1\rangle_{\text{b}}|1\rangle_{\text{a}}|^{2}$:

$\displaystyle|\langle\text{HOM}|\;U_{\text{BS}}\;|1\rangle_{\text{b}}|1\rangle_{\text{a}}|^{2}=|\frac{1}{\sqrt{2}}\sum\limits_{n=0}^{\infty}\frac{(-1)^{n+1}\Theta^{2n+1}(2\sqrt{2})^{n+1}(\sqrt{2})^{\text{n}}}{(2n+1)!}|^{2}=4\mathsf{T}(1-\mathsf{T}),\qquad(\cos^{2}\Theta=\mathsf{T}).$
(123)

If we restrict ourselves by the limited input subspace, including only
$|1\rangle_{\text{a}}|1\rangle_{\text{b}},|0\rangle_{\text{a}}|1\rangle_{\text{b}},|1\rangle_{\text{a}}|0\rangle_{\text{b}}$
and $|0\rangle_{\text{a}}|0\rangle_{\text{b}}$, the formulas can be
simplified:

$\displaystyle\langle\text{HOM}|\;U_{\text{BS}}=\sum\limits_{n=0}^{\infty}\frac{(-1)^{n+1}\left(2\Theta\right)^{2n+1}}{(2n+1)!}\langle
1|_{a}\langle 1|_{b}=-\sin\left(2\Theta\right)\;\langle 1|_{\text{a}}\langle
1|_{\text{b}},$ (124)
$\displaystyle\langle\text{HOM}|\;U_{\mathsf{G}}=\sum\limits_{k=0}^{\infty}\frac{(-1)^{k}(\mathsf{G}/2)^{2k+1}(2k+1)!(1+k)}{(k!)^{2}2^{2k-1}}\;\langle
1|_{a}\langle
1|_{b}+\frac{(-1)^{k}(\mathsf{G}/2)^{2k+2}(2k+2)!}{(k!)(k+1)!2^{2k+1}}\;\langle
0|_{\text{a}}\langle 0|_{\text{b}}=$
$\displaystyle=\frac{4\mathsf{G}(8-\mathsf{G}^{2})}{(4+\mathsf{G}^{2})^{5/2}}\;\langle
1|_{\text{a}}\langle
1|_{\text{b}}+\frac{2\mathsf{G}^{2}}{(4+\mathsf{G}^{2})^{3/2}}\;\langle
0|_{\text{a}}\langle 0|_{\text{b}}.$ (125)

Table 1: Matrix elements of the output state corresponding the different cases
of the input for the BS and QND gate transformations

input | $\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$
---|---
|  BS |  QND gate
$|1\rangle_{\text{a}}|1\rangle_{\text{b}}$ |  $4\mathsf{T}(1-\mathsf{T})$ |  $\frac{16\mathsf{G}^{2}(\mathsf{G}^{2}-8)^{2}}{(4+\mathsf{G}^{2})^{5}}$
$|1\rangle_{\text{a}}|0\rangle_{\text{b}}$ |  $0$ |  $0$
$|0\rangle_{\text{a}}|1\rangle_{\text{b}}$ |  $0$ |  $0$
$|0\rangle_{\text{a}}|0\rangle_{\text{b}}$ |  $0$ |  $\frac{4\mathsf{G}^{4}}{(4+\mathsf{G}^{2})^{3}}$

In the case of the general input BS transformation provides non-zero matrix
elements exclusively for
$|\varphi\rangle_{\text{in}}=|1\rangle_{\text{a}}|1\rangle_{\text{b}},|0\rangle_{\text{a}}|2\rangle_{\text{b}},|2\rangle_{\text{a}}|0\rangle_{\text{b}}$.
However, the gate transformation generates and annihilates excitations in
pairs. That means that, for instance, if at the gate input there are states
with an even number of bosons (as
$|0\rangle_{\text{a}}|0\rangle_{\text{b}},|1\rangle_{\text{a}}|1\rangle_{\text{b}},|1\rangle_{\text{a}}|3\rangle_{\text{b}}$
etc.), then $\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$ will be
nonzero (some functions of the gain). If at the input of the gate there is a
state with an odd number of excitations (as
$|0\rangle_{\text{a}}|1\rangle_{\text{b}},|1\rangle_{\text{a}}|0\rangle_{\text{b}},|1\rangle_{\text{a}}|2\rangle_{\text{b}},|2\rangle_{\text{a}}|1\rangle_{\text{b}}$
etc.), then the matrix elements of the output state will be zero.

There is a significant difference between these two transformations. The BS
transformation is passive, it neither creates nor annihilates excitations in a
system of two harmonic oscillators. If initially there are exactly one
excitation in each of the oscillators, at the output of a BS they can appear
bunched in a single mode via the Hong-Ou-Mandel (HOM) effect. It is not
sufficient to have just one excitation in a single mode to observe bunching.
Unlike BS, the QND transformation is active, which means it can possibly
change the total number of excitations in the system (the energy of the
system). The creation of quanta by QND interaction can produce effects that
resemble the HOM interference. Although such case can be confused with the HOM
interference effect, it is still possible to analyze whether the QND
interaction is capable of generating the non-classical two-quanta
superpositions going beyond any classical states serving as input to the QND
interaction. However, such analysis requires a general approach to the HOM
interference beyond the simple case with the passive BS interaction.

In this generalized description, the matrix element of the output state
$\rho_{\text{out}}$, describing the success probability of detection of two-
photon HOM-entangled states (the HOM element), can be introduced as
$|\langle\text{HOM}|U|\varphi\rangle_{\text{in}}|^{2}$, where
$|\varphi\rangle_{\text{in}}$ is an initial state, and $U$ is a unitary
transformation. Here, the HOM-state is determined as
$\langle\text{HOM}|=\left(\langle 2|_{\text{b}}\langle 0|_{\text{a}}-\langle
0|_{\text{b}}\langle 2|_{\text{a}}\right)/\sqrt{2}$. It is well known that a
BS provides an ideal photon bunching (the HOM effect). This means that the
success probability
$|\langle\text{HOM}|U_{\text{BS}}|\varphi\rangle_{\text{in}}|^{2}$ at the
output of the BS equals one. This effect occurs when two identical quanta
enter a balanced beam splitter ($\mathsf{T}=0.5$), one in each input port (the
input state
$|\varphi\rangle_{\text{in}}=|1\rangle_{\text{a}}|1\rangle_{\text{b}}$).

In order to compare the two transformations in the context of the HOM effect,
let us look at the matrix elements of the output state of each of the
transformations. For simplicity, first let us restrict the subspace of the
input and assume that $|\varphi\rangle_{\text{in}}$ belongs to the space of
coherent superpositions of vacuum and one excitation of each mode, that is an
arbitrary pure superposition of
$|1\rangle_{\text{a}}|1\rangle_{\text{b}},|0\rangle_{\text{a}}|1\rangle_{\text{b}},|1\rangle_{\text{a}}|0\rangle_{\text{b}}$
and $|0\rangle_{\text{a}}|0\rangle_{\text{b}}$. Then, to obtain the desired
matrix element we need:

$\displaystyle\langle\text{HOM}|\;U_{\text{BS}}=-\sin\left(2\Theta\right)\;\langle
1|_{\text{a}}\langle
1|_{\text{b}},\;\langle\text{HOM}|\;U_{\mathsf{G}}=\frac{4\mathsf{G}(8-\mathsf{G}^{2})}{(4+\mathsf{G}^{2})^{5/2}}\;\langle
1|_{\text{a}}\langle
1|_{\text{b}}+\frac{2\mathsf{G}^{2}}{(4+\mathsf{G}^{2})^{3/2}}\;\langle
0|_{\text{a}}\langle 0|_{\text{b}}.$ (126)

It clearly shows that the HOM matrix element provided by the inputs
$|0\rangle_{\text{a}}|1\rangle_{\text{b}}$ and
$|1\rangle_{\text{a}}|0\rangle_{\text{b}}$ is equal to zero for both BS and
QND transformations. Matrix element provided by
$|0\rangle_{\text{a}}|0\rangle_{\text{b}}$ input is equal to zero in the case
of a beamsplitter. However, for the QND gate this element is a function of the
gain $\mathsf{G}$ and equals zero only in the trivial case with
$\mathsf{G}=0$. That is, by varying the gain of the QND gate, it is impossible
to make the contribution of the input vacua
$|0\rangle_{\text{a}}|0\rangle_{\text{b}}$ vanish, in order to render these
two transformations fully analogous.

For the case of a QND gate with $|1\rangle_{\text{a}}|1\rangle_{\text{b}}$ at
the input, one can observe that for a certain region of the parameter
$\mathsf{G}$, the probability of bunching of both excitations in one subsystem
is higher than the probability of equal redistribution of the excitations
between the subsystems. Visually, it is characterized by the presence of the
maximum of $\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$ (see Fig.
4(b) )approximately equal to $0.26$ for
$\mathsf{G}=\sqrt{11-\sqrt{105}}\approx 0.87$ (as compared to $1$ for the BS
with $\Theta=\pi/4$). However, we should keep in mind that this correspondence
to the case of a BS is not complete due to the non-zero contribution from the
vacuum input for the gate case that does not exist in the case of a BS.

###  B.2 BS and QND transformations using quadratures.

Both Quantum Non-Demolition (QND) gate and beam splitter BS, transform
quadratures of two quantum oscillators a and b, correlating the quadratures of
the two oscillators. The initial quadratures
$\mathbf{r}^{\text{in}}=(\hat{X}_{\text{a}}(0),\hat{P}_{\text{a}}(0),\hat{X}_{\text{b}}(0),\hat{P}_{\text{b}}(0))^{T}$
after transformation should relate to the final quadratures
$\mathbf{r}^{\text{out}}\equiv\mathbf{r}=(x_{\text{a}},p_{\text{a}},x_{\text{b}},p_{\text{b}})^{T}$
as

$\displaystyle\mathbf{r}^{\text{out}}=T_{\text{G,BS}}\mathbf{r}^{\text{in}}+N,\;\text{where}$
(127) $\displaystyle
T_{\text{G}}=\left(\begin{array}[]{cccc}\mathfrak{T}_{\text{a}}&0&\mathfrak{G}_{\text{a}}&0\\\
0&\mathfrak{T}_{\text{a}}&0&0\\\ 0&0&\mathfrak{T}_{\text{b}}&0\\\
0&\mathfrak{G}_{\text{b}}&0&\mathfrak{T}_{\text{b}}\end{array}\right),\;T_{\text{BS}}=\left(\begin{array}[]{cccc}\sqrt{\mathsf{T}}&0&\sqrt{1-\mathsf{T}}&0\\\
0&\sqrt{\mathsf{T}}&0&\sqrt{1-\mathsf{T}}\\\
-\sqrt{1-\mathsf{T}}&0&\sqrt{\mathsf{T}}&0\\\
0&-\sqrt{1-\mathsf{T}}&0&\sqrt{\mathsf{T}}\\\ \end{array}\right).$ (136)

where $\mathsf{T}$ is the transmittance coefficient (BS transformation);
$\mathfrak{G}_{\text{a,b}}$ are the controllable gains of the built QND gate,
$\mathfrak{T}_{\text{a,b}}$ are the transfer factors (QND transformation).
$N\equiv(\mathfrak{N}_{\text{X}_{\text{a}}},\mathfrak{N}_{\text{P}_{\text{a}}},\mathfrak{N}_{\text{X}_{\text{b}}},\mathfrak{N}_{\text{P}_{\text{b}}})^{T}$
describe the excess noises, thus for the ideal transformation they should be
negligible
($\mathfrak{N}_{\text{X}_{\text{a}},\text{P}_{\text{a}},\text{X}_{\text{b}},\text{P}_{\text{b}}}\rightarrow
0$) (and for QND transfer factors equal one ($\mathfrak{T}_{\text{a,b}}=1$),
while the gains are of the same magnitude but opposite sign
($\mathsf{G}=\mathfrak{G}_{\text{a}}=-\mathfrak{G}_{\text{b}}$). Thus, gain
$\mathsf{G}$ is the only parameter characterizing the ideal QND gate
transformation, while the BS transformation also is characterized by a single
parameter $\mathsf{T}$.

###  B.3 Wigner function of the state

The Wigner functions (WF) can be used to calculate the HOM element. Let us
demonstrate how to obtain the WF corresponding to an arbitrary operator.
First, let us remind that the wave-function of the $n$-th excited level (Fock
state $\left|n\right\rangle$) can be derived using Hermite polynomials as:

$\displaystyle\psi_{n}(x)=\langle
x|n\rangle=\frac{1}{\sqrt{2^{n+1/2}\pi^{1/2}n!}}H_{n}\left(\frac{x}{\sqrt{2}}\right)\exp\left(-\frac{x^{2}}{4}\right),$
(137) $\displaystyle H_{0}(z)=1,$ (138) $\displaystyle H_{1}(z)=2z,$ (139)
$\displaystyle H_{2}(z)=4z^{2}-2=2(2z^{2}-1),...$ (140)

The WF of the operator $|\varphi\rangle\langle\psi|$ in the $s$-dimensional
space is as follows:

$\displaystyle
W_{|\varphi\rangle\langle\psi|}(\mathbf{x},\mathbf{p})=\frac{1}{(4\pi)^{s}}\int\limits_{-\infty}^{\infty}d\mathbf{y}\;e^{-\frac{i\mathbf{p}\mathbf{y}}{2}}\langle\mathbf{x}+\frac{\mathbf{y}}{2}|\varphi\rangle\langle\psi|\mathbf{x}-\frac{\mathbf{y}}{2}\rangle,\;\;\text{where}\;\mathbf{q}=(q_{1},q_{2},...q_{s})\;\;\text{for}\;\;\mathbf{q=x,y,p}.$
(141)

Here,

$\displaystyle\langle
x+\frac{y}{2}|2\rangle=\psi_{|2\rangle}(x+\frac{y}{2})=\frac{1}{\sqrt{2}\sqrt[4]{2\pi}}\left(\left(x+\frac{y}{2}\right)^{2}-1\right)\exp\left(-\frac{(x+\frac{y}{2})^{2}}{4}\right)$
(142) $\displaystyle\langle
x+\frac{y}{2}|1\rangle=\psi_{|1\rangle}(x+\frac{y}{2})=\frac{1}{\sqrt[4]{2\pi}}\left(x+\frac{y}{2}\right)\exp\left(-\frac{(x+\frac{y}{2})^{2}}{4}\right)$
(143) $\displaystyle\langle
x+\frac{y}{2}|0\rangle=\psi_{|0\rangle}(x+\frac{y}{2})=\frac{1}{\sqrt[4]{2\pi}}\exp\left(-\frac{(x+\frac{y}{2})^{2}}{4}\right)$
(144) $\displaystyle\langle
0|x-\frac{y}{2}\rangle=\psi^{\ast}_{|0\rangle}(x-\frac{y}{2})=\frac{1}{\sqrt[4]{2\pi}}\exp\left(-\frac{(x-\frac{y}{2})^{2}}{4}\right)\qquad\text{etc}.$
(145)

In the article, we work in the $4$-dimensional space. Thus, for the operator
$A=|\text{HOM}\rangle\langle\text{HOM}|$, associated with the state
$|\text{HOM}\rangle=\left(|0_{2}\;2_{1}\rangle-|2_{2}\;0_{1}\rangle\right)/\sqrt{2}$,
the WF can be calculated as:

$\displaystyle
W_{\text{HOM}}(x_{1},p_{1},x_{2},p_{2})=\frac{1}{(4\pi)^{2}}\iint
dy_{1}dy_{2}\;e^{-\frac{i(p_{1}y_{1}+p_{2}y_{2})}{2}}f_{\text{HOM}}(x_{1},y_{1},x_{2},y_{2})$
$\displaystyle=\frac{1}{16\pi^{2}}\exp\left(-\frac{p_{1}^{2}+p_{2}^{2}+x_{1}^{2}+x_{2}^{2}}{2}\right)((p_{1}-p_{2})^{2}+(x_{1}-x_{2})^{2}-2)((p_{1}+p_{2})^{2}+(x_{1}+x_{2})^{2}-2).$
(146)

where

$\displaystyle f_{\text{HOM}}(x_{1},y_{1},x_{2},y_{2})=\langle
x_{1}+\frac{y_{1}}{2}|\langle
x_{2}+\frac{y_{2}}{2}|A\;|x_{2}-\frac{y_{2}}{2}\rangle|x_{1}-\frac{y_{1}}{2}\rangle=$
$\displaystyle=\frac{1}{2}\left(\langle
x_{1}+\frac{y_{1}}{2}||2_{1}\rangle\langle
x_{2}+\frac{y_{2}}{2}||0_{2}\rangle-\langle
x_{2}+\frac{y_{2}}{2}||2_{2}\rangle\;\langle
x_{1}+\frac{y_{1}}{2}||0_{1}\rangle\right)\otimes$
$\displaystyle\otimes\left(\langle 0_{1}||x_{1}-\frac{y_{1}}{2}\rangle\langle
2_{2}||x_{2}-\frac{y_{2}}{2}\rangle-\langle
2_{1}||x_{1}-\frac{y_{1}}{2}\rangle\langle
0_{2}||x_{2}-\frac{y_{2}}{2}\rangle\right)=$
$\displaystyle=\frac{1}{8\pi}\exp\left(-\frac{(x^{2}_{2}+\frac{y^{2}_{2}}{4})+(x^{2}_{1}+\frac{y^{2}_{1}}{4})}{2}\right)\left(\left(x_{1}+\frac{y_{1}}{2}\right)^{2}-\left(x_{2}+\frac{y_{2}}{2}\right)^{2}\right)\left(\left(x_{1}-\frac{y_{1}}{2}\right)^{2}-\left(x_{2}-\frac{y_{2}}{2}\right)^{2}\right).$
(147)

Analogically, for the operator $A^{\prime}=|1_{2}1_{1}\rangle\langle
1_{1}1_{2}|$:

$\displaystyle f_{|11\rangle\langle
11|}(x_{1},y_{1},x_{2},y_{2})=\frac{1}{2\pi}\left(x_{1}^{2}-\left(\frac{y_{1}}{2}\right)^{2}\right)\left(x_{2}^{2}-\left(\frac{y_{2}}{2}\right)^{2}\right)\exp\left(-\frac{x_{1}^{2}+(\frac{y_{1}}{2})^{2}+x_{2}^{2}+(\frac{y_{2}}{2})^{2}}{2}\right),$
(148) $\displaystyle W_{|11\rangle\langle
11|}(x_{1},p_{1},x_{2},p_{2})=\frac{1}{4\pi^{2}}\exp\left(-\frac{p_{1}^{2}+p_{2}^{2}+x_{1}^{2}+x_{2}^{2}}{2}\right)(-1+p_{1}^{2}+x_{1}^{2})(-1+p_{2}^{2}+x_{2}^{2})).$
(149)

###  B.4 Calculation of the HOM matrix element

In order to calculate the matrix elements for an ideal (without additional
noise) transformation, it is enough to know the form of the unitary
transformation $U$. Thus the matrix element of the output can be calculated
as:

$\mathsf{M}_{\text{HOM}}=\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle=|\langle\text{HOM}|U|\varphi\rangle_{\text{in}}|^{2}.$
(150)

To evaluate the robustness of the QND gate against photon loss, we examine an
incoherent mixture of vacuum and single-photon states at each input port of
the gate

$\rho_{\text{in}}^{(p)}=(p_{\text{a}}\;|1\rangle\langle
1|+(1-p_{\text{a}})\;|0\rangle\langle
0|)_{\text{a}}\cdot(p_{\text{b}}\;|1\rangle\langle
1|+(1-p_{\text{b}})\;|0\rangle\langle 0|)_{\text{b}},$ (151)

where the parameter $p_{\text{a,b}}$ characterizes how much vacuum has been
admixed to the single-photon state at the input ports, and calculate matrix
elements for the output state of the gate. Using Eq. (17), we can obtain the
HOM matrix element of the output state of the gate:

$\langle\text{HOM}|\rho_{\text{out}}^{(p)}|\text{HOM}\rangle=p_{\text{a}}p_{\text{b}}\frac{16\mathsf{G}^{2}(\mathsf{G}^{2}-8)^{2}}{(4+\mathsf{G}^{2})^{5}}+(1-p_{\text{a}})(1-p_{\text{b}})\frac{4\mathsf{G}^{4}}{(4+\mathsf{G}^{2})^{3}}.$
(152)

This matrix element is symmetrical with respect to $p_{\text{a}}$ and
$p_{\text{b}}$. Surprisingly, the independent coherent superpositions
$(\sqrt{p_{\text{a}}}\;|1\rangle+\sqrt{1-p_{\text{a}}}\;|0\rangle)_{\text{a}}\cdot(\sqrt{p_{\text{b}}}\;|1\rangle+\sqrt{1-p_{\text{b}}}\;|0\rangle)_{\text{b}}$
at the input give rise to the same matrix element Eq. (23) as the mixture Eq.
(22).

This approach is good to calculate the ME when the input state is a
$|n\rangle_{\text{a}}|m\rangle_{\text{b}}$-boson state or some combination of
them. Sometimes it is more convenient to take a different approach described
below. We can use the Wigner function (WF) and the matrix elements
$\langle\psi|\rho_{\text{out}}|\varphi\rangle$ of the output state can be
defined as:

$\mathsf{M}_{|\varphi\rangle\langle\psi|}=\langle\psi|\rho_{\text{out}}|\varphi\rangle=(4\pi)^{2}\iiiint
d\mathbf{r}\;W_{|\varphi\rangle\langle\psi|}(\mathbf{r})\cdot
W_{\text{out}}(\mathbf{r}),$ (153)

where $W_{\text{out}}(\mathbf{r})$ is the WF of the output state,
$\mathbf{r}=\mathbf{r}^{\text{out}}=(x_{\text{a}},p_{\text{a}},x_{\text{b}},p_{\text{b}})^{T}$.
$W_{|\varphi\rangle\langle\psi|}(\mathbf{r})$ corresponds to the projector
$|\varphi\rangle\langle\psi|$, e.g. to calculate $\mathsf{M}_{\text{HOM}}$ we
need the WF of the HOM operator:

$W_{\text{HOM}}(\mathbf{r})=\frac{1}{16\pi^{2}}\exp\left(-\frac{p_{\text{a}}^{2}+p_{\text{b}}^{2}+x_{\text{a}}^{2}+x_{\text{b}}^{2}}{2}\right)((p_{\text{a}}+p_{\text{b}})^{2}+(x_{\text{a}}+x_{\text{b}})^{2}-2)((p_{\text{a}}-p_{\text{b}})^{2}+(x_{\text{a}}-x_{\text{b}})^{2}-2).$
(154)

The WF approach is convenient when we know the exact WF of the input state and
how the transformation changes the quadratures of the oscillators. This
approach is also suitable for the case of
$|n\rangle_{\text{a}}|m\rangle_{\text{b}}$ input. We demonstrate it for the
gate transformation with $|1\rangle_{\text{a}}|1\rangle_{\text{b}}$ input. Let
us assume that both quantum oscillators were initially in single-boson states.
Since they are independent, the exact Wigner function of the initial state of
the system can be obtained by multiplying of the two single-photon state
Wigner functions:

$\displaystyle
W_{\text{in}}(\mathbf{r})=W_{\text{a}}|_{|1\rangle}(\mathbf{r}_{\text{a}})W_{\text{b}}|_{|1\rangle}(\mathbf{r}_{\text{b}}),\;\;\text{where}\;\;W_{\text{a,b
}}|_{|1\rangle}(\mathbf{r}_{\text{a,b
}})=\frac{1}{2\pi}\left(\mathbf{r}_{\text{a,b
}}^{2}-1\right)\exp\left(-\frac{\mathbf{r}_{\text{a,b }}^{2}}{2}\right),$
(155) $\displaystyle\mathbf{r}_{\text{a,b }}=(x_{\text{a,b }},p_{\text{a,b
}})^{T}.$ (156)

The WF of the output state of the system for the ideal gate:

$\displaystyle
W_{\text{out}}(\mathbf{r})=\frac{1}{4\pi^{2}}((p_{\text{b}}+\mathsf{G}p_{\text{a}})^{2}+x_{\text{b}}^{2}-1)(p_{\text{a}}^{2}+(x_{\text{a}}-\mathsf{G}x_{\text{b}})^{2}-1)\times$
$\displaystyle\times\exp\left(-\frac{(p_{\text{b}}+\mathsf{G}p_{\text{a}})^{2}+x_{\text{b}}^{2}+p_{\text{a}}^{2}+(x_{\text{a}}-\mathsf{G}x_{\text{b}})^{2}}{2}\right).$
(157)

Using Eq.(153), we will obtain the same result as in Eq.( 23) for
$p_{\text{a}}=p_{\text{b}}=1$. Both described approaches are identical and,
being applied to any ideal transformations, they give the same results.

For the non-ideal case we need to take to account the noises, which is
possible using the language of covariance matrices. A Gaussian quantum state
(such as vacuum, coherent, squeezed or thermal states) can be fully described
by the first and second statistical moments, that is, a vector of means and a
covariance matrix. The Wigner function of an arbitrary Gaussian state in
$\mathbb{R}_{s}$-dimension can be represented using the covariance matrices
$V$ as:

$\displaystyle W(\mathbf{r})=\text{W}(V,\mathbf{r})=\frac{1}{2\pi\sqrt{\det
V}}\exp\left(-\frac{1}{2}(\mathbf{r}-\bar{\mathbf{r}})^{T}V^{-1}(\mathbf{r}-\bar{\mathbf{r}})\right),$
(158)
$\displaystyle\mathbf{r}=\left(x_{1},p_{1},...,x_{s},p_{s}\right),\;\bar{\mathbf{r}}=\left(\langle
X_{1}\rangle,\langle P_{1}\rangle,...,\langle X_{s}\rangle,\langle
P_{s}\rangle\right).$

After a nonideal transformation $T$ of the vector $\mathbf{r}$, the evolution
of the covariance matrices and the WF are as follows:

$\displaystyle
W_{\text{in,out}}(\mathbf{r})=\text{W}(\mathbf{r},V_{\text{in,out}}),\qquad\text{where}\;V^{\text{out}}=T_{\mathsf{G}}V^{\text{in}}T_{\mathsf{G}}^{T}+V_{\text{N}},$
(159) $\displaystyle[V_{N}]_{ij}=\frac{1}{2}\langle
N_{i}N_{j}+N_{j}N_{i}\rangle\quad(i,j=1,...,s).$ (160)

Thus, we can use Eq.(153) and obtain the HOM matrix element of the output
state for the nonideal transformation.

However, neither single-photon nor HOM states are not Gaussian. Nevertheless,
the single-photon state can be approximated by the superposition of the
thermal $W_{\text{th}}(\mathbf{r})$ and vacuum $W_{|0\rangle}(\mathbf{r})$
states as:

$\displaystyle
W|_{|1\rangle}(\mathbf{r})\approx\frac{1}{n}\left((n+1)W_{\text{th}}(\mathbf{r})-W_{|0\rangle}(\mathbf{r})\right),\qquad
n\ll 1.$ (161)

Here $W_{\text{th}}$ is the Wigner function of a thermal state with mean
occupation $n$, that is a Gaussian state with zero means and covariance matrix
equal to $V_{\text{th}}=(2n+1)\mathbb{I}_{2\times 2}.$

Thermal and vacuum states are Gaussian, and, in turn, they can be represented
using Eq. (158) with $\bar{\mathbf{r}}=0$ and $s=4$. Thus, to describe the
initial state at the input of the gate between a-oscillator and b-oscillator
we can use the following approximated WF:

$\displaystyle
W_{\text{in}}(\mathbf{r})=W_{\text{a}}|_{|1\rangle}(\mathbf{r}_{\text{a}})W_{\text{b}}|_{|1\rangle}(\mathbf{r}_{\text{b}})\approx\frac{1}{n}\left((n+1)W_{\text{th}}(\mathbf{r}_{\text{a}})-W_{|0\rangle}(\mathbf{r}_{\text{a}})\right)\frac{1}{n}\left((n+1)W_{\text{th}}(\mathbf{r}_{\text{b}})-W_{|0\rangle}(\mathbf{r}_{\text{b}})\right)=$
$\displaystyle=\frac{1}{n^{2}}\left((n+1)^{2}\text{W}(V^{\text{in}}_{\text{th}_{\text{a}},\text{th}_{\text{b}}},\mathbf{r})-(n+1)\text{W}(V^{\text{in}}_{\text{vac}_{\text{a}},\text{th}_{\text{b}}},\mathbf{r})-(n+1)\text{W}(V^{\text{in}}_{\text{th}_{\text{a}},\text{vac}_{\text{b}}},\mathbf{r})+\text{W}(V^{\text{in}}_{\text{vac}_{\text{a}},\text{vac}_{\text{b}}},\mathbf{r})\right)=$
$\displaystyle=\frac{1}{n^{2}}\sum\limits_{k,l=0,1}(-(n+1))^{k+l}\text{W}(V^{\text{in}}_{\text{k,l}},\mathbf{r}).$
(162)

Here, the corresponding covariance matrices are

$\displaystyle
V^{\text{in}}_{1,1}=(2n+1)\mathbb{I}_{4},\;V^{\text{in}}_{0,0}=\mathbb{I}_{4},\;V^{\text{in}}_{1,0}=\begin{pmatrix}2n+1&0&0&0\\\
0&2n+1&0&0\\\ 0&0&1&0\\\
0&0&0&1\end{pmatrix},\;V^{\text{in}}_{0,1}=\left(\begin{array}[]{cccc}1&0&0&0\\\
0&1&0&0\\\ 0&0&2n+1&0\\\ 0&0&0&2n+1\end{array}\right),$ (167) $\displaystyle
V^{\text{out}}_{k,l}=T_{\mathsf{G}}V^{\text{in}}_{k,l}T_{\mathsf{G}}^{T}+V_{\text{N}},\qquad[V_{\text{N}}]_{ij}=\frac{1}{2}\langle
N_{i}N_{j}+N_{j}N_{i}\rangle,$ (168)

where $\mathbb{I}_{4}$ is the identity matrix of size $4$, $k=0,1;\;l=0,1$ and
$i,j=1,...,4$.

After the gate transformation of the vector $\mathbf{r}$, the WF of the state
at the output of the gate:

$\displaystyle
W_{\text{out}}(\mathbf{r})=\frac{1}{n^{2}}\sum\limits_{k,l=0,1}(-(n+1))^{k+l}\text{W}(V^{\text{out}}_{k,l},\mathbf{r}).$
(169)

Let us derive the approximated WF of the HOM state. Since
$W_{\text{HOM}}(\mathbf{r})$ is the WF of the state that would be at the
output of the 1:1 beamsplitter if at the input there were two single-photon
states, we use the same approach as we used to obtain
$W_{\text{out}}(\mathbf{r})$ and get

$\displaystyle
W_{\text{HOM}}(\mathbf{r})=\frac{1}{n^{2}}\sum\limits_{k,l=0,1}(-(n+1))^{k+l}W(V^{\prime}_{k,l},\mathbf{r}),$
(170) $\displaystyle
V^{\prime}_{k,l}=T_{\text{bs}}V^{\text{in}}_{k,l}T_{\text{bs}}^{T},\qquad
T_{\text{bs}}=\frac{\sqrt{2}}{2}\left(\begin{array}[]{cccc}1&0&1&0\\\
0&1&0&1\\\ -1&0&1&0\\\ 0&-1&0&1\end{array}\right).$ (175)

To calculate the HOM matrix element let us use the rule $\int
d\mathbf{r}\;\exp\\{-\frac{1}{2}\mathbf{r}^{T}V^{-1}\mathbf{r}\\}=4\pi^{2}\sqrt{\text{det}[V]}$
that works for every symmetric positive-definite matrix $V$. Then, keeping in
mind
$\mathbf{r}^{T}(A+B)\mathbf{r}=\mathbf{r}^{T}A\mathbf{r}+\mathbf{r}^{T}B\mathbf{r}$,
we can calculate the HOM matrix element as:

$\displaystyle\mathsf{M}_{\text{HOM}}=\frac{4}{n^{4}}\sum\limits_{k,l,m,d=0,1}(-(n+1))^{k+l+m+d}\cdot
F(V^{\prime}_{k,l},V^{\text{out}}_{m,d}),\qquad
F(A,B)=\left(\sqrt{\det[A+B]}\right)^{-1}.$ (176)

Using the same approach we can derive the HOM element for the case of the the
mixture state input:

$\displaystyle\mathsf{M}_{\text{HOM}}=\frac{4}{n^{4}}p^{2}\sum\limits_{k,l,m,d=0,1}(-(n+1))^{k+l+m+d}\cdot
F(V^{\prime}_{k,l},V^{\text{out}}_{m,d})+$
$\displaystyle+\frac{4}{n^{3}}p(1-p)\sum\limits_{k,l,m=0,1}(-(n+1))^{k+l+m}\cdot
F(V^{\prime}_{k,l},V^{\text{out}}_{m,0})+$
$\displaystyle+\frac{4}{n^{3}}p(1-p)\sum\limits_{k,l,d=0,1}(-(n+1))^{k+l+d}\cdot
F(V^{\prime}_{k,l},V^{\text{out}}_{0,d})+$
$\displaystyle+\frac{4}{n^{2}}(1-p)^{2}\sum\limits_{k,l=0,1}(-(n+1))^{k+l}\cdot
F(V^{\prime}_{k,l},V^{\text{out}}_{0,0}).$ (178)

###  B.5 Input and output thresholds (nonclassicality borders)

We define two nonclassicality thresholds by evaluating the maximum of
$\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$ over (i) all
superpositions of the coherent states at the output of the QND interaction:
$\rho_{\text{out}}=\rho_{\text{coh}}=\left|\alpha_{\text{a}}\beta_{\text{b}}\middle\rangle\middle\langle\alpha_{\text{a}}\beta_{\text{b}}\right|$
and (ii) before the QND interaction:
$\rho_{\text{out}}=U_{\text{QND}}\rho_{\text{coh}}U_{\text{QND}}^{\dagger}$.
For the BS interaction, such thresholds coincide.

If the output state $\rho_{\text{out}}$ is classical, i.e. is a mixture of
coherent states, then it turns out that
$0\leq\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle\leq 1/e^{2}$.
Thus, $1/e^{2}$ is the output threshold for the HOM interference. That is,
when measuring the HOM element, if we get a value greater than $1/e^{2}$, then
the measured state for sure is a non-classical one. This threshold is shown by
a thin gray dashed line in the Fig. 4(a).

To derive the _input threshold_ , let us use two random coherent states as the
input states of the gate and calculate the HOM element for the output state
$\rho^{\text{coh}}_{\text{out}}=U\rho^{\text{coh}}U^{\dagger}$. We assume
phase-randomized input state which means that phases of the input coherent
states are averaged. If we examine the dependence of the
$\langle\text{HOM}|\rho^{\text{coh}}_{\text{out}}|\text{HOM}\rangle$ on the
gain $\mathsf{G}$ for all the coherent states (see Fig. 4(b)) we will obtain
the area restricted by an input threshold (blue curve) that has a specific
complex shape.

To calculate the input threshold (phase randomized) let us take two random
coherent states as the inputs of the QND gate. The WF of the input state (two
independent coherent states) is the following

$\displaystyle
W_{in}(\mathbf{r},\mathbf{R})=\frac{1}{4\pi^{2}\sqrt{\text{det}V_{in}}}\exp\left(-\frac{1}{2}(\mathbf{r}-\mathbf{R})^{T}V_{in}^{-1}(\mathbf{r}-\mathbf{R})\right).$
(179)

Here
$\mathbf{r}=\left(x_{\text{a}},p_{\text{a}},x_{\text{b}},p_{\text{b}}\right)^{T}$,
$\mathbf{R}\equiv\left(X_{\text{a}},Y_{\text{a}},X_{\text{b}},Y_{\text{b}}\right)^{T}$
is the colomn-vector of means and $V_{in}=\mathbb{I}_{4\times 4}$ is the
covariance matrix of the initial state (coherent). The vector of means changes
as $\mathbf{R}\rightarrow T_{\mathsf{G}}\mathbf{R}$. Thus, the WF of the
output state is as follows:

$\displaystyle
W_{out}(\mathbf{r},\mathbf{R},\mathsf{G})=\frac{1}{4\pi^{2}\sqrt{\text{det}V_{\text{out}}}}\exp\left(-\frac{1}{2}(\mathbf{r}-T_{\mathsf{G}}\mathbf{R})^{T}V_{\text{out}}^{-1}(\mathbf{r}-T_{\mathsf{G}}\mathbf{R})\right).$
(180)

Using this WF we can obtain the HOM element as

$\displaystyle\mathsf{M}_{\text{HOM}}(\mathbf{R},\mathsf{G})=\langle\text{HOM}|\rho_{\text{out}}\text{HOM}\rangle\equiv(4\pi)^{2}\iiiint
d\mathbf{r}\;W_{\text{HOM}}(\mathbf{r})\cdot
W_{\text{out}}(\mathbf{r},\mathbf{R},\mathsf{G}).$ (181)

If, calculating the matrix elements for $\rho^{coh}_{out}$ over all the
possible coherent states averaged over phases, then the range of possible
values for the matrix elements will significantly change. Assuming
$\mathbf{R}=\left(X_{\text{a}},Y_{\text{a}},X_{\text{b}},Y_{\text{b}}\right)^{T}=\left(R_{\text{a}}\cos[\varphi_{\text{a}}],R_{\text{a}}\sin[\varphi_{\text{a}}],R_{\text{b}}\cos[\varphi_{\text{b}}],R_{\text{b}}\sin[\varphi_{\text{b}}]\right)^{T}$
we can obtain

$\displaystyle\mathsf{M}^{\text{av}}_{\text{HOM}}(R_{\text{a}},R_{\text{b}},\mathsf{G})=\frac{1}{4\pi^{2}}\iint
d\varphi_{\text{a}}\;d\varphi_{\text{b}}\;\;\mathsf{M}_{\text{HOM}}(\mathbf{R},\mathsf{G})$
(182)

and investigate it over all possible $R_{\text{a}},R_{\text{b}}$ for the
certain $\mathsf{G}$.

a

b

Figure 8:  a) Shape of the input coherent threshold. b)
$\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$ matrix element of the
output state for the ideal QND gate as a function of the gain $\mathsf{G}$
calculated for the different cases of the input: quantum input
$|1\rangle_{\text{a}}|1\rangle_{\text{b}}$ (solid black curves), mixture input
Eq. (22) (dashed black curves, dashing scale indicates parameter $p$). Dashed
gray line is the output threshold. Blue curve is the input thresholds (phase
randomized) restricting area that covers all the possible values of the matrix
elements of the output state of the gate in the case of the random coherent
input with averaged phases.

Figure 8(b) demonstrates
$\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$ depending on the gain,
assuming $p_{\text{a}}=p_{\text{b}}=p$, compared with the case of the pure
input $|1\rangle_{a}|1\rangle_{b}$. Expectedly, as the parameter $p$
decreases, the contribution from $|1\rangle_{\text{a}}|1\rangle_{\text{b}}$
term decreases, while the contribution from
$|0\rangle_{\text{a}}|0\rangle_{\text{b}}$ term increases. Visually, it is
reflected in the gradual change of the curves’ shape – for relatively high $p$
maximum first decreases, then smoothly shifts to the right. Thus, the maximum
of the HOM element decreases from $0.26$ at $p=1$ to $1/e^{2}$ at $p\approx
0.7$, which corresponds to the output threshold. At $p\approx 0.48$ it already
crosses the input threshold, so for $p<0.48$ the HOM element lies below the
input threshold at the gain $\mathsf{G}\approx 0.87$. For $p<0.40$ the HOM
element lies below the input threshold for any gain.

Figure 9:  $\langle\text{HOM}|\rho_{\text{out}}|\text{HOM}\rangle$ matrix
element of the output state for the ideal QND gate as a function of the gain
$\mathsf{G}$ calculated for the different cases of the input: the vacuum input
$|0\rangle_{\text{a}}|0\rangle_{\text{b}}$(black dashed), the state with one-
boson at the first input and coherent at the second
$|1\rangle_{\text{a}}|\beta\rangle_{\text{b}}$ (red), two independent coherent
states at the input $|\alpha\rangle_{\text{a}}|\beta\rangle_{\text{b}}$ (blue
dot-dashed), phase averaged coherent states at the input
$|\alpha\rangle_{\text{a}}|\beta\rangle_{\text{b}}$ (blue), two-mode squeezed
state TMS (green). The curve of the vacuum input depends only on the gain, all
other inputs depend on many parameters (phases, displacements, squeezings
etc), so their curves are the highest possible values (thresholds) of the HOM
element calculated over all their parameters.