# Metrological detection of entanglement generated by non-Gaussian operations

David Barral, Mathieu Isoard, Giacomo Sorelli, Manuel Gessner, Nicolas

Treps, Mattia Walschaers

► To cite this version:

David Barral, Mathieu Isoard, Giacomo Sorelli, Manuel Gessner, Nicolas Treps, et al. Metrological detection of entanglement generated by non-Gaussian operations. New Journal of Physics, 2024, 26, 10.1088/1367-2630/ad6475. hal-05523769

HAL Id: hal-05523769

https://hal.science/hal-05523769v1

Submitted on 23 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

PAPER • OPEN ACCESS

# Metrological detection of entanglement generated by non-Gaussian operations

To cite this article: David Barral et al 2024 *New J. Phys.* **26** 083012

View the [article online](#) for updates and enhancements.

You may also like

- *Multipartite Gaussian entanglement measure with applications to graph states and bosonic field theory*
Matteo Gori, Matthieu Sarkis and Alexandre Tkatchenko

- *Dynamic entanglement for continuous variables in an electric field background*
Shu-Min Wu, Jin-Xuan Li, Xiao-Wei Fan et al.

- *From squeezing to Gaussian entanglement via beamsplitters*
Shuangshuang Fu, Shunlong Luo and Yue Zhang
<figure><img src="image_4.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

4 April 2024

REVISED

2 July 2024

ACCEPTED FOR PUBLICATION

17 July 2024

PUBLISHED

12 August 2024

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_5.png" /></figure>

PAPER

# Metrological detection of entanglement generated by non-Gaussian operations

David Barral¹,*, Mathieu Isoard¹, Giacomo Sorelli¹,² ID, Manuel Gessner³ ID, Nicolas Treps¹ and Mattia Walschaers¹ ID

¹ Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-Université PSL, Collège de France, 4 place Jussieu, Paris F-75252, France

² Fraunhofer IOSB, Ettlingen, Fraunhofer Institute of Optronics, System Technologies and Image Exploitation, Gutleuthausstr. 1, 76275 Ettlingen, Germany

³ Departamento de Física Teórica, IFIC, Universidad de Valencia-CSIC, C/ Dr Moliner 50, Burjassot, Valencia 46100, Spain

* Author to whom any correspondence should be addressed.

E-mail: david.barral@lkb.ens.fr and mattia.walschaers@lkb.upmc.fr

**Keywords:** entanglement, Fisher information, non-Gaussian, continuous variables, metrology, homodyne detection

Supplementary material for this article is available online

## Abstract

Entanglement and non-Gaussianity are physical resources that are essential for a large number of quantum-optics protocols. Non-Gaussian entanglement is indispensable for quantum-computing advantage and outperforms its Gaussian counterparts in a number of quantum-information protocols. The characterization of non-Gaussian entanglement is a critical matter as it is in general highly demanding in terms of resources. We propose a simple protocol based on the Fisher information for witnessing entanglement in an important class of non-Gaussian entangled states: photon-subtracted states. We demonstrate that our protocol is relevant for the detection of non-Gaussian entanglement generated by multiple photon-subtraction and that it is experimentally feasible through homodyne detection.

## 1. Introduction

Entanglement is considered one of the most striking breakthroughs of 20th century science. The gedanken experiment proposed by Einstein, Podolsky and Rosen in 1935 [1] pointed out the notion of inseparability of a state composed by two quantum particles spatially distanced with maximally correlated momenta and maximally anti-correlated positions. Nowadays, entanglement stands as a physical resource underpinning most of current development in quantum technologies [2]. The efficient detection and measurement of entanglement is a very active area of quantum physics [3], being far from simple especially for continuous variable (CV) systems which involve physical quantities with a continuous spectrum of values [4].

Multimode squeezed states of light are the cornerstone of CV quantum networks [5, 6]. They exhibit Gaussian statistics and their entanglement properties are completely specified by their covariance matrix. Criteria and witnesses for this Gaussian entanglement have been proposed and tested for decades [7–11]. However, Gaussian entanglement can always be undone with passive linear optics, a phenomenon generally refereed to as passive separability [12]. It was recently found that one requires states that are not passively separable as a resource for a quantum computational advantage [13]. Because all Gaussian states are passively separable, we can always find mode bases in which the covariance matrix of the state will not show any direct signature of entanglement. Yet, this is not the case for the class of non-passively-separable states, where those changes of basis do not disentangle the state. Because this entanglement is generated by non-Gaussian operations and hidden in the non-Gaussian features of the state, we will here refer to it as *non-Gaussian entanglement* [14, 15]. The goal of this work is to find a practical way to detect this non-Gaussian entanglement.

In order to characterize non-Gaussian entanglement, a number of criteria based on high-order moments and on uncertainty relations of different classes of operators have been proposed [16–23]. Nevertheless, these criteria are far from being feasible with current experimental methods. Other more experimentally-friendly
criteria are based on the Shannon entropy and the fidelity of teleportation in quantum channels [24, 25]. Here we tackle the problem from an operational point of view: non-Gaussian quantum correlations are also known to improve metrological sensitivity, the performance of quantum key distribution and quantum teleportation protocols [26–29]. The advantage of relying on the improvement of quantum protocols is that the detected entanglement is useful by design. In this article, we will focus specifically on metrological protocols, where quantum estimation tools have been devised to witness entanglement [30–32]. These witnesses are based on the fact that metrological sensitivity determines precision of measurements and this sensitivity is limited for separable states. This can be used to detect entanglement. Two powerful assets of these sensitivity-based witnesses are i) they do not make assumptions about the quantum state—Gaussianity, purity, etc and ii) they contain information about all high-order moments.

By using the approach of [31, 32] taking into account the limitations of CV quantum optics, we propose a general protocol based solely on homodyne detection, using both the variance of the measurement outcomes and the joint measurement statistics. Our protocol is efficient in terms of resources as the parameter estimation is done in postprocessing using solely the data collected by homodyne detection. We show its relevance analyzing an important class of non-Gaussian entangled states: photon-subtracted states. We demonstrate that our protocol is pertinent for the detection of non-Gaussian entanglement and that it is experimentally feasible.

The article is organized as follows: we first present our protocol to detect entanglement through homodyne detection and postprocessing of the joint probability distribution based on the metrological witness introduced in [31, 32] in section 2. We then present in section 3 the probe states that we will use to test our non-Gaussian entanglement witness. In section 4 we analyze which parameter is best suited to measure entanglement in our metrological protocol and calculate entanglement in an ideal case. In section 5 we study a realistic case taking into account unbalanced input squeezing, losses and discretization of the measurement outcomes. We finally discuss possible experimental implementations of our scheme, their limitations and feasibility in section 6 and we present our conclusions in section 7.

## 2. Entanglement detection via local homodyne detection and postprocessing

We consider here the following problem: two experimenters, Alice and Bob, who share an optical quantum state $\hat{\rho}_{AB}$, want to elucidate if their shared state is entangled or not, while minimizing the amount of experimental resources. If the input state is Gaussian, they just need to measure the variances of linear combinations of optical-field quadratures and apply second-order moment-based criteria like for instance those of Duan et al, Simon or Giovanetti et al [7–9]. This can be easily implemented experimentally using homodyne detection. However, the larger class of non-Gaussian states do not always present entanglement that can be revealed by second-order moment-based criteria. Particularly, the majority of entanglement criteria for quantum states with purely non-Gaussian correlations are based on either carrying out full quantum-state tomography [33] or measuring high-order moment correlations, protocols which are very demanding experimentally.

Here, we apply a metrological protocol to detect entanglement. Alice and Bob share information in order to estimate jointly a parameter $\theta$ generated by a Hamiltonian $\hat{H} = \hat{H}_A + \hat{H}_B$ that acts locally on both subsystems such that $\hat{\rho}_{AB}^{\theta} = e^{-i\theta\hat{H}}\hat{\rho}_{AB}e^{i\theta\hat{H}}$ (see figure 1). We define the joint probability density, i.e. the conditional probability to obtain a set of local measurement outcomes ($\xi_A, \xi_B$) given the parameter $\theta$, as

$$ \mathcal{P}(\xi_A, \xi_B | \theta) = \mathrm{Tr} \left[ \hat{\rho}_{AB}^{\theta} \hat{\Pi} \right], $$

where $\hat{\Pi}$ is a positive-operator valued measure (POVM) such that $\int \hat{\Pi} d^2\xi = 1$, where $d^2\xi = d\xi_A d\xi_B$. The metrological protocol consists of measuring the Fisher information (FI) defined as

$$ F(\mathcal{P}(\xi_A, \xi_B | \theta)) = \int_{\mathbb{R}^2} \mathcal{P}(\xi_A, \xi_B | \theta) \left( \frac{\partial \mathcal{L}(\xi_A, \xi_B | \theta)}{\partial \theta} \right)^2 d^2\xi, \quad (1) $$

where $\mathcal{L}(\xi_A, \xi_B | \theta) = \log(\mathcal{P}(\xi_A, \xi_B | \theta))$ represents the logarithmic likelihood related to the probability density $\mathcal{P}(\xi_A, \xi_B | \theta)$. In our case, as illustrated in figure 1, the observables will correspond to local homodyne measurements $\hat{\xi}_A = \cos \phi_A \hat{x}_A + \sin \phi_A \hat{p}_A$ and $\hat{\xi}_B = \cos \phi_B \hat{x}_B + \sin \phi_B \hat{p}_B$, where $\phi_A, \phi_B$ are two angles, and $\hat{x}_A, \hat{x}_B, \hat{p}_A, \hat{p}_B$ are the quadratures operators defined from the annihilation operators as

$$ \hat{a}_I = \frac{\hat{x}_I + i\hat{p}_I}{2}, \quad I \in A, B. \quad (2) $$

The quadrature operators thus satisfy the commutation rules $[\hat{x}_I, \hat{p}_J] = 2i\delta_{IJ}$, $I, J \in \{A, B\}$.
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Sketch of the proposed metrological protocol for entanglement detection. Alice and Bob share a quantum state $\hat{\rho}_{AB}$. They jointly estimate a parameter $\theta$ generated by two local Hamiltonians $\hat{H}_{A/B}$. Using two homodyne detectors with a common phase reference, Alice and Bob can retrieve the parameter-dependent joint probability distribution $\mathcal{P}(x_A, x_B|\theta)$, and thus the Fisher information related to this parameter estimation, and the local variances of the Hamiltonians $\hat{H}_{A/B}$. With this in hand, Alice and Bob can jointly compute the metrological entanglement witness of equation (4).</figcaption></figure>

Then, if $\hat{\rho}_{AB}$ is separable, the FI of equation (1) for a state $\hat{\rho}_{AB}^{\theta}$ generated by $\hat{H}$ is upper bounded by [31, 32]

$$F(\mathcal{P}(\xi_A, \xi_B|\theta)) \le 4\mathrm{Var}[\hat{\rho}_A, \hat{H}_A] + 4\mathrm{Var}[\hat{\rho}_B, \hat{H}_B], \quad (3)$$

where $\hat{\rho}_{A/B}$ are the reduced density matrices for systems A and B, respectively. Because this is a necessary condition for separability, its violation is a sufficient criterion for entanglement.

Therefore, we can introduce the following metrological witness of entanglement

$$E = F(\mathcal{P}(\xi_A, \xi_B|\theta)) - 4(\mathrm{Var}[\hat{\rho}_A, \hat{H}_A] + \mathrm{Var}[\hat{\rho}_B, \hat{H}_B]) > 0. \quad (4)$$

This inequality can reveal entanglement but not its origin—Gaussian or non-Gaussian. However, as we will discuss below, this witness can reveal entanglement generated by non-Gaussian operations where Gaussian criteria fail to detect it. Therefore, this non-Gaussian entanglement leads to a better sensitivity to estimate a given parameter $\theta$ than one would have obtained with a separable state and becomes a useful resource for metrology.

The other interest of the witness (4) is that it also holds for any state, pure or mixed, and its major asset is the practicability of its computation. Homodyne measurements in each mode with a common phase reference allow us to access experimentally i) the joint probability distribution $\mathcal{P}(\xi_A, \xi_B|\theta)$, and thus the FI, and ii) the variances associated to the local generators, enabling to test the entanglement witness given by equation (4). Moreover, in some cases (see section 4) the parameter-dependence of the joint probability distribution $\mathcal{P}(\xi_A, \xi_B|\theta)$ can be generated in postprocessing applying appropriate transformations directly to the joint probability distribution as $\mathcal{P}(\xi_A, \xi_B|\theta) = \mathcal{P}(U_{\theta}(\xi_A), U_{\theta}(\xi_B))$, with $U_{\theta}(\xi_{A/B})$ the transformation related to the Hamiltonian $\hat{H}_{A/B}$ in the quadrature space $\xi_{A/B}$ [36, 37]. This important feature avoids to apply impractical inline transformations to the state simplifying greatly the detection of entanglement.

In equation (4), the FI is always bounded by the Quantum Fisher information (QFI) defined as

$$F_Q[\hat{\rho}_{AB}, \hat{H}] = \max_{\hat{\Pi}} F(\mathrm{Tr}[\hat{\rho}_{AB}^{\theta} \hat{\Pi}]). \quad (5)$$

This bound is reached when the measurement observable is optimized over all possible POVM $\hat{\Pi}$ [35]. In this case, the entanglement witness (4) is maximized and we denote this quantity by $E_Q$:

$$\begin{aligned} E_Q &\equiv \max_{\hat{\Pi}} E \\ &= F_Q[\hat{\rho}_{AB}, \hat{H}] - 4(\mathrm{Var}[\hat{\rho}_A, \hat{H}_A] + \mathrm{Var}[\hat{\rho}_B, \hat{H}_B]). \end{aligned} \quad (6)$$

Note that $E_Q$ might not be attained when we restrict ourselves to homodyne measurements since it might not correspond to the optimal POVM to saturate the QFI.

For pure states we can easily obtain the QFI from the variance of the generator $\hat{H}$ of the parameter $\theta$ as

$$F_Q[\hat{\rho}_{AB}, \hat{H}] = 4\mathrm{Var}[\hat{\rho}_{AB}, \hat{H}].$$
Applying this identity into equation (4) we obtain the following simple condition for entanglement

$$E_Q = 8 \operatorname{Cov}[\rho_{AB}; \hat{H}_A, \hat{H}_B] > 0, \qquad (7)$$

where $\operatorname{Cov}[\rho_{AB}; \hat{H}_A, \hat{H}_B] = \langle \hat{H}_A \hat{H}_B \rangle - \langle \hat{H}_A \rangle \langle \hat{H}_B \rangle$. As a consequence any correlation that is seen in a bipartite pure state is a signature of entanglement.

As a final remark, we note that there are other limits on the metrological sensitivity that have been derived for states with a classical P function [38], but non-classicality and entanglement are not necessarily directly related since non-classicality can exist in a single mode system. However, in some cases quantum correlations and single-mode non-classical features are intertwined. For instance, photon subtraction can induce Wigner-negativity in an entangled mode if that entangled mode can perform Einstein–Podolsky–Rosen steering in the mode of subtraction [39].

### 3. Application to photon-subtracted states

The protocol described in the previous section is valid for any CV system, regardless of the nature of the state under consideration, as long as one has access to the probability distributions of each subsystem. In this section we introduce the states that we will use as a probe of our non-Gaussian entanglement criterion, namely, photon-subtracted states. In particular we will analyze bipartite states without Gaussian correlations in order to focus on their non-Gaussian features.

We consider two-mode photon subtracted states. This class of states has been demonstrated in optical systems using different degrees of freedom, such as polarization or frequency modes [33, 40]. In section 6 we will explain in detail different experimental methods for their production. Let us consider two independent single-mode squeezed states respectively related to Alice and Bob

$$|\Psi_0\rangle = \hat{S}_A(r_A) \hat{S}_B(r_B) |00\rangle, \qquad (8)$$

where $\hat{S}_I(r_I) = \exp\{(-r_I/2)(\hat{a}_I^2 - \hat{a}_I^{\dagger 2})\}$ is the single-mode squeezing operator, and $r_I \in \mathbb{R}$ is the squeezing parameter for each mode $I = A, B$. The amount of squeezing in decibels is given by $s_I = 10 \log_{10}(e^{-2r_I})$.

In what follows we analyze two cases: in-phase squeezing ($r_A > 0$ and $r_B > 0$) and in-quadrature squeezing ($r_A > 0$ and $r_B < 0$).

Equation (8) corresponds to a Gaussian state and all its information is encoded in the covariance matrix $V_0 = \operatorname{diag}(e^{-2r_A}, e^{2r_A}, e^{-2r_B}, e^{2r_B})$, written with respect to the vector of amplitude and phase quadratures in each mode $\xi = (x_A, p_A, x_B, p_B)^T$. Note that $V_0$ does not present off-diagonal terms, thus the input Gaussian state is fully separable.

Next, we perform a delocalized subtraction of n photons in the same or different modes on this state. This operation produces in general a superposition of squeezed Fock states. For multiphoton subtraction the resulting state is

$$|\Psi\rangle \propto \prod_{j=1}^{n} (\cos(\phi_j) \hat{a}_A + \sin(\phi_j) \hat{a}_B) |\Psi_0\rangle,$$

where the parameters $\phi_j$ control the probability of subtraction in each mode for each subtraction $j = 1, \dots, n$. For $n = 1$ we have [41]

$$|\Psi\rangle \propto (\cos(\phi) \hat{a}_A + \sin(\phi) \hat{a}_B) |\Psi_0\rangle = \hat{S}_A(r_A) \hat{S}_B(r_B) (\cos(\phi) \sinh(r_A) |10\rangle + \sin(\phi) \sinh(r_B) |01\rangle),$$

where we have used the Bogolyubov transformation $\hat{S}_I^\dagger(r_I) \hat{a}_I \hat{S}_I(r_I) = \cosh(r_I) \hat{a}_I + \sinh(r_I) \hat{a}_I^\dagger$. A sketch of this operation is shown in figure 2. The wavefunction of a single-photon subtracted state in the amplitude quadratures of the optical field is given by

$$\Psi(x_A, x_B) \equiv \langle x_A, x_B | \Psi \rangle \propto e^{-\frac{c^2 r_A x_A^2 + c^2 r_B x_B^2}{4}} [(e^{2r_A} - 1) \cos(\phi) x_A + (e^{2r_B} - 1) \sin(\phi) x_B]. \qquad (9)$$

Examples of joint probability distributions $P(x_A, x_B) = |\Psi(x_A, x_B)|^2$ for a photon subtracted state given by equation (9) with $\phi = \pi/4$ and $r_A = r_B = 0.2$, $r_A = -r_B = 0.2$, are respectively shown in figures 3(a) and (b).

For $n = 2$ we can consider for instance subtraction in the same ($\phi_1 = \phi_2 = \phi$) or orthogonal modes ($\phi_1 = \phi, \phi_2 = \phi + \pi/2$) as
<figure><img src="image_7.png" /><figcaption>Figure 2. Sketch of an optical setup for delocalized single-photon subtraction. Alice and Bob prepare two squeezed states in given optical modes. A small fraction of each mode power is diverted to a common beam splitter with a transmittivity controlled by a parameter $\phi$. An event measured by a single-photon detector heralds the subtraction of a photon delocalized between the two modes.</figcaption></figure>

$$|\Psi_{\parallel}\rangle \propto (\cos(\phi) \hat{a}_{A} + \sin(\phi) \hat{a}_{B})^2 |\Psi_{0}\rangle,$$

$$|\Psi_{\perp}\rangle \propto (\cos(\phi) \hat{a}_{A} + \sin(\phi) \hat{a}_{B}) (\sin(\phi) \hat{a}_{A} - \cos(\phi) \hat{a}_{B}) |\Psi_{0}\rangle.$$

Taking in-phase squeezing, $r_A = r_B \equiv r$ and subtraction along $\phi = \pi/4$, we obtain

$$|\Psi_{\parallel}\rangle \propto \hat{S}_{A}(r) \hat{S}_{B}(r) (\cosh(r)|00\rangle + \sinh(r)|11\rangle + (1/\sqrt{2})(|20\rangle + |02\rangle)), \quad (10)$$

$$|\Psi_{\perp}\rangle \propto \hat{S}_{A}(r) \hat{S}_{B}(r) (|20\rangle - |02\rangle).$$

Thus, subtracting photons along different axes produce states with completely different entanglement features. Examples of joint probability distributions for two-photon subtracted states are gathered in the supplementary material.

Coming back to the case of single-photon subtracted states, the entanglement present in these states is not grasped by Gaussian entanglement witnesses: this can be generally understood from the covariance matrix of a photon-subtracted state. [12] shows that this covariance matrix for a single-photon subtraction can generally be written as

$$V = V_0 + 2 \frac{(V_0 - 1) P (V_0 - 1)}{\mathrm{Tr}[(V_0 - 1) P]}, \quad (11)$$

where $V_0$ is the initial Gaussian state's covariance matrix and $P$ is a matrix that projects on the phase space axes associated with the mode of photon subtracted states. In our present case, we find that

$$P = \begin{pmatrix} \cos^2(\phi) & 0 & \frac{1}{2} \sin(2\phi) & 0 \\ 0 & \cos^2(\phi) & 0 & \frac{1}{2} \sin(2\phi) \\ \frac{1}{2} \sin(2\phi) & 0 & \sin^2(\phi) & 0 \\ 0 & \frac{1}{2} \sin(2\phi) & 0 & \sin^2(\phi) \end{pmatrix}.$$

Since $P$ is a positive matrix and $V_0 - 1$ a symmetric matrix, $(V_0 - 1)P(V_0 - 1)$ is also positive. Thus, as we see from equation (11), on the level of the covariance matrix the photon subtraction only adds Gaussian noise. This implies that no additional entanglement can be witnessed by purely looking at the covariance matrix [34]. As a consequence, since we set $V_0 = \mathrm{diag}(e^{-2r_A}, e^{2r_A}, e^{-2r_B}, e^{2r_B})$, we find that $V$ should not display any entanglement.

On the contrary, for two-photon subtracted states, Gaussian witnesses can reveal entanglement in some cases. This results from constructive interferences between both photon subtractions which induce quantum correlations at the level of the covariance matrix. In particular, when the two subtractions happen in the same mode along $\phi = \pi/4$ and when the squeezing is low, we can see directly from expression (10) that Gaussian correlations appear. Indeed, in that case, a part of $|\Psi_{\parallel}\rangle$ is proportional to $|00\rangle + \tanh(r)|11\rangle$, which is the expansion of a two-mode squeezed state at order $O(r)$, and thus induces Gaussian correlations that can be witnessed at the level of the covariance matrix.

However, Gaussian witnesses are very sensitive to the parameters of the subtracted states: for instance, when the squeezing parameters of the initial Gaussian state are too high or when the subtractions happen in
<figure><img src="image_8.png" /><figcaption>Figure 3. Contour plots of the joint probability distribution for a single-photon subtracted state given by equation (9) with $\phi = \pi/4$ and (a) $r_A = r_B = 0.2$, (b) $r_A = -r_B = 0.2$ ($|s_{A/B}| = 1.74$ dB).</figcaption></figure>

orthogonal modes, Gaussian witnesses cease to work. The metrological protocol proposed in this paper circumvents this issue since it hinges on the Fisher information which can reveal correlations overlooked by witnesses only based on second order moments. As discussed in the following sections, the metrological criterion will always be able to detect entanglement for two-photon subtracted states regardless to their properties. We discuss in particular at the end of section 5 an interesting case where Gaussian witnesses fail to detect entanglement while the metrological criterion works and can be used directly with a limited amount of homodyne data.

## 4. Ideal detection of non-Gaussian entanglement

In order to decide which Hamiltonian $\hat{H}$ is best suited to witness entanglement we can calculate theoretically $E_Q$ through equation (7). This can guide us deciding which parameter is best suited to detect entanglement of a given quantum state in a realistic scenario. Below, we use the estimation of parameters related to the four single-mode Gaussian gates in CV quantum optics. Namely: displacement, phase-shift, shearing and squeezing [42]. We analyze in which cases the joint estimation of these parameters reveals the entanglement of the single-photon subtracted state given by equation (9). For the sake of simplicity, we focus here on the case $\phi = \pi/4$. A generalization to any angle $\phi$ and an equivalent analysis for two photon-subtracted states can be found in the supplemental material.

### 4.1. Displacement

A displacement of $\theta$ along the axis $x_A = \pm x_B$ is produced by the following operator

$$ \hat{D}(\theta) = e^{-i\theta(\hat{p}_A \pm \hat{p}_B)/2}. $$
<figure><img src="image_9.png" /><figcaption>Figure 4. Displacement-estimation entanglement witness $E_Q$ for one-photon subtracted states (9) given by equation (12) versus squeezing of the input squeezed states $s_A$ and $s_B$. (a) Displacement along $x_A = x_B$ optimizes $E_Q$ for states with in-phase input squeezing. (b) Displacement along $x_A = -x_B$ optimizes $E_Q$ for states with in-quadrature input squeezing. The black dashed lines correspond to the maximum value of $E_Q$.</figcaption></figure>

The Hamiltonian related to this displacement operator is $H_{\pm} = (\hat{p}_A \pm \hat{p}_B)/2$. The optimal entanglement witness $E_Q$ obtained with displacement operators along the axis $x_A = \pm x_B$ for a pure photon-subtracted state given by equation (9) with $\phi = \pi/4$ is

$$E_Q = \pm 2 (\langle \hat{p}_A \hat{p}_B \rangle - \langle \hat{p}_A \rangle \langle \hat{p}_B \rangle) = \pm 2 e^{r_A + r_B} \cos(\epsilon), \quad (12)$$

with

$$\cos(\epsilon) = \frac{2 \sinh(r_A) \sinh(r_B)}{\sinh^2(r_A) + \sinh^2(r_B)},$$

where we have used the condition of entanglement for pure states given by equation (7).

Displacement along either $x_A = x_B$ or $x_A = -x_B$ detects entanglement respectively for in-phase squeezing ($r_A, r_B > 0$) and in-quadrature—orthogonal—squeezing ($r_A > 0, r_B < 0$). Figures 4(a) and (b) show contour plots of optimal entanglement witness $E_Q$ in the two cases. States with in-phase input squeezing show always a larger violation of $E_Q$ due to the argument $r_A + r_B$ in equation (12). One can also optimize the input squeezing parameters $r_A$ and $r_B$ to maximize the witness $E_Q$. These values correspond to the black dashed lines in figures 4(a) and (b). For in-phase squeezing the witness reaches its maximum at $r_A = r_B$ and is given by $E_Q = 2e^{2r_A}$. Likewise, for in-quadrature squeezing the maximum value of $E_Q$ is not along the diagonal (see the dashed line in figure 4(b)), but below it. For a given value of $r_A$, $E_Q$ reaches its maximum for $r_B = \log(1/(1 + 2 \sinh(r_A)))^{1/2}$ and is given by

$$E_Q = \frac{2 e^{r_A}}{1 + \sinh(r_A)}.$$

The shapes of figures 4(a) and (b) can be explained in terms of the symmetries of the two functions that compose equation (12): $\pm \cos(\epsilon)$ is a symmetric function with respect to the diagonal $s_A = s_B$ for every input
<figure><img src="image_10.png" /><figcaption>Figure 5. Phase-estimation entanglement witness $E_Q$ for one-photon subtracted states (9) given by equation (13) versus squeezing of the input squeezed states $s_A$ and $|s_B|$. The black dashed line corresponds to the maximum value of $E_Q$.</figcaption></figure>

squeezing, whereas $2e^{r_A+r_B}$ is symmetric with respect to the diagonal (antidiagonal, in this case along $s_A - s_B = 6$ dB) for in-phase (in-quadrature) squeezing.

Importantly, we obtain the same result calculating the entanglement witness through equation (4), $E = E_Q$, indicating that the FI saturates the QFI. The result of equation (12) is particularly interesting because, following equation (7), second order moments of the distribution reveal entanglement with a non-Gaussian origin.

Recently, Tian et al analyzed the multipartite entanglement in a nondegenerate triple photon state using a metrological criterion [43]. They claimed there that non-Gaussian entanglement cannot be sufficiently captured by linear quadratures, i.e. displacements. While this is the case for triple photon states, we have shown that it does not hold in general: displacements can detect non-Gaussian entanglement of photon-subtracted states.

## 4.2. Phase shift

The phase-shift operator

$$ \hat{R}(\theta) = e^{-i\theta(\hat{N}_A \pm \hat{N}_B)} \propto e^{-i\theta(\hat{x}_A^2 + \hat{p}_A^2 \pm \hat{x}_B^2 \pm \hat{p}_B^2)/4} $$

rotates the state by a phase $\theta$ in local phase subspaces in either clockwise-clockwise (+) or clockwise-counterclockwise (-) directions. The related Hamiltonian is $H_{\pm} = \hat{N}_A \pm \hat{N}_B$. The optimal entanglement witness is in this case

$$ E_Q = \pm 8 (\langle \hat{N}_A \hat{N}_B \rangle - \langle \hat{N}_A \rangle \langle \hat{N}_B \rangle) = \mp 2 \cosh(2r_A) \cosh(2r_B) \cos^2(\epsilon). \quad (13) $$

Entanglement is always detected for clockwise-counterclockwise (-) phase shifts, but not for clockwise-clockwise (+) as it is just a global phase shift. Figure 5 shows contour plots of the entanglement witness $E_Q$ for different values of squeezing. Notably, the detection of entanglement does not depend on the sign of the input squeezed states as $E_Q$ is invariant under change of sign of the squeezing parameters $r_{A/B}$. The detected entanglement witness is maximum for $r_A = |r_B|$ (dashed line along the diagonal in figure 5) being $E_Q = 2 \cosh^2(2r_A)$.

One can wonder if the Fisher information in equation (4) reaches the QFI in this case. While measuring the joint probability distribution in the $(x_1, x_2)$-plane was enough to obtain the maximal value of the FI and saturates the QFI for the displacement estimation, here the situation is a bit more complicated. For simplicity, we consider the case $r_A = r_B$ in what follows. The FI can be optimized by finding the set of angles $(\phi_A, \phi_B)$ of the measurement outcomes $\xi_A = \cos \phi_A x_A - \sin \phi_A p_A$, $\xi_B = \cos \phi_B x_B - \sin \phi_B p_B$ for which the joint probability distribution $P(\xi_A, \xi_B|\theta)$ leads to the best value of the FI (see the supplemental material). However, we find that such local rotations are not enough to saturate the QFI, and that only a mixing of modes A and B before the homodyne detectors can lead to a saturation of the QFI. It is indeed possible to prove that a non-local rotation of $-\pi/4$ between modes A and B and measuring the joint probability distribution $P(x'_A, p'_B|\theta)$, with $x'_A = (x_A - x_B)/\sqrt{2}$ and $p'_B = (p_A + p_B)/\sqrt{2}$ is needed to saturate the QFI.
### 4.3. Shearing

The shearing—also known as phase-gate—operator

$$ \hat{S}(\theta) = e^{-i\theta(\hat{x}_A^2 \pm \hat{x}_B^2)/4} $$

shears the state with respect to the axes $x_A$ and $\pm x_B$ by a gradient of $\theta$. The transformation in phase space is given by $(x_I, p_I) \to (x_I, p_I + \theta x_I)$ [42]. The related Hamiltonian is $H_{\pm} = (\hat{x}_A^2 \pm \hat{x}_B^2)/4$. The optimal entanglement witness is in this case

$$ E_Q = \pm \frac{\langle \hat{x}_A^2 \hat{x}_B^2 \rangle - \langle \hat{x}_A^2 \rangle \langle \hat{x}_B^2 \rangle}{2} = \mp \frac{e^{-2(r_A + r_B)}}{2} \cos^2(\epsilon). \quad (14) $$

Thus, shearing with respect to $x_A$ and $x_B$ does not detect entanglement. However, shearing with respect to $x_A$ and $-x_B$ captures it. Note that in this case the entanglement witness $E_Q$ is maximized for $r_{A/B} < 0$, i.e. squeezing along the quadratures $p_{A/B}$, unlike displacement and phase estimation where $E_Q$ is maximized for squeezing along $x_{A/B}$.

Figures 6(a) and (b) show contour plots of $E_Q$ in the cases of input squeezing along the same quadrature (a) or along different quadratures (b). For input squeezing along the same quadratures the detected entanglement witness is maximum again for $r_A = r_B$ (dashed line along the diagonal in figure 6(a)) and given by $E_Q = e^{-4r_A}/2$. However, for input squeezing along different quadratures, the maximum value for $E_Q$ is below the diagonal, as it also happens for displacement. For a given value of $r_A$, the maximum $E_Q$ is obtained for $r_B = (-r_A + \log(1 + e^{r_A} - e^{2r_A}))/2$ (dashed line in figure 6(b)) and is given by

$$ E_Q = \frac{e^{-2r_A}}{2(-1 + \sinh(r_A))^2}. $$

The shape of figures 6(a) and (b) is explained in the same way as for displacement.

Here again, we optimize the FI to see if it is possible to reach the bound $E = E_Q$. The same analysis as in the case of the phase-shift operator (by performing local rotations before the homodyne detection) is summarized in the supplemental material. The same conclusion follows: the FI never reaches the QFI, and only a non-local rotation of $-\pi/4$ between modes A and B leads to a saturation of the QFI from quadrature measurements.

### 4.4. Squeezing

The squeezing operator

$$ \hat{S}(\theta) = e^{-i\theta(\hat{x}_A\hat{p}_A + \hat{p}_A\hat{x}_A \pm \hat{x}_B\hat{p}_B \pm \hat{p}_B\hat{x}_B)/4} $$

squeezes the position quadratures of modes A and B by a factor of $e^{\theta}$ (+) or squeezes the position quadratures of A by $e^{\theta}$ and stretches those of B by $e^{-\theta}$ (-). The related Hamiltonian is $H_{\pm} = (\hat{x}_A\hat{p}_A + \hat{p}_A\hat{x}_A \pm \hat{x}_B\hat{p}_B \pm \hat{p}_B\hat{x}_B)/4$. The optimal entanglement witness is here

$$ E_Q = \pm \frac{\langle \{\hat{x}_A, \hat{p}_A\} \{\hat{x}_B, \hat{p}_B\} \rangle - \langle \{\hat{x}_A, \hat{p}_A\} \rangle \langle \{\hat{x}_B, \hat{p}_B\} \rangle}{2} = 0, $$

where $\{\hat{x}, \hat{p}\} = \hat{x}\hat{p} + \hat{p}\hat{x}$. Interestingly, the joint estimation of the squeezing parameter does not detect entanglement in any of the above two cases.

### 4.5. Comparison and resource evaluation

In order to decide which parameter-estimation strategy is best suited to detect entanglement we show in figure 7 the evolution of the maximal value of $E_Q$ when optimizing input squeezing parameters $r_A$ and $r_B$ to follow the black dashed curves in figures 4–6 versus amount of squeezing in dB in Alice’s mode in the case of in-phase and in-quadrature input squeezing. As discussed in sections 4.1–4.3, for in-phase input squeezing the maximum of $E_Q$ is obtained for $r_A = r_B$ (dashed line along the diagonal in figures 4(a), 5 and 6(a)). For in-quadrature input squeezing the maximum is obtained for $r_A = -r_B$ for phase-shift (dashed line along the diagonal in figure 5), and for $r_B = \log(1/(1 + 2\sinh(r_A))^{1/2})$ and $r_B = (-r_A + \log(1 + e^{r_A} - e^{2r_A}))/2$ for displacement and shearing, respectively (dashed line below the diagonal in figures 4(b) and 6(b)). Remarkably, for values of input squeezing lower than ≈5 dB, the best strategy is to jointly estimate the displacement (solid, blue). For larger values of input squeezing, phase shift and shearing estimation offer a greater sensitivity to entanglement (green and orange, respectively). On the contrary, as we saw above the joint estimation of the squeezing parameter does not offer any information on the entanglement of this state (solid, gray).
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Shearing-estimation entanglement witness $E_Q$ for one-photon subtracted states (9) given by equation (14) versus squeezing of the input squeezed states $s_A$ and $s_B$. (a) $E_Q$ for states with input squeezing along the same quadrature. (b) $E_Q$ for states with input squeezing along different quadratures. The black dashed lines correspond to the maximum value of $E_Q$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 7.** Maximal value of the entanglement witness $E_Q$, i.e. when input squeezing parameters $r_A$ and $r_B$ are optimized to follow the dashed curves in figures 4–6, versus squeezing in Alice's mode: displacement estimation (blue), shearing estimation (orange), phase shift estimation (green), and squeezing estimation (gray). In-phase (in-quadrature) input squeezing in solid (dashed). For phase shift estimation (green) the curve is the same in both cases. $E_Q > 0$ witnesses entanglement.</figcaption></figure>

In terms of resources, displacement estimation is also advantageous. Both probability distributions and quadrature variances corresponding to Alice and Bob can be directly measured with homodyne detection. Likewise, shearing estimation can be performed with homodyne detection, but fourth-order moments of the distributions (kurtosis) are necessary, which implies in general larger data sets. In the case of phase estimation, photon-number variances are necessary, which implies adding complexity to the detection.

Another great advantage of displacement estimation is that the displacement operation can be applied in post-processing: once the probability distribution $P(x_A, x_B|0)$ is measured, the displaced probability
distribution [under the Hamiltonian $\hat{H}_{\pm} = (\hat{p}_{A} \pm \hat{p}_{B})/2$] is directly given by $P(x_{A}, x_{B}|\theta) = P(x_{A} + \theta, x_{B} \pm \theta|0)$ [37], from which one can compute the classical FI (see section 5.3)—which we know saturates the QFI in this case, leading to the best possible estimation. On the contrary, the shearing and phase-shift operations can not be implemented in post-processing using just the probability distribution as full information about the state is needed for such operations. Thus, shearing and phase-shift unitaries have to be implemented at the level of the experimental setup or by post-processing after measuring the full quantum state by Wigner tomography [33] or Husimi Q function sampling [44]. In addition to this complication, contrary to the displacement operation, as we pointed out in sections 4.2 and 4.3, the FI can be optimized with local rotations of the measured quadratures, but only saturates the QFI when one mixes modes A and B.

# 5. Realistic detection of non-Gaussian entanglement

In this section we study the measurement of entanglement in a realistic scenario. As we found above, estimating displacement is the best strategy for ideal detection at moderate values of squeezing. Moreover, it is the simplest one, as the variances of the generators—field quadratures—are directly measured with homodyne detection. We thus focus on this option in the following. A similar analysis could be carried out for shearing and phase-shift estimation. Below we analyze the effect of unbalancing the sensitivity in the displacement estimation, the effect of losses on the detection of entanglement and the discretization of the sampled data to build a joint probability distribution and calculate the Fisher information.

## 5.1. Optimization of displacement axis for entanglement witness

In the previous section, we analyzed the detection of entanglement through displacement estimation when displacing the input state along the axes $x_A = \pm x_B$. However, we can optimize the entanglement witness displacing the input state along an axis different to $x_A = \pm x_B$ or, in other words, unbalancing the sensitivity related to Alice and Bob in the joint parameter estimation. The idea is the following: instead of displacing the same amount $(1, \pm 1)$ in both amplitude quadratures, we displace $(\sqrt{2} \cos(\delta + \pi/4), \pm \sqrt{2} \sin(\delta + \pi/4))$ along $x_A$ and $x_B$, respectively, where $\delta \in [0, \pi]$ is an angle that we can optimize for each pair of values of $r_A$ and $r_B$. This leads to a new Hamiltonian $\hat{H}_{\pm} = (\cos(\delta + \pi/4)\hat{p}_A \pm \sin(\delta + \pi/4)\hat{p}_B)/\sqrt{2}$. Calculating the optimal entanglement witness $E_Q$ of equation (7) we find now

$$E_Q^\delta = E_Q \cos(2\delta). \qquad (15)$$

Therefore, displacing along $x_A = \pm x_B$ ($\delta = 0, \pi$) is indeed the optimal strategy and displacement along any other axis can only degrade the detection of entanglement since $|\cos(2\delta)| \le 1$.

## 5.2. Optical losses

The effect of optical losses can be entirely absorbed by the covariance matrix when it is the same in both modes [45]. The covariance matrix of the input squeezed state $V_0$ is modified in the following way $V_\eta = (1 - \eta)V_0 + \eta I$, where $\eta$ represents the amount of losses. For instance, the probability distribution related to the quantum state given by equation (9) with $\phi = \pi/4$ and $r_A = r_B \equiv r$ is now

$$P_{\eta}(x_A, x_B) \propto e^{-\frac{x_A^2+x_B^2}{2\sigma^2}} (2\eta e^{2r}\sigma^2 + (1-\eta)(x_A+x_B)^2),$$

with $\sigma^2 = (1 - \eta)e^{-2r} + \eta$. A similar but less straightforward result is obtained for general values of $\phi$, $r_A$ and $r_B$.

Figure 8 shows the effect of losses on the detection of entanglement for single photon-subtracted states with $\phi = \pi/4$ and different squeezing parameters: figure 8(a) corresponds to in-phase squeezing $s_A = s_B < 0$ ranging from -1 to -6 dB; in figure 8(b) $s_A$ is positive, while $s_B$ is negative and equals the optimal in-quadrature squeezing—given by $r_B = \log(1/(1 + 2 \sinh(r_A)))^{1/2}$ (see the dashed line in figure 4(b)). For in-quadrature input squeezing, the metrological detection of entanglement is resilient up to ≈5%, whereas for in-phase input squeezing, the metrological detection of entanglement is resilient up to ≈15% for input values of squeezing between -4 and -6 dB. Moreover, entanglement is more resilient to losses in comparison with quantum steering, where the losses threshold is about 7% for the same states [46].

For comparison, we show in figure 9 the effect of losses on the detection of entanglement for a two photon-subtracted state with $\phi = \pi/4$ in two cases: in-phase squeezing $s_A = s_B$ [figure 9(a)], and optimal in-quadrature squeezing [figure 9(b)]. Contrary to the one-photon subtracted case where we found the optimal value of $s_B$ for $\eta = 0$ and used the same value for $\eta > 0$, here we compute for each squeezing parameter $s_A > 0$ and losses $\eta$ the optimal in-quadrature squeezing $s_B^{opt} < 0$. More details on the procedure to obtain these optimal values can be found in the supplemental material. Table 1 gathers the values we found
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Effect of losses $\eta$ on displacement-estimation-based entanglement witness for a single photon-subtracted quantum state with $\phi = \pi/4$ and (a) optimal in-phase input squeezing $s_A = s_B < 0$ (see legend), (b) optimal in-quadrature input squeezing $s_B(s_A)$ (see legend) as given by the black dashed line in figure 4(b). $E > 0$ witnesses entanglement.</figcaption></figure>

<table><caption>**Table 1.** Optimal in-quadrature input squeezing parameters $s_B^{opt}$ for two-photon subtracted states. The third row corresponds to the optimal value $s_B^{opt} = 0$ when $\eta = 0$ (no losses).</caption><thead><tr><th>s<sub>A</sub> (dB)</th><th>1</th><th>1.5</th><th>2</th><th>2.5</th><th>3</th></tr></thead><tbody><tr><th>s<sub>B</sub><sup>opt</sup> (dB)</th><td>-1.4</td><td>-2.6</td><td>-5</td><td>-6.1</td><td>-7.3</td></tr><tr><th>s<sub>B</sub><sup>opt</sup> (dB)</th><td>-0.9</td><td>-1.3</td><td>-1.7</td><td>-2</td><td>-2.3</td></tr></tbody></table>

for the five input squeezing parameters $s_A$ of figure 9(b). As we can notice from this table, the squeezing parameter $s_B^{opt}$ can be very different and always larger in absolute value than the one found for $\eta = 0$ (see third row of table 1); it mainly depends on how the Fisher information is affected by losses. Large anti-squeezing leads to a better resilience of the Fisher information against losses and this explains why the optimal value for $|s_B|$ shifts to larger values. An explanation of the above statement based on the shape of probability distributions from which the Fisher information is computed can be found in the supplemental material.

A similar optimization can also be done for one-photon subtracted states, but it leads to too large squeezing values difficult to reach experimentally with the currently available technology.

As seen from figure 9, the metrological detection of entanglement for in-phase squeezing is resilient up to ≈99% of losses for low input squeezing, and larger than ≈25% for 1.5 dB of squeezing. For in-quadrature squeezing, one can expect a resilience up to 30% of losses for $s_A = 2$ dB. Thus we can conclude that the effect of losses on the metrological detection of entanglement is both parameter and state-dependent, and that two-photon subtracted states are better suited to detect entanglement using displacement estimation in a realistic experiment.

It must be emphasised that our entanglement witness detects only entanglement related to the metrological sensitivity of the state [30]. The losses produce quantum decoherence and impair the metrological power of the quantum state. We have checked that other entanglement witness, such as the logarithmic negativity, detect entanglement in regions where our metrological witness cannot. However, while sensitive to losses, the metrological witness turns out to be more robust than usual Gaussian witnesses for certain states. It is difficult to predict which states lead to better resilience of the metrological witness, though, given the number of parameters to optimize. In some configurations, reducing the level of squeezing might help to mitigate losses as observed in figure 9(a) and as already pointed out in [46] in the case of steering. However, this only reflects the fact that these states are less entangled and thus more noise tolerant, a behavior that translates into a better resilience of the Fisher information for these states.
<figure><img src="image_14.png" /><figcaption>Figure 9. Effect of losses $\eta$ on displacement-estimation-based entanglement witness for a two photon-subtracted quantum state with $\phi = \pi/4$ and in-phase input squeezing $s_A = s_B > 0$ (see legend), (b) optimal in-quadrature input squeezing $s_B(s_A)$ (see supplemental material). $E > 0$ witnesses entanglement.</figcaption></figure>

### 5.3. Discretization of sampled data

The FI and the variances associated to the local displacement generators Var($\hat{p}_{A/B}$) have to be measured experimentally in order to compute the entanglement witness $E$ of equation (4). The variances are directly obtained measuring the phase quadratures with homodyne detection. Estimating the FI experimentally from discrete outcomes, in contrast with the theoretical computation that assumes a continuum of outcomes, relies on the computation of a statistical distance—the Hellinger distance—between a reference probability distribution and the parameter-dependent one [26]. The squared Hellinger distance between a parameter-dependent probability distribution $P(x_A, x_B|\theta)$ and a reference $P(x_A, x_B|0)$ is defined as

$$d_{H,P}^2(\theta) = \frac{1}{2} \iint_R \left( \sqrt{P(x_A, x_B|\theta)} - \sqrt{P(x_A, x_B|0)} \right)^2 dx_A dx_B.$$

The Taylor expansion of the squared Hellinger distance to second order yields [26]

$$d_{H,P}^2(\theta) = \frac{F}{8}\theta^2 + \mathcal{O}(\theta^3),$$

with $F \equiv F(\hat{p}_{AB}, \hat{H})$ the FI. Thus, a quadratic fit is enough to calculate the FI.

However, in an experimental implementation we do not have exact probability distributions $P(x_A, x_B|\theta)$, but relative frequency distributions $\{\mathcal{F}(x_A, x_B|\theta)\}$, which approach the probability distributions for infinitely many independent measurements. In this case, due to statistical fluctuations $\delta\mathcal{F}$, the squared Hellinger distance varies when repeating the measurement. Taking the outcome frequencies from a sample of $M$ experimental realizations, the sample average of the squared Hellinger distance $d_{H,\mathcal{F}}^2(\theta)$ between two relative frequencies $\{\mathcal{F}(x_A, x_B|\theta)\}$ and $\{\mathcal{F}(x_A, x_B|0)\}$ is approximately [26]

$$\langle d_{H,\mathcal{F}}^2(\theta) \rangle = c_0 + \left(\frac{F}{8} + c_2\right) \theta^2 + \mathcal{O}(\theta^3, \delta\mathcal{F}^3), \qquad (16)$$

with $c_0 = (n-1)/4M$, $c_2 \approx F(1+n)/32M$ and $n$ the number of discrete bin pairs $\{x_A, x_B\}$ with measured values. In our case the values of $c_0$ and $c_2$ are obtained through the fitting procedure of simulated data with a quadratic function. Note that $\langle d_{H,\mathcal{F}}^2(\theta) \rangle$ converges asymptotically to $d_{H,P}^2(\theta)$ as $M \to \infty$ and hence the estimation of $F$ is asymptotically unbiased with the bias decreasing as $M^{-1}$.

In the following we study the protocol by simulating homodyne detection with rejection sampling of the theoretical probability distributions obtained from equation (9). We partition the real line corresponding to
<figure><img src="image_15.png" /><figcaption>**Figure 10.** Sampled joint relative frequency distributions $\mathcal{F}(x_A, x_B)$ obtained through rejection sampling of photon-subtracted states. (a) single-photon subtraction with $r_A = r_B = -0.5$ ($|s_{A/B}| = -4.3$ dB) and (b) two-photon subtraction with $r_A = r_B = 0.2$ ($|s_{A/B}| = 1.73$ dB) for $5 \times 10^5$ samples. Bin size $\Delta = 0.2$ (in units of $x_{A/B}$).</figcaption></figure>

the outcomes of the quadrature measured by Alice and Bob in a series of bins with a given bin size $\Delta$. We consider an even number of bins as the mean value of the field is zero for our non-Gaussian probe state. Figure 10 shows two examples of sampled joint relative frequency distributions $\{\mathcal{F}(x_A, x_B)\}$ obtained through rejection sampling of a single photon-subtracted state probability distribution $\mathcal{P}(x_A, x_B)$ given by equation (9) for $r_A = r_B = -0.5$ and of a two photon-subtracted state probability distribution with $r_A = r_B = 0.2$ (see supplementary material). The number of samples is $M = 5 \times 10^5$ and the bin size $\Delta = 0.2$ (in the units of $x_{A/B}$).

We list below the steps to follow in order to calculate the FI:

1. we take the two sets of sampled data corresponding to Alice $\vec{x}_A$ and to Bob $\vec{x}_B$ and split the sampled data ($\vec{x}_A, \vec{x}_B$) of total size $M$ in two equal sets.

2. we bin the data in areas of given size and compute the relative frequencies $\{\mathcal{F}(x_A, x_B|0)\} \approx \mathcal{P}(x_A, x_B|0)$ of the first set that is used as a reference.

3. we displace the data of the second set by an amount $\theta$-the displacement parameter-, bin the data and compute the relative frequencies $\{\mathcal{F}(x_A, x_B|\theta)\} \approx \mathcal{P}(x_A, x_B|\theta)$ that are used as a probe.

4. we calculate the square root of each relative frequency for the reference and the displaced data, take the difference and square it.

5. we calculate the sample average of the squared Hellinger distance $\langle d_{H,\mathcal{F}}^2(\theta) \rangle$ for a value of $\theta$.

6. we repeat this process for different values of $\theta$ and fit the results to a parabola, obtaining the FI with its statistical error through equation (16).

Using this value of FI and the sum of the variances of the phase quadratures we calculate the entanglement witness E through equation (4).

We show in figures 11 and 12 the effect of data discretization and number of samples in the detection of entanglement for lossless and lossy cases. Figure 11 corresponds to the case of a one-photon subtracted state given by equation (9) with $r_A = r_B = -0.5$, either in the lossless case, i.e. $\eta = 0$ [see figure 11(a)], or in the lossy case [$\eta = 0.1$, see figure 11(b)]. Results for a two-photon subtracted state with $r_A = r_B = 0.1$ are depicted in figures 12(a) (lossless case), (b) (lossy case with $\eta = 0.1$), and (c) (lossy case with $\eta = 0.25$). We displace our second data set of size $M/2$ between $\theta \in \{-0.05, 0.05\}$ in steps of $5 \times 10^{-3}$, resulting in 20 data points that we fit with a parabola using equation (16). We partition the outcome quadratures measured by Alice and Bob in a series of bins of size $\Delta$. We perform 100 simulations for each value of bin size and total number of samples to obtain statistical averages and errors. We show the value of entanglement witness E obtained using a continuous probability distribution in solid gray, and the values and errors obtained for
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Effect of bin size $\Delta$, number of total samples $M$ and losses $\eta$ on the entanglement witness $E$ for a single-photon-subtracted quantum state with $\phi = \pi/4$ and $r_A = r_B = -0.5$ ($|S_{AB}| = -4.3$ dB), and (a) $\eta = 0$ or (b) $\eta = 0.1$. Averages and errors are calculated over 100 simulations. $E > 0$ witnesses entanglement. The non-zero horizontal gray dashed lines correspond to the theoretical value.</figcaption></figure>

different bin size $\Delta$ and total samples $M$ in color. The colors represent different number of samples: $M = 10^6$ (blue), $M = 2 \times 10^6$ (orange), $M = 4 \times 10^6$ (green) and $M = 10^7$ (red).

We find that the distance between the computed value from the simulated data and the theoretical value decreases as the bin size shrinks. For large bin size, the number of samples does not affect significantly the accuracy of the measurement. However, for smaller bin size, the accuracy of the discretized estimation raises as the number of samples increases. In general, the statistical error obtained from the fit is lower as the bin size increases. Note that a discretization with an insufficient number of points can lead to an overestimation of the entanglement witness $E$. We find that for both one and two-photon subtracted states in the lossless case the estimation is in good agreement with the theoretical value for $M \ge 2 \times 10^6$ and $\Delta \le 0.1$. In the lossy case, more samples are necessary for the same value of bin size $\Delta$ and overestimation is more significant. To not overestimate the entanglement we should use $M \ge 2 \times 10^6$ and $\Delta > 0.1$. Notably, in both cases we detect entanglement even using a coarse-grained bin size $\Delta = 0.4$ and a relatively low number of samples $M = 10^6$.

In figure 13 we show that the above method is working well for a wide range of two-photon subtracted states with various squeezing parameters. The square with error bars are directly computed from the sampled data with $\Delta = 0.15$ and $M = 2 \times 10^6$ and all agree with the theoretical entanglement witness $E$. It makes it possible to detect non-Gaussian entanglement up to 45% of losses for squeezing parameters up to 1 dB, and up to 30% of losses if we include the 1.5 dB case.

In the above paragraphs we studied cases where the displacement-based metrological criterion works well, but for which usual Gaussian witnesses are also able to detect entanglement for two-photon subtracted states and for certain values of input squeezing parameters. Therefore, in the following, we would like to emphasize a case where Gaussian witnesses do not work, but the metrological criterion does. As discussed in section 3, when we subtract two photons in orthogonal modes, quantum correlations cease to appear at the level of the covariance matrix, making second-order-moment-based entanglement criteria inefficient. Likewise, it turns out that for these type of states the first-order metrological criterion (i.e. using displacement generators) will also fail; thus, it means that we need to consider second order generators. In this case, we can use squeezing operators for which $E_Q \ne 0$ for two-photon subtracted states. The interest of using these generators is two-fold here: (1) the squeezing operation onto the state amounts to squeeze the position quadratures of modes A and B by a factor $e^\theta$ or $e^{-\theta}$ (see section 4.4). Therefore, once the joint probability distribution $P(x_A, x_B|0)$ is measured experimentally, the transformation can be easily implemented in post-processing to obtain $P(x_A, x_B|\theta)$; (2) the FI saturates the QFI when considering the joint probability distributions in the basis $(x_A, x_B)$.
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Effect of bin size $\Delta$, number of total samples *M* and losses $\eta$ on the entanglement witness *E* for a two-photon-subtracted quantum state with $\phi = \pi/4$, $r_A = r_B = 0.1$ ($|s_{A/B}| = 0.87$ dB) and (a) $\eta = 0$, (b) $\eta = 0.1$, (c) $\eta = 0.25$. Averages and errors are calculated over 100 simulations. $E > 0$ witnesses entanglement. The non-zero horizontal gray dashed lines correspond to the theoretical value.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 13.** Displacement-estimation-based entanglement witness *E* for a two-photon subtracted state with $\phi = \pi/4$, $s_A = s_B$, $s_A$ ranging from 0.3 dB to 1.5 dB (see legend). The square with error bars correspond to simulated sampled data calculated over 100 simulations with bin size $\Delta = 0.15$ and total samples $M = 2$ millions. The solid curves correspond to the expected theoretical results.</figcaption></figure>

For all the reasons explained above, we applied the squeezing-estimation-based entanglement on a two photon subtracted states with $r_A = -0.1$ ($s_A = -0.87$ dB) and $r_B = -0.2$ ($s_B = -1.74$ dB), $\phi_A = \pi/4$ and $\phi_B = -\pi/4$, where $\phi_A$ and $\phi_B$ control the modes where we subtract the photons. As shown in figure 14, with a bin size $\Delta = 0.2$ and a number of samples $M = 2 \times 10^6$, there is a good agreement between the theoretical
<figure><img src="image_19.png" /><figcaption>**Figure 14.** Squeezing-estimation-based entanglement witness *E* for a two-photon subtracted state with $\phi_A = -\phi_B = \pi/4$, $r_A = -0.1$ ($s_A = -0.87$ dB), $r_B = -0.2$ ($s_B = -1.74$ dB). The square with error bars correspond to simulated sampled data calculated over 100 simulations with bin size $\Delta = 0.2$ and total samples $M = 2$ millions. The solid curves correspond to the expected theoretical results.</figcaption></figure>

entanglement value *E* and the entanglement estimation based on simulated sampled data calculated over 100 simulations. To obtain these results we also needed an estimation of the local variances of the squeezing generators that can be directly obtained from the statistics of the sampling data using the kurtosis in different quadrature bases (see section 4 of the Supplemental material for more details). One sees that in this case we can detect entanglement up to almost 30% of losses.

## 6. Experimental implementations

Let us discuss possible practical implementations of this protocol. There are few approaches depending on the degree of freedom—or mode—selected to encode the quantum information: path, polarization, frequency and so on. The shared feature of the input modes is that they are independent and excited in squeezed states. An event measured by a single-photon detector fed by a small fraction of power from Alice and Bob's modes where which-mode information is erased, heralds the subtraction of a photon delocalized between the two modes [33]. Two balanced homodyne detectors with a common local oscillator LO retrieve then the joint probability distribution. The sketch of figure 1 is pretty accurate for path-encoded modes where a common beam splitter erases the which-path information.

In the case of spectral modes where the number of modes is usually larger than two—for instance in a multimode frequency-comb Gaussian resource [48, 49]—mode-selective photon-subtraction is accomplished by sum-frequency generation [40]. The detection of an up-converted photon heralds the subtraction of a photon from a multimode input state in one or various spectral modes selected by a pump suitably tailored in frequency. The joint probability distribution of photon-subtracted spectral modes can be retrieved by spectrally-resolved homodyne detection [50]. This approach allows to measure simultaneously the quadratures of the electric field in a number of frequency-band modes. Then, applying a change of basis between the photon-subtracted spectral modes and these frequency-band modes one retrieves the quadrature traces in the modes of interest and hence, the joint probability distribution.

Moreover, we outline that in an experiment, in order to prove that the entanglement results entirely from the non-local photon subtraction, one would use the data from the unconditioned state to test our entanglement witness and demonstrate the independence of the two input squeezed states.

Finally, comparing our simulations with the values measured by Ra *et al* [40] for one-photon subtracted states, where the squeezing of the first and second spectral modes is $s_A = -2.3$ dB and $s_B = -1.7$ dB, respectively, with purities above 90% and detection losses of the order of 12%, and with the values obtained in Takanashi *et al* [51] for two-photon subtracted states, with -3.0 dB of squeezing and detection losses of 15%, we conclude that with a reasonable number of samples ($\approx 10^6$) it is possible to witness non-Gaussian entanglement using exclusively homodyne detection with an experimentally feasible protocol.

## 7. Conclusions and outlook

We proposed a protocol based on Fisher information for witnessing entanglement in an important class of non-Gaussian states: photon-subtracted CV states. The protocol is based on the metrological entanglement criterion proposed in [31], and its strength comes from its simplicity, as it relies solely on homodyne detection. Our approach witnesses entanglement not detected by Gaussian criteria, like for instance Duan *et al* criterion, using the same resources, i.e. quadrature measurements.
We characterized the optimal metrologically-useful entanglement of single- and two-photon-subtracted states analyzing their metrological power in estimation of parameters generated by all single-mode Gaussian gates, namely: displacement, phase shift, shearing and squeezing. We analyzed displacement estimation in details since it gives the largest sensitivity for currently experimentally-relevant values of squeezing ($\le 5$ dB) and it can be applied in postprocessing, thus minimizing the resources necessary in non-Gaussian entanglement characterization and outperforming other protocols where quantum-state tomography is needed.

We demonstrated that our protocol is relevant and experimentally feasible using data from a simulated experiment where the effect of losses, data discretization, and number of samples were taken into account. Our results show that non-Gaussian entanglement can be detected with a feasible number of measurements and data binning. It is well known that losses impair the metrological power of quantum states. However, we found that our metrology-based entanglement detection is resilient up to 30% losses for some purely non-Gaussian entangled states.

The general setup of figure 1 is versatile and can be used to both implement Gaussian entanglement detection protocols based on the covariance matrix and our metrological protocol for non-Gaussian entanglement. For certain classes of states, we believe that this should be sufficient to be able to detect entanglement in any mode basis. However, to determine whether or not a state is passively separable, as would be required for the sampling protocols in [13], one would still need to certify the presence of entanglement in every possible mode basis. While our work certainly offers us a useful experimental tool, we also hope that it will be a step towards finding new techniques that allow us to certify entanglement in every possible mode basis. After all, non-Gaussian entangled states encompass a huge state space and we have just started to scratch its surface. In order to gain insight about general features of this exotic quantum feature, in future work we will analyze entangled states based on other non-Gaussian operations and connect our entanglement criterion with others based on higher-order covariance matrices [52, 53].

## Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files).

## Acknowledgments

This work received funding from the ANR JCJC Project NoRdiC (ANR-21-CE47-0005), the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 899587, and the QuantERA II Project SPARQL that has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 101017733. This work was also funded by MCIN/AEI/10.13039/501100011033 and the European Union 'NextGenerationEU' PRTR fund [RYC2021-031094-I], by the Ministry of Economic Affairs and Digital Transformation of the Spanish Government through the QUANTUM ENIA Project call—QUANTUM SPAIN Project, by the European Union through the Recovery, Transformation and Resilience Plan—NextGenerationEU within the framework of the Digital Spain 2026 Agenda, and by the CSIC Interdisciplinary Thematic Platform (PTI+) on Quantum Technologies (PTI-QTEP+). This work was supported in part by the Valencian Government grant with Reference Number CIAICO/2021/184. It was carried out during the tenure of an ERCIM 'Alain Bensoussan' Fellowship Programme.

## ORCID iDs

David Barral https://orcid.org/0000-0001-6658-4684

Giacomo Sorelli https://orcid.org/0000-0003-2197-2768

Manuel Gessner https://orcid.org/0000-0003-4203-0366

Mattia Walschaers https://orcid.org/0000-0003-4649-0044

## References

[1] Einstein A, Podolsky B and Rosen N 1935 Can quantum-mechanical description of physical reality be considered complete? *Phys. Rev.* **47** 777

[2] Acín A *et al* 2018 The quantum technologies roadmap: a European community view *New J. Phys.* **20** 080201

[3] Horodecki R, Horodecki P, Horodecki M and Horodecki K 2009 Quantum entanglement *Rev. Mod. Phys.* **81** 865

[4] Braunstein S L and van Loock P 2005 Quantum information with continuous variables *Rev. Mod. Phys.* **77** 513

[5] Larsen M V, Guo X, Breum C R, Neergaard-Nielsen J S and Andersen U L 2019 Deterministic generation of a two-dimensional cluster state *Science* **366** 369–72

[6] Asavanant *W et al* 2019 Generation of time-domain-multiplexed two-dimensional cluster state *Science* **366** 373–6
[7] Duan L-M, Giedke G, Cirac J I and Zoller P 2000 Inseparability criterion for continuous variable systems Phys. Rev. Lett. 84 2722

[8] Simon R 2000 Peres-Horodecki separability criterion for continuous variable systems Phys. Rev. Lett. 84 2726

[9] Giovannetti V, Mancini S, Vitali D and Tombesi P 2003 Characterizing the entanglement of bipartite quantum systems Phys. Rev. A 67 022320

[10] van Loock P and Furusawa A 2003 Detecting genuine multipartite continuous-variable entanglement Phys. Rev. A 67 052315

[11] Abiuso P, Bauml S, Cavalcanti D and Acín A 2021 Measurement-device-independent entanglement detection for continuous-variable systems Phys. Rev. Lett. 126 190502

[12] Walschaers M, Fabre C, Parigi V and Treps N 2017 Entanglement and Wigner function negativity of multimode non-Gaussian states Phys. Rev. Lett 119 183601

[13] Chabaud U and Walschaers M 2023 Resources for bosonic quantum computational advantage Phys. Rev. Lett. 130 090602

[14] Namiki R 2012 Photonic families of non-Gaussian entangled states and entanglement criteria for continuous-variable systems Phys. Rev. A 85 062307

[15] Walschaers M 2021 Non-Gaussian states and where to find them PRX Quantum 2 030204

[16] Schukin E and Vogel W 2005 Inseparability criteria for continuous bipartite quantum states Phys. Rev. Lett. 95 230502

[17] Miranowicz A and Piani M 2006 Comment on “Inseparability criteria for continuous bipartite quantum states” Phys. Rev. Lett. 97 058901

[18] Schukin E and van Loock P 2016 Higher-order Einstein-Podolsky-Rosen correlations and inseparability conditions for continuous variables Phys. Rev. A 93 032114

[19] Agarwal G S and Biswas A 2005 Inseparability inequalities for higher order moments for bipartite systems New J. Phys. 7 211

[20] Hillery M and Zubairy M S 2006 Entanglement conditions for two-mode states Phys. Rev. Lett. 96 050503

[21] Shen Y, Assad S M, Grosse N B, Li X Y, Reid M D and Lam P K 2015 Nonlinear entanglement and its application to generating cat states Phys. Rev. Lett. 114 100403

[22] Kogias I, Skrzypczyk P, Cavalcanti D, Acín A and Adesso G 2015 Hierarchy of steering criteria based on moments for all bipartite quantum systems Phys. Rev. Lett. 115 210401

[23] Tian M, Zou Z, Zhang D, Barral D, Bencheikh K, He Q, Sun F-X and Xiang Y 2024 Quantum Sci. Technol. 9 015021

[24] Walborn S P, Taketani B G, Salles A, Toscano F and de Matos Filho R L 2009 Entropic entanglement criteria for continuous variables Phys. Rev. Lett 103 160505

[25] Nha H, Lee S-Y, Ji S-W and Kim M S 2012 Efficient entanglement criteria beyond Gaussian limits using Gaussian measurements Phys. Rev. Lett. 108 030503

[26] Strobel H, Muessel W, Linnemann D, Zibold T, Hume D B, Pezze L, Smerzi A and Oberthaler M K 2014 Fisher information and entanglement of non-Gaussian spin states Science 345 424

[27] Hu L, Al-amri M, Liao Z and Zubairy M S 2020 Continuous-variable quantum key distribution with non-Gaussian operations Phys. Rev. A 102 012608

[28] Opatrny T, Kurizki G and Welsch D-G 2000 Improvement on teleportation of continuous variables by photon subtraction via conditional measurement Phys. Rev. A 61 032302

[29] Filippov S and Ziman M 2014 Entanglement sensitivity to signal attenuation and amplification Phys. Rev. A 90 010301(R)

[30] Pezze L and Smerzi A 2009 Entanglement, nonlinear dynamics and the Heisenberg limit Phys. Rev. Lett. 102 100401

[31] Gessner M, Pezze L and Smerzi A 2016 Efficient entanglement criteria for discrete, continuous and hybrid variables Phys. Rev. A 94 020101(R)

[32] Gessner M, Pezze L and Smerzi A 2017 Entanglement and squeezing in continuous-variable systems Quantum 1 17

[33] Ourjoumtsev A, Ferreyrol F, Tualle-Brouri R and Grangier P 2009 Preparation of non-local superpositions of quasi-classical light states Nat. Phys. 5 189–92

[34] Hyllus P and Eisert J 2006 Optimal entanglement witnesses for continuous-variable systems New J. Phys. 8 51

[35] Braunstein S L and Caves C M 1994 Statistical distance and the geometry of quantum states Phys. Rev. Lett. 72 3439

[36] Nieto M M 1997 Displaced and squeezed number states Phys. Lett. A 229 135–43

[37] We can use the following expressions for operators acting on a wavefunction Ψ(x) [36]

$$ \exp[\theta \partial_x] \Psi(x) = \Psi(x - \theta), $$

$$ \exp[\theta(x\partial_x)]\Psi(x) = \Psi(xe^{\theta}), $$

$$ \exp[\theta(\partial_x)^2]\Psi(x) = \frac{1}{\sqrt{4\pi\theta}} \int_{-\infty}^{\infty} \exp[-\frac{(y-x)^2}{4\theta}] \Psi(y)dy, $$

where ∂_x is related to the phase quadrature ˆp through the functional relation ˆp = −2i∂_x.

[38] Rivas A and Luis A 2010 Precision quantum metrology and nonclassicality in linear and nonlinear detection schemes Phys. Rev. Lett. 105 010403

[39] Walschaers M and Treps N 2020 Remote generation of Wigner negativity through Einstein-Podolsky-Rosen steering Phys. Rev. Lett. 124 150501

[40] Ra Y S, Dufour A, Walschaers M, Jacquard C, Michel T, Fabre C and Treps N 2020 Non-Gaussian quantum states of a multimode light field Nat. Phys. 16 144–7

[41] Agarwal G S 2013 Quantum Optics (Cambridge University Press)

[42] Gu M, Weedbrook C, Menicucci N C, Ralph T C and van Loock P 2009 Quantum computing with continuous variables Phys. Rev. A 79 062318

[43] Tian M, Xiang Y, Sun F-X, Fadel M and He Q 2022 Characterizing multipartite non-Gaussian entanglement for a three-mode spontaneous parametric down-conversion process Phys. Rev. Appl. 18 024065

[44] Chabaud U, Roeland G, Walschaers M, Grosshans F, Parigi V, Markham D and Treps N 2021 Certification of non-Gaussian States with Operational Measurements PRX Quantum 2 020333

[45] Walschaers M, Ra Y-S and Treps N 2019 Mode-dependent-loss model for multimode photon-subtracted states Phys. Rev. A 100 023828

[46] Lopetegui C E, Gessner M, Fadel M, Treps N and Walschaers M 2022 Homodyne detection of non-Gaussian quantum steering PRX Quantum 3 030347

[47] Qin Z, Gessner M, Ren Z, Deng X, Han D, Li W, Su X, Smerzi A and Peng K 2019 Characterizing the multipartite continuous-variable entanglement structure from squeezing coefficients and the fisher information npj Quantum Inf. 5 3

[48] Roslund J, de Araujo R M, Jiang S, Fabre C and Treps N 2014 Wavelength-multiplexed quantum networks with ultrafast frequency combs Nat. Photon. 8 109–12
[49] Cai Y, Roslund J, Ferrini G, Arzani F, Xu X, Fabre C and Treps N 2017 Multimode entanglement in reconfigurable graph states using optical frequency combs *Nat. Commun.* **8** 15645

[50] Ansquer M, Thiel V, De S, Argence B, Gredat G, Bretenaker F and Treps N 2021 Unveiling the dynamics of optical frequency combs from phase-amplitude correlations *Phys. Rev. Res.* **3** 033092

[51] Takanashi H, Wakui K, Suzuki S, Takeoka M, Hayasaka K, Furusawa A and Sasaki M 2008 Generation of large-amplitude coherent state superposition via ancilla-assisted photon subtraction *Phys. Rev. Lett.* **101** 233605

[52] Zhang D, Barral D, Cai Y, Zhang Y, Xiao M and Bencheikh K 2021 Hierarchy of nonlinear entanglement dynamics for continuous variables *Phys. Rev. Lett.* **127** 150502

[53] Zhang D, Barral D, Zhang Y, Xiao M and Bencheikh K 2023 Genuine tripartite non-Gaussian entanglement *Phys. Rev. Lett.* **130** 093602