Dual-path methods for propagating quantum microwaves

In the quantum regime, a beam splitter is a linear device that superposes two input signals (see figure 1(a)). Here, we consider only 50:50 beam splitters with input–output relations

$$\begin{equation} \left(\begin{array}{@{}l@{}} \hat{a}_1\\ \hat{a}_2 \end{array}\right)= \frac{1}{\sqrt{2}} \left(\begin{array}{@{}ll@{}} 1 & 1\\ -1 & 1 \end{array}\right) \left(\begin{array}{@{}l@{}} \hat{a}\\ \hat{v} \end{array}\right), \end{equation} \tag{ 1 }$$

where $\hat{a}$, $\hat{v}$ are the input signal operators and $\hat{a}_1$, $\hat{a}_2$ are the output signal operators, obeying the bosonic commutation relations, e.g. $[\hat{a},\hat{a}^{\dagger}]=1$. In the case of propagating signals in the frequency range of 1–10 GHz, hybrid ring structures [26, 27] have already been used as beam splitters in several experiments.

Zoom In Zoom Out Reset image size

In our scheme, we use devices that amplify all quadratures in the same way, the so-called phase-insensitive amplifiers (represented in figure 1(b)). The input–output relations for these devices are [28]

$$\begin{equation} \hat{a}'=\sqrt{g}\,\hat{a}+\sqrt{g-1}\,\hat{h}_{\rm{amp}}^{\dagger}, \end{equation} \tag{ 2 }$$

where $\hat{h}_{\rm {amp}}$ is a bosonic operator needed in order to make $\hat{a}'$ to fulfil the bosonic commutation relations and g > 1 is the power gain. An example for devices working in our frequency range and providing sufficient gain are high-electron-mobility transistor (HEMT) amplifiers. In these devices, the amount of noise added depends on the design, material system and operation temperature, but $\hat{h}_{\rm {amp}}$ is typically associated with a state with 10–20 photons.

In quantum optics, a quadrature measurement can be implemented via a homodyne technique using two photodetectors. This method consists in superimposing the signal with a strong coherent state on a 50:50 beamsplitter, measuring the number of photons at the outputs and calculating their difference. It involves photon-number measurements, which can be implemented in quantum optics using photodetectors, and it allows for single shot quadrature measurements. In the microwave regime, efficient, low-noise photodetectors are not readily available. However, if the signal is strong enough, we can measure its quadratures using an in-phase-quadrature (IQ) mixer [29]. The IQ mixer allows us to measure simultaneously both quadratures at the cost of adding noise. It can be modelled using a beam splitter and an ancilla $\hat{v}_{\rm {IQ}}$ in a thermal state [30], e.g. the vacuum (see figure 1(c)). From the input–output relations (1), we see that the quadrature operators $\skew3\hat{X}_1\equiv \sqrt {2}\hat{x}_{1}$ and $\skew3\hat{P}_2\equiv \sqrt {2}\skew3\hat{p}_{2}$ (where 1 and 2 label the two outputs of the beam splitter) have the same first moment of $\hat{x}_a$ and $\skew3\hat{p}_a$ respectively, but different variances

$$\begin{equation} \langle\skew3\hat{X}_1\rangle=\langle\hat{x}_a\rangle,\quad \Delta \skew3\hat{X}_1^2=\Delta \hat{x}_{a}^2+\Delta \hat{x}_{v_{\rm{IQ}}}^2, \end{equation} \tag{ 3 }$$

$$\begin{equation} \langle\skew3\hat{P}_2\rangle=\langle\skew3\hat{p}_a\rangle,\quad \Delta \skew3\hat{P}_2^2=\Delta \skew3\hat{p}_{a}^2+\Delta \skew3\hat{p}_{v_{\rm{IQ}}}^2. \end{equation} \tag{ 4 }$$

Here, $\Delta \skew3\hat{A}^2\equiv \langle \skew3\hat{A}^2\rangle -\langle \skew3\hat{A}\rangle ^2$ indicates the variance of the observable $\skew3\hat{A}$, and we have defined the quadratures as $\hat{x}\equiv (\hat{a}+\hat{a}^{\dagger })/\sqrt {2}$ and $\skew3\hat{p}\equiv -\mathrm {i}(\hat{a}-\hat{a}^{\dagger })/\sqrt {2}$.

