Distributed quantum sensing in a continuous variable entangled network

TL;DR

This paper demonstrates that a four-node entangled continuous variable quantum network can achieve phase sensing precision beyond separable probes, with potential applications in biological imaging and atomic clock networks.

Contribution

It experimentally shows how entangled CV states enable enhanced distributed quantum phase sensing across multiple nodes.

Findings

Achieved quantum phase sensing with precision surpassing separable probes.

Utilized a four-mode entangled continuous variable state.

Demonstrated deterministic sensing in a multi-node network.

Abstract

Networking plays a ubiquitous role in quantum technology. It is an integral part of quantum communication and has significant potential for upscaling quantum computer technologies that are otherwise not scalable. Recently, it was realized that sensing of multiple spatially distributed parameters may also benefit from an entangled quantum network. Here we experimentally demonstrate how sensing of an averaged phase shift among four distributed nodes benefits from an entangled quantum network. Using a four-mode entangled continuous variable (CV) state, we demonstrate deterministic quantum phase sensing with a precision beyond what is attainable with separable probes. The techniques behind this result can have direct applications in a number of primitives ranging from biological imaging to quantum networks of atomic clocks.

Tables

Table 1. Table 1: The phase calibration result

Squeezing
k1	-3.96	k2	-3.97	k3	-3.95	k4	-3.96
b1	-0.13	b2	-0.59	b3	-0.64	b4	0.37
Anti-squeezing
k1	-3.99	k2	-3.99	k3	-4.06	k4	-4.06
b1	-89.50	b2	-89.97	b3	-89.87	b4	-88.90

Equations

$$\sigma =\frac{\sqrt{⟨\Delta {\hat{P}}_{\mathrm{avg}}^2⟩}}{|\partial ⟨{\hat{P}}_{\mathrm{avg}}⟩/\partial {\varphi }_{\mathrm{avg}}|},$$

$${\sigma }_s$$

$${\sigma }_e$$

$${\sigma }_s^{\mathrm{opt}}$$

$${\sigma }_e^{\mathrm{opt}}$$

$${\hat{p}}_j=\left(\sqrt{\eta }{\hat{x}}_{s,j}+\sqrt{1-\eta }{\hat{x}}_{\mathrm{vac},j}\right)\sin {\varphi }_j+\left(\sqrt{\eta }{\hat{p}}_{s,j}+\sqrt{1-\eta }{\hat{p}}_{\mathrm{vac},j}\right)\cos {\varphi }_j,$$

$$⟨{\hat{p}}_j⟩=\sqrt{\eta }⟨{\hat{x}}_{s,j}⟩\sin {\varphi }_j=\sqrt{2\eta }{\alpha }_s\sin {\varphi }_j\approx \sqrt{2\eta }{\alpha }_s{\varphi }_j.$$

$$⟨{\hat{P}}_{\mathrm{avg}}⟩\approx \sqrt{2\eta }{\alpha }_s{\varphi }_{\mathrm{avg}}.$$

$${\sigma }_s=\frac{\sqrt{⟨\Delta {\hat{P}}_{\mathrm{avg}}^2⟩}}{\left|\partial ⟨{\hat{P}}_{\mathrm{avg}}⟩/\partial {\varphi }_{\mathrm{avg}}\right|}.$$

$$\partial ⟨{\hat{P}}_{\mathrm{avg}}⟩/\partial {\varphi }_{\mathrm{avg}}\approx \sqrt{2\eta }{\alpha }_s,$$

$$⟨\Delta {\hat{P}}_{\mathrm{avg}}^2⟩$$

$$=\frac{1}{M^2}\sum _{j=1}^M\left({\sin }^2{\varphi }_j(\eta ⟨\Delta {\hat{x}}_{s,j}^2⟩+(1-\eta )⟨\Delta {\hat{x}}_{\mathrm{vac},j}^2⟩)+{\cos }^2{\varphi }_j(\eta ⟨\Delta {\hat{p}}_{s,j}^2⟩+(1-\eta )⟨\Delta {\hat{p}}_{\mathrm{vac},j}^2⟩)\right)$$

$$=\frac{1}{M^2}\sum _{j=1}^M\left(\frac{\eta e^{2r_s}}{2}{\sin }^2{\varphi }_j+\frac{\eta e^{-2r_s}}{2}{\cos }^2{\varphi }_j+\frac{1-\eta }{2}\right).$$

$$⟨\Delta {\hat{P}}_{\mathrm{avg}}^2⟩\approx \frac{1}{M}\left(\eta ⟨\Delta {\hat{p}}_{s,j}^2⟩+\frac{1-\eta }{2}\right)=\frac{\eta e^{-2r_s}+1-\eta }{2M}.$$

$${\sigma }_s=\frac{\sqrt{e^{-2r_s}+1/\eta -1}}{2{\alpha }_s\sqrt{M}}.$$

$$N_s=N_{s,\mathrm{coh}}+N_{s,\mathrm{sqz}}=\eta ({\alpha }_s^2+{\sinh }^2{r}_s).$$

$${\hat{p}}_j=\left(\sqrt{\eta }{\hat{x}}_j^{\prime }+\sqrt{1-\eta }{\hat{x}}_{\mathrm{vac},j}\right)\sin {\varphi }_j+\left(\sqrt{\eta }{\hat{p}}_j^{\prime }+\sqrt{1-\eta }{\hat{p}}_{\mathrm{vac},j}\right)\cos {\varphi }_j,$$

$${\hat{x}}_e=\frac{1}{\sqrt{M}}\sum _{j=1}^M{\hat{x}}_j^{\prime },{\hat{p}}_e=\frac{1}{\sqrt{M}}\sum _{j=1}^M{\hat{p}}_j^{\prime }.$$

$$⟨{\hat{P}}_{\mathrm{avg}}⟩=\frac{1}{M}\sum _{j=1}^M⟨{\hat{p}}_j⟩\approx \frac{\sqrt{\eta }}{M}\sum _{j=1}^M⟨{\hat{x}}_j^{\prime }⟩{\varphi }_j=\sqrt{\frac{2\eta }{M}}{\alpha }_e{\varphi }_{\mathrm{avg}}.$$

$$⟨\Delta {\hat{P}}_{\mathrm{avg}}^2⟩$$

$$\approx \frac{\eta }{M^2}⟨\Delta {\left(\sum _{j=1}^M{\hat{x}}_j^{\prime }{\varphi }_j\right)}^2⟩+\frac{\eta }{M^2}⟨\Delta {\left(\sum _{j=1}^M{\hat{p}}_j^{\prime }\right)}^2⟩+\frac{1-\eta }{M^2}⟨\Delta {\left(\sum _{j=1}^M{\hat{p}}_{\mathrm{vac},j}\right)}^2⟩$$

$$\approx \frac{\eta }{M^2}⟨\Delta {\left(\sum _{j=1}^M{\hat{p}}_j^{\prime }\right)}^2⟩+\frac{1-\eta }{M^2}\sum _{j=1}^M⟨{\hat{p}}_{\mathrm{vac},j}^2⟩$$

$$=\frac{\eta }{M^2}M⟨\Delta {\hat{p}}_e^2⟩+\frac{1-\eta }{M^2}\sum _{j=1}^M⟨{\hat{p}}_{\mathrm{vac},j}^2⟩$$

$$=\frac{\eta e^{-2r_e}+1-\eta }{2M}.$$

$${\tilde{\varphi }}^2⟨\Delta {\left(\sum _{j=1}^M{\hat{x}}_j^{\prime }\right)}^2⟩$$

$$\Rightarrow {\tilde{\varphi }}^2$$

$$\Rightarrow \tilde{\varphi }$$

$${\sigma }_e=\frac{\sqrt{⟨\Delta {\hat{P}}_{\mathrm{avg}}^2⟩}}{\left|\partial ⟨{\hat{P}}_{\mathrm{avg}}⟩/\partial {\varphi }_{\mathrm{avg}}\right|}=\frac{\sqrt{e^{-2r_e}+1/\eta -1}}{2{\alpha }_e}.$$

$$N_e=N_{e,\mathrm{coh}}+N_{e,\mathrm{sqz}}=\frac{\eta }{M}({\alpha }_e^2+{\sinh }^2{r}_e).$$

$${\mathcal{L}}_e({\alpha }_e,r_e,\lambda )$$

$$=\frac{\sqrt{e^{-2r_e}+1/\eta -1}}{2{\alpha }_e}+\lambda \frac{\eta }{M}({\alpha }_e^2+{\sinh }^2{r}_e)-\lambda N,$$

$$0={\nabla }_{{\alpha }_e}{\mathcal{L}}_e$$

$$0={\nabla }_{r_e}{\mathcal{L}}_e$$

$$0={\nabla }_{\lambda }{\mathcal{L}}_e$$

$$e^{2r_e}=\frac{{\Lambda }_M-\eta }{1-\eta },$$

$${\alpha }_e^2=N_{\mathrm{tot}}-{\sinh }^2{r}_e=N_{\mathrm{tot}}-\frac{e^{2r_e}+e^{-2r_e}-2}{4}=\frac{MN}{\eta }-\frac{({\Lambda }_M-1{)}^2}{4(1-\eta )({\Lambda }_M-\eta )},$$

$${\mu }_e=\frac{N_{e,\mathrm{sqz}}}{N}=\frac{{\sinh }^2{r}_e}{N_{\mathrm{tot}}}=\frac{\eta ({\Lambda }_M-1{)}^2}{4MN(1-\eta )({\Lambda }_M-\eta )},$$

Full text

Distributed quantum sensing in a continuous variable entangled network

Xueshi Guo

xguo@fysik.dtu.dk

  

Casper R. Breum

Center for Macroscopic Quantum States (bigQ), Department of Physics, Technical University of Denmark, Fysikvej, 2800 Kgs. Lyngby, Denmark

  

Johannes Borregaard

QMATH, Department of Mathematical Sciences, University of Copenhagen, Universitetsparken 5, 2100 Copenhagen, Denmark

  

Shuro Izumi

  

Mikkel V. Larsen

Center for Macroscopic Quantum States (bigQ), Department of Physics, Technical University of Denmark, Fysikvej, 2800 Kgs. Lyngby, Denmark

  

Tobias Gehring

Center for Macroscopic Quantum States (bigQ), Department of Physics, Technical University of Denmark, Fysikvej, 2800 Kgs. Lyngby, Denmark

  

Matthias Christandl

QMATH, Department of Mathematical Sciences, University of Copenhagen, Universitetsparken 5, 2100 Copenhagen, Denmark

  

Jonas S. Neergaard-Nielsen

jsne@fysik.dtu.dk

  

Ulrik L. Andersen

ulrik.andersen@fysik.dtu.dk

Center for Macroscopic Quantum States (bigQ), Department of Physics, Technical University of Denmark, Fysikvej, 2800 Kgs. Lyngby, Denmark

Networking plays a ubiquitous role in quantum technology Wehner et al. (2018). It is an integral part of quantum communication and has significant potential for upscaling quantum computer technologies Nickerson et al. (2013). Recently, it was realized that sensing of multiple spatially distributed parameters may also benefit from an entangled quantum network Kómár et al. (2014); Eldredge et al. (2018); Proctor et al. (2018); Ge et al. (2018); Zhuang et al. (2018); Humphreys et al. (2013); Gagatsos et al. (2016); Polino et al. (2019). Here we experimentally demonstrate how sensing of an averaged phase shift among four distributed nodes benefits from an entangled quantum network. Using a four-mode entangled continuous variable (CV) state, we demonstrate deterministic quantum phase sensing with a precision beyond what is attainable with separable probes. The techniques behind this result can have direct applications in a number of primitives ranging from molecular tracking to quantum networks of atomic clocks.

Quantum noise associated with quantum states of light and matter ultimately limits the precision by which measurements can be carried out Giovannetti et al. (2006); Escher et al. (2011); Giovannetti et al. (2011). However, by carefully designing the coherence of this quantum noise to exhibit properties such as entanglement and squeezing, it is possible to measure various physical parameters with significantly improved sensitivity compared to classical sensing schemes Caves (1981). Numerous realizations of quantum sensing utilizing non-classical states of light Yonezawa et al. (2012); Berni et al. (2015); Slussarenko et al. (2017) and matter Muessel et al. (2014) have been reported, while only a few applications have been explored. Examples are quantum-enhanced gravitational wave detection The LIGO Scientific Collaboration (2011), detection of magnetic fields Wolfgramm et al. (2010); Li et al. (2018); Jones et al. (2009) and sensing of the viscous-elasticity parameter of yeast cells Taylor et al. (2013). All these implementations are, however, restricted to the sensing of a single parameter at a single location.

Spatially distributed sensing of parameters at multiple locations in a network is relevant for applications from local beam tracking Qi et al. (2018) to global scale clock synchronization Kómár et al. (2014). The development of quantum networks enables new strategies for enhanced performance in such scenarios. Theoretical works Humphreys et al. (2013); Knott et al. (2016); Baumgratz and Datta (2016); Pezzè et al. (2017); Eldredge et al. (2018); Proctor et al. (2018); Ge et al. (2018); Zhuang et al. (2018) have shown that entanglement can improve sensing capabilities in a network using either twin-photons or Greenberger-Horne-Zeilinger (GHZ) states combined with photon number resolving detectors Proctor et al. (2018); Ge et al. (2018) or using CV entanglement for the detection of distributed phase space displacements Zhuang et al. (2018). In this Letter, we experimentally demonstrate an entangled CV network for sensing the average of multiple phase shifts inspired by the theoretical proposal of Ref. Zhuang et al. (2018). We focus on the task of estimating small variations around a known phase in contrast to ab initio phase estimation. For the first time in any system, we demonstrate deterministic distributed sensing in a network of four nodes with a sensitivity beyond that achievable with a separable approach using similar quantum states.

We start by introducing a theoretical analysis of the networked sensing scheme assuming the existence of an external phase reference. Consider a network of $M$ nodes with optical inputs that undergo individual phase shifts, $\phi_{j}\ (j=1,\dots,M)$. The goal is to estimate the averaged phase shift, $\phi_{\mathrm{avg}}=\sum_{j=1}^{M}\phi_{j}/M$, among all nodes with as high precision as possible. Two different sensing setups are considered: A separable system where the nodes are interrogated with independent quantum states (Figure 1a) and an entangled system where they are interrogated with a joint quantum state (Figure 1b). We assume the squeezers give out pure single-mode Gaussian quantum states described by the state vectors $\hat{D}(\alpha)\hat{S}(r)|0\rangle$, where $\hat{D}$ and $\hat{S}$ are the displacement and squeezing operators, respectively, $\alpha$ is the displacement amplitude and $r$ is the squeezing factor. We assume that each probe state undergoes loss in a channel with transmission $\eta$. We furthermore restrict the estimator to be the joint phase quadrature, $\hat{P}_{\mathrm{avg}}=\sum_{j=1}^{M}\hat{p}_{j}/M$ (where $\hat{p}_{j}$ are the phase quadratures of the individual modes), practically corresponding to the averaged outcome of $M$ individual homodyne detectors. These states and detectors are of particular interest due to their experimental feasibility, inherent deterministic nature, high efficiency, and robustness to noise.

Using the separable approach, $M$ identical Gaussian probe states are prepared and individually detected, while in the entangled approach, a single squeezed Gaussian state is distributed evenly to the $M$ nodes via a beam splitter array and likewise measured individually with homodyne detectors at the nodes. If one wanted to estimate different linear combinations of the phase shifts than the simple average, other beam splitter divisions would be required Eldredge et al. (2018); Proctor et al. (2018). The sensitivity of the measurement can be defined as the standard deviation of the measurement which, by error propagation, is Giovannetti et al. (2011)

[TABLE]

where $\langle\Delta\hat{P}^{2}_{\mathrm{avg}}\rangle=\langle\hat{P}^{2}_{\mathrm{avg}}\rangle-\langle\hat{P}_{\mathrm{avg}}\rangle^{2}$ is the variance of the estimator. We are only interested in the sensitivity for small phase shifts, since one can always use an initial rough phase estimation to adjust the homodyne detector (the local oscillator phase) to the maximum sensitivity setting Berni et al. (2015). For small phase shifts, we obtain the sensitivities for the separable ($\sigma_{s}$) and entangled ($\sigma_{e}$) approaches (see Supplementary Material Sec. I):

[TABLE]

We now constrain the average number of photons, $N$, hitting each sample. The photons can be separated into those originating from coherent displacement and those originating from squeezing: $N=N_{s,\mathrm{coh}}+N_{s,\mathrm{sqz}}=\eta(\alpha_{s}^{2}+\sinh^{2}r_{s})$ for the separable case and $N=N_{e,\mathrm{coh}}+N_{e,\mathrm{sqz}}=\eta(\alpha_{e}^{2}+\sinh^{2}r_{e})/M$ for the entangled case. The ratio between photon numbers, parametrized as $\mu_{s(e)}=N_{s(e),\mathrm{sqz}}/N$ can be tuned to give the optimal sensitivities

[TABLE]

For perfect efficiency ($\eta=1$), it is clear that the sensitivity of the entangled system yields Heisenberg scaling both in the number of nodes $(1/M)$ and the number of photons per mode ($1/N$) whereas the separable system only achieves the latter and a classical $1/\sqrt{M}$-scaling with the number of modes. The gain in sensitivity of the entangled network relative to the separable network (denoted $G=\sigma^{\mathrm{opt}}_{s}/\sigma^{\mathrm{opt}}_{e}$) is thus $G=\sqrt{M}$.

For non-ideal efficiency, the Heisenberg scaling ceases to exist in accordance with previous work on single parameter estimation Knysh et al. (2011). In fact, for $\eta\rightarrow 0$, both sensitivities approach $1/2\sqrt{MN}$. Still, it is important to note that the entangled network exhibits superior behavior for any value of $\eta$, $M$ and $N$ for optimized $\mu_{s},\mu_{e}$. Some examples for the sensitivity gain are illustrated in Figures 1c and d. From Fig 1d where a network of $M=4$ nodes is considered, it is clear that the highest gain in sensitivity is attained at a finite photon number. We also note that for large photon numbers, the gain tends to unity for non-zero loss meaning no enhanced sensitivity when using the entangled approach. However, there is still a practical advantage for the entangled approach: Only one squeezed state is needed compared to the M squeezed states with similar squeezing levels for the separable approach (see Supplementary Material Sec. I).

Next, we demonstrate experimentally the superiority of using an entangled network for distributed sensing. The entangled network is realized by dividing equally a displaced single mode squeezed state into four spatial modes by means of three balanced beam splitters (Fig. 2a, see Supplementary Material Sec. III for more details). These modes are then sent to the four nodes of the network at which they each undergo a phase shift $\phi_{j}$ and are finally measured with high-efficiency homodyne detectors (HD) that are set to measure the phase quadrature, $\hat{p}_{j}$. The external phase reference is set by the local oscillator which co-propagates with the probes through the setup but in a different polarization mode. This ensures that the relative phases between the probes and the local oscillator can be controlled. The resulting photo-currents from the four detectors are further processed and subsequently combined to produce the averaged phase shift. For demonstration purpose, we set all $\phi_{j}$ to the same value $\phi_{j}=\phi_{avg}$, but in principle they could be different.

We choose to define our quantum states within a narrow spectral mode at the 3 MHz sideband frequency. There are no fundamental restrictions in the scheme on the optical modes employed. In any practical setting, they would be chosen based on the nature of both the squeezing source and the samples being probed. Here, the 3 MHz sideband is chosen to maximize the squeezing from our source, an optical parametric oscillator (OPO) operating below threshold: At higher frequencies, the squeezing reduces due to the limited bandwidth of the OPO, while at lower frequencies, it is degraded by technical noise. A displaced squeezed state is obtained by injecting into the OPO a coherent state produced by phase modulating the injected beam at 3 MHz. The maximum squeezing measured through the joint measurement of 4 HDs is $\sim$5 dB at 3 MHz. More details on the probe generation are in Supplementary Material Sec. III.

An experimental run is shown in Fig. 2b. In this particular run, a displaced squeezed state with an average photon number of $N=2.48\pm 0.12$ in each mode is prepared of which $N_{e,\mathrm{sqz}}=0.30\pm 0.01$ photons are from the squeezing operation and $N_{e,\mathrm{coh}}=2.19\pm 0.11$ are from the phase modulation as this distribution is near-optimal for the entangled case. We then impose 12 different $\phi_{\mathrm{avg}}$ values by phase shifts at each node while recording the Fourier transformed homodyne detector outputs; the spectra around the 3 MHz sideband for six of the $\phi_{\mathrm{avg}}$ values are shown in Fig. 2b (see Supplementary Material Sec. V for more details). These outputs yield poor estimates of the individual phase shifts (because the squeezing in each mode is only $\sim$0.8 dB) but the averaged phase shift obtained by summing the photo-currents produces an entanglement-enhanced estimate with significantly lower noise. The spectra for the averaged photo-currents are shown in Figure 2c. For comparison, we also simulate the separable approach by directing the entire displaced squeezed state (with properly optimized parameters) to a single node. We then perform the phase estimation at that node and scale the obtained sensitivity by $\sqrt{4}$ to get the projected performance for an average over four identical sites. An example is shown in Fig 2b for $N=2.63\pm 0.11$, with $N_{s,\mathrm{sqz}}=0.31\pm 0.01$ and $N_{s,\mathrm{coh}}=2.32\pm 0.10$.

We quantify the performance of the sensing network by estimating the sensitivities of the two approaches based on the averaged homodyne measurement outcomes, $P_{\mathrm{avg}}$. By extracting the rate of change with respect to a phase rotation, $|\partial\langle\hat{P}_{\mathrm{avg}}\rangle/\partial\phi_{\mathrm{avg}}|$, as well as the variance, $\langle\Delta\hat{P}_{\mathrm{avg}}^{2}\rangle$, of $P_{\mathrm{avg}}$, we deduce the sensitivity using Eq. (1). For the experimental runs described above, we obtain sensitivities of $\sigma_{e}=0.099\pm 0.003$ and $\sigma_{s}=0.118\pm 0.002$ for the entangled and separable approach, respectively. This corresponds to single shot resolvable distributed phase shifts (that is, phase shifts for which the signal-to-noise ratio is unity) of $5.66^{\circ}\pm 0.18^{\circ}$ for the entangled case and $6.76^{\circ}\pm 0.11^{\circ}$ for the separable case with $\sim 2.5$ photons. Using a coherent state in replacement of the squeezed state, the minimal resolvable phase for $2.5$ photons is $9.06^{\circ}\pm 0.07^{\circ}$ corresponding to the standard quantum limit. Note that these angles are larger than our small phase shift approximation (which requires $\phi_{avg}$ to be much smaller than $\sim 7^{\circ}$ for the conditions in this experimental run, see Supplementary Sec. I). In practice this means that it is necessary to probe the sample more than once to resolve the small phases implemented in the experiment. Sampling the phases $K$ times will result in $\sqrt{K}$ times smaller resolvable phase shift angles. The entangled strategy will still benefit from the enhanced sensitivity per probe.

We find the sensitivities for different total average photon numbers both for the entangled and separable network, and plot the results in Figure 3a. For every selection of the total photon number, we adjust $\mu$ to a near-optimal value for optimized sensitivity (Figure 3b,c). It is clear in Figure 3a that both realizations beat the standard quantum limit (reachable by coherent states of light), and most importantly, we see that the entangled network outperforms the separable network. The ultimate sensitivity of our entangled approach is not reached in our implementation. However, homodyne detection will not even in principle saturate this bound and non-Gaussian measurements are in general needed (see Supplementary Material Sec. II).

Our results experimentally demonstrate how mode entanglement, here in the form of squeezing of a collective quadrature of a multi-mode light field, can enhance the sensitivity in a distributed sensing scenario. Importantly, we show this enhancement in an experimentally feasible setting where the sensitivity of standard coherent probes are enhanced through quadrature squeezing. This approach allows for easily tunable probe powers in order to adapt the setup to the specific application. Furthermore, because the entanglement is generated from a simple beam-splitter network, it is straight-forward to scale to more modes where the sensitivity gain may be even larger, cf. Fig. 1c. The main limitation will be the channel efficiency which will eventually limit the gain. Consequently, we believe that techniques demonstrated here have direct applications in a number of areas. Specifically, beam tracking relevant for molecular tracking Taylor et al. (2013); Qi et al. (2018) could directly benefit from these techniques. Such applications impose limits on the allowed probe power to prevent photon damage and heating of the systems. Mode-entanglement can thus be used to increase sensitivity without increasing the probe power. Using squeezed coherent light for quantum non-demolition (QND) measurement has also been exploited for generation of spin squeezing in atomic ensembles Hammerer et al. (2010) and optical magnetometry Wolfgramm et al. (2010). While this is usually considered for single ensembles, the generalization to multiple ensembles can provide enhanced sensitivity and new primitives for quantum information processing. Combining several ensembles for magnetometry and utilizing mode-entanglement would further reduce the shot-noise and increase sensitivity of a collective optical measurement. Performing a collective optical QND measurement of several atomic ensembles can prepare a distributed spin-squeezed state for quantum network applications. In particular, squeezing of multiple optical lattice clocks could be used for collective enhancement of clock stability Kómár et al. (2014); Polzik and Ye (2016). In Ref. Polzik and Ye (2016), this was obtained by letting a single probe interact with all ensembles in a sequential manner. However, utilizing mode-entanglement, this can be performed in a parallel fashion with no quantum signal being transmitted between the ensembles.

**DATA AVAILABILITY

**Experimental data and analysis code is available on request.

**Acknowledgments

**M.C. and J.B. acknowledge support from VILLUM FONDEN via the QMATH Centre of Excellence (Grant no. 10059), the European Research Council (ERC Grant Agreements no 337603), and from the QuantERA ERA-NET Cofund in Quantum Technologies implemented within the European Union’s Horizon 2020 Programme (QuantAlgo project) via the Innovation Fund Denmark. X.G., C.B., S.I., M.L., T.G., J.N. and U.A. acknowledge support from Center for Macroscopic Quantum States (bigQ DNRF142). X.G., S.I. and J.N. acknowledge support from VILLUM FONDEN via the Young Investigator Programme (Grant no. 10119).

**AUTHOR CONTRIBUTIONS

**J.B., U.A., J.N., T.G. and X.G. conceived the experiment. X.G., C.B. and M.L. performed the experiment and analyzed the data. J.B., X.G., S.I., M.C. and J.N. worked on the theoretical analysis. X.G. wrote the paper with contributions from J.B., C.B., S.I., J.N. and U.A. J.N. and U.A. supervised the project.

**ADDITIONAL INFORMATION

Competing interests:** The authors declare that there are no competing interests.

I Averaged phase shift sensing with

$\hat{P}_{\mathrm{avg}}$ estimator

Our distributed phase sensing scenario is as follows. At each of $M$ spatially separated locations, an optical phase shift $\phi_{j}$ occurs. We are interested in estimating the average phase shift $\phi_{\mathrm{avg}}=\frac{1}{M}\sum_{j=1}^{M}\phi_{j}$. It is straight-forward to generalize to other linear combinations of the phase shifts, but for the sake of demonstrating the power of the entangled approach it suffices to consider the simple average, where the gain is maximum Proctor et al. (2018). We consider two different approaches: The separable scheme where each phase shift is probed individually by squeezed coherent states, and the entangled scheme where the $M$ locations are part of an optical network endowed with a single squeezed coherent state that is distributed among the nodes to serve as an entangled probe. In either case, the phase shifted probes are measured by homodyne detection of their phase quadratures and the results are communicated classically to establish the average.

We furthermore make the following assumptions to simplify the analysis:

All the phase shifts are small, giving the small-angle approximation $\sin\phi\approx\phi$. 2. 2.

All probes in the separable approach are identical, having real-valued displacement amplitude $\alpha_{s}$ and squeezing in the phase quadrature with squeezing parameter $r_{s}$. That is, the $M$ probes are each in the state $|\psi^{(s)}\rangle=\hat{D}(\alpha_{s})\hat{S}(r_{s})|0\rangle$, where $\hat{D}(\alpha)=\exp(\alpha\hat{a}^{\dagger}-\alpha^{\ast}\hat{a})$ is the displacement operator and $\hat{S}(r)=\exp(\frac{r}{2}(\hat{a}^{{\dagger}2}-\hat{a}^{2})$) is the squeezing operator. 3. 3.

In the entangled approach, the single initial resource state has real-valued displacement amplitude $\alpha_{e}$ and phase squeezing with squeezing parameter $r_{e}$, that is, it is in the state $|\psi^{(e)}\rangle=\hat{D}(\alpha_{e})\hat{S}(r_{e})|0\rangle$. This resource is divided evenly through the network to the $M$ nodes. 4. 4.

The channel losses, quantified by the efficiency parameter $\eta$, are identical for the $M$ channels and they occur entirely prior to the probes reaching the phase samples. In other words, we assume the phase samples themselves and the detection to be lossless. While this assumption is not quite realistic, even in our experiment, it mostly has consequences when keeping track of the number of photons hitting the sample but does not influence the sensitivity as such. In a truly distributed setting, most losses would also happen in the distribution of the resources.

For high-sensitivity estimation of larger phase shifts, these assumptions can still be fulfilled, as long as the local oscillator in the homodyne detector is pre-adjusted to be roughly 90∘ out of phase with the shifted probe. This rough estimation can be done with just a few initial probings Berni et al. (2015).

I.1 General sensitivity for small phase shift

I.1.1 Separable scheme

With probe states given as described above, we use the notation defined in Figure S1 to analyse the separable scheme. The phase quadrature of a single mode after channel loss and the phase shift $\phi_{j}$ is

[TABLE]

where $\hat{x}_{s,j},\hat{p}_{s,j}$ are the quadrature operators of the initial squeezed states with mean values $\langle\hat{x}_{s,j}\rangle=\sqrt{2}\alpha_{s}$, $\langle\hat{p}_{s,j}\rangle=0$ and variances $\langle\Delta\hat{x}_{s,j}^{2}\rangle=\frac{1}{2}e^{2r_{s}}$, $\langle\Delta\hat{p}_{s,j}^{2}\rangle=\frac{1}{2}e^{-2r_{s}}$, while $\hat{x}_{\mathrm{vac},j}$, $\hat{p}_{\mathrm{vac},j}$ are vacuum mode operators admixed through the losses. The expectation value of the rotated phase quadrature is

[TABLE]

The phase shift can thus be directly estimated from the measured $\hat{p}_{j}$ values. The average phase shift of $M$ modes, $\phi_{\mathrm{avg}}=\frac{1}{M}\sum_{j=1}^{M}\phi_{j}$, can then be estimated with the estimator $\hat{P}_{\mathrm{avg}}=\frac{1}{M}\sum_{j=1}^{M}\hat{p}_{j}$:

[TABLE]

The sensitivity of the estimation is defined as the standard deviation which, from standard error propagation analysis, is given by

[TABLE]

The slope of $\hat{P}_{\mathrm{avg}}$ versus $\phi_{\mathrm{avg}}$ is

[TABLE]

and its variance is

[TABLE]

The second equality comes from the fact that in the separable approach there are no correlations between the modes. Under a stronger bound on the magnitude of the phase shifts, $\phi_{j}\ll\sqrt{\langle\Delta\hat{p}^{2}_{s,j}\rangle/\langle\Delta\hat{x}^{2}_{s,j}\rangle}$, this expression reduces to

[TABLE]

Hence, the sensitivity is

[TABLE]

The average number of photons hitting each sample is

[TABLE]

I.1.2 Entangled scheme

With entangled probes (the notation used in our analysis is summarized in Figure S2), we use the same estimator, $\hat{P}_{\mathrm{avg}}$. The individual modes that combine to form the average are, however, now related through the distributed single initial resource $\hat{p}_{e}$:

[TABLE]

where the primed mode operators are obtained after symmetric distribution in the beam-splitter network, that is,

[TABLE]

The mean value of the estimator is

[TABLE]

The variance is

[TABLE]

In the second line, we made use of the fact that there are no correlations between $\hat{x}$ and $\hat{p}$ quadratures for the given probe state in our entangled scheme as well as the small angle approximation $\cos(\phi_{j})\approx 1$, $\sin(\phi_{j})\approx\phi_{j}$. In the third line, we further tightened the small angle approximation by taking a $\tilde{\phi}$ such that for all $j$, $|\phi_{j}|<\tilde{\phi}$ and assuming

[TABLE]

This approximation gives a sensitivity for the entangled approach of

[TABLE]

Note that this, in contrast with the separable approach, does not depend on the number of modes $M$. The sensitivity is therefore the same as the sensitivity for a single mode with the same resource state - but in the single mode case the sample would of course be exposed to $M$ times as many photons. The average number of photons hitting each sample in the distributed, entangled scheme is

[TABLE]

I.2 Optimized parameters and sensitivities

I.2.1 Entangled scheme

With the sensitivities given by eqs. (S10) and (S17), we wish to find the values for the displacement amplitudes and squeezing parameters that optimize the sensitivity for a fixed photon number on the sample. This problem can be solved with Lagrangian multipliers, using the constraint $N_{s,e}-N=0$, where $N$ is the photon number to be held fixed during optimization. The total photon number of the resource state(s) before loss is then $N_{\mathrm{tot}}=MN/\eta$.

The Lagrange function for the entangled scheme is

[TABLE]

and the equations for the stationary point of the Lagrangian become

[TABLE]

After some manipulation, the solutions can be expressed as

[TABLE]

with $\Lambda_{M}=\sqrt{1+4MN(1-\eta)}$. The optimal photon number ratio is

[TABLE]

and the optimal sensitivity obtained with these parameters becomes

[TABLE]

which for $\eta=1$ reduces to $\sigma_{e}^{\mathrm{opt}}(\eta=1)=\frac{1}{2MN}\sqrt{\frac{MN}{MN+1}}$. This sensitivity exhibits Heisenberg scaling in both photon number (due to the squeezing) and mode number (due to the entanglement).

I.2.2 Separable scheme

Doing the same derivation for the separable scheme, that is, starting from the Lagrange function $\mathcal{L}_{s}(\alpha_{s},r_{s},\tilde{\lambda})=\sigma_{s}+\tilde{\lambda}(N_{s}-N)$, results in the following optimal parameters for squeezing and displacement:

[TABLE]

with $\Lambda_{1}=\sqrt{1+4N(1-\eta)}$, and a corresponding photon number ratio

[TABLE]

Finally, the optimal sensitivity becomes

[TABLE]

which for $\eta=1$ reduces to $\sigma_{s}^{\mathrm{opt}}(\eta=1)=\frac{1}{2\sqrt{M}N}\sqrt{\frac{N}{N+1}}$, thus no longer showing Heisenberg scaling in the mode number. The result enable us to obtain the simulation result in Fig 1c and 1d in the main text, and the optimal $\mu$ and corresponding squeezing rated need to get the optimal $\mu$ is shown in Fig. S3.

II Quantum Cramér-Rao Bound for distributed phase sensing

In this section, we discuss the Quantum Cramér-Rao Bound (QCRB) for $\phi_{avg}$ sensing, which defines the ultimate sensitivity limit for a given probe with multi-mode Gaussian state. We also compare the QCRB for $\phi_{avg}$ sensing to the counterpart of single parameter phase sensing by using a single-mode Gaussian state as probe, which was discussed in detail in Ref. Pinel et al. (2013). Throughout the analysis, we assume the sensing channel has a constant efficiency $\eta$.

For a general sensing problem, the Quantum Cramér-Rao Bound sets a lower bound, minimized over all possible measurements, on the uncertainty with which a parameter $\phi$ can be estimated through an unbiased estimator $\tilde{\phi}$, given a probe in a certain quantum state: $\langle\Delta\tilde{\phi}^{2}\rangle\geq F_{\phi}^{-1}$, where $F_{\phi}$ is the quantum Fisher information. The ultimate sensitivity limit for sensing of a single phase shift is thus $\sigma_{\phi}^{CR}=1/\sqrt{F_{\phi}}$. The quantum Fisher information for single mode phase sensing using a Gaussian probe with initial displacement $\alpha$ and squeezing $r$ is given by Pinel et al. (2013)

[TABLE]

with an effective thermalization photon number due to loss

[TABLE]

and and effective squeezing parameter

[TABLE]

The QCRB for sensing of an average $\phi_{\mathrm{avg}}$ of multiple phase shifts with coherent probes, $\sigma^{CR}_{\mathrm{coh}}$, and with our separable approach, $\sigma^{CR}_{\mathrm{sep}}$, can also be found from Eq. (S32) divided by $\sqrt{M}$ to account for the $M$ independent phase estimations.

For non-trivial estimation involving multiple parameters, the quantum Fisher information matrix (QFIM) is needed. The variance of an unbiased estimator $\tilde{q}$ of an arbitrary linear combination of $M$ parameters, $q=\sum_{i=1}^{M}w_{i}\phi_{i}$, is $\langle\Delta\tilde{q}^{2}\rangle=\mathbf{w}^{T}\mathbf{K}\mathbf{w}$ with the weight coefficients $\mathbf{w}^{T}=(w_{1},\ldots,w_{M})$ and parameter covariance matrix $\mathbf{K}$ with $K_{ij}=\langle(\tilde{\phi}_{i}-\phi_{i})(\tilde{\phi}_{k}-\phi_{j})\rangle$. Given a quantum Fisher information matrix $\mathbf{F}$, the QCRB is expressed as

[TABLE]

The QCRB for $\phi_{\mathrm{avg}}$ with our entangled scheme where the weights are $\mathbf{w}_{\mathrm{avg}}^{T}=(1/4,1/4,1/4,1/4)$ is then

[TABLE]

How to calculate the QFIM for arbitrary multi-mode Gaussian states, as well as the existence (or not) of a measurement that reaches the bound, is discussed in Šafránek (2018); Nichols et al. (2018). We use Eqs. (16)-(21) in Ref. Nichols et al. (2018) to numerically calculate $\mathbf{F}_{\mathrm{ent}}$. As discussed in Ref. Šafránek (2018), when any of the symplectic eigenvalues of the quadrature covariance matrix of the Gaussian state has unity value, the process in Nichols et al. (2018) gives a singular result. This applies to our entangled scheme since it has 3 vacuum input modes. We solve this numerical problem pragmatically by replacing the 3 vacuum states with very weak thermal states (10*-6* mean photon number in each).

The optimal QCRBs for different scenarios together with the optimal sensitivity of our measurement schemes $\sigma^{opt}_{s}$ and $\sigma^{opt}_{e}$ are shown in Fig. S4. The QCRBs are optimized over $\alpha$ and $r$ for a fixed mean photon number. It is interesting to note, that for the $\eta=0.735$ detection efficiency, the states optimizing the QCRBs are all squeezed vacuum states. From the result in Fig. S4, we find at our detection efficiency:

[TABLE]

The sensitivities we obtain do not reach the corresponding ultimate limits. Furthermore, these relations show that in principle it should be possible to reach a better sensitivity with a separable scheme than what we obtain with the entangled scheme. However, the difference between them is small and—to the best of our knowledge—no efficient way of experimentally implementing a measurement to reach the ultimate limit $\sigma^{CR}_{sep}$ is known. In Oh et al. (2019), it is discussed that an $\hat{X}\hat{P}+\hat{P}\hat{X}$ type of measurement is needed. This is non-Gaussian and can not be realized by only Gaussian operations such as squeezing, beam splitters, phase shifts and homodyne/heterodyne detection. For sensing involving multiple parameters, it is unclear how to reach the optimal bound $\sigma^{CR}_{\mathrm{ent}}$—or whether it is at all possible. A joint, non-Gaussian measurement may be optimal.

Finally, note that the QCRB of the entangled scheme overlaps with that of the single parameter estimation. This makes sense intuitively: Splitting a resource state equally in four and using these to probe the average of four phase shifts should result in the same sensitivity as that of a single phase shift probed by the same un-split resource.

III Preparation of entangled probes

The entangled probes are prepared in two steps. First, we generate a squeezed coherent state, denoted as the squeezed probe (SP), by an optical parametric oscillator (OPO). Second, we send the SP through a beam-splitter network (BSN) to generate 4 entangled probes. We define the mode of the SP to be a narrow sideband at 3 MHz, since this is region where we have high squeezing quality.

III.1 Generation of squeezed probes with OPO

The laser source for the experiment is an amplified NKT Photonics X-15 fibre laser operating at 1550 nm. Most of the light is used for pumping a second harmonic generation (SHG) cavity (same design as the OPO described below) to produce 775 nm light to act as the OPO pump. The rest is used for the local oscillator and the probe and lock beams. As shown in Fig. S5, we use a bow-tie shaped OPO with a periodically poled potassium titanyl phosphate (PPKTP) crystal to generate the SP by type-0 parametric down conversion. The bandwidth of the cavity is 8.0 MHz half width half maximum (HWHM) and the OPO pump power threshold is 850 mW. The 775 nm pump, which for the measurements presented here varied between 150 mW and 350 mW, is coupled through the dichroic curved cavity mirrors and dumped after passing the crystal. A 3.6 mW coherent beam at 1550 nm, weakly phase-modulated by an electro-optic modulator (EOM) at 3 MHz and 28.7 MHz, is coupled into the OPO in the counter-propagating direction through a high reflectivity mirror (HR) with a transmittance of about 100 ppm. This beam is used to lock the cavity by the Pound–Drever–Hall technique with the 28.7 MHz side band and the resonant detector D1. All cavity and phase locks in the experiment are handled by Red Pitaya FPGA boards running the PyRPL lockbox software Neuhaus et al. (2017).

The reflection from the HR mirror is re-coupled into the forward-propagating mode of the OPO with a 0∘ mirror to serve as the carrier of the sideband mode that defines our probe state. A variable attenuator (Att.) is inserted to control the optical power. In the OPO, the forward-propagating beam is squeezed by the parametric process and coupled out through a 10$\%$ transmittive out coupling mirror (light gray in Fig. S5). A half-wave plate ($\lambda/2$) and a polarization beam-splitter (PBS) is used to tap around $1\%$ of the OPO output towards detector D2 to lock the phase between the carrier and the pump for de-amplification. As a result, the carrier is squeezed in the amplitude quadrature, leading to squeezing of the phase quadrature of the probe in the 3 MHz modulated sideband frequency mode since the sideband is encoded by phase modulation.

III.2 Generation and detection of entangled probes

The detailed experimental setup is shown in Fig. S6. It is essentially a multi-port version of a squeezed-light-enhanced polarization interferometer Grangier et al. (1987). We create four entangled probes by sending the squeezed probe, SP, through a BSN consisting of three 50:50 beam-splitters. Prior to this, the SP is spatially combined on a PBS with a strong beam (LO) which will act as the local oscillator for all four modes. The LO phase is locked to either the $\hat{p}$ or $\hat{x}$ quadrature of the SP by tapping $\sim 1\%$ towards a polarization-based homodyne detection setup, the output of which is used to control a piezo-mounted mirror in the LO path. In each of the four modes, the phase between LO and SP can be further controlled by a $\lambda/4$ and a $\lambda/2$ wave plate. The $\lambda/4$ plates change the LO and SP into left-hand and right-hand circular polarization, respectively. The $\lambda/2$ plates introduce phase shifts between SP and LO and play two roles: First, they are used to synchronize the phases for the entangled probes by compensating the phase difference induced by 50:50 beam splitters. Second, they are used to simulate the phase samples, that is, the imposed phases $\phi_{1},\dots,\phi_{4}$. For details, see section IV.1. Finally, the four outputs are measured on homodyne detectors.

III.2.1 Homodyne detection and data acquisition

All five homodyne detection setups use the same scheme, illustrated in Fig. S7a. The circularly polarized SP and LO interfere after the PBS. The optical power of the LO is about 3 mW on each HD and it detects a SP of about 10 nW. The output of the detector is electronically split into AC and DC parts with a bias-tee of about 100 kHz. The AC signal, $V_{ac}$, is used for phase sensing. It includes the 3 MHz side-band, but filters out the carrier at DC and the side-band for cavity locking at 28.7 MHz with a low pass filter at around 14 MHz. The DC signal, $V_{dc}$, detects the carrier. It is used for phase locking and phase calibration (see section IV.2).

The $V_{ac}$ outputs of HD1 to HD4 are sent to a 4-channel oscilloscope (LeCroy HDO6034), which acquires time-voltage traces of 200 $\mu$s with a 50 MHz sampling rate. The power spectral densities (PSDs) of the individual HD outputs and the averaged output is obtained by Fast Fourier Transform (FFT) on a computer. Fig. S7b shows PSDs for the averaged voltage of the 4 HDs in different experimental conditions with no modulation from the EOM. All the PSDs in Fig. S7b are the averaged result of 400 oscilloscope measurements. To show the signal-to-noise ratio of the data acquisition system, we measure the PSDs of the shot noise level (SNL, measured when SP is blocked) and electronic noise (EL, measured when both SP and LO are blocked). The result is shown in green and dark grey traces in Fig. S7b. We see that the electronic noise clearance is about 23 dB at the 3 MHz side band, which corresponds to about $0.5\%$ effective loss in detection efficiency. We will discuss the other PSDs shown in Fig. S7b in the following subsection.

III.2.2 Input-output relations of the BSN

The BSN we use in the experiment is shown in Figure S6. The only non-vacuum input mode $\hat{a}_{1}$ is the SP, whose mode operator is $\hat{a}_{1}=\hat{S}^{\dagger}(r)\hat{a}\hat{S}(r)+\alpha$ in the Heisenberg picture, with $\hat{a}$ being the annihilation operator of the OPO input at 3 MHz and the real-valued $\alpha$ being the effective coherent excitation of the mode after modulation by the EOM and de-amplification in the OPO. All the other input modes $\hat{a}_{2}$, $\hat{a}_{3}$ and $\hat{a}_{4}$ are vacuum modes. The output modes of the BSN, $\hat{b}_{j}$ can be explicitly written as:

[TABLE]

Here we have introduced an identical overall efficiency $\eta$ and vacuum mode operator $\hat{a}_{vac,j}$ for $j=1\dots 4$. Although the various inefficiencies occur at different points in the experiment, for simplicity we have assumed (as in section I) that they all occur after the distribution of the probes in the BSN and that they are identical for the four channels. Experimentally, we use eight variable irises before the PDs of all four HDs to equalize the overall detection efficiency.

III.2.3 Overall detection efficiency estimation

The loss budget of our experiment setup is as follows: the escape efficiency of the OPO $\sim 95\%$; the quantum efficiency of the photo diodes in HD $\sim 98\%$; the imperfection of the mode matching between SP and LO $\sim 90\%$; the electronic noise of the homodyne detection $\sim 99\%$; the efficiency introduced by tapping for phase locking $\sim 97\%$ and the efficiency of all optics between OPO output and the PD of the HD $\sim 92\%$. The loss budget of the experiment system gives an estimation of the overall detection efficiency of $\eta\sim 74\%$.

We also estimate the overall detection efficiency by measuring the squeezing/anti-squeezing degrees (notated with $v^{2}_{sq}$ and $v^{2}_{asq}$) for the entangled approach at 3 MHz. Since

[TABLE]

we can calculate $\eta$ and $r_{e}$ with measured $v^{2}_{sq}$ and $v^{2}_{asq}$. The overall efficiency estimated with 5 different pump powers to the OPO is $\eta=73.5\%\pm 1.5\%$. This result coincide with the loss budget estimation, and we use this result to theoretically calculate the sensitivity.

For the separable approach, where the BSN is removed, the overall efficiency is $\sim 1.5\%$ higher. However, we compensate this by tapping more to the lock detector D2 in Figure S5 so the separable approach has similar efficiency to that of the entangled approach.

III.2.4 Entanglement characterization of the probes

The squeezing degree for each individual output mode will not be better than 3/4 shot noise due to the splitting of the SP in the BSN. However, the squeezing of SP is converted into entanglement between all the probes. By joint measurement of the 4 probes (simply averaging the voltage from the four HDs), we can recover the squeezing degree of the SP: From Eq. (III.2.2), the joint measurement recovering the squeezing of SP is simply the sum of the four HD${}_{\text{1--4}}$ outputs. The recovered squeezed and anti-squeezed quadratures are shown as SQ (blue) and ASQ (red) in Fig. S7b. We see the joint measurement gives about 4.8 dB of squeezing at the 3 MHz side band frequency. The additional noise seen below 2 MHz is due to technical noise from our laser. As a calibration of the noise of the probe before the parametric process, we measure the PSDs of $\hat{X}$ and $\hat{P}$ quadrature by blocking the pump of our OPO, and the result is shown with PX (light blue) and PP (light red) in Fig. S7b (noting that here we refer to the amplitude/phase quadrature of the carrier of the SP since there is no side-band). We see the technical noise of both $\hat{X}$ and $\hat{P}$ quadrature decreases as the frequency increases and overlap with the SNL when the frequency is above 1.8 MHz. Therefore, in our estimation of the overall detection efficiency at the side band frequency (3 MHz), we ignore contributions from technical noise of the laser.

With the measurement described in Fig. S7b, we can get the squeezing/anti-squeezing degree in SNL units. Fig. S8a and b shows the squeezing and anti-squeezing of an individual channel (HD1) and that from the joint measurement, respectively. The dashed lines show the squeezing and anti-squeezing predicted by Polzik et al. (1992)

[TABLE]

where $S^{-}(f)$ and $S^{+}(f)$ denotes the squeezing and anti-squeezing spectrum, $\eta=0.735$ is the estimated overall detection efficiency, $f_{cav}=8.0$ MHz is the HWHM of the OPO cavity and $P_{th}=850$ mW is the threshold of the OPO. In this measurement, $P=300$ mW pump power is used. Here both $f_{cav}$ and $P_{th}$ are obtained from independent measurements.

We quantitatively verify the entanglement of the probes by reconstructing the covariance matrix of the 4 modes. As we do not expect correlations between $\hat{x}$ and $\hat{p}$ quadratures, we only experimentally reconstruct $M_{x}=\mathrm{Cov}(\hat{x}_{j},\hat{x}_{k})$ and $M_{p}=\mathrm{Cov}(\hat{p}_{j},\hat{p}_{k})$ for $j,k=1$ to $4$ at around 3 MHz. After balancing the length of cables from HD${}_{\text{1--4}}$ to the oscilloscope, we digitally filter the recorded traces by a 50 kHz band pass filter centered around 3 MHz, and measure $M_{x}$ and $M_{p}$, respectively. The covariance matrices in shot noise units from the average of 400 oscilloscope measurements are:

[TABLE]

where symmetric elements are not shown. We show the entanglement property of the probes by calculating the logarithmic negativity $\mathcal{N}(\hat{\rho})$ between them, where $\mathcal{N}(\hat{\rho})>0$ is a sufficient condition for entanglement Vidal and Werner (2002). For a Gaussian state this can be obtained through the symplectic eigenvalues of the partially transposed covariance matrix, so that $\mathcal{N}(\rho)=\sum_{k}f(\tilde{v}_{k})$, where $\tilde{v}_{k}$ are the symplectic eigenvalues and $f(x)=-\mathrm{log}_{2}(x)$ for $x<1$ and 0 otherwise. By constructing the full $M_{x,p}$ covariance matrix from $M_{x}$ and $M_{p}$, we find that for any two, three or four modes the value of $\mathcal{N}(\hat{\rho})$ is within the range of $0.20\pm 0.02$, $0.33\pm 0.02$ and $0.51\pm 0.02$ respectively, confirming the presence of quadrature entanglement across all mode combinations.

IV Phase control and calibration

In this section we first calculate the interference at the two photo diodes of the HD in Fig. S7. The result shows that the phase between SP and LO can be controlled by rotating the $\lambda/2$ wave plates. After that, we describe the phase calibration procedure and result in our experiment. With the phase calibration result, we can control the phase $\phi_{j}$ for $j=1\dots 4$ (and therefore $\phi_{\mathrm{avg}}$) by rotating the $\lambda/2$ wave plates to a specific position.

IV.1 Phase control with $\lambda/2$ wave plates

The LO with p polarization and OPO output with s polarization are combined by the PBS in Fig. S6, and the Jones vector after the PBS is

[TABLE]

where $E_{LO}\cdot e^{-i\phi_{LO}}$ is the LO and $E_{SP}\cdot e^{-i\phi_{SP}}$ is the OPO output (squeezed probe). The Jones Matrix for a wave plate is Collett (2005)

[TABLE]

where $\theta$ is the angle between the fast axis of the wave plate and p polarization (the direction of LO), and $\phi$ is the retardance of the wave-plate ($\phi=\pi$ or $\phi=\pi/2$ for an ideal $\lambda/2$ or $\lambda/4$ wave-plate, respectively). We fix the $\lambda/4$ wave plate at $\theta=45^{\circ}$ and put the $\lambda/2$ wave plate at a variable angle $\theta_{v}$, resulting in the output Jones vector

[TABLE]

with

[TABLE]

Therefore, the interference between the two polarization modes observed at the two diodes of the HD after the second PBS is:

[TABLE]

where $\phi_{d}=\phi_{SP}-\phi_{LO}$ is the initial phase difference between OPO output and LO mode before being overlapped at the first PBS. The result show that if we rotate the $\lambda/2$ wave plate by an angle of 1∘, the phase between LO and OPO output will change 4∘. The form of Eq. (S44) also shows the visibility of the HD is not affected by the polarization transformation since it doesn’t have any constant term. However, if the wave plates or PBS are not perfect, which means that the wave plates have either more or less retardance or that the PBS has a finite extinction ratio between s and p polarization, a similar calculation shows the rotation of $\lambda/2$ wave plate by 1∘ will result in a phase shift slightly deviating from 4∘, and that the visibility of the interference at HD can be reduced. We experimentally measure these imperfections as shown in the following subsection.

IV.2 Phase calibration

During the experiment we lock $\phi_{d}$ to be either $0^{\circ}$ or $90^{\circ}$ with HDL, and use the rotation of the $\lambda/2$ wave plate before each HD to control the phase of each mode. In order to account for potential imperfections in our experiment, we first measured the visibility reduction from imperfect polarization components. We find a worst-case reduction of the HD visibility from $98.5\%$ to $95.2\%$. We also perform a phase calibration by scanning the phase between LO and SP carrier with a ramp at 27 Hz while the interference fringe measured from $V_{dc}$ of HDL and HD1,2,3,4 is recorded. The phase between LO and signal in each path is inferred from sine curve fitting. We calibrate the phase with 40 repeated measurements at each $\lambda/2$ wave plate position, and the result is shown in Fig. S9. The SQ (blue dots) shows the result when we lock $\phi_{i}=0^{\circ}$ and the HDs measure the squeezed quadrature, and the ASQ (red dots) shows the result when we lock $\phi_{i}=90^{\circ}$. For both SQ and ASQ, we rotate the $\lambda/2$ wave-plate position in each channel by an actuator in the wave-plate mount, allowing us to faithfully use the calibration result in the experiment.

From the calibration results, we see that the phase is linear within the whole range of the actuator (8∘) on the wave plate mounts. The result of the linear fitting to HD channel $j=1$ to $4$ with the equation

[TABLE]

is summarized in Table 1. With these fitted parameters, we can control both the phase in each channel $\phi_{j}$ or the averaged phase $\phi_{\mathrm{avg}}$ accurately. Particularly, if we lock $\phi_{i}$ to $0^{\circ}$, we can change the $\phi_{\mathrm{avg}}$ by a slope of $3.96^{\circ}\pm 0.02^{\circ}$; if we lock $\phi_{i}$ to $90^{\circ}$, can change the $\phi_{\mathrm{avg}}$ by a slope of $4.02^{\circ}\pm 0.02^{\circ}$.

V Data analysis

In this section we introduce the details of our data analysis procedure, which includes measuring the sensitivity by fitting and counting how many photons in average is used in the SP.

As the estimator of $\phi_{\mathrm{avg}}$, $\hat{P}_{\mathrm{avg}}$ is experimentally estimated from the PSD of the averaged output of the four HDs in each mode. Fig. S10 shows the PSD results measured for different $\phi_{\mathrm{avg}}$. Each PSD is obtained from the FFT of an average of 2000 oscilloscope traces. The spectrum peak at 3 MHz $S_{pk}$ gives the value of

[TABLE]

where $V_{sn}$ is the 4-mode shot noise limit (SNL) voltage from HDs decided by LO power, electronic gain and the digital filtering. The constant 2 in Eq. (S46) comes from the commutation relationship we choose $[\hat{X},\hat{P}]=i$. We start our data analysis by separating the peak into two voltage parts

[TABLE]

where $V_{s}=\sqrt{2}V_{sn}|\langle\hat{P}_{\mathrm{avg}}\rangle|$ is the signal part induced by the coherent photons of the side band, and $V_{n}=\sqrt{2}V_{sn}\sqrt{\langle\Delta\hat{P}_{\mathrm{avg}}^{2}\rangle}$ is the part induced by the fluctuation of the light. Except at the 3 MHz peak, the spectra in Fig. S10 vary slowly with frequency. This enables us to extract $V_{n}$ from the adjacent frequencies of the 3 MHz peak. The procedure of $V_{n}$ estimation is illustrated with the anti-squeezing quadrature (ASQ, $\phi_{\mathrm{avg}}$ = -89.5 $\pm 0.8^{\circ}$) PSD in Fig. S10 as an example. We first do a linear fit with the frequency range indicated by the red dots, which is slightly away from 3 MHz. This fitting gives the black dashed line labeled as ”Fitting for ASQ”. $V_{n}$ is then inferred by the square root value of the fitted line at 3 MHz. Since our side band line width is obviously smaller than the 5 kHz resolution of the FFT, only one peak point is observed in the PSDs in Fig. S10. Therefore, $V_{s}$ can simply be calculated by the difference between the blue dot at 3 MHz and the fitting result.

In our experiment we always introduce equal positive phase shift in all channels. In this case we know that $\langle\hat{P}_{\mathrm{avg}}\rangle>0$, and $V_{s}$ and $V_{n}$ relate to the averaged phase $\phi_{\mathrm{avg}}$ by

[TABLE]

Here $\alpha$ is the real coherent amplitude from modulation, M is the mode number, $\theta_{1}$ and $\theta_{2}$ are parameters indicating the imperfections of the experimental setup (ideally they should be 0), where $\theta_{1}$ parametrizes the residual amplitude modulation of the phase modulating EOM and $\theta_{2}$ parametrizes the phase locking offset of the squeezing measurement. $v^{2}_{sq}=\eta e^{-2r_{e}}+(1-\eta)$ and $v^{2}_{asq}=\eta e^{2r_{e}}+(1-\eta)$ are squeezing and anti-squeezing degrees in SNL units.

Note that the form of Eq. (V) rely on two assumptions: First, we assume that the modulation signal on the EOM is perfectly coherent so $V_{s}(\phi_{\mathrm{avg}})$ doesn’t have an offset term. This assumption is consolidated by the fact that we drive the EOM with a sine wave generated from a function generator with phase noise less than -65 dBc. Second, we ignore the phase fluctuations of the phase locking. This assumption is consolidated by the high ($\sim$32 dB) signal-to-noise ratio of the locking detector HDL, though this signal-to-noise ratio is not a direct measurement of the phase fluctuation.

V.1 Sensitivity fitting

With $V_{s}$ and $V_{n}$ extracted for a range of $\phi_{\mathrm{avg}}$ settings, we can estimate the sensitivity. By comparing the definition of $\sigma$ in Eq. (S4) with Eq. (V), the sensitivity to a small phase shift at a given $\phi_{\mathrm{avg}}$ offset is

[TABLE]

where $V^{\prime}_{s}=\partial V_{s}/\partial\phi_{\mathrm{avg}}$ is the partial derivative of $V_{s}$ with respect to $\phi_{\mathrm{avg}}$ and the $\sigma$ estimation is independent of SNL measurement since dividing $V_{n}$ with $V^{\prime}_{s}$ can cancel $V_{sn}$ out.

In the experiment we give an identical local phase shift to all 4 modes so that $\phi_{j}=\phi_{\mathrm{avg}}$ for all $j=1$ to $4$, and change the value of $\phi_{\mathrm{avg}}$ around both the squeezing $\phi_{\mathrm{avg}}=0$ and the anti-squeezing $|\phi_{\mathrm{avg}}|=90^{\circ}$. The $\phi_{\mathrm{avg}}$ we choose to induce as well as the fitting to Eq. (V) with measured $V_{s}$ and $V_{n}$ from Figure S10 is shown in Figure S11. In the $V_{s}(\phi_{\mathrm{avg}})$ fitting, the parameters to fit are the slope $k=\sqrt{\frac{2}{M}}V_{sn}\alpha_{e}$ and $\theta_{1}$. In the $V_{n}$ fitting, the squeezing noise voltage scaled by SNL, $k_{sq}=V_{sn}\cdot v_{sq}$, the anti-squeezing noise voltage scaled by SNL, $k_{asq}=V_{sn}\cdot v_{asq}$, and $\theta_{2}$ are fitting parameters. With the fitting result, we estimate the small angle sensitivity of our system $\sigma_{min}$ by

[TABLE]

We do fitting for 5 different pump power of OPO and find the fitted values of $\theta_{1}$ and $\theta_{2}$ are $3.4^{\circ}\pm 0.2^{\circ}$ and $1.6^{\circ}\pm 0.6^{\circ}$, respectively. These values are reasonably small, and in principle could be further reduced by better locking and phase modulation techniques. For the most sensitive case (maximum squeezing rate) in our result, the fitted $\theta_{1}$ and $\theta_{2}$ indicate $\sigma_{min}$ could have been further improved by $\sim 0.2\%$ and $\sim 0.9\%$, respectively. The $\sigma_{min}$ extracted in Eq. (S50) is shown in our experiment results for the entangled approach in main text Fig 3 as $\sigma_{e}$. The uncertainty of $V^{\prime}_{s}$ is obtained from the fitting, and the uncertainty of $V_{n}$ is obtained from the standard deviation of 2000 measurements. The error bars for $\sigma$ in Fig. 3 are calculated by error propagation of Eq. (S50).

A similar analysis method is used for the separable approach, but the PSDs used in the separable approach are from only one HD instead of averaged HD outputs. By removing the BSN, our setup gives the separable approach of $M=1$. To compare with the entangled approach of $M=4$, we rescale our result with $1/\sqrt{M}$ as a result of classical averaging. The scaled sensitivity is quoted as our experiment result for the separable approach in main text Fig. 3 as $\sigma_{s}$.

V.2 Resource counting

In this section, we show how to experimentally measure the average photon number per mode that we use in the phase sensing.

For the entangled approach, we estimate $MN_{e}=MN_{e,coh}+MN_{e,sqz}$ by comparing the joint measurement PSDs for squeezing and anti-squeezing quadrature to that for SNL, where $M=4$ is the mode number. With the notation defined above, the average number of squeezed photons for all modes in the entangled approach are obtained by comparing $V_{n}$ to $V_{sn}$ with

[TABLE]

Similarly, the average number of coherent photons are obtained by comparing $V_{s}$ to $V_{sn}$ with

[TABLE]

With Eq. (S51) and (S52), $N_{e}=N_{e,coh}+N_{e,sqz}$ gives the $N$ values in main text Fig. 3 for the entangled approach $\sigma_{e}$. The error bars for $N$ in Fig. 3, entangled approach are calculated by error propagation of Eqs. (S51-S52).

For the separable approach, we use a very similar technique. However, the PSD is from a single HD instead of joint measurement. Explicitly, the photon number per mode for the separable approach is $N_{s}=N_{s,coh}+N_{s,sqz}$, with

[TABLE]

and

[TABLE]

where $\phi$ is the phase shift of the single mode, and we use $V_{sn^{\prime}}$ to denote the 1-mode SNL, which is about $1/4$ of the 4-mode SNL $V_{sn}$ used in the entangled approach. The photon number per mode $N_{s}=N_{s,coh}+N_{s,sqz}$ gives the $N$ values in main text Fig. 3 for the separable approach $\sigma_{s}$. The error bars for $N$ in Fig. 3, separable approach are calculated by error propagation of Eqs. (60-61).