# Detection of mode-intrinsic quantum entanglement

Carlos E Lopetegui, Mathieu Isoard, Nicolas Treps, Mattia Walschaers

► To cite this version:

Carlos E Lopetegui, Mathieu Isoard, Nicolas Treps, Mattia Walschaers. Detection of mode-intrinsic quantum entanglement. Optica Quantum, 2025, 3 (4), pp.312-328. 10.1364/OPTICAQ.560594. hal-04668245

HAL Id: hal-04668245

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

Submitted on 30 Jul 2025

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Detection of mode-intrinsic quantum entanglement

CARLOS E. LOPETEGUI,\* ID MATHIEU ISOARD, NICOLAS TREPS, ID AND MATTIA WALSCHAERS

Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-Université PSL, Collège de France, 4 place Jussieu, F-75252 Paris, France
\*carlos-ernesto.lopetegui-gonzalez@lkb.upmc.fr

Received 3 March 2025; accepted 8 May 2025; published 27 June 2025

Quantum correlations are at the core of quantum information processing and are necessary to achieve a quantum computational advantage. In the context of continuous-variable quantum systems, another necessary resource for quantum advantages is non-Gaussianity. In this work, we propose a witness based on previously known relations between metrological power and quantum correlations, to detect a strong form of entanglement that only non-Gaussian states possess and that cannot be undone by passive optical operations, i.e., entanglement in all mode bases. The strength of our witness is two-fold: it only requires measurements in one basis to check entanglement in any arbitrary mode basis; it can be made applicable experimentally using homodyne measurements and without requiring a full tomography of the state.

Published by Optica Publishing Group under the terms of the Creative Commons Attribution 4.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

https://doi.org/10.1364/OPTICAQ.560594

## 1. INTRODUCTION

Major efforts are being made worldwide to harness the power of quantum mechanics for technological applications. The use of these systems for performing computations beyond the capabilities of classical processors is one of the most encouraging perspectives. A wide variety of physical platforms are considered for this purpose [1–8], among which optical systems have proven to be a promising framework given their potential for scalability [9,10].

The full power of optical quantum technologies is attained when the collective state of many modes can be coherently controlled. These collective states manifest entanglement, a unique feature of quantum systems that has challenged the common conceptions of physicists over the last century [11–13]. Entanglement, and quantum correlations in general, reflect in measurement statistics that are inconsistent with classical correlations. In the context of quantum information processing and quantum computing, entanglement has been revealed to be an essential resource [14,15]. However, in continuous variables (CV) quantum computing, entanglement does not suffice to go beyond the reach of classical simulations if the state is fully Gaussian [16], or if its associated Wigner function – which describes all the statistical properties of the quantum state – is positive. Indeed, it has been shown that negativity in the Wigner function is a necessary resource for realizing sampling protocols that cannot be efficiently simulated with classical resources [17]. This result significantly boosted the already ongoing efforts aimed at producing non-Gaussian states [18–22], with remarkable results achieved using microwave cavities [23], superconducting circuits [24], or motional states of ions [25,26] and mechanical oscillators [27].

However, recent results have pointed out that not every kind of quantum entanglement (even if the state is Wigner negative)

is sufficient to render computations with CV quantum systems difficult to simulate with a classical device. In particular, in Ref. [28], it was shown that in a certain family of sampling protocols, a strong form of quantum entanglement is required: not passive separability, i.e., entanglement that cannot be undone with optical passive transformations (beam splitters and phase shifters). This kind of entanglement was previously discussed in Refs. [29–31], from a different perspective. In Ref. [29], they refer to it as mode-independent quantum entanglement and point out its potential usefulness for communication protocols. We refer to this kind of quantum correlations as *mode-intrinsic entanglement* [32]. Finding and characterizing non-Gaussian states that exhibit this refined quantum feature and that can be produced in the lab are thus important steps in the domain of CV quantum information.

Many witnesses and criteria have been derived for detecting entanglement in Gaussian states [33–37]. All of them rely in some way on the symplectic structure of the covariance matrix. Once measured in one basis, such a matrix can be obtained in any mode basis by means of an orthogonal transformation. For this reason, these methods are quite versatile for studying entanglement in every mode basis, which is desirable for the study of mode-intrinsic entanglement. However, Williamson’s decomposition [38,39] implies that there is always a basis on which the covariance matrix is consistent with a separable Gaussian state, regardless of the nature of the quantum state itself. Therefore, these methods are not suitable for demonstrating mode-intrinsic entanglement. Furthermore, it shows that mode-intrinsic entanglement is inherently a non-Gaussian feature.

Over the years, many protocols have been proposed to detect entanglement in non-Gaussian states. If one has access to the density matrix of a given state in the Fock space, it is possible to
rely on techniques from discrete variables settings, like the positive partial transpose (PPT) criterion [40], to detect quantum correlations. However, the reconstruction of the density matrix requires a full tomography of the state [41], a very demanding and time-consuming experimental task. Efforts to circumvent this issue for non-Gaussian states have lead to the introduction of different methods to witness quantum correlations [42–50].

Among those, there is a subset [46,47] that rely on connections to quantum metrology [51–53]. These witnesses are based on the Fisher information [54], a quantity that reflects the highest sensitivity with which we can estimate a parameter of interest using measurements of an observable. These methods rely on the idea that quantum correlations can lead to metrological sensitivities that could not be achieved with only classical correlations. CV-specialized protocols based on these witnesses [55,56] have been shown to succeed in detecting quantum correlations in relevant non-Gaussian states where usual Gaussian criteria fail. This suggests that non-Gaussian entanglement is not only crucial in the framework of quantum computation, but can also potentially offer advantages in sensing applications. This idea relates naturally to recent results [57] that link the negativity of some quasiprobability distributions to enhancements in quantum phase estimation in quantum optics.

All the methods we have mentioned so far to detect entanglement in non-Gaussian states are developed to detect entanglement once we have fixed a mode basis of interest. However, to study entanglement in every mode basis, a method that retains the versatility of the Gaussian protocols, but potentially succeeds where the former fails, is necessary. In this paper, we propose such a protocol based on metrological witnesses. We show that our protocol works for relevant non-Gaussian states that can be produced with state-of-the-art experiments, namely photon subtracted states [18,19]. Likewise, our witness demonstrates effectiveness in multimode systems such as (non-Gaussian) cluster states – that can be obtained by degaussifying readily available large CV cluster states [9,10,18], which have been widely studied for their scalability in the pursuit of large-scale quantum computing. While sensitive to losses, the mode-intrinsic-entanglement witness that we suggest here can be applied experimentally through homodyne detection and the estimation of the corresponding classical Fisher information – using the Hellinger distance as proposed in Ref. [56]; it thus avoids full tomography of the state, making metrological witnesses a very efficient tool in the detection of mode-intrinsic entanglement and non-Gaussian features. Moreover, it proves useful for the less specific task of studying the entanglement properties of a given state beyond the basis in which measurements are performed, with higher accuracy than Gaussian methods.

The paper is organized as follows. In Sec. 2, we introduce the metrological entanglement witness derived in Ref. [46] and extend it to the detection of mode-intrinsic entanglement in Sec. 3. Then, in Sec. 4, we apply our witness to experimentally relevant multimode states, namely photon subtracted states. In the last section (Sec. 5), we consider a lower bound of the mode-intrinsic entanglement witness derived in the previous sections; we show that this witness is experimentally accessible through homodyne measurements; in particular, we test this lower bound witness on a two-photon-subtracted state using simulated experimental data.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Schematic of the metrological protocol in which the entanglement witness on which we present is based. Alice and Bob share a state $\hat{\rho}$ and try to estimate a parameter $\kappa$ encoded in a unitary evolution under the Hamiltonian $\hat{H}_A + \hat{H}_B$, given by a sum of local Hamiltonians, which implies that the evolution can actually be factorized. To perform the estimation, they can perform any measurement associated with a positive operator valued measure $\Pi_{AB}$, with elements $\{\hat{\Pi}_{M_A, M_B}\}$ whose associated outcomes are $(M_A, M_B)$. The line joining the measurement devices is made to highlight that they can perform any set of non-local compatible measurements.</figcaption></figure>

## 2. METROLOGY-BASED ENTANGLEMENT WITNESSES

In this section, we introduce metrological witnesses of entanglement on which the protocol that we present is based [46]. These witnesses rely on the analysis of a specific metrological task. We consider two parties, Alice and Bob, who share a state $\hat{\rho}$.

They are interested in jointly estimating a parameter $\kappa$ that they encode through a unitary transformation induced by the sum of local generators $\hat{H} = \hat{H}_A + \hat{H}_B$ (see Fig. 1). In this case, the density matrix of the state transforms as $\hat{\rho}_{\kappa} = \exp(-i\kappa\hat{H})\hat{\rho}\exp(i\kappa\hat{H})$. The sensitivity with which they can estimate the parameter, given that they can measure some (possibly non-local) observables $\hat{M}_{A(B)}$, is upper bounded by the Fisher information associated to the probability distribution of outcomes of these observables. The Fisher information is given by

$$F(p(M_A, M_B|\kappa)) = \int_{M_A, M_B} dM_A dM_B \left( \frac{\partial \mathcal{L}(p(M_A, M_B|\kappa))}{\partial \kappa} \right)^2 \times p(M_A, M_B|\kappa), \quad (1)$$

where $\mathcal{L}(p(M_A, M_B|\kappa)) = \log(p(M_A, M_B|\kappa))$ stands for the logarithmic likelihood, and $p(M_A, M_B|\kappa) = \text{Tr}(\hat{\rho}_{\kappa} \hat{\Pi}_{M_A, M_B})$ is the joint probability distribution conditioned by the parameter $\kappa$, of outcomes $(M_A, M_B)$.

The ultimate sensitivity achievable in the estimation of $\kappa$ is given by the quantum Fisher information (QFI) of the state, defined as the maximum Fisher information, optimized over all possible measurements, i.e.,

$$F_Q(\hat{\rho}, \hat{H}) = \max_{\Pi_{AB}} F \left[ \text{Tr} \left( \hat{\rho}_{\kappa} \hat{\Pi}_{M_A, M_B} \right) \right]. \quad (2)$$

It can be shown that for pure states,

$$F_Q(|\psi\rangle\langle\psi|, \hat{H}) = 4\text{Var}(|\psi\rangle\langle\psi|, \hat{H}), \quad (3)$$
where $\text{Var}(|\psi\rangle\langle\psi|, \hat{H})$ is the variance of the generator $\hat{H}$ in the state $|\psi\rangle$. In the case of mixed states, there is also an analytical expression for the quantum Fisher information, but for the moment, it suffices to say that $F_Q(\hat{\rho}, \hat{H}) \le 4\text{Var}(\hat{\rho}, \hat{H})$. Another useful property for our purpose is the additivity of the QFI, namely,

$$F_Q(\hat{\rho}_A \otimes \hat{\rho}_B, \hat{H}_A + \hat{H}_B) = F_Q(\hat{\rho}_A, \hat{H}_A) + F_Q(\hat{\rho}_B, \hat{H}_B), \quad (4)$$

where $\hat{\rho}_{A(B)} = \text{Tr}_{B(A)}(\hat{\rho})$ is the reduced density matrix. Combining these two properties and the convexity of the quantum Fisher information, one can show that if Alice and Bob share a separable state $\hat{\rho}$, it leads to the inequality [46]

$$F_Q(\hat{\rho}, \hat{H} = \hat{H}_A + \hat{H}_B) \le 4 (\text{Var}(\hat{\rho}_A, \hat{H}_A) + \text{Var}(\hat{\rho}_B, \hat{H}_B)). \quad (5)$$

This inequality can only be violated by entangled states and can be thus used to witness it. In practical terms, this implies that if they share an entangled state, they better join forces than trying to estimate independently the parameter. One would be tempted to say that this implies that entanglement offers a metrological advantage in this setup; yet the situation is a bit more subtle than that, because we are comparing the estimation they can do jointly with respect to what they can estimate with their noisy reduced states. However, this does not imply that there is no better pure single mode state they could independently use and each get better Fisher information.

The full power of this kind of inequalities can be exploited when we optimize over a set of allowed choices of generators. For any $\hat{H}_{A(B)}$, the violation of the metrological inequality provides a sufficient condition for entanglement. Yet, if a violation is not attained, it might easily be that there is a different choice of generators that produces one. If Alice and Bob each have access to a set of local generators $\mathcal{A}_{A(B)} = \{\hat{H}_{A(B)}^{(1)}, \dots, \hat{H}_{A(B)}^{(N_{A(B)})}\}$, where $N_{A(B)}$ is the number of generators in each set, they can choose to implement generators given by $\hat{H}_{A(B)} = \mathbf{c}_{A(B)}^T \mathcal{A}_{A(B)}$. We group $\mathbf{c}_{A(B)}$ in a single vector $\mathbf{c} = \{\mathbf{c}_A, \mathbf{c}_B\}$, and similarly $\mathcal{A} = \{\mathcal{A}_A, \mathcal{A}_B\}$. This allows to write the quantum Fisher information in a matrix form as [46]

$$F_Q(\hat{\rho}, \mathbf{c}^T \mathcal{A}) = \mathbf{c}^T \hat{Q}_{\hat{\rho}}^{\mathcal{A}} \mathbf{c}, \quad (6)$$

where $\hat{Q}_{\hat{\rho}}^{\mathcal{A}}$ is the quantum Fisher information matrix and its elements are given, provided the spectral decomposition $\hat{\rho} = \sum_k p_k |\psi_k\rangle\langle\psi_k|$, by

$$(\hat{Q}_{\hat{\rho}}^{\mathcal{A}})_{i,j}^{(m,n)} = \sum_{k,l} \frac{(p_k - p_l)^2}{p_k + p_l} \langle\psi_k|\hat{H}_i^{(m)}|\psi_l\rangle \langle\psi_l|\hat{H}_j^{(n)}|\psi_k\rangle, \quad (7)$$

where the sum extends over all pairs $k, l$ such that $p_k + p_l \ne 0$. The state $\hat{\rho}$ is typically considered to be the one after the implementation of the evolution. Yet here, we consider only small values of the parameters, so that the spectral decomposition considered is that of the probe state, and all higher orders on the value of the parameters are omitted.

The right-hand side of Eq. (5), however, corresponds to the variance of the generator in the product state of the local reduced states, i.e., $\Pi(\hat{\rho}) = \hat{\rho}_A \otimes \hat{\rho}_B$. This term can be expressed also in a matrix form by means of the covariance matrix of the generators defined by

$$\Gamma_{i,j}^{(m,n)} = \text{Cov}(\hat{H}_i^{(m)}, \hat{H}_j^{(n)})_{\Pi(\hat{\rho})}. \quad (8)$$

The variance of the generator for a particular choice of $\mathbf{c}$ is given by

$$\text{Var}(\mathbf{c}_A^T \mathcal{A}_A) + \text{Var}(\mathbf{c}_B^T \mathcal{A}_B) = \mathbf{c}^T \hat{\Gamma}_{\Pi(\hat{\rho})}^{\mathcal{A}} \mathbf{c}. \quad (9)$$

Combining Eqs. (6) and (9) into Eq. (5), we get that, for any separable state,

$$\mathbf{c}^T (\hat{Q}_{\hat{\rho}}^{\mathcal{A}} - 4 \Gamma_{\Pi(\hat{\rho})}^{\mathcal{A}}) \mathbf{c} \le 0, \quad \forall \mathbf{c}. \quad (10)$$

Given that this condition is satisfied for all choices of $\mathbf{c}$, it can be rewritten as

$$\lambda_{\max} \left( \hat{Q}_{\hat{\rho}}^{\mathcal{A}} - 4 \Gamma_{\Pi(\hat{\rho})}^{\mathcal{A}} \right) \le 0, \quad (11)$$

where $\lambda_{\max}(M)$ stands for the maximum eigenvalue of $M$. This inequality provides us a witness of entanglement

$$E_{\hat{\rho}}^{\mathcal{A}} \equiv \lambda_{\max} \left( \hat{Q}_{\hat{\rho}}^{\mathcal{A}} - 4 \Gamma_{\Pi(\hat{\rho})}^{\mathcal{A}} \right). \quad (12)$$

If this witness takes positive values, we can certify the presence of entanglement. If the state is pure and each set $\mathcal{A}_{A(B)}$ forms a complete set of observables that spans the full local Hilbert spaces, the condition in Eq. (11) becomes not just a sufficient but also a necessary condition [46].

All the above results can be straightforwardly generalized to multimode systems and to multipartite entanglement detection. In this case, the $m$ parties have access to a set of local generators $\mathcal{A}_m = \{\hat{H}_m^{(1)}, \dots, \hat{H}_m^{(N_m)}\}$. The elements of the quantum Fisher information matrix and the covariance matrix can be then computed from Eqs. (7) and (8), with $\Pi(\hat{\rho}) = \hat{\rho}_{p_1} \otimes \dots \otimes \hat{\rho}_{p_p}$, where we split the $m$-multimode system into $p$ partitions $\mathcal{P}_i$ to detect $p$-partite entanglement.

# 3. WITNESS OF MODE-INTRINSIC ENTANGLEMENT

## 3.1. Multi-Mode Quantum Optics

We start this section by providing a quick introduction to multimode quantum optics [58]. In quantum optics, modes are the normalized solutions of Maxwell equations. A mode basis is a set of orthonormal solutions to Maxwell equations. Classically, this implies that we can write the electromagnetic field as weighted superpositions of the elements of a given mode basis. In quantum optics, we write the complex operator corresponding to the electromagnetic field as

$$\hat{E}^{(+)}(\mathbf{r}, t) = \sum_k \epsilon \hat{a}_k \mathbf{u}_k(\mathbf{r}, t), \quad (13)$$

where $\mathbf{u}_k(\mathbf{r}, t)$ are all the modes in the basis of interest, $\epsilon$ is a scaling term (independent of $k$ within the narrow band approximation), and $\hat{a}_k$ are the annihilation operators associated to modes $k$. Notice that the quantization step does not modify the mode structure but rather the field fluctuations associated with each of them, which is encoded through the statistics associated with the creation and annihilation operators of the field. Those being complex operators, they are written as

$$\hat{a}_k = \frac{\hat{q}_k + i\hat{p}_k}{2}, \quad \hat{a}_k^\dagger = \frac{\hat{q}_k - i\hat{p}_k}{2}, \quad \forall k \in \{1, \dots, m\}, \quad (14)$$

where $\hat{a}_k^\dagger$ is the creation operator corresponding to mode $k$. Operators $\hat{q}_k$ and $\hat{p}_k$ are the amplitude and phase quadratures of the field, and they are the fundamental observables in CV quantum optics. Quadrature operators satisfy the canonical commutation relation $[\hat{q}_k, \hat{p}_l] = 2i \delta_{k,l}$, just like the position and momentum operators of a quantum harmonic oscillator.
Different mode bases are related to each other by unitary transformations U, i.e., $\mathbf{v}_k(\mathbf{r}, t) = \sum_{l=1}^{M} U_{kl} \mathbf{u}_l(\mathbf{r}, t)$. This implies that the creation and annihilation operators in different bases are related by a similar transformation, so that the annihilation operators $\hat{b}_k$ associated with modes $\mathbf{v}_k(\mathbf{r}, t)$ are given by

$$ \hat{b}_k = \sum_l U_{kl}^* \hat{a}_l, \qquad (15) $$

which can be rewritten in matrix form as $\vec{\hat{b}} = U^* \vec{\hat{a}}$, where $\vec{\hat{a}}$ is the column vector of annihilation operators associated with the mode basis $\{\mathbf{u}_k(\mathbf{r}, t)\}$. These unitary transformations correspond to orthogonal real transformation at the level of the quadratures of the field. If we group the quadrature operators as $\vec{\mathbf{x}}^{(u)} = \{\hat{x}_1, \dots, \hat{x}_M, \hat{p}_1, \dots, \hat{p}_M\}$, and similarly for the basis $\mathbf{v}$, they are related to each other by

$$ \vec{\mathbf{x}}^{(v)} = \mathbf{O} \vec{\mathbf{x}}^{(u)}, \qquad (16) $$

where

$$ \mathbf{O} = \begin{pmatrix} \mathrm{Re}(U) & \mathrm{Im}(U) \\ -\mathrm{Im}(U) & \mathrm{Re}(U) \end{pmatrix}. $$

A mode basis change thus corresponds to orthogonal transformations.

## 3.2. Mode-Intrinsic Entanglement

A key aspect of non-Gaussian entanglement lies in its robustness against mode basis change. In contrast to Gaussian states whose entanglement can always be destroyed through a passive linear optics transformation, for some non-Gaussian states, no such transformation exists; these states are said to be *not passively separable* and we call this type of entanglement *mode-intrinsic entanglement*.

Note that not all non-Gaussian states that manifest entanglement are mode intrinsic entangled. In particular, if it is possible to obtain them from a separable non-Gaussian state by applying a passive transformation, then they can be disentangled by applying the inverse (passive) transformation. For example, any state in which a photon is shared (split) between two modes, i.e., $|\psi\rangle = \cos\theta|01\rangle + \sin\theta|10\rangle$, is non-Gaussian, and is entangled for every $\theta \in (0, \pi/2)$. Nevertheless, it is passive separable, as it suffices to apply a beam splitter (mode basis change), with angle $-\theta$, to undo the entanglement and bring a state to the form $|01\rangle$. A similar example is provided by the Hong-Ou-Mandel state $|^{(02)-[20]}_\sqrt{2}$, which is obtained by applying a balanced beam splitter on two separable photons $|11\rangle$. However, there are well-known examples of mode-intrinsic entangled states, i.e., states that cannot be constructed from a separable state by applying passive linear optics, for example, $|\psi\rangle = |^{(02)+[20]+[11]}_\sqrt{3}$, any NOON state $|\mathrm{NOON}\rangle = |^{(0N)+[N0]}_\sqrt{2}$, for $N \ge 3$, and the two-photon Bell state $|\psi\rangle = |^{[00]+[11]}_\sqrt{2}$.

To detect such entanglement, we propose to use the metrological witness Eq. (11) introduced in Sec. 2 as follows: a state will be mode-intrinsic entangled if in every mode basis, $E_\rho^A > 0$, which is equivalent to the nonexistence of a passive linear optics operation that disentangles the state. By analyzing the transformation properties of the Fisher information matrix and the covariance matrix under mode basis change, we show that the check for mode-intrinsic entanglement does not require the experimental reconstruction of those quantities in each mode basis. Thus, we propose a method to do this check without the burden of

measuring in each mode basis, but rather post-processing the measurements performed in a single one.

Quadrature operators and their higher order combinations are a judicious choice of generators for evaluating the witness, as their statistics can be recovered experimentally through homodyne measurements [59]. Moreover, the effect of the unitaries they implement, up to order two, are either displacements in the quadratures themselves, or Gaussian operations, which can either be implemented in post-processing, or lie within the reach of experimental implementations in the near term.

To that aim, we introduce the following generators:

$$ \hat{H}_i^S = S \left( \hat{q}_1^{(i)} \dots \hat{q}_m^{(i)} \hat{p}_1^{(i+1)} \dots \hat{p}_m^{(i)} \right), \quad \text{with } \sum_{j=1}^{2m} k_j^{(i)} \le N, \qquad (17) $$

where $S(\hat{H})$ is the symmetrized version of $\hat{H}$ such that $\hat{H}_i^S$ is Hermitian. The generators $\hat{H}_i^S$ contain products of $N$ quadrature operators or less, and are entirely described by the sequences $(k_1^{(i)}, \dots, k_{2m}^{(i)})$, where each $k_j^{(i)}$ corresponds to the power of the associated quadrature operator in the decomposition Eq. (17). These generators can be either *local* or *non-local*. Here, locality is defined as follows: *local* generators are composed of quadratures operators belonging to only one subsystem, while *non-local* generators involve quadrature operators of multiple subsystems. For instance, the quadrature operator $\hat{p}_1^2$ is a local generator, while $\hat{p}_1 \hat{p}_2$ is a non-local generator shared between subsystems 1 and 2.

To evaluate the witness Eq. (12), in a specific mode basis, in which the quantum state of our system is $\hat{\rho}$, we need to evaluate the quantum Fisher information matrix $\mathcal{Q}_{\hat{\rho}}^{\mathcal{A}_{N,m}}$ associated with the given set of local generators, using Eq. (7). We denote by $\mathcal{A}_{N,m}^{(L)}$ the set of local generators of the form in Eq. (17), up to order $N$, in an $m$-mode system. For instance, at order 2, for 2 modes, this set is given by

$$ \mathcal{A}_{2,2}^{L} = \left\{ \hat{q}_1, \hat{p}_1, \hat{q}_2, \hat{p}_2, \hat{q}_1^2, \frac{1}{2}(\hat{q}_1 \hat{p}_1 + \hat{p}_1 \hat{q}_1), \hat{p}_1^2, \hat{q}_2^2, \frac{1}{2}(\hat{q}_2 \hat{p}_2 + \hat{p}_2 \hat{q}_2) \right\}. $$

To further proceed with the witnessing of mode-intrinsic entanglement, a similar evaluation would be required for each $\hat{O}\hat{\rho}\hat{O}^\dagger$. Yet, by looking at Eq. (7), we can see that going to the Heisenberg picture allows one to compute the elements of the quantum Fisher information matrix in this new basis from those associated with $\hat{\rho}$.

There is a subtlety in this step: when you apply the change of basis $\hat{O}$, on a local generator, contributions from non-local generators pop up. For example, the local generator $\hat{q}_1^2$, will be transformed, under a real basis change, onto $\hat{q}_1^2 \cos^2 \theta + \hat{q}_2^2 \sin^2 \theta + \sin(2\theta)\hat{q}_1\hat{q}_2$, which contains the non-local generator $\hat{q}_1\hat{q}_2$. The quantum Fisher information matrix elements associated with such a generator cannot be considered in the evaluation of the witness Eq. (12) in a specific mode basis. Nevertheless, the example above shows that, to compute the matrix associated with the local generators in an arbitrary basis, the elements of the quantum Fisher information matrix associated with non-local generators have to be known in the initial (reference) basis. Thus, we trade the burden of measuring in each basis, with measuring in a single basis and computing a larger quantum Fisher information matrix that includes contribution from non-local generators. The corresponding set of non-local generators is denoted $\mathcal{A}_{N,m}^{(NL)}$, and include all elements of the form in Eq. (17), up to order $N$. As an illustration, for $m = 2$ and $N = 2$, this set
of operators reads

$$ \mathcal{A}_{2,2}^{(\text{NL})} = \left\{ \hat{q}_1, \hat{p}_1, \hat{q}_2, \hat{p}_2, \hat{q}_1^2, \frac{1}{2}(\hat{q}_1\hat{p}_1 + \hat{p}_1\hat{q}_1), \hat{p}_1^2, \hat{q}_2^2, \right. \\ \left. \frac{1}{2}(\hat{q}_2\hat{p}_2 + \hat{p}_2\hat{q}_2), \hat{p}_2^2, \hat{q}_1\hat{q}_2, \hat{p}_1\hat{p}_2, \hat{q}_1\hat{p}_2, \hat{p}_1\hat{q}_2 \right\}. \quad (18) $$

Generally, for *m* modes and generators of order *N*, the set $\mathcal{A}_{N,m}^{(\text{NL})}$ contains

$$ \ell(N, m) = \sum_{k=1}^{N} \left( \binom{2m + k - 1}{k} \right) \quad (19) $$

elements.

The overall procedure to witness mode intrinsic entanglement reads as follows.

1. Compute the quantum Fisher information matrix $Q_{\hat{\rho}}^{\mathcal{A}_{N,m}^{(\text{NL})}}$ and the higher order covariance matrix $\Gamma_{\hat{\rho}}^{\mathcal{A}_{N,m}^{(\text{NL})}}$, associated with the set of non-local generators $\mathcal{A}_{N,m}^{(\text{NL})}$, for the state $\hat{\rho}$, on the reference basis. When the state is pure, $Q_{\hat{\rho}}^{\mathcal{A}_{N,m}^{(\text{NL})}} = 4\Gamma_{\hat{\rho}}^{\mathcal{A}_{N,m}^{(\text{NL})}}$.

2. Obtain $Q_{\hat{\rho}(\boldsymbol{\theta})}^{\mathcal{A}_{N,m}^{(\text{NL})}}$ and $\Gamma_{\hat{\rho}(\boldsymbol{\theta})}^{\mathcal{A}_{N,m}^{(\text{NL})}}$, as

$$ Q_{\hat{\rho}(\boldsymbol{\theta})}^{\mathcal{A}_{N,m}^{(\text{NL})}} = O(\boldsymbol{\theta}) Q_{\hat{\rho}}^{\mathcal{A}_{N,m}^{(\text{NL})}} O^T(\boldsymbol{\theta}), \quad (20) $$

$$ \Gamma_{\hat{\rho}(\boldsymbol{\theta})}^{\mathcal{A}_{N,m}^{(\text{NL})}} = O(\boldsymbol{\theta}) \Gamma_{\hat{\rho}}^{\mathcal{A}_{N,m}^{(\text{NL})}} O^T(\boldsymbol{\theta}), $$

where $\boldsymbol{\theta}$, parameterize the set of mode basis changes $\hat{O}$, and $O(\boldsymbol{\theta})$ stands for the matrix that represents the action of $\hat{O}$ over the quantum Fisher information matrix and the higher order covariance matrix. The exact form of this matrix is provided in Eq. (S2) of Supplement 1.

3. Remove the contribution of non-local operators. The final matrices denoted by $Q_{\hat{\rho}(\boldsymbol{\theta})}^{\mathcal{A}_{N,m}^{(\text{L})}}$ and $\Gamma_{\hat{\rho}(\boldsymbol{\theta})}^{\mathcal{A}_{N,m}^{(\text{L})}}$ will have dimension $\ell_{\text{loc}}(N, m) \times \ell_{\text{loc}}(N, m)$, with $\ell_{\text{loc}}(N, m) = mN(N+3)/2$. This step is performed by applying projectors onto the local subspaces (rectangular matrices of dimension $(\ell_{\text{loc}}(N, m) \times \ell(N, m))$, on the extended matrices.

4. Compute $\Gamma_{\Pi(\hat{\rho}(\boldsymbol{\theta}))}^{\mathcal{A}_{N}^{(\text{L})}} = \Gamma_{\hat{\rho}_A(\boldsymbol{\theta})}^{\mathcal{A}_{N}^{(\text{L})}} \oplus \Gamma_{\hat{\rho}_B(\boldsymbol{\theta})}^{\mathcal{A}_{N}^{(\text{L})}}$. This is done by projecting independently onto the covariance matrices of the reduced states $\hat{\rho}_A$ and $\hat{\rho}_B$, and construct their direct sum, which is equivalent to keeping only the diagonal blocks of the full covariance matrix.

5. Compute $E_{\hat{\rho}}^{\mathcal{A}_{N,m}}(\boldsymbol{\theta}) = \lambda_{\max} \left[ Q_{\hat{\rho}(\boldsymbol{\theta})}^{\mathcal{A}_{N,m}} - 4 \Gamma_{\Pi(\hat{\rho}(\boldsymbol{\theta}))}^{\mathcal{A}_{N,m}} \right]$, i.e., the entanglement criterion in the new basis. Then, if

$$ W_Q(\hat{\rho}, \mathcal{A}_{N,m}) = \min_{\boldsymbol{\theta}} E_{\hat{\rho}}^{\mathcal{A}_{N,m}}(\boldsymbol{\theta}) > 0, \quad (21) $$

the state will not be passively separable.

Equation (21) provides thus a witness of mode intrinsic quantum entanglement that can be evaluated through numerical optimization over all mode bases. For a small number of modes, it can be rather instructive to evaluate the witness over a grid covering all possible mode bases. In this way, the procedure described above not only provides a way to witness mode intrinsic entanglement, but also provides a way to characterize the entanglement properties of a given resource state over different mode bases.

For a given number of modes *m*, the most general orthogonal (passive), transformation *O* that we are concerned with can be

represented in a compact way, using, for example, Clements' decomposition [60], by the combination of $m(m-1)/2$ real beam splitters and $m(m-1)/2$ local rotations (phase shifts), as represented in Fig. 2. The transformation depends on $m(m-1)$ parameters that we gather in the vector $\boldsymbol{\theta} = (\theta, \varphi)$, where $\theta$ contains the $m(m-1)/2$ angles of the beam splitters and $\varphi$ contains the $m(m-1)/2$ phases of the phase shifters. To reproduce the full set of passive transformations, an additional set of *m* local phase shifters should be included; nevertheless, it is irrelevant for our task, given that these do not alter the entanglement properties of the states we consider. The complexity of the protocol increases with both the number of modes and the order of the generators, since the number of generators contributing grows exponentially on both *N* and *m*. The most demanding steps are the first and second. The first involves independently computing all the upper triangular part of both the covariance matrix and the quantum Fisher information matrix. Each such computation can, in itself, be quite costly, depending on the representation of the quantum state at hand. The quantum Fisher information matrix can only be computed from the density matrix of states, which involves a number of elements exponential in the number of modes, making it by far the most demanding step, even at low order in the generators. The second step involves the computation of the transformation matrix and then matrix multiplication. Importantly, steps 1–4 have to be performed only once, as we obtain a parameterized representation of the matrices involved in step 5, which is then performed until the minimum is found. A full exploration over the space of modes requires $(N_{\text{grid}})^{m^2}$ evaluations, where $N_{\text{grid}}$ is the number of bins in which we partition the set of possible values of the phases and beam splitter angles. Nevertheless, a more clever optimization drastically reduces the number of such evaluation, depending on the roughness of the landscape of the function.

Typically, we do not require to go on the degree of the generators further than order two. At order two, the complexity of the second step is $O(m^6)$, and for the first step, it is $O(N_{\text{Fock}}^m m^4)$, where $N_{\text{Fock}}$ is the cutoff in Fock space, thus the dimensionality of the density matrix. This implies that the dominant burden in complexity is given by the computation using the density matrix. In latter sections of this manuscript, we explore ways to reduce this complexity by relying on lower bounds of the quantum Fisher information matrix, given by the classical Fisher information matrix associated with a given measurement setting.

## 4. APPLICATION OF THE WITNESS TO EXPERIMENTALLY RELEVANT STATES

To explore the capability of the witness derived in the previous section to detect mode intrinsic entanglement for bipartite and multimode systems, we use it on a specific class of non-Gaussian states: photon subtracted states. These states are produced by subtracting one or more photons from a Gaussian state. The removal of photons is a non-Gaussian operation, thus rendering the state non-Gaussian. Therefore, it makes it possible to explore how quantum properties are affected by the non-Gaussianity using a class of states which falls within a well-established theoretical framework [61,62], and has the advantage to be already produced in the lab with state-of-the-art experiments using different spatial and temporal modes of indistinguishable sources [19–21], or in the temporal-frequency domain [18]. The general
<figure><img src="image_4.png" /><figcaption>Fig. 2. Clements’ design for the representation of an arbitrary m-mode passive interferometer [60]. Red blocks represent local transformations parameterized by an angle $\phi \in [0, 2\pi)$, and blue blocks involving two neighboring modes represent real beam splitters, parameterized by an angle $\theta \in [0, 2\pi)$. A set of m local rotations was omitted after the last set of beam splitters as they do not modify the entanglement properties of the probe state.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Graphical representation of the family of states that we consider. On the left, we depict the cluster state as generated by applying an orthogonal transformation on a set of squeezed states. On the right, we depict the cluster state by representing the links between the modes that interact directly along the transformation $\hat{O}$. The final state to build our probe states is to de-Gaussify these clusters by performing photon subtraction in specific modes of our choice. In particular, we will consider only up to two-photon subtractions.</figcaption></figure>

expression of the states that we consider is given by

$$|\psi\rangle \propto \left( \prod_{i=1}^{N_S} \hat{a}_i \right) \hat{O} \left( \prod_{j=1}^{M} \hat{S}(r_j, \bar{\phi}_j) \right) |0\rangle, \qquad (22)$$

up to a normalization factor, where $N_S$ is the number of photons subtracted and $M$ is the number of modes in the cluster. Here, $\hat{S}(r_j, \bar{\phi}_j)$ is a single-mode squeezing operation that squeezes the vacuum state of mode $j$ along the quadrature $\hat{\xi}_{\phi_j} = \hat{q}_j \cos \bar{\phi}_j + \hat{p}_j \sin \bar{\phi}_j$, with a squeezing level $r_j$. The operation $\hat{O}$ is a general passive transformation, and $\hat{a}_i$ are the annihilation operators corresponding to mode $i$. In Fig. 3, we illustrate the case $N_i = 2$. We also point out that photon subtractions can happen in a superposition of different modes of the cluster state. This scenario will only be considered for two-mode systems, as in any case, it can be cast into the form above by including the mixing operation in the unitary $\hat{O}$ (see Sec. 4.2).

We split the following discussion in two subsections. We start by analyzing multimode cluster states and show that the witness can be applied and successfully detect mode-intrinsic entanglement. In particular, we apply the metrological witness to three-mode, four-mode, and five-mode cluster states. Then, we move to two-mode photon subtracted states. In some configurations discussed below, this leads to mode-intrinsic entangled

states that the witness can detect. For such a two-mode system, the unitary transformation $O(\vartheta)$ to go from one basis to another (as represented in Fig. 2) only depends on two parameters $\vartheta = (\theta, \varphi)$. Therefore, one can easily visualize the metrological witness $E_{\varphi}^{N,2}(\theta, \varphi)$ defined in Eq. (12) in all possible bases in a two-dimensional parameter space. In Sec. 5, we exploit the classical Fisher information associated with specific measurements, to construct lower bounds for the metrological witness, that can be obtained experimentally. Thus, we explore the potential of the method to detect mode-intrinsic entanglement in state-of-the-art experiments.

## 4.1. Multimode CV-Cluster States

A CV-cluster state is made of a set of modes which are coupled together by applying entangling gates $C_Z = \exp(i\hat{q}_i \otimes \hat{q}_j)$. Overall, they form a graph structure ($V, E$), where $V$ is the set of nodes, i.e., the set of modes in the multimode state, while $E$ represents the set of edges between the nodes, i.e., the quantum correlations that connect the modes together after applying the $C_Z$ gates. The graph formed by the cluster state can be described by the so-called adjacency matrix $V$ whose coefficients $V_{ij}$ are 1 if an edge exists between nodes $i$ and $j$ (in other words, if modes $i$ and $j$ are entangled to each other), 0 otherwise. We show in
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Photon subtracted (a) three-mode, (b) four-mode, (c) five-mode cluster states under consideration. Nodes represent modes and black edges represent entanglement links between modes. The red nodes indicate the mode where we subtract a photon.</figcaption></figure>

Fig. 4 the three cluster states that we will consider in this section. As an example, we give the adjacency matrix of the three-mode cluster state:

$$V_3\text{-modes} = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 0 \end{pmatrix}.$$

The red nodes in Fig. 4 indicate in which mode we subtract a photon to create our photon subtracted cluster state.

A cluster state associated with a given adjacency matrix $V$ is obtained from a set of single-mode squeezed states by applying a unitary matrix

$$U_v = (1 + iV)(V^2 + 1)^{-1/2} O, \quad (23)$$

where $1$ is the unitary matrix and $O$ can be any real orthogonal matrix [63–65]. The choice of $O$ can be optimized to target a specific cluster configuration as discussed in Ref. [64]. Indeed, a perfect cluster state where correlations between modes are exactly given by the adjacency matrix can only be obtained in the limit of infinite squeezing. In this case, the variances of certain combinations of quadratures, called nullifiers, are exactly zero; these nullifiers are defined as

$$\hat{\delta}_i \equiv \hat{p}_i - \sum_{k=1}^{m} V_{ik} \hat{q}_k \quad \forall i \in \{1, \dots, m\}, \quad (24)$$

where $m$ is the number of modes in the cluster. In the case of finite squeezing, the variances of the nullifiers are not zero but one can minimize them, by optimizing the orthogonal matrix $O$ in Eq. (23). We thus follow the procedure described in Ref. [64] using a genetic algorithm to minimize the variances of the nullifiers and obtain the best cluster states given realistic input squeezing for the three graph structures in Fig. 4. To consider realistic squeezing parameters, we take those measured in the 4-mode cluster states of Ref. [66]: –7 dB, –6 dB, –4 dB, 0 dB. For the other cluster states considered here, we use similar squeezing parameters: –7 dB, –4 dB, 0 dB (3-mode cluster state); –7 dB, –6 dB, –4 dB, –2 dB, 0 dB (5-mode cluster state).

Before we subtract a photon, the cluster states are fully Gaussian. Therefore, we do not expect to witness any mode-intrinsic entanglement since there will always exist a basis where the states are separable. The subtraction of a photon in the first mode (the red nodes in Fig. 4) renders the state non-Gaussian and will affect the modes in its vicinity; we thus anticipate to detect mode-intrinsic entanglement.

In the case of an $m$-mode system, we have the freedom to choose how we split it into $p$ different partitions $P_1|\dots|P_p$ and eventually detect $p$-partite mode-intrinsic entanglement. It is

<table><caption>**Table 1.** Values of the Mode-intrinsic Entanglement Witness Derived in Eq. (21) for the 3-mode Cluster State Considered in Fig. 4(a) and for All Possible Partitions (up to Permutations)</caption><thead><tr><th>Partitions</th><th>1|2|3</th><th>12|3</th></tr></thead><tbody><tr><td>$W_Q(\hat{\rho}, \mathcal{A}_{N=1,3})$</td><td>0.94</td><td>0</td></tr></tbody></table>

<table><caption>**Table 2.** Values of the Mode-intrinsic Entanglement Witness Derived in Eq. (21) for the 4-mode Cluster State Considered in Fig. 4(b) and for All Possible Partitions (up to Permutations)</caption><thead><tr><th>Partitions</th><th>1|2|3|4</th><th>12|3|4</th><th>12|34</th><th>1|234</th></tr></thead><tbody><tr><td>$W_Q(\hat{\rho}, \mathcal{A}_{N=1,4})$</td><td>0.87</td><td>0.33</td><td>0.33</td><td>0</td></tr></tbody></table>

<table><caption>**Table 3.** Values of the Mode-intrinsic Entanglement Witness Derived in Eq. (21) for the 5-mode Cluster State Considered in Fig. 4(c) and for all Possible Partitions (up to Permutations)</caption><thead><tr><th>Partitions</th><th>1|2|3|4|5</th><th>12|3|4|5</th><th>12|34|5</th><th>123|4|5</th><th>123|45</th><th>1234|5</th></tr></thead><tbody><tr><td>$W_Q(\hat{\rho}, \mathcal{A}_{N=1,5})$</td><td>0.92</td><td>0.31</td><td>0.31</td><td>0.17</td><td>0.17</td><td>0</td></tr></tbody></table>

important to note that the labeling of the modes is irrelevant in the study of mode-intrinsic entanglement; only the number of modes in each partition matters. Indeed, any two different partitions, differing only by the labeling of the modes, are related by a swapping of the modes, which is a passive linear optics operation. Therefore, any two such partitionings are equivalent regarding mode-intrinsic entanglement.

Tables 1, 2, 3 gather the values of the mode-intrinsic entanglement witness $W_Q(\hat{\rho}, \mathcal{A}_{N=1,m})$ – using first-order generators only – obtained for the three clusters and for all possible mode splittings (up to a reordering of the mode labels as discussed above). Given the number of parameters that we need to optimize to find the minimum of the entanglement witness $E_{\rho}^{\mathcal{A}_N}(\theta)$ (see Fig. 2, where each blue and red block of the unitary transformation adds one degree of freedom), we also use here a genetic algorithm to solve this optimization problem. As we can clearly see from the results, for many mode splittings, the subtraction of the photon renders the cluster states mode-intrinsic entangled. However, it is very instructive to notice that only in the case where we split the cluster states into a bipartition where one of the subsystems contain $m-1$ modes, the value of the witness becomes zero. Notice that in all clusters considered, the last mode was in the vacuum state (zero squeezing), meaning that the non-classicality was constrained to $m-1$ modes. This offers an interesting insight on the interplay between non-Gaussianity and Gaussian non-classicality, that only together can lead to features like mode-intrinsic entanglement.

Although it is difficult to draw any systematic conclusion from the results above, this shows that the witness we propose, already at first-order of the generators, is able to detect mode-intrinsic entanglement for experimentally relevant cluster states. It also opens a new toolbox to explore and better understand the entanglement structure of large multimode non-Gaussian states such as those considered in Ref. [67]. A deeper exploration of
this tool for the study of large clusters is left for future works. In the rest of this paper, we rather focus on the study of the smallest version of these systems, i.e., two-mode, two-photon-subtracted states. In this more controlled setting, we can fully explore the capabilities and limitations of the method we develop, while also restricting to a scenario with more immediate potential to be tested experimentally.

## 4.2. Two-Mode Photon Subtracted States

As explained above, we now turn our attention to two-mode photon subtracted states of the form

$$|\psi\rangle = \prod_{i=1}^{N_S} (\cos \Theta_i \hat{a}_i + \sin \Theta_i \hat{a}_2) \hat{S}(r_1) \hat{S}(r_2) |00\rangle, \quad (25)$$

where $\hat{S}(r_1)$ and $\hat{S}(r_2)$ are single-mode squeezing operators, equivalent to those in Eq. (22), where we removed the phase dependence $\hat{\phi}_i$ and allow $r_i$ to be negative. Angles $\Theta_i$ control mode superposition in which photons are subtracted. In general, we could have considered the subtraction to happen in an already entangled Gaussian state; nevertheless, those can always be brought to the form in Eq. (25) by the action of some orthogonal transformation. In this subsection, we consider two cases: $N_S = 1$ and $N_S = 2$. In both cases, we consider the single-mode squeezed states to be either squeezed in position ($r_i>0$) or in momentum ($r_i<0$). One-photon-subtracted pure states ($N_S = 1$) are not mode-intrinsic entangled if: (i) $r_1 = r_2$ for all $\Theta_1$; or (ii) $\Theta_1 = n\pi/2, n \in \mathbb{Z}$ for all $r_1, r_2$. Similarly, two-photon-subtracted states ($N_S = 2$) are not mode-intrinsic entangled if: (i) $r_1 = r_2, \Theta_1$ and $\Theta_2$ corresponding to orthogonal modes; (ii) $\Theta_1 = n\pi/2$ and $\Theta_2 = m\pi/2, m \in \mathbb{Z}$ for all $r_1, r_2$. In all other cases, one can prove that these states are entangled in all mode bases. To illustrate the capability of the metrological witness Eq. (21) to detect mode intrinsic entanglement, we now consider two specific examples: $N_S = 1, r_1 = -r_2 = 0.2, \Theta_1 = \pi/4; N_S = 2, r_1 = -r_2 = 0.2, \Theta_1 = -\Theta_2 = \pi/4$.

In the case of a two-mode state, the mode basis change transformation at the level of quadrature operators using Clement's decomposition can be written as

$$O(\theta, \varphi) = O_{BS}(\theta) O_{PS}(\varphi), \quad (26)$$

where

$$O_{BS}(\theta) = \begin{pmatrix} \cos \theta & -\sin \theta & 0 & 0 \\ \sin \theta & \cos \theta & 0 & 0 \\ 0 & 0 & \cos \theta & -\sin \theta \\ 0 & 0 & \sin \theta & \cos \theta \end{pmatrix},$$

$$O_{PS}(\varphi) = \begin{pmatrix} \cos \varphi & 0 & -\sin \varphi & 0 \\ 0 & 1 & 0 & 0 \\ \sin \varphi & 0 & \cos \varphi & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix},$$

correspond respectively to a beam splitter and a phase shift transformation.

We can then determine the action of the mode basis transformation on the generators, $O(\theta, \varphi)$, following the indications in Appendix 2, and apply the metrological witness $E_{\rho}^{\mathcal{A}_N}(\theta, \varphi)$ [see Eq. (21)] using first-order and second-order generators $\mathcal{A}_{N=1}^{(L)}$ and $\mathcal{A}_{N=2}^{(L)}$, respectively. Notice that we have omitted the subscript $m=2$ from the definition of the set of generators, as from this point on, we will only be dealing with two-mode states.

The results are shown in Fig. 5 for $N_S = 1$ and Fig. 6 for $N_S = 2$. In both figures, the top graphs, labeled panel (a), correspond to the results using generators in the set $\mathcal{A}_{N=1}$, while we used the set $\mathcal{A}_{N=2}$ to produce the bottom plots, labeled panel (b). The black dots indicate the minimal value of $E_{\rho}^{\mathcal{A}_N}(\theta, \varphi)$ over the parameter space $(\theta, \varphi)$. For the one-photon-subtracted state, the metrological witness is able to detect mode intrinsic entanglement both at first order and second order on the generators: the minimal value over the parameter space is always strictly greater than zero, thus demonstrating entanglement in all mode bases. The robustness of the witness against optical losses is addressed below (see in particular Figs. 7 and 8). In the case of the two-photon-subtracted state that we analyze, the situation is different; as shown in Fig. 6(a), at first order, we do not detect entanglement in every mode-basis: the minimum of $E_{\rho}^{\mathcal{A}_{N=1}}(\theta, \varphi)$ reaches zero for a range of values $\theta = 0 (\forall \varphi), \theta = \pi/2 (\forall \varphi)$, and $\varphi = \pi/2 (\forall \theta)$ (see the black lines at the borders of the figure). Therefore, the first-order protocol cannot be used to assess the mode intrinsic entanglement of this state. In contrast, as shown in Fig. 6(b), once we use second-order generators, the minimal value over all mode-bases becomes strictly positive (it equals 0.98 in this precise example), thus meaning that the state is mode-intrinsic entangled. This illustrates the interest to use higher order generators and to not restrict the metrological witness to first-order generators. More generally, one can show that for two-photon-subtracted states ($N_S = 2$), the criteria based on first-order generators always fails when modes corresponding to $\Theta_1$ and $\Theta_2$ are orthogonal to each other.

So far, we did not include optical losses, yet they are inherent to any experimental setup. Therefore, in what follows, we apply the metrological witness to lossy one- and two-photon-subtracted states. To simulate such losses, we apply two single-mode loss channels $\hat{L}_1(\eta)$ and $\hat{L}_2(\eta)$ to each mode as described in Ref. [68], where $\eta$ represents the amount of losses in the channel, i.e., $\eta = 1$ implies 100% losses.

Figures 7(a) and 7(b) show the value of $W_Q(\hat{\rho}, \mathcal{A}_N) = \min_{\varphi, \alpha} E_{\rho}^{\mathcal{A}_N}(\theta, \varphi)$ for first- and second-order generators ($N = 1$ and $N = 2$) as a function of the optical losses $\eta$ and of the angle $\Theta_1$ for a one-photon-subtracted state with $r_1 = -r_2 = 0.2$. We recall that $\Theta_1$ parameterizes the basis where we subtract the photon. Figures 8(a) and 8(b) also show the value of $W_Q(\hat{\rho}, \mathcal{A}_N)$ for $N = 1$ and $N = 2$ as a function of $\eta$, and of the angle of the second photon subtraction $\Theta_2$, for a two-photon-subtracted state with $r_1 = 0.17, r_2 = -0.3$, and $\Theta_1 = \pi/4$. This choice of squeezing parameters is motivated from the results obtained in Ref. [56], where it has been observed that such squeezing parameters lead to a better resilience to losses. The black dashed curve delimits the regions of positive (red area) and negative (blue area) values. The red area thus corresponds to the region where the state is detected to be not passively separable. As one can see from the figures, the metrological witness is able to detect mode-intrinsic entanglement even in the presence of losses. However, using second-order generators results in a much better resilience to losses compared with first order both for one- and two-photon-subtracted states; it thrives to detect entanglement in all mode bases for a wide range of losses (up to 95 % of losses in the case of a two-photon-subtracted state when $\Theta_2 \approx -0.25$).

The maximal amount of losses to detect intrinsic mode entanglement shown in Figs. 7 and 8 (black dashed lines) is a theoretical bound and represents more a theoretical case study to demonstrate the robustness of the witness against losses, as we compute it using the quantum Fisher information matrix $Q_{\hat{\rho}}^{\mathcal{A}_N}$
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Entanglement witness (a) $E_{\rho}^{\mathcal{A}_{N=1}}(\theta, \varphi)$ for the set of generators $\mathcal{A}_{N=1}$ (first order) and (b) $E_{\rho}^{\mathcal{A}_{N=2}}(\theta, \varphi)$ for generators up to second order ($\mathcal{A}_{N=2}$) as a function of the mode basis change parameters $(\theta, \varphi)$. Here, we consider a one-photon-subtracted state with $\Theta_1 = \pi/4$, and $r_1 = -r_2 = 0.2$, which corresponds to 1.7 dB of squeezing. The black dots indicate the minimal value over the parameter space.</figcaption></figure>

to evaluate the metrological criterion in Eq. (21). Such a quantity is difficult, if not impossible, to measure experimentally; what can be measured though is the Fisher information matrix (upper bounded by the quantum Fisher information), resulting in a weaker version of the metrological-witness of mode-intrinsic entanglement. However, even in this case, we show in Sec. 5 that the witness is able to detect entanglement in all mode bases while considering quantities that are experimentally accessible.

## 5. RELAXATION OF THE WITNESS

So far, we have considered a protocol that relies on the quantum Fisher information matrix associated with a given set of

generators, for a quantum state. This is a strong requirement if we want to design a protocol that aims to be applicable in experimental settings since it requires full state tomography, an unfeasible task in a multimode setting. Nevertheless, for a given measurement setting, we know that the Fisher information is a lower bound of the quantum Fisher information [see Eq. (2)]. Therefore, one can derive a lower bound – experimentally accessible – of the mode-intrinsic entanglement witness introduced in the previous sections by relying on the (classical) Fisher information matrix – instead of the quantum Fisher information matrix – whose elements will be defined below in Eq. (27).
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Entanglement witness (a) $E_{\rho}^{A_{N=1}}(\theta, \varphi)$ for the first-order set of generators $\mathcal{A}_{N=1}$ and (b) $E_{\rho}^{A_{N=2}}(\theta, \varphi)$ for generators up to second order ($\mathcal{A}_{N=2}$) as a function of the mode basis change parameters $(\theta, \varphi)$. Here, we consider a two-photon-subtracted state with $r_1 = -r_2 = 0.2$, $\Theta_1 = -\Theta_2 = \pi/4$. The thick black lines for panel (a) and the black dot for panel (b) indicate the minimal value over the parameter space.</figcaption></figure>

To proceed, as illustrated in Fig. 9, we consider a setting in which we implement, either physically or, if possible, in post-processing, a transformation $\exp(-i \sum_{k=1}^{l(N,m)} \kappa_k \hat{H}_k^S)$ on a quantum state $\hat{\rho}$. We recall that the generators $\hat{H}_k^S$ are the symmetrized operators defined in Eq. (17) and $l(N, m)$ is the cardinality of the set containing all generators of order $N$ for an $m$-mode system [see Eq. (19)]. Each generator implements a parameter $\kappa_j$ that we gather in a vector $\mathbf{\kappa} = (\kappa_1, \dots, \kappa_{l(N,m)})$. Then, we perform homodyne detection in each mode $i$ of the system in quadratures $\hat{\xi}_{\phi_i} = \cos \phi_i \hat{q}_i + \sin \phi_i \hat{p}_i$ with $\phi_i \in [0, \pi)$, i.e., we consider measurement projectors $\hat{\Pi}_i = |\xi_{\phi_i}\rangle\langle\xi_{\phi_i}|$, where $|\xi_{\phi_i}\rangle$ is

an eigenstate of the quadrature operator $\hat{\xi}_{\phi_i}$. The elements of the classical Fisher information matrix for a given state $\hat{\rho}$ and the measurement setting described above are computed as

$$ F(\hat{\rho}, \mathcal{A})_{i,j} = \int_{\mathbb{R}^m} (\Pi_{i=1}^m d\xi_{\phi_i}) p(\xi_{\phi_1}, \dots, \xi_{\phi_m} | \mathbf{\kappa} = \mathbf{0}) \frac{\partial \log p(\xi_{\phi_1}, \dots, \xi_{\phi_m} | \mathbf{\kappa}_1, \dots, \mathbf{\kappa}_{l(N,m)})}{\partial \mathbf{\kappa}_i} \bigg|_{\mathbf{\kappa}=\mathbf{0}} \quad (27) \\ \frac{\partial \log p(\xi_{\phi_1}, \dots, \xi_{\phi_m} | \mathbf{\kappa}_1, \dots, \mathbf{\kappa}_{l(N,m)})}{\partial \mathbf{\kappa}_j} \bigg|_{\mathbf{\kappa}=\mathbf{0}}, $$
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Metrological witness $W_Q$ [see Eq. (21)] for a single-photon-subtracted state ($N_s = 1$) with $r_1 = -r_2 = 0.2$. The witness is computed (a) using the set of first-order generators $\mathcal{A}_{N=1}$ and (b) using generators up to second order ($\mathcal{A}_{N=2}$), as a function of the optical losses $\eta$ and the angle of photon subtraction $\Theta_1$ ranging from $-\pi/4$ to $\pi/4$. The black dashed lines delimit the region of negative values (blue area) from the region of positive values (red area). The red area corresponds to the region where the state is detected to be not passively separable.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Metrological witness $W_Q$ [see Eq. (21)] for a two-photon-subtracted state ($N_s = 2$) with $\Theta_1 = \pi/4$, and $r_1 = 0.17$, $r_2 = -0.3$, which respectively correspond to 1.5 dB, -2.6 dB. The witness is computed (a) using the set of first-order generators $\mathcal{A}_{N=1}$ and (b) using generators up to second order ($\mathcal{A}_{N=2}$), as a function of the optical losses $\eta$ and the angle of the second photon subtraction $\Theta_2$ ranging from $-\pi/4$ to $\pi/4$. The black dashed lines delimit the region of negative values (blue area) from the region of positive values (red area). The red area corresponds to the region where the state is detected to be mode-intrinsic entangled.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 9. Measurement setting considered for the estimation of a lower bound on the quantum Fisher information matrix. The joint measurement of quadratures in the different modes allows to build a probability distribution of outcomes $p(\xi_{\phi_1}, \dots, \xi_{\phi_N})$. From the influence of the different generators in this object, we can estimate the associated classical Fisher information matrix.</figcaption></figure>

where $p(\xi_{\phi_1}, \dots, \xi_{\phi_m} | \kappa_1, \dots, \kappa_{l(N,m)})$ corresponds to the probability distribution of homodyne outcomes after all the parameters $\kappa$ are implemented by the corresponding generators $\mathcal{A}$. To compute this quantity, we need to know how each generator in our set acts on the probability distribution. Both theoretically and experimentally, the problem is feasible only up to generators of order two. For two modes, this set is given by $\mathcal{A}_{N=2}^{(\text{NL})}$ in Eq. (18). All of them correspond to symplectic transformations or displacements that can be implemented theoretically by a change of coordinates in the Wigner function,

$$W(\vec{\xi} | \kappa_j) = W(S_j^{-1} \vec{\xi} - \vec{\alpha}_j), \qquad (28)$$

where $S_j$ and $\vec{\alpha}_j$ correspond to the symplectic transformation and displacements, respectively, associated with the generator $\hat{H}_j^S$ contained in $\mathcal{A}_{N=2}^{(\text{NL})}$, while $\vec{\xi} = (q_1, \dots, q_m, p_1, \dots, p_m)$ gather all quadratures. Once the Wigner function is transformed, to obtain the probability distribution of outcomes for the measurement setting that we consider, we can, at the level of the Wigner function, rotate each quadrature by $\phi_i$, $i \in \{1, \dots, m\}$, i.e.,

$$W(\vec{\xi}_{\phi} | \kappa_j) = W[\mathcal{R}(\phi) (S_j^{-1} \vec{\xi} - \vec{\alpha}_j)], \qquad (29)$$

where $\mathcal{R}(\phi)$ is a rotation operator that transforms the set of quadratures $(q_1, \dots, q_m, p_1, \dots, p_m)$ to $\vec{\xi}_{\phi} \equiv (\xi_{\phi_1}, \dots, \xi_{\phi_m}, \vec{\xi}_{\phi_1}, \dots, \vec{\xi}_{\phi_m})$, with $\vec{\xi}_{\phi_i} = -q_i \sin \phi_i + p_i \cos \phi_i$. As a final step, we just need to marginalize to the $m$ quadratures of interest to obtain the corresponding parameterized probability distribution

$$p(\xi_{\phi_1}, \dots, \xi_{\phi_m} | \kappa_j) = \int_{\mathbb{R}^m} \Pi_{i=1}^m d\vec{\xi}_{\phi_i} W(\vec{\xi}_{\phi} | \kappa_j). \qquad (30)$$

From an experimental point of view, symplectic operations are expected to be feasible [69–71], even though currently challenging to implement. However, they are a necessary requirement for universal continuous variable quantum computation [8,71]. Even though we are not currently capable to experimentally implement all of these unitaries, we should keep in mind that some of them can also be implemented in post-processing of the measurement data [72]. The most simple example is the case of displacements, generated by the quadrature operators (i.e., at order one), for example,

$$p(\xi_{\phi_1}, \dots, \xi_{\phi_m} | \kappa_{m+1}) = p(\xi_{\phi_1} - \kappa_{m+1} \cos(\phi_1), \dots, \xi_{\phi_m}), \qquad (31)$$

where $\kappa_{m+1}$ is the parameter generated by the $(m+1)$th generator contained in $\mathcal{A}_{N=1}^{(\text{NL})}$, i.e., $\hat{p}_1$. At order 2, the only generators that

cannot be treated in post-processing are the shearing, i.e., those of the form $\hat{q}_i^2$, whose effect is given by

$$W(\vec{\xi} | \kappa_{\hat{q}_i^2}) = W(q_1, \dots, q_m, p_1, \dots, p_i - 2\kappa_{\hat{q}_i^2} q_i, \dots, p_m), \qquad (32)$$

which mixes quadratures that cannot be simultaneously measured.

We can then construct the Fisher information matrix $F(\hat{\rho}, \mathcal{A})$ by inserting the resulting probability distributions in Eq. (30) in a given measurement setting $\vec{\xi}_{\phi}$ in Eq. (27). Each choice of measurements gives rise to a different Fisher information matrix. To illustrate this point, we refer to Fig. 10 which shows how this matrix is affected by two different choices of measurement settings for a one-photon-subtracted state. For the sake of simplicity, we restrict ourselves to two modes and to displacement operators, i.e., first-order generators contained in the set $\mathcal{A}_{N=1} = (\hat{q}_1, \hat{q}_2, \hat{p}_1, \hat{p}_2)$ (denoted as the vector operator $\hat{\xi}$ in Fig. 10). In the first scenario, we measure position quadratures $q_1$ and $q_2$ [see Fig. 10(a)], which results in the probability distributions $p(q_1, q_2|0)$ and $p(q_1, q_2|\kappa)$ [shown in panel (b)] before and after we apply the transformation on the state. Note that for displacement operators, the transformation can be applied in post-processing, i.e., directly from the probability distribution $p(q_1, q_2|0)$ reconstructed from experimental data. When measuring the quadratures $(q_1, q_2)$, the probability distribution of outcomes is only affected by displacement operators $\hat{p}_1$ and $\hat{p}_2$, resulting in non-zero elements in the Fisher information matrix only for these generators [as shown in Fig. 10(c)]. In the second scenario, we measure momentum quadratures $p_1$ and $p_2$ [see Fig. 10(d)]. In this case, only generators $\hat{q}_1$ and $\hat{q}_2$ displace the distribution in the $(p_1, p_2)$-setting resulting in a different Fisher information matrix. We thus obtain a set of Fisher information matrices $F(\hat{\rho}, \mathcal{A}_N, \phi)$ depending on the state $\hat{\rho}$, on the set of generators $\mathcal{A}_N$, and on the measurement setting parameterized by $\phi$.

In all possible measurement settings, the classical Fisher information matrix is a lower bound on the quantum Fisher information matrix, i.e., $F(\hat{\rho}(\vartheta), \mathcal{A}_N, \phi) \ll Q(\hat{\rho}(\vartheta), \mathcal{A}_N)$, leading to the following measurement-setting-dependent entanglement witness:

$$E_{\rho(\vartheta), \text{hom}}^{\mathcal{A}_N}(\phi) \equiv \lambda_{\text{max}} \left( F(\hat{\rho}, \mathcal{A}_N, \phi) - 4 \Gamma_{\Pi(\hat{\rho}(\vartheta))}^{\mathcal{A}_N} \right). \qquad (33)$$

Here, we remark the dependence on the mode basis $\vartheta$ in which we probe the entanglement of the state $\hat{\rho}$. The objective is to then lift the criteria to probe in different basis, in the same way as we
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Schematic procedure to obtain the Fisher information matrix for different measurement settings in the case of a one-photon-subtracted state and first-order generators. (a) Schematic representation of the metrological parameter estimation procedure generated by operators contained in the vector $\xi = (\hat{q}_1, \hat{q}_2, \hat{p}_1, \hat{p}_2)$, i.e., $\mathcal{A}_{N=1}$, and whose parameters are gathered in $\kappa = (\kappa_1, \kappa_2, \kappa_3, \kappa_4)$. The measurement setting consists of homodyne detection along the position quadratures $q_1$ and $q_2$. (b) Effect of the displacement transformation $\exp(-i\kappa \cdot \xi)$ on the probability distribution $p(q_1, q_2)$. (c) Resulting Fisher information matrix. (d) Same schematic representation of the metrological protocol in the measurement basis $(p_1, p_2)$. (e) Effect of the displacement transformation $\exp(-i \sum_{k=1}^4 \kappa_k \xi_k)$ on the probability distribution $p(p_1, p_2)$. (f) Resulting Fisher information matrix.</figcaption></figure>

did for the quantum Fisher information matrix. This witness is a lower bound of the witness given in Eq. (12)

$$E_{\rho, \vartheta_0, \text{hom}}^{\mathcal{A}_N}(\phi) \le E_{\rho, \vartheta_0}^{\mathcal{A}_N}(\phi).$$

We can then transform the classical Fisher information measured in a mode basis $\vartheta_0$, to any other basis $\vartheta$ using the same transformation as for the QFI, in Eq. (20). The resulting matrix, though not necessarily a Fisher information matrix corresponding to a specific measurement, is still a lower bound of the quantum information matrix in the new basis, i.e.,

$$O(\vartheta)F(\hat{\rho}(\vartheta_0), \mathcal{A}_N, \phi)O(\vartheta)^\top \le O(\vartheta)Q(\hat{\rho}(\vartheta_0), \mathcal{A}_N)O(\vartheta)^\top.$$

This allows us to use the classical Fisher information matrix obtained in a specific basis $\vartheta_0$ to probe any other basis $\vartheta$. It is worth noticing that once we probe in this way a new basis, the witness we obtain is dependent on $\vartheta_0$, and is not the same as that in Eq. (33). We define it as

$$E_{\hat{\rho}, \text{hom}}^{\mathcal{A}_N}(\phi, \vartheta | \vartheta_0) \\ \equiv \lambda_{\max} \left( O(\vartheta) F(\hat{\rho}, \mathcal{A}_N, \phi) O(\vartheta)^{\top} - 4 \left( O(\vartheta) \Gamma^{\mathcal{A}_N} O(\vartheta)^{\top} \right)_{\text{\textbf{Pi}}(\hat{\rho})} \right). \quad (34)$$

Using this witness, we can, similarly as in Eq. (21), define a (measurement-setting dependent) mode-intrinsic entanglement witness given by

$$\min_{\vartheta} E_{\hat{\rho}, \text{hom}}^{\mathcal{A}_N}(\phi, \vartheta | \vartheta_0) > 0.$$

We can create stronger versions of this witness by considering for each mode basis $\vartheta$ the best value of the witness $E_{\hat{\rho}, \text{hom}}^{\mathcal{A}_N}(\phi^{(i)}, \vartheta | \vartheta_0)$ for different homodyne settings $\{\phi^{(1)}, \dots, \phi^{(n)}\}$, and then minimizing over the mode-basis parameter space $\vartheta$:

$$W_{\text{hom}}(\hat{\rho}, \mathcal{A}_N | \vartheta_0) = \min_{\vartheta} \max_{\{\phi^{(1)}, \dots, \phi^{(n)}\}} E_{\hat{\rho}, \text{hom}}^{\mathcal{A}_N}(\phi^{(i)}, \vartheta | \vartheta_0). \quad (35)$$

Usually, we will consider $\vartheta_0 = 0$, which corresponds to measure directly in the basis in which the state is prepared. As it turns out, this seems to be the best choice for the kind of states we consider, i.e., states obtained by photon subtractions in a superposition of

two independently squeezed modes. For other kind of states, an exploration of different initial basis for the measure is advised to potentially increase the performance of the witness.

In what follows, we apply this procedure on two-mode two-photon-subtracted states that we have considered in previous sections. In all cases, the plots were constructed using the Fisher information matrices corresponding to the measurement settings $\{q_1, q_2\}$ and $\{p_1, p_2\}$ (as in Fig. 10), which are the maximally informative homodyne settings for the states that we consider. We apply the criteria on both of these Fisher information matrices as explained above and for each mode basis $\vartheta = (\theta, \varphi)$, we pick the highest value between the two of them. This precludes us from doing optimization for the determination of the mode-intrinsic entanglement witness Eq. (35), but in general, at order two of the generators and for a low number of modes, it is still an efficient procedure.

We can observe in Fig. 11 an analysis of the evaluation of the criteria at second order of the generators, on two-photon-subtracted states under the effect of losses. In particular, we consider, as in Fig. 8 photon subtraction from two squeezed modes with squeezing levels $s_1 = 1.5$ dB and $s_2 = -2.6$ dB, for which some loss resilience is expected [56]. In Figs. 11(a)–11(c), we consider the case in which the two photons are subtracted in different modes parameterized by $\Theta_1 = \pi/4$ and $\Theta_2 = -\pi/4$, respectively. We can observe that the witness is positive in all bases for losses lower than 10%. In Fig. 11(d)–11(f), we show the behavior of the witness when both photons are subtracted in the same mode, parameterized by $\Theta = \pi/4$. In this case, the witness becomes negative around 8% loss rate. In both cases, the resilience to losses of the mode-intrinsic entanglement witness is quite low. Nevertheless, it is worth analyzing the values of the witness at different mode basis in the presence of losses. Figures 11(b) and 11(d) show the contour plot of the entanglement witness over a grid of relevant mode basis change at a level of losses of 5%. In this case, both plots have positive values in all mode bases, with very large values in some regions. However, in Figs. 11(c) and 11(f), we show the behavior for a level of losses of 20%. At this point, we can observe that both plots have regions in which the witness is non-positive, while still having positive values in a large set of bases. Even though we are not
<figure><img src="image_13.png" /><figcaption>Fig. 11. Behavior of the metrological witness for two-mode two-photon-subtracted states, under the effect of losses. We consider that the two photons are subtracted from two single-mode squeezed states with squeezing levels $r_1 = 1.5$ dB and $r_2 = -2.6$ dB. In panels (a)–(c), we consider the two photon subtractions to happen in complementary modes ($\Theta_1 = \pi/4$, $\Theta_2 = -\pi/4$), while in panels (d)–(f), we consider the subtractions to happen in the same mode, given by $\Theta_1 = \Theta_2 = \pi/4$. Panels (a) and (d) show the metrological witness derived in Eq. (35), $W_{\text{hom}}(\hat{\rho}, \mathcal{A}_{N=2}|\vartheta_0)$. Panels (b), (d) and (c), (f) show the behavior of the mode-specific entanglement witness, $\max_{\{\phi^{(1)}, \phi^{(2)}\}} E_{\hat{\rho}_{\text{hom}}}^{\mathcal{A}_{N=2}}(\phi^{(i)}, \vartheta|\vartheta_0)$ for the measurement settings $\phi^{(1)} = \{q_1, q_2\}$ and $\phi^{(2)} = \{p_1, p_2\}$ and for all relevant mode basis at levels of losses of 5% and 20%, respectively.</figcaption></figure>

witnessing mode-intrinsic entanglement in this case, which is the central concern of this paper, these plots show the potential of the method we have developed for the less specific task of entanglement detection in different mode bases, even under high optical losses. We can think of it as a way to unveil hidden entanglement measuring in only one basis. In the cases that we considered, the highest values of the witness were actually obtained for bases in which we did not perform any measurement on the state.

The results shown in Fig. 11 show that this method has potential experimental utility, nevertheless, the fact that this is done at second order in the generators makes it far from the reach of current experimental capabilities. In the next subsection, we show a more experimentally friendly approach, by evaluating the criteria at first order in the same kind of states. We will compute all the quantities involved in the criteria from homodyne data, simulated to resemble what we would obtain from an actual experimental realization.

Despite the fact that they cannot detect mode intrinsic entanglement, Gaussian methods are also a natural fit for this more general task of mapping the entanglement properties of a given state over the whole range of possible bases beyond the one in which measurements are performed. In Sec. 1 of Supplement 1, we perform a detailed comparison, for the same states analyzed in Fig. 11, between the results that are obtained using the method we have developed, and those obtained using a Gaussian witness, namely the Duan et al. criterion [34]. We find in this comparison that our method generally outperforms Gaussian witnesses for this task.

## 5.1. Implementation with Simulated Homodyne Data

We consider the simulation of an experimental evaluation of the witness Eq. (35). We should keep in mind that this witness

is a lower bound of that in Eq. (21). We consider in particular the evaluation of the witness relying on two homodyne settings $\phi_1 = \{0, 0\}$ and $\phi_2 = \{\pi/2, \pi/2\}$, which correspond to Alice and Bob measuring simultaneously in the same quadratures $\hat{q}_A, \hat{q}_B$ and $\hat{p}_A, \hat{p}_B$, respectively. Our goal is to explore the relevance of the protocol for currently feasible experiments. For that reason, we restrict ourselves to the evaluation of the criteria for first-order generators, i.e., $\mathcal{A}^{(1)} = \{\hat{q}_A, \hat{q}_B, \hat{p}_A, \hat{p}_B\}$, whose effect can be accounted for in post-processing of the measurement data.

To evaluate the criteria, we sample $10^6$ points from the probability marginals associated with the different measurement settings $\phi^{(i)} = \{\phi_1^{(i)}, \phi_2^{(i)}\}$, to obtain a histogram of relative frequencies $F(\xi_{\phi_1^{(i)}}, \xi_{\phi_2^{(i)}})$ that in the limit of a large number of samples should converge to the exact distribution $p(\xi_{\phi_1^{(i)}}, \xi_{\phi_2^{(i)}})$, obtained through Eq. (30). Once we have all the frequency histograms, we use them to compute the classical Fisher information matrices associated with each of them, which can be done through the use of the Hellinger distance [73], defined as

$$d_H^2(\varphi) = \int_{\xi_{\phi_1}^{(i)}, \xi_{\phi_2}^{(i)}} \left( \sqrt{p(\xi_{\phi_1}^{(i)}, \xi_{\phi_2}^{(i)} | \mathbf{\kappa})} - \sqrt{p(\xi_{\phi_1}^{(i)}, \xi_{\phi_2}^{(i)} | 0)} \right)^2 d\xi_{\phi_1}^{(i)} d\xi_{\phi_2}^{(i)}, \quad (36)$$

where $p(\xi_{\phi_1}, \xi_{\phi_2} | \mathbf{\kappa})$ is the probability distribution after the implementation of the parameter $\mathbf{\kappa}$, which in our case, corresponds to a displacement that can be trivially evaluated following Eq. (31). In the computation of the Hellinger distance, we replace the exact probability distribution by the corresponding histograms and consequently the integral is approximated by a discrete sum. For small values of the parameter the Hellinger distance can be
<figure><img src="image_14.png" /><figcaption>Fig. 12. Simulation of an experimental evaluation of the witness Eq. (35) averaged over 100 samples, each of them containing 1 000 000 data points obtained from homodyne detection. We consider here a two-photon-subtracted state with squeezing levels $r_1 = 1.5$ dB and $r_2 = -2.6$ dB, and subtractions in the same mode $\Theta_1 = \Theta_2 = \pi/4$. Evaluation of the witness $E_{\rho}^{A_{N=1}}$ in the measurement basis (a) $(q_1, q_2)$ [see Fig. 10(a)], (b) $(p_1, p_2)$ [see Fig. 10(d)]. (c) maximal value between both measurement settings $(q_1, q_2)$ and $(p_1, p_2)$ for each mode basis change $(\theta, \varphi)$ of the transformation $U(\theta)$. Panels (a), (b), (c) correspond to zero loss, while panels (d), (e), (f) are obtained from the same measurement settings as the top row but for 10% of losses. The gray areas shown for each contour line correspond to the uncertainties resulting from the 100 samples.</figcaption></figure>

approximated by the quadratic form

$$d_H^2(\kappa) \approx \kappa^T \frac{F}{8} \kappa + O(|\kappa|^3), \qquad (37)$$

where $F$ is the Fisher information matrix associated with the histograms $p(\cdot|\kappa)$.

In Fig. 12, we show the behavior of the entanglement witness on different bases for the same states considered in Figs. 11(d)–11(f). The quantities in each plot were computed over 100 realizations of the sampling process, which allowed to estimate the uncertainties on the values of the witness. These uncertainties are represented through the gray areas surrounding each contour. Figures 12(a)–12(c) correspond to the pure state case. In Figs. 12(d)–12(f), we consider the state after a lossy channel with 10% loss rate. Figures 12(a) and 12(d) show, for both instances, the evaluation of the witness Eq. (33) from measuring simultaneously the momentum quadratures in both modes. In Figs. 12(b) and 12(e), we show the evaluation from the measurement of the position quadratures. In the lossy case, only the position quadrature based witness gives positive results, which is related to the inhomogeneous distribution of squeezing in the Gaussian state from which the photons were subtracted. While we are not able to detect mode-intrinsic entanglement at this level of losses, these results show the potential of the method we have developed to detect entanglement in a wide range of bases, in an experimentally feasible way and from a rather restricted set of measurement settings. This last remark also emphasizes and strengthens the ability to detect entanglement features from a limited set of probability marginals, thereby justifying the high effectiveness of neural network methods that rely on the same type of correlation patterns to identify entanglement [49].

## 6. CONCLUSION

In this paper, we derive a witness capable of detecting mode-intrinsic entanglement, a specific form of entanglement, only exhibited by non-Gaussian states, and a relevant resource in the analysis of the quantum computational advantage for a family of sampling protocols in the context of CV systems [28]. As we demonstrate in this work, our mode-intrinsic entanglement witness based on metrological tools, and more specifically on the quantum Fisher information, proves to be efficient when applied to experimentally relevant states such as multimode photon-subtracted states.

Moreover, we underline the capability of the witness to be experimentally accessible. First, it only requires measurements in one mode basis. Once the (quantum) Fisher information has been measured in one basis, one can check entanglement in any other basis in post-processing. Second, while the witness relies on the quantum Fisher information, a quantity very difficult, if not impossible, to access experimentally, one can obtain a lower bound of the witness using the classical Fisher information associated with homodyne measurements. As we explain in Sec. 5, the witness can be evaluated through different homodyne measurement settings and produce good results. We highlight and check its experimental applicability with simulations of experimental data obtained from probability marginals.

Overall, what we propose is a new and promising tool for exploring the specificity of quantum correlations exhibited by non-Gaussian states. Even if the witness is sensitive to losses, thus limiting its capability to certify mode-intrinsic entanglement, it works on a wide range of mode bases, surpassing the efficiency of Gaussian witnesses in most cases. This makes it a tool of relevance to unveil, in post-processing, entanglement on multiplexed optical systems beyond their preferred mode basis. This is particularly relevant in systems where some mode
bases are much easier to access than others. Think, for example, of spectral modes where frequency bands are much easier to measure than other mode bases.

In future studies, it would be desirable to find methods to mitigate the effect of losses; one could, for example, complement our witness with other entanglement criteria in a specific basis where it fails to detect entanglement. With this new tool at hand, we envisage to further explore the effect of non-Gaussian operations, like photon subtraction, on the entanglement properties of large cluster states, complementing the results presented in Sec. 4.1. Moreover, it remains to better understand the interplay between entanglement and the metrological advantage that non-Gaussian states offer; or, in other words, to which extent states that strongly violate the metrological witness can be useful for some applications in metrology.

**Funding.** Agence Nationale de la Recherche (ANR-21-CE47-0005, ANR-22-PETQ-0013); HORIZON EUROPE European Innovation Council (101114899); Horizon 2020 Framework Programme (101017733).

**Acknowledgment.** We thank M. Frigerio for inspiring discussions.

**Disclosures.** The authors declare no conflicts of interest.

**Data availability.** Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

**Supplemental document.** See Supplement 1 for supporting content.

## REFERENCES

1. M. Steffen, D. P. DiVincenzo, J. M. Chow, et al., "Quantum computing: An IBM perspective," IBM J. Res. Dev. 55, 13:1–13:11 (2011).

2. R. Acharya, "Suppressing quantum errors by scaling a surface code logical qubit," Nature 614, 676–681 (2023).

3. L. Henriet, L. Beguin, A. Signoles, et al., "Quantum computing with neutral atoms," Quantum 4, 327 (2020).

4. S. J. Evered, D. Bluvstein, M. Kalinowski, et al., "High-fidelity parallel entangling gates on a neutral-atom quantum computer," Nature 622, 268–272 (2023).

5. D. Kielpinski, C. Monroe, and D. J. Wineland, "Architecture for a large-scale ion-trap quantum computer," Nature 417, 709–711 (2002).

6. E. Gouzien, D. Ruiz, F.-M. Le Régent, et al., "Performance analysis of a repetition cat code architecture: Computing 256-bit elliptic curve logarithm in 9 hours with 126 133 cat qubits," Phys. Rev. Lett. 131, 040602 (2023).

7. S. Author Bartolucci, P. Birchall, H. Bombín, et al., "Fusion-based quantum computation," Nat. Commun. 14, 912 (2023).

8. J. E. Bourassa, R. N. Alexander, M. Vasmer, et al., "Blueprint for a scalable photonic fault-tolerant quantum computer," Quantum 5, 392 (2021).

9. M. V. Larsen, X. Guo, C. R. Breum, et al., "Deterministic generation of a two-dimensional cluster state," Science 366, 369–372 (2019).

10. W. Asavanant, Y. Shiozawa, S. Yokoyama, et al., "Generation of time-domain-multiplexed two-dimensional cluster state," Science 366, 373–376 (2019).

11. A. Einstein, B. Podolsky, and N. Rosen, "Can quantum-mechanical description of physical reality be considered complete?" Phys. Rev. 47, 777–780 (1935).

12. E. Schrödinger, "Discussion of probability relations between separated systems," Math. Proc. Camb. Phil. Soc. 31, 555–563 (1935).

13. R. Horodecki, P. Horodecki, M. Horodecki, et al., "Quantum entanglement," Rev. Mod. Phys. 81, 865–942 (2009).

14. A. Acín, I. Bloch, H. Buhrman, et al., "The quantum technologies roadmap: a European community view," New J. Phys. 20, 080201 (2018).

15. R. Jozsa and N. Linden, "On the role of entanglement in quantum-computational speed-up," Proc. R. Soc. Lond. A. 459, 2011–2032 (2003).

16. S. D. Bartlett, B. C. Sanders, S. L. Braunstein, et al., "Efficient classical simulation of continuous variable quantum information processes," Phys. Rev. Lett. 88, 097904 (2002).

17. A. Mari and J. Eisert, "Positive wigner functions render classical simulation of quantum computation efficient," Phys. Rev. Lett. 109, 230503 (2012).

18. Y.-S. Ra, A. Dufour, M. Walschaers, et al., "Non-Gaussian quantum states of a multimode light field," Nat. Phys. 16, 144–147 (2020).

19. M. Endo, R. He, T. Sonoyama, et al., "Non-Gaussian quantum state generation by multi-photon subtraction at the telecommunication wavelength," Opt. Express 31, 12865–12879 (2023).

20. J. S. Neergaard-Nielsen, B. M. Nielsen, C. Hettich, et al., "Generation of a superposition of odd photon number states for quantum information networks," Phys. Rev. Lett. 97, 083604 (2006).

21. K. Wakui, H. Takahashi, A. Furusawa, et al., "Photon subtracted squeezed states generated with periodically poled KTiOPO₄," Opt. Express 15, 3568–3574 (2007).

22. S. Konno, W. Asavanant, F. Hanamura, et al., "Logical states for fault-tolerant quantum computation with propagating light," Science 383, 289–293 (2024).

23. M. Kudra, M. Kervinen, I. Strandberg, et al., "Robust preparation of Wigner-negative states with optimized snap-displacement sequences," PRX Quantum 3, 030301 (2022).

24. S. Touzard, A. Grimm, Z. Leghtas, et al., "Coherent oscillations inside a quantum manifold stabilized by dissipation," Phys. Rev. X 8, 021005 (2018).

25. C. Flühmann, T. L. Nguyen, M. Marinelli, et al., "Encoding a qubit in a trapped-ion mechanical oscillator," Nature 566, 513–517 (2019).

26. L. Podhora, L. Lachman, T. Pham, et al., "Quantum non-Gaussianity of multiphonon states of a single atom," Phys. Rev. Lett. 129, 013602 (2022).

27. A. P. Reed, K. H. Mayer, J. D. Teufel, et al., "Faithful conversion of propagating quantum information to mechanical motion," Nat. Phys. 13, 1163–1167 (2017).

28. U. Chabaud and M. Walschaers, "Resources for bosonic quantum computational advantage," Phys. Rev. Lett. 130, 090602 (2023).

29. J. Sperling, A. Perez-Leija, K. Busch, et al., "Mode-independent quantum entanglement for light," Phys. Rev. A 100, 062129 (2019).

30. U. Chabaud and M. Mehraban, "Holomorphic representation of quantum computations," Quantum 6, 831 (2022).

31. M. Walschaers, C. Fabre, V. Parigi, et al., "Statistical signatures of multimode single-photon-added and subtracted states of light," Phys. Rev. A 96, 053835 (2017).

32. We prefer this terminology over mode-independent quantum entanglement, used in Ref. [29], because we consider the latter could mislead the reader to think that the amount of entanglement in all basis is the same, which is generally not the case.

33. R. Simon, "Peres-Horodecki separability criterion for continuous variable systems," Phys. Rev. Lett. 84, 2726–2729 (2000).

34. L.-M. Duan, G. Giedke, J. I. Cirac, et al., "Inseparability criterion for continuous variable systems," Phys. Rev. Lett. 84, 2722–2725 (2000).

35. V. Giovannetti, S. Mancini, D. Vitali, et al., "Characterizing the entanglement of bipartite quantum systems," Phys. Rev. A 67, 022320 (2003).

36. P. van Loock and A. Furusawa, "Detecting genuine multipartite continuous-variable entanglement," Phys. Rev. A 67, 052315 (2003).

37. P. Abiuso, S. Bäuml, D. Cavalcanti, et al., "Measurement-device-independent entanglement detection for continuous-variable systems," Phys. Rev. Lett. 126, 190502 (2021).

38. J. Williamson, "On the algebraic problem concerning the normal forms of linear dynamical systems," American Journal of Mathematics 58, 141–163 (1936).

39. S. L. Braunstein, "Squeezing as an irreducible resource," Phys. Rev. A 71, 055801 (2005).

40. A. Peres, "Separability criterion for density matrices," Phys. Rev. Lett. 77, 1413–1415 (1996).

41. A. I. Lvovsky and M. G. Raymer, "Continuous-variable optical quantum-state tomography," Rev. Mod. Phys. 81, 299–332 (2009).

42. E. Shchukin and W. Vogel, "Inseparability criteria for continuous bipartite quantum states," Phys. Rev. Lett. 95, 230502 (2005).
43. C. Rodó, G. Adesso, and A. Sanpera, "Operational quantification of continuous-variable correlations," Phys. Rev. Lett. **100**, 110505 (2008).

44. S. P. Walborn, A. Salles, R. M. Gomes, et al., "Revealing hidden Einstein-Podolsky-Rosen nonlocality," Phys. Rev. Lett. **106**, 130402 (2011).

45. A. Saboia, F. Toscano, and S. P. Walborn, "Family of continuous-variable entanglement criteria using general entropy functions," Phys. Rev. A **83**, 032307 (2011).

46. M. Gessner, L. Pezzè, and A. Smerzi, "Efficient entanglement criteria for discrete, continuous, and hybrid variables," Phys. Rev. A **94**, 020101 (2016).

47. B. Yadin, M. Fadel, and M. Gessner, "Metrological complementarity reveals the Einstein-Podolsky-Rosen paradox," Nat. Commun. **12**, 2410 (2021).

48. M. Gårttner, T. Haas, and J. Noll, "General class of continuous variable entanglement criteria," Phys. Rev. Lett. **131**, 150201 (2023).

49. X. Gao, M. Isoard, F. Sun, et al., "Correlation-pattern-based continuous variable entanglement detection through neural networks," Phys. Rev. Lett. **132**, 220202 (2024).

50. A. A. Valido, F. Levi, and F. Mintert, "Hierarchies of multipartite entanglement for continuous-variable states," Phys. Rev. A **90**, 052321 (2014).

51. V. Giovannetti, S. Lloyd, and L. Maccone, "Advances in quantum metrology," Nat. Photonics **5**, 222–229 (2011).

52. L. Pezzè, A. Smerzi, M. K. Oberthaler, et al., "Quantum metrology with nonclassical states of atomic ensembles," Rev. Mod. Phys. **90**, 035005 (2018).

53. G. Tóth and I. Apellaniz, "Quantum metrology from a quantum information science perspective," J. Phys. A: Math. Theor. **47**, 424006 (2014).

54. J. Liu, H. Yuan, X.-M. Lu, et al., "Quantum Fisher information matrix and multiparameter estimation," J. Phys. A: Math. Theor. **53**, 023001 (2019).

55. C. E. Lopetegui, M. Gessner, M. Fadel, et al., "Homodyne detection of non-Gaussian quantum steering," PRX Quantum **3**, 030347 (2022).

56. D. Barral, M. Isoard, G. Sorelli, et al., "Metrological detection of purely-non-Gaussian entanglement," arXiv (2023).

57. N. Lupu-Gladstein, Y. B. Yilmaz, D. R. M. Arvidsson-Shukur, et al., "Negative quasiprobabilities enhance phase estimation in quantum-optics experiment," Phys. Rev. Lett. **128**, 220504 (2022).

58. C. Fabre and N. Treps, "Modes and states in quantum optics," Rev. Mod. Phys. **92**, 035005 (2020).

59. U. Leonhardt and H. Paul, "Measuring the quantum state of light," Progress in Quantum Electronics **19**, 89–130 (1995).

60. W. R. Clements, P. C. Humphreys, B. J. Metcalf, et al., "Optimal design for universal multiport interferometers," Optica **3**, 1460–1465 (2016).

61. M. Walschaers, "Non-Gaussian quantum states and where to find them," PRX Quantum **2**, 030204 (2021).

62. M. Walschaers, Y.-S. Ra, and N. Treps, "Mode-dependent-loss model for multimode photon-subtracted states," Phys. Rev. A **100**, 023828 (2019).

63. P. van Loock, C. Weedbrook, and M. Gu, "Building Gaussian cluster states by linear optics," Phys. Rev. A **76**, 032321 (2007).

64. G. Ferrini, J. Roslund, F. Arzani, et al., "Optimization of networks for measurement-based quantum computation," Phys. Rev. A **91**, 032314 (2015).

65. F. Sansavini and V. Parigi, "Continuous variables graph states shaped as complex networks: Optimization and manipulation," Entropy **22**, 1 (2020).

66. R. Medeiros de Araújo, J. Roslund, Y. Cai, et al., "Full characterization of a highly multimode entangled state embedded in an optical frequency comb using pulse shaping," Phys. Rev. A **89**, 053828 (2014).

67. M. Walschaers, B. Sundar, N. Treps, et al., "Emergent complex quantum networks in continuous-variables non-Gaussian states," Quantum Sci. Technol. **8**, 035009 (2023).

68. M. Eaton, C. González-Arciniegas, R. N. Alexander, et al., "Measurement-based generation and preservation of cat and grid states within a continuous-variable cluster state," Quantum **6**, 769 (2022).

69. X. Sun, Y. Wang, Y. Tian, et al., "Deterministic and universal quantum squeezing gate with a teleportation-like protocol," Laser Photonics Rev. **16**, 2100329 (2022).

70. J. Zhao, K. Liu, H. Jeng, et al., "A high-fidelity heralded quantum squeezing gate," Nat. Photonics **14**, 306–309 (2020).

71. M. V. Larsen, X. Guo, C. R. Breum, et al., "Deterministic multi-mode gates on a scalable photonic quantum computing platform," Nat. Phys. **17**, 1018–1023 (2021).

72. G. Ferrini, J. P. Gazeau, T. Coudreau, et al., "Compact Gaussian quantum computation by multi-pixel homodyne detection," New J. Phys. **15**, 093015 (2013).

73. S. L. Braunstein and C. M. Caves, "Statistical distance and the geometry of quantum states," Phys. Rev. Lett. **72**, 3439–3443 (1994).