Suppose we have a microwave signal $\hat{a}$ and we want to know its state. As the interesting signals are too weak to be detected, we need to amplify before performing a measurement, as in equation (2) (see figure 2). In the SPM, we first send a known state to characterize $\hat{h}_{\rm {amp}}$, and then we measure with $\hat{a}$ as input. If we send, for instance, a known coherent state |α〉 as input, we can retrieve the moments of $\hat{h}_{\rm {amp}}$ from the moments of the outcoming signals in a recursive way

$$\begin{eqnarray} &&\fl \langle\hat{h}_{\rm{amp}}^{ l}\hat{h}_{\rm{amp}}^{\dagger m}\rangle =\frac{\langle\hat{a}'^{\dagger l}\hat{a}'^{m}\rangle}{(g-1)^{\frac{l+m}{2}}}-\sum_{i_1=0}^l\sum_{i_2=0}^{m-1}\left({l}\atop{i_1}\right)\,\left({m}\atop{i_2}\right)\,\left(\frac{g}{g-1}\right)^{\frac{l+m-(i_1+i_2)}{2}}(\alpha^{*})^{l-i_1}\alpha^{m-i_2}\langle\hat{h}_{\rm{amp}}^{i_1}\hat{h}_{\rm{amp}}^{\dagger i_2}\rangle\nonumber\\ &&-\sum_{i_1=0}^{l-1}\left({l}\atop{i_1}\right)\,\left(\frac{g}{g-1}\right)^{\frac{l-i_1}{2}}(\alpha^{*})^{l-i_1}\langle\hat{h}_{\rm{amp}}^{ i_1}\hat{h}_{\rm{amp}}^{\dagger m}\rangle. \end{eqnarray} \tag{ 5 }$$

The gain g can, in principle, be found by comparing the first moments ($g=|\langle \hat{a}'\rangle /\alpha |^2$, assuming α ≠ 0). However, in reality this is difficult since in a typical experiment situation α is not well known. Hence, g is determined in a calibration experiment with thermal states [14, 26]. Once the calibration is done, usually the vacuum (α = 0) is chosen as reference state. This leads to the formula

$$\begin{equation} \langle\hat{h}_{\rm{amp}}^{ l}\hat{h}_{\rm{amp}}^{\dagger m}\rangle=\frac{\langle\hat{a}'^{\dagger l}\hat{a}'^{m}\rangle}{(g-1)^{\frac{l+m}{2}}}. \end{equation} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

Once we have characterized the amplifier noise, we can perform tomography of a general quantum state $\hat{a}$. In fact, by inverting equation (2), we obtain

$$\begin{eqnarray} \langle\hat{a}^{\dagger l}\hat{a}^{ m}\rangle &=&\frac{\langle\hat{a}'^{\dagger l}\hat{a}'^{m}\rangle}{g^{\frac{l+m}{2}}}-\sum_{i_1=0}^l\sum_{i_2=0}^{m-1}\left({l}\atop{i_1}\right)\left({m}\atop{i_2}\right)\left(\frac{g-1}{g}\right)^{\frac{l+m-(i_1+i_2)}{2}}\langle\hat{a}^{\dagger i_1}\hat{a}^{i_2}\rangle\langle\hat{h}_{\rm{amp}}^{ l-i_1}\hat{h}_{\rm{amp}}^{\dagger m-i_2}\rangle\nonumber \\ \quad& &-\sum_{i_1=0}^{l-1}\left({l}\atop{i_1}\right)\left(\frac{g-1}{g}\right)^{\frac{l-i_1}{2}}\langle\hat{a}^{\dagger i_1}\hat{a}^m\rangle\langle\hat{h}_{\rm{amp}}^{ l-i_1}\rangle, \end{eqnarray} \tag{ 7 }$$

where $\langle \hat{h}_{\rm {ampl}}^{l}\skew3\hat h_{\rm {ampl}}^{\dagger m}\rangle$ is derived in equation (5).

In an experiment, the channel is lossy. We can model the losses with a beam splitter

$$\begin{equation} \hat{b}_{\rm{out}}=\sqrt{\eta} \,\hat{b}_{\rm{in}} +\sqrt{1-\eta}\,\hat{h}_{\rm{loss}}, \end{equation} \tag{ 8 }$$

where 0 < η < 1 is the power loss, and $\hat{h}_{\rm {loss}}$ is the noise added by the environment, strictly related to its temperature. If we have losses and amplification in series as in figure 2, we can consider the channel as amplifying or lossy, depending on the value of gη. For the detection of quantum microwaves, a large value of the effective gain gη on the order of 10⁴ is required, so that the amplified signal amplitudes are well above the noise added by the next components. Firstly, we consider the case that the losses occur before the signal is amplified (see figure 2(a)) and rewrite equation (2) as

$$\begin{equation} \hat{a}'=\sqrt{g\eta}\hat{a}+\sqrt{g\eta-1}\hat{V}^{\dagger} \end{equation} \tag{ 9 }$$

with gη > 1 and $\hat{V}=\sqrt {\frac {1-\eta }{\eta -\frac {1}{g}}}\hat{h}_{\rm {loss}}^{\dagger }+\sqrt {\frac {1-\frac {1}{g}}{\eta -\frac {1}{g}}}\hat{h}_{\rm {amp}}\simeq \sqrt {\frac {1}{\eta }-1}\hat{h}_{\rm {loss}}^{\dagger }+\sqrt {\frac {1}{\eta }}\hat{h}_{\rm {amp}}\equiv \hat{V}_{\rm {a}}$ in the limit g ≫ 1, obtaining an equivalent model, but with an effective gain and noise mode $\hat{V}$. Therefore, we can still use equation (5) to reconstruct the moments of $\hat{V}$ and equation (7) to retrieve the moments of $\hat{a}$. Secondly, if the losses occur after the amplification (see figure 2(b)), we get the same effective gain gη, but the noise mode $\hat{V}=\sqrt {\frac {1-\frac {1}{g}}{1-\frac {1}{g\eta }}}\hat{h}_{\rm {amp}}+\sqrt {\frac {1-\eta }{g\eta -1}}\hat{h}_{\rm {loss}}^{\dagger }\simeq \hat{h}_{\rm {amp}}\equiv \hat{V}_{\rm {b}}$ in the limit gη ≫ 1. We conclude that the lossy models are equivalent to the ideal one, but with different amplifier gains. Furthermore, we note that if the gain is large enough, the effect of the losses occurring after the amplification are negligible. In classical network theory, this is well-known [31]. In fact, assuming $\hat{h}_{\rm {amp}}$ and $\hat{h}_{\rm {loss}}$ in thermal states and g ≫ 1, we have that $\langle \hat{V}_{\rm {b}}^{\dagger }\hat{V}_{\rm {b}}\rangle =\langle \hat{h}_{\rm {amp}}^{\dagger }\hat{h}_{\rm {amp}}\rangle <(\frac {1}{\eta }-1)+(\frac {1}{\eta }-1)\langle \hat{h}_{\rm {loss}}^{\dagger }\hat{h}_{\rm {loss}}\rangle +\frac {1}{\eta }\langle \hat{h}_{\rm {amp}}^{\dagger }\hat{h}_{\rm {amp}}\rangle =\langle \hat{V}_{\rm {a}}^{\dagger } \hat{V}_{\rm {a}}\rangle$, ∀η∈(0,1), so the model in figure 2(b) contains less noise than the one in figure 2(a). Regarding the optimization of experimental setups, minimizing the losses occurring before amplification is crucial under the assumption that g is large enough, which usually holds. In the remainder of this work, we use the ideal amplifier model of (2), treating g and $\hat{h}_{\rm {amp}}$ as effective parameters.

We now describe the quantum formalism of the DPM, including explicit formulae for the reconstruction of the moments. Our method assumes minimum information on an ancilla, i.e. the knowledge of its first two moments, and it allows to reconstruct the moments of the input signal and of the amplifier noise fields at the same time. The scheme is depicted in figure 3. We send a signal $\hat{a}$ to a 50:50 beam splitter with an ancilla $\hat{v}$, and we direct the output fields ($\skew3\hat c_1$ and $\skew3\hat c_2$) of the beam splitter to two amplifiers, obtaining

$$\begin{equation} \skew3\hat c'_1=\frac{\sqrt{g}}{\sqrt{2}}(\hat{a}+\hat{v})+\sqrt{g-1}\hat{h}_{1}^{\dagger}, \end{equation} \tag{ 10 }$$

$$\begin{equation} \skew3\hat c'_2=\frac{\sqrt{g}}{\sqrt{2}}(-\hat{a}+\hat{v})+\sqrt{g-1}\hat{h}_{2}^{\dagger}, \end{equation} \tag{ 11 }$$

where we have assumed for simplicity the same gain g for both channels. Lastly, we measure the quadratures $\hat{x}$ and $\skew3\hat{p}$ of both channels with IQ mixers. It is useful to introduce the complex envelopes [30]

$$\begin{equation} \skew3\hat{S}_{1,2}\equiv \frac{\hat{x}_{1,2}+\mathrm{i}\skew3\hat{p}_{1,2}}{\sqrt{g}}=\frac{\skew3\hat c'_{1,2}+\hat{v}_{1,2}^{\dagger}}{\sqrt{g}}, \end{equation} \tag{ 12 }$$

where $\hat{x}_{1,2}$ and $\skew3\hat{p}_{1,2}$ are the quadrature operators for the channels 1 and 2, $\hat{v}_{1,2}$ are the ancillas in the IQ mixer model, and the last equality can be easily checked from the quadrature definition and the beam splitter relation [30]. Finally, we obtain

$$\begin{equation} \skew3\hat{S}_1=\frac{\hat{a}+\hat{v}}{\sqrt{2}}+\hat{V}_1^{\dagger},\quad \skew3\hat{S}_2=\frac{-\hat{a}+\hat{v}}{\sqrt{2}}+\hat{V}_2^{\dagger} \end{equation} \tag{ 13 }$$

with $\hat{V}_{1,2}=\sqrt {1-\frac {1}{g}}\hat{h}_{1,2}+\sqrt {\frac {1}{g}}\hat{v}_{1,2}\simeq \hat{h}_{1,2}$ in the limit g ≫ 1.

Zoom In Zoom Out Reset image size

We now iteratively find the relation between the moments of the detected signals and the moments of the input signal. For the first moment, we obtain

$$\begin{equation} \langle\skew3\hat{S}_1\rangle=\frac{\langle\hat{a}\rangle+\langle\hat{v}\rangle}{\sqrt{2}},\quad\langle\skew3\hat{S}_2\rangle=\frac{-\langle\hat{a}\rangle+\langle\hat{v}\rangle}{\sqrt{2}}, \end{equation} \tag{ 14 }$$

where we have assumed $\langle \hat{V}_1^{\dagger }\rangle =\langle \hat{V}_2^{\dagger }\rangle =0$. This is a safe assumption for any kind of noise. From (14) we derive

$$\begin{equation} \langle\hat{a}\rangle=\frac{1}{\sqrt{2}}(\langle\skew3\hat{S}_1\rangle-\langle\skew3\hat{S}_2\rangle), \end{equation} \tag{ 15 }$$

$$\begin{equation} \langle\hat{v}\rangle=\frac{1}{\sqrt{2}}(\langle\skew3\hat{S}_1\rangle+\langle\skew3\hat{S}_2\rangle). \end{equation} \tag{ 16 }$$

We assume the knowledge of the first two moments of the ancilla $\hat{v}$, so we leave (16) as a consistency check. From the second moments

$$\begin{equation} \langle\skew3\hat{S}_1^2\rangle=\frac{\langle\hat{a}^2\rangle}{2}+\frac{\langle\hat{v}^2\rangle}{2}+\langle\hat{a}\rangle\langle\hat{v}\rangle+\langle\hat{V}_1^{\dagger 2}\rangle, \end{equation} \tag{ 17 }$$

$$\begin{equation} \langle\skew3\hat{S}_2^2\rangle=\frac{\langle\hat{a}^2\rangle}{2}+\frac{\langle\hat{v}^2\rangle}{2}-\langle\hat{a}\rangle\langle\hat{v}\rangle+\langle\hat{V}_2^{\dagger 2}\rangle, \end{equation} \tag{ 18 }$$

$$\begin{equation} \langle\skew3\hat{S}_1\skew3\hat{S}_2\rangle=-\frac{\langle\hat{a}^2\rangle}{2}+\frac{\langle\hat{v}^2\rangle}{2}, \end{equation} \tag{ 19 }$$

we get by linear inversion

$$\begin{equation} \langle\hat{a}^2\rangle=-2\langle\skew3\hat{S}_1\skew3\hat{S}_2\rangle+\langle\hat{v}^2\rangle, \end{equation} \tag{ 20 }$$

$$\begin{equation} \langle\hat{V}_1^{\dagger 2}\rangle=\langle\skew3\hat{S}_1^2\rangle-\frac{\langle\hat{a}^2\rangle}{2}-\frac{\langle\hat{v}^2\rangle}{2}-\langle\hat{a}\rangle\langle\hat{v}\rangle, \end{equation} \tag{ 21 }$$

$$\begin{equation} \langle\hat{V}_2^{\dagger 2}\rangle=\langle\skew3\hat{S}_2^2\rangle-\frac{\langle\hat{a}^2\rangle}{2}-\frac{\langle\hat{v}^2\rangle}{2}+\langle\hat{a}\rangle\langle\hat{v}\rangle. \end{equation} \tag{ 22 }$$

Similar formulae hold for $\langle \hat{a}^{\dagger } \hat{a}\rangle$ and $\langle \hat{V}_{1,2}\hat{V}_{1,2}^{\dagger }\rangle$ (see the appendix). In the derivation of the above formulae, the mutual statistical independence of $\hat{a}$, $\hat{v}$, $\hat{V}_1$ and $\hat{V}_2$ is crucial. This is a reasonable assumption if we consider, for instance, that the two amplifiers are spatially well separated and, hence, have different environments.

We note that we have derived the moments of the input field via the cross-correlation terms, and, using this result, we have derived the moments of the noise via the auto-correlation terms. This calculation can be iterated over every order, and in this way, we can retrieve the moments of the input state and the noise. We also note that from the third moment on, we do not need further assumptions on the ancilla $\hat{v}$. For instance, the third moments of the measurement outcomes are

$$\begin{equation} \langle \skew3\hat{S}_1^3\rangle=\frac{\langle\hat{a}^3\rangle}{2\sqrt{2}}+\frac{\langle\hat{v}^3\rangle}{2\sqrt{2}}+\langle\hat{V}_1^{\dagger3}\rangle+\langle f_{\skew3\hat{S}_1^3 }(\hat{a},\hat{v}, \hat{V}_1)\rangle, \end{equation} \tag{ 23 }$$

$$\begin{equation} \langle\skew3\hat{S}_1^2\skew3\hat{S}_2\rangle=-\frac{\langle\hat{a}^3\rangle}{2\sqrt{2}}+\frac{\langle\hat{v}^3\rangle}{2\sqrt{2}}+\langle f_{\skew3\hat{S}_1^2\skew3\hat{S}_2 }(\hat{a},\hat{v}, \hat{V}_1,\hat{V}_2)\rangle, \end{equation} \tag{ 24 }$$

$$\begin{equation} \langle\skew3\hat{S}_1\skew3\hat{S}_2^2\rangle=\frac{\langle\hat{a}^3\rangle}{2\sqrt{2}}+\frac{\langle\hat{v}^3\rangle}{2\sqrt{2}}+\langle f_{\skew3\hat{S}_1\skew3\hat{S}_2^2 }(\hat{a},\hat{v}, \hat{V}_1,\hat{V}_2)\rangle, \end{equation} \tag{ 25 }$$

$$\begin{equation} \langle \skew3\hat{S}_2^3\rangle=-\frac{\langle\hat{a}^3\rangle}{2\sqrt{2}}+\frac{\langle\hat{v}^3\rangle}{2\sqrt{2}}+\langle\hat{V}_2^{\dagger3}\rangle+\langle f_{\skew3\hat{S}_2^3 }(\hat{a},\hat{v}, \hat{V}_2)\rangle, \end{equation} \tag{ 26 }$$

where the $f_{\ldots }$ are functions of the operators, and $\langle f_{\ldots }\rangle$ contain moments up to the second order. Making the right combinations, one derives

$$\begin{equation} \langle\hat{a}^3\rangle=\sqrt{2}(\langle\skew3\hat{S}_1\skew3\hat{S}_2^2\rangle-\langle\skew3\hat{S}_1^2\skew3\hat{S}_2\rangle+\langle f_{\skew3\hat{S}_1^2\skew3\hat{S}_2}\rangle-\langle f_{\skew3\hat{S}_1\skew3\hat{S}_2^2 }\rangle), \end{equation} \tag{ 27 }$$

$$\begin{equation} \langle\hat{v}^3\rangle=\sqrt{2}(\langle\skew3\hat{S}_1\skew3\hat{S}_2^2\rangle+\langle\skew3\hat{S}_1^2\skew3\hat{S}_2\rangle-\langle f_{\skew3\hat{S}_1^2\skew3\hat{S}_2}\rangle-\langle f_{\skew3\hat{S}_1\skew3\hat{S}_2^2 }\rangle), \end{equation} \tag{ 28 }$$

$$\begin{equation} \langle\hat{V}_1^{\dagger3}\rangle=\langle \skew3\hat{S}_1^3\rangle-\frac{\langle\hat{a}^3\rangle}{2\sqrt{2}}-\frac{\langle\hat{v}^3\rangle}{2\sqrt{2}}-\langle f_{\skew3\hat{S}_1^3 }\rangle, \end{equation} \tag{ 29 }$$

$$\begin{equation} \langle\hat{V}_2^{\dagger3}\rangle=\langle \skew3\hat{S}_2^3\rangle+\frac{\langle\hat{a}^3\rangle}{2\sqrt{2}}-\frac{\langle\hat{v}^3\rangle}{2\sqrt{2}}-\langle f_{\skew3\hat{S}_2^3 }\rangle. \end{equation} \tag{ 30 }$$

This calculation can be generalized to any higher-order moments (see the appendix). Summarizing, we present a method to simultaneously reconstruct the moments of an input state and the moments of the noise added in the channels, assuming a minimum information on an ancilla.

In this section, we compare the DPM to the SPM with respect to the reconstruction of the quantum state moments of an input signal. We numerically test the methods with a squeezed state with squeezing parameter ξ = 0.5i [32] as input state, corresponding to a photon number of 0.27. We model the amplifier noise fields and the ancilla states with thermal states with photon numbers n_amp and n_anc and restrict our analysis to the first four moments. We have simulated several experiments to study the dependence of the statistical uncertainty of the reconstructed signal moments on the number of measurements N and on the photon number of the amplifier noise n_amp. For fixed n_amp, we compare the two methods in the following way. For the SPM we perform N measurements allowing the reconstruction of the amplifier noise by using a vacuum state as input. Then, we simulate N measurement results for the reconstruction of the squeezed state chosen as input. For the DPM, instead, we simulate N measurements of each channel, that allow us to calculate all the needed auto- and cross-correlations to reconstruct the input state. To obtain the individual measurements at the outputs of the IQ mixers, we use random numbers obeying a multivariate Gaussian statistics. The latter is determined by a covariance matrix, which is calculated using the input–output relations of all components. For a certain N, we estimate the average values of the simulated quadrature moments by using the sample mean. For the DPM the sample mean is

$$\begin{equation} \langle {x'}_1^{k_1}{p'}_1^{l_1}{x'}_2^{k_2}{p'}_2^{l_2}\rangle_{\rm{est}}=\frac{1}{N}\sum_{j=1}^N{x'}_{1j}^{k_1}{p'}_{1j}^{l_1}{x'}_{2j}^{k_2}{p'}_{2j}^{l_2}, \end{equation} \tag{ 31 }$$

where j labels the single outcomes. For the SPM we apply an analogous formula. We use the reconstruction formulae of the SPM and DPM to retrieve the moments of the input states. In order to access the statistical properties of the two methods, the described procedure is repeated 5000 times (for each parameter combination). We divide the data into 20 blocks with 250 samples each. The 250 samples are used to calculate the standard deviation σ of the reconstructed signal moments for the DPM and SPM. Finally, we compare the two methods by calculating the ratio σ_DP/σ_SP. Thereby, a higher standard deviation means worse performance. The 20 blocks allow us to estimate the uncertainty of the obtained mean value of the ratio σ _DP/σ_SP.

In experiments at millikelvin temperatures n_anc is negligible small. Therefore, we set n_anc = 0 in all presented simulations. The results are depicted in figures 4–6, where we assume n_amp = 10 for all amplifiers. Additionally, we have assumed that the first moment of the noise is zero in the SPM. This is not needed in the single-path reconstruction, but it is necessary if we want to have a fair comparison to the DPM, as in the latter this assumption is fundamental. In figures 4–6, we see that the DPM and the SPM have the same performance regarding the first and second moments, while for the third and the fourth moments the dual-path has an edge. The result is not dependent on the number of measurements (see figure 4). For an individual moment, the standard deviation scales as $1/\sqrt {N}$ (data not shown). In figure 5 we note that the moment $\langle \hat{a}^{\dagger}{a}^{2}\rangle$ performs better than $\langle \hat{a}^{3}\rangle$ even if they are of the same order. This is because $\langle \hat{a}^{\dagger}a^2\rangle$ can be estimated in more ways than $\langle \hat{a}^{3}\rangle$, and we are considering the average over all the possible combinations (see the appendix). Figure 6 shows that for both methods the standard deviation scales exponentially with the moment order, but the rate of the DPM is lower.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As second numerical test, we fix the number of measurements to N = 10⁷, and we simulate the two methods with different n_amp. Throughout this work, we assume that n_amp is the same for all amplifiers. In figures 7(c), (d) and 8(c), (d), we observe that the DPM performs better than the SPM for the third and fourth moments for high n_amp, while for n_amp = 0 both methods have approximately the same performance. From the good agreement between the polynomial fit curves and the data in figure 8, we conclude that the standard deviation depends on the amplifier noise as σ∝n^k/2_amp for large n_amp, where k is the order of the considered moment. The exponent k/2 has also been observed in simulations of an early version of the DPM for sinusoidal signals [34].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We want to check and estimate the entanglement between the output modes of the beamsplitter $\skew3\hat{s}_1$ and $\skew3\hat{s}_2$, by minimal assumptions on the beam splitter. To achieve this goal, we can perform a reference-state method adapted for two channels. First, we send a known state to characterize the correlations between the two amplifier noise fields $\langle \hat{V}_1^{k_1}\hat{V}_1^{\dagger j_1}\hat{V}_2^{k_2}\hat{V}_2^{\dagger j_2}\rangle$ by making use of the known moments of the reference state $\langle \skew3\hat{s}_1^{\dagger l_1}\skew3\hat{s}_1^{m_1}\skew3\hat{s}_2^{\dagger l_2}\skew3\hat{s}_2^{m_2}\rangle$, and algebraically inverting the following equation:

$$\begin{eqnarray} \langle S_1^{\dagger l_1}S_1^{m_1} S_2^{\dagger l_2}S_2^{m_2}\rangle &=& \sum_{k_1=0}^{l_1}\sum_{k_2=0}^{l_2} \sum_{j_1=0}^{m_1}\sum_{j_2=0}^{m_2} {l_1\choose k_1}{l_2\choose k_2} {m_1\choose j_1}{m_2\choose j_2}\nonumber\\ &&\times\langle \hat{s}_1^{\dagger l_1-k_1}\hat{s}_1^{m_1-j_1} \hat{s}_2^{\dagger l_2-k_2}\hat{s}_2^{m_2-j_2} \rangle \langle \hat{V}_1^{k_1} \hat{V}_1^{\dagger j_1} \hat{V}_2^{k_2} \hat{V}_2^{\dagger j_2} \rangle\,. \end{eqnarray} \tag{ 32 }$$

Here, we have used the complex envelopes defined in (12)

$$\begin{equation} \skew3\hat{S}_1=\hat{s}_1+\hat{V}_1^{\dagger}, \end{equation} \tag{ 33 }$$

$$\begin{equation} \skew3\hat{S}_2=\hat{s}_2+\hat{V}_2^{\dagger}. \end{equation} \tag{ 34 }$$

We note that by measuring correlations, we can also evaluate whether the two amplifier noise contributions are correlated or not, while in the DPM it was assumed that they are uncorrelated. In the case of vacuum as inputs of the beam splitter, which is the most simple choice in the case of microwave experiments at millikelvin temperatures, we have simply

$$\begin{equation} \langle\hat{s}_1^{\dagger l_1-k_1}\hat{s}_1^{m_1-j_1}\hat{s}_2^{\dagger l_2-k_2}\hat{s}_2^{m_2-j_2}\rangle_{\rm{known}}= \delta_{l_1,k_1}\delta_{m_2,j_2}\delta_{m_1,j_1}\delta_{l_2,k_2}. \end{equation} \tag{ 35 }$$

Then, we send the desired input signal and measure the correlations between $\skew3\hat{S}_1$ and $\skew3\hat{S}_2$. Inverting (32) with respect to $\langle \skew3\hat{s}_1^{\dagger l_1-k_1}\skew3\hat{s}_1^{m_1-j_1}\skew3\hat{s}_2^{\dagger l_2-k_2}\skew3\hat{s}_2^{m_2-j_2}\rangle$, we find auto- and cross-correlations between the two modes $\skew3\hat{s}_1$ and $\skew3\hat{s}_2$. We note that we have not used any details of the beam splitter, such as input–output relations, in the derivation of equations (33)–(35) and thus treat it as a black box. We only assume that if the input states are vacuum states, also the output is a vacuum state. This assumption is well justified since a cold beam splitter is a passive device containing no energy sources [27, 29].

Once we have measured the correlations of the two output modes, we can use criteria based on moments to check if there is entanglement. These kind of criteria have been well investigated in [23–25]. As example, we consider an entanglement witness based on the entanglement of Gaussian states. Another issue is how to quantify entanglement starting from the measured moments. There are different entanglement measures based on different features, but all of them obey some basic properties [35]. The estimation of an entanglement measure for an arbitrary state is still an open problem. Analytical solutions have been proposed for some classes of states, e.g. Gaussian states [36] (see section 5.2). Furthermore, techniques have been proposed [37, 38] to determine a lower bound for the degree of entanglement even in the case when only a finite number of moments is available (incomplete tomography). Usually, they require the solution of a convex optimization problem, that can be done efficiently if the dimension of the Hilbert space is not too large. For instance, let us define the negativity [39]

$$\begin{equation} \mathcal{N}(\rho)\equiv\frac{\|\rho^{T_1}\|_1-1}{2}, \end{equation} \tag{ 36 }$$

where ρ is the density matrix of the two-mode state after the beam splitter and ρ^T₁ denotes the partial transpose with respect to the system 1. A lower bound on $\mathcal {N}(\rho )$ is given by solving the minimization problem

$$\begin{eqnarray} {\rm minimize} \quad&\mathcal{N}(\sigma)\nonumber\\ {\rm subject \ to}\quad & \langle\hat{s}_1^{\dagger k_1}\hat{s}_1^{k_2}\hat{s}_2^{\dagger k_3}\hat{s}_2^{k_4}\rangle_{\sigma}=\langle\hat{s}_1^{\dagger k_1}\hat{s}_1^{k_2}\hat{s}_2^{\dagger k_3}\hat{s}_2^{k_4}\rangle_{\rho},\nonumber \\ \\ \quad&\sigma\succeq0,\nonumber\\ \quad& \rm{Tr}\,\sigma=1, \nonumber \end{eqnarray} \tag{ 37 }$$

that can be reshaped as a semidefinite program [37].

Finally, we emphasize that the proposed detection method respects some basic criteria for a reliable experimental entanglement verification [40]. In particular, it does not make any assumption on the input state and it is independent of the state generation process.

In the experiment, we send a squeezed vacuum state and a vacuum state to a beam splitter. The squeezed state is generated by sending the vacuum signal produced by a 50 Ω resistor to a JPA operated in the degenerate mode. The operation point of the JPA is characterized by a non-degenerate signal gain of 5.1 dB. The two outputs of the beam splitter are amplified via two HEMT amplifiers working at 4 K. The measurements are performed using two IQ mixers working at room temperature. The data processing is realized using an FPGA logic [14]. The number of the measurements registered is ∼5 × 10⁹. In figure 9, we represent the reconstructed moments for the input signal $\hat{a}$. We compare these moments with the ones of the Gaussian state defined by the first two reconstructed moments. We see that the moments correspond very well to the ones expected for Gaussian states. We note that the reconstructed squeezed state is not pure, because we have $\Delta \hat{x}_{\theta }^2=0.179\pm 0.001$ and $\Delta \hat{x}_{\theta +\frac {\pi }{2}}^2=5.326\pm 0.006$, where $\hat{x}_{\theta }$ is the quadrature with minimum variance, and these values lead to $\Delta \hat{x}_{\theta } \Delta \hat{x}_{\theta +\frac {\pi }{2}}= 0.978\pm 0.006>0.5$. For this reason, four moments are not enough for full-tomography, so we simply check the Gaussianity of the state by comparing the third and fourth moment with their theoretical values (see figure 9). The squeezing obtained is 4.45 ± 0.03 dB with respect to the variance of the vacuum state [14], and this is a witness of non-classicality of the state [42] (we note that the squeezing depends only on the first two moments, which are very well estimated). The moments of the amplifier noise operators agree with a thermal state with an average photon number 12.239 ± 0.002 and 12.805 ± 0.002, similar to the value used for the simulations in section 3.3.

Zoom In Zoom Out Reset image size

In the next example, we use the extended reference-state method to detect the entanglement of the output of the beam splitter, in the same experimental conditions as in the previous example. As we want to prove that the beam splitter is an entangler, we do not assume its input–output relations. In this sense, we need first to send vacuum as input to reconstruct the HEMT amplifier noise, and, then, we measure the output moments sending a squeezed state as input. We can reconstruct the moments of the beam splitter output by using equation (32), and by using witnesses based on moments, we can assert whether there is entanglement. The witness that we have used is based on the analytical solution of the negativity for a two-mode Gaussian state, given by [36]

$$\begin{equation} \mathcal{N}(\rho_{\rm{G}})=\max\left\{0,\frac{1-\nu}{2\nu}\right\}\equiv\max\{0,\widetilde{ \mathcal{N}_{\rm K}}\}, \end{equation} \tag{ 38 }$$

where ρ_G indicates a Gaussian state, and ν is the smallest symplectic value of the two-mode Gaussian state [36] (see supplementary material of [14] for an explicit expression). Here, the negativity kernel $\widetilde {\mathcal {N}_{\mathrm { K}}}$ [14] is an entanglement witness, given that if a non-Gaussian state has the same first two moments as an entangled Gaussian state, then it is entangled [43]. If we assume the Gaussianity of the state, as it is in our case, then $\widetilde {\mathcal {N}_{\mathrm { K}}}$ is directly related to an entanglement measure via equation (38). By using the same sample data of the previous example, we obtain a negativity of 0.489 ± 0.004.