Detecting entanglement of unknown continuous variable states with random measurements

The most valuable characteristics of quantum systems are quantum correlations such as entanglement, which represents a useful resource for applications unattainable in the framework of classical theory, such as quantum teleportation, quantum cryptography and dense coding [1, 2]. In the last two decades, the use of continuous variable (CV) systems became a very powerful approach to quantum information processing and communication [3–5], opening the way to various protocols and tasks, like quantum cryptography [6], quantum teleportation [7], quantum dense coding [8]. CV quantum systems provide the quantum description of the propagating electromagnetic field, and therefore manifest particular relevance for quantum communication and quantum techniques like detection, imaging and sensing [9].

A significant problem in quantum information theory and any application is to efficiently reveal the properties of an unknown quantum state, in particular to certify the presence of entanglement in a given unknown state. Usual entanglement criteria for CV systems consist in certain operations on the second moments, or uncertainties, of quantum states, such as the positive partial transpose (PPT) criterion [10, 11]. Therefore, first a full tomography is required for completely unknown states, in order to reconstruct the entire covariance matrix (CM). However, this method may be a very resource-consuming and demanding experimental procedure, especially for quantum states with a high number of modes. In addition, the full information about the second moments of the state can be excessive for the characterisation of entanglement present in the state. Instead, one can choose to measure certain fixed combinations of second moments, giving rise to a specific test, which detects entanglement in some states and does not detect it in others [12–14].

Entanglement witnesses (EWs) represent another commonly used entanglement test, being directly accessible in experiments through measurable observables [15]. A Hermitian observable W is an EW if for all separable states ρ_s , Tr[Wρ_s ] ⩾ 0 holds, while for some entangled state ρ we have that Tr[Wρ] < 0 [15]. For CV systems a special instance of EWs can be defined, which embodies the entanglement criterion in terms of the variances of the canonical observables of the state [16–19]. This is possible because all the relevant information on the quantum entanglement of a CV system is encoded in its CM, specifically in its symplectic spectrum [20]. Typically, EWs are employed when certain knowledge about the state is available.

Given an unknown quantum state, however, the complexity of the state and the absence of any information about it deprive us of a specific experimental strategy in tackling the problem of efficient entanglement detection. Therefore, the best strategy in this case would be to perform random measurements, serving as building blocks for the construction of an EW by means of a semidefinite optimization algorithm. In this context, one has the opportunity to represent the performance of entanglement detection in terms of the number of measurements required for this task.

We introduce new linear constraints for EWs and use them to build a semidefinite program (SDP). Our method based on EWs is applicable to arbitrary unknown CV states, and can be easily adjusted to detect any type of entanglement encoded in the second moments, including bound entanglement. Important is that we find a high probability of success in detecting entanglement with fewer measurements than in a tomographically complete set. Our method presents a clear advantage over entanglement criteria subjected to specific measurement directions on the CM, or to full information about the CM [11, 12]. As this is of high interest for experimental applications, we also provide a statistical analysis showing the robustness of this method to statistical errors. This idea is inspired by an analogous method for the discrete-variable case, which was developed in [21].

The paper is organized as follows. In section 2 we present the theoretical framework of CV states, mainly based on the second moments description of the state. In section 3 we introduce the EWs based on the CM of the state, as presented in references [16, 17], and propose a set of stronger linear semidefinite constraints in order to characterize the EWs. Then, we simulate random homodyne measurements for two-mode CMs in section 4 and formulate a SDP optimizing the witness constructible from given experimental data. We present the results of the efficiency of entanglement detection for random two-mode CMs and, in particular, for the class of squeezed vacuum states in section 5, and illustrate an example of bipartite bound EW. A statistical analysis of our method is provided in section 6, and a summary and conclusions are presented in section 7.

A CV system of N canonical bosonic modes, like the quantized electromagnetic field with a Hamiltonian of a system of N harmonic oscillators (modes), is defined in a Hilbert space $\mathcal{H}=\underset{k=1}{\overset{N}{\oplus }}{\mathcal{H}}_{k}$, each one with an infinite-dimensional space ${\mathcal{H}}_{k}={L}^{2}\left(\mathbb{R}\right)$ and two canonical observables ${\hat{x}}_{k}$ and ${\hat{p}}_{k}$, with the corresponding phase space variables of position x_k and momentum p_k [3–5]. One can define a vector of quadrature operators ${\hat{R}}^{\text{T}}\equiv \left({\hat{R}}_{1},\dots ,{\hat{R}}_{2N}\right)=\left({\hat{x}}_{1},{\hat{p}}_{1},\dots ,{\hat{x}}_{N},{\hat{p}}_{N}\right)$ satisfying the bosonic commutation relations

$$\begin{equation}\left[{\hat{R}}_{i},{\hat{R}}_{j}\right]=\mathrm{i}{{\Omega}}_{ij}\hat{\text{I}},\quad i,j=1,\dots ,2N,\end{equation} \tag{ 1 }$$

where $\hat{\text{I}}$ is the identity matrix, and Ω_ij are the elements of the fundamental symplectic matrix (we assume ℏ = 1)

$$\begin{equation}{{\Omega}}_{N}=\underset{1}{\overset{N}{\oplus }}\left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 0\hfill \end{matrix}\right).\end{equation} \tag{ 2 }$$

The primary role in this study is played by the statistical moments of the quadrature operators, that characterize the state with density operator ρ [3, 4], up to the second order: the displacement vector, which is the real vector d of first order moments ${d}_{i}={\langle {\hat{R}}_{i}\rangle }_{\rho }=\mathrm{Tr}\left[{\hat{R}}_{i}\rho \right]$, and the CM, which is the real, symmetric matrix γ whose entries are the second order moments in symmetrized form (the variances) of the quadrature operators, defined as [5, 22]:

$$\begin{equation}{\gamma }_{ij}={\langle {\left\{{\hat{R}}_{i}-\langle {\hat{R}}_{i}\rangle ,{\hat{R}}_{j}-\langle {\hat{R}}_{j}\rangle \right\}}_{+}\rangle }_{\rho },\end{equation} \tag{ 3 }$$

where {,}₊ represents the anticommutator. The Robertson–Schrödinger uncertainty relation in terms of the CM reads

$$\begin{equation}\gamma +\mathrm{i}{{\Omega}}_{N}{\geqslant}0,\end{equation} \tag{ 4 }$$

assuring that it is a CM of a physical quantum state. Gaussian states represent the class of CV states which are completely characterized by their first and second moments. The entanglement criteria discussed in this paper can also detect entanglement of non-Gaussian states.

A quantum state of a bipartite system is entangled if it cannot be prepared by means of operations acting locally on the subsystems. In the case of separable states correlations are attributed to possible classical communication between subsystems, and hence are of classical origin. This reasoning carries over to CV systems, where a separability criterion can be defined in terms of CMs. If a CM γ of a state of N modes is fully separable, then there exist CMs γ_i , i = 1, ..., N, corresponding to N subsystems, such that [23]

$$\begin{equation}\gamma {\geqslant}\underset{i=1}{\overset{N}{\oplus }}{\gamma }_{i}.\end{equation} \tag{ 5 }$$

Conversely, if this holds, then Gaussian states with CM γ are separable. Therefore, if this criterion is violated, then the corresponding state is entangled, irrespective of whether it is Gaussian or not. If it is not violated, then a Gaussian state is separable, while a non-Gaussian state might be entangled.

In the following, we will refer to the situation of a k-partition of an N-mode system as the splitting or distribution of an N-mode system into k subsystems, where every subsystem j (j = 1, ..., k) is composed of N_j modes, such that ${\sum }_{j=1}^{k}{N}_{j}=N$.

2.1. Symplectic transformations

Unitary operators acting on the quantum state space are equivalent to symplectic transformations which preserve the commutation relations of canonical variables. The real symplectic group is defined by [24]:

$$\begin{equation}Sp\left(2N,\mathbb{R}\right)=\left\{S\in \mathcal{M}\left(2N,\mathbb{R}\right):S{{\Omega}}_{N}{S}^{\text{T}}={{\Omega}}_{N}\right\},\end{equation} \tag{ 6 }$$

where S is a symplectic transformation acting in phase space as $\hat{R}\to {\hat{R}}^{\prime }=S\hat{R}$, and $\mathcal{M}\left(2N,\mathbb{R}\right)$ denotes the set of 2N × 2N real matrices. Symplectic transformations act by congruence on CMs: γ' = SγS^T.

Every symplectic transformation can be decomposed using the Euler decomposition, which represents the singular value decomposition for real symplectic matrices [24]:

$$\begin{equation}S=K\left[\underset{i=1}{\overset{N}{\oplus }}S\left({r}_{i}\right)\right]L,\end{equation} \tag{ 7 }$$

where K, L are symplectic and orthogonal matrices, while

$$\begin{equation}S\left({r}_{i}\right)=\left(\begin{matrix}\hfill {\mathrm{e}}^{-{r}_{i}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\mathrm{e}}^{{r}_{i}}\hfill \end{matrix}\right)\end{equation} \tag{ 8 }$$

are one-mode squeezing matrices (symplectic and nonorthogonal) with r_i the squeezing parameter. The symplectic and orthogonal matrices form the maximal compact subgroup K(N) within the noncompact group $Sp\left(2N,\mathbb{R}\right)$ [24]. The group K(N) is isomorphic to the group U(N) of N × N complex unitary matrices:

$$\begin{equation}K\left(N\right)=\left\{S\left(X,Y\right)\vert X-\mathrm{i}Y\in U\left(N\right)\right\},\end{equation} \tag{ 9 }$$

where the corresponding symplectic matrices are of the following form:

$$\begin{equation}S\left(X,Y\right)=\left(\begin{matrix}\hfill X\hfill & \hfill Y\hfill \\ \hfill -Y\hfill & \hfill X\hfill \end{matrix}\right)\in Sp\left(2N,\mathbb{R}\right).\end{equation} \tag{ 10 }$$

Such transformations describe multiport interferometers and are called passive canonical unitaries, which preserve the photon number [4]. The active canonical unitaries correspond to nonorthogonal symplectic transformations, such as one-mode squeezers.

In the following we will use the theorem by Williamson [20], according to which every matrix $M\in \mathcal{M}\left(2N,\mathbb{R}\right)$, M ⩾ 0 can be brought to a diagonal form through symplectic transformations:

$$\begin{equation}SM{S}^{\mathrm{T}}=\mathrm{diag}\left({s}_{1},{s}_{1},\dots ,{s}_{N},{s}_{N}\right),\end{equation} \tag{ 11 }$$

where s₁, ..., s_N ⩾ 0 are called symplectic eigenvalues of M. By

$$\begin{equation}\text{str}\left[M\right]{:=}\sum _{i=1}^{N}{s}_{i}\end{equation} \tag{ 12 }$$

we will denote the symplectic trace of M.

An EW based on CMs is characterised by a real symmetric matrix Z ⩾ 0 such that [17]

$$\begin{equation}\mathrm{Tr}\left[Z{\gamma }_{\mathrm{s}}\right]{\geqslant}1,\quad \text{for}\enspace \text{all}\;\text{separable}\;\text{CMs}\quad {\gamma }_{\mathrm{s}},\end{equation} \tag{ 13a }$$

$$\begin{equation}\mathrm{Tr}\left[Z\gamma \right]{< }1,\quad \text{for}\enspace \text{some}\;\text{entangled}\;\text{CM}\quad \gamma .\end{equation} \tag{ 13b }$$

The EWs based on second moments defined in equations (13a) and (13b) represent hyperplanes in the space of CMs that separate some entangled states from the set of separable CMs. If there exists Z which fulfills conditions (13b), then the state with CM γ is entangled, irrespective of whether it is Gaussian or not, while if this test does not detect entanglement in a given non-Gaussian state, then the result is inconclusive. The following theorem fully characterises the EWs for multimode CV states defined in equations (13a) and (13b) for different entanglement classes.

Theorem . (Taken from [16, 17]) A CM γ of a k-partite system with ${\sum }_{j=1}^{k}{N}_{j}=N$ modes is entangled with respect to this partition iff there exists Z such that

$$\begin{equation}\mathrm{Tr}\left[Z\gamma \right]{< }1,\end{equation} \tag{ 14 }$$

where Z is a real, symmetric 2N × 2N matrix satisfying

$$\begin{align}\hfill & Z{\geqslant}0,\quad \sum _{j=1}^{k}\mathrm{str}\left[{Z}_{j}\right]{\geqslant}\frac{1}{2},\hfill \end{align} \tag{ 15 }$$

where Z_j is the block matrix on the diagonal of Z acting on the subsystem j. Matrices Z are called EWs based on second moments.

Due to the convexity of the set of separable CMs there always exists an EW Z giving the result of equation (14) for an entangled CM γ. In reference [17] the authors formulated the above theorem in a slightly different way: in addition to equation (15) it is stated that such an EW has to satisfy also $\text{str}\left[Z\right]{< }\frac{1}{2}$, instead of condition Tr[Zγ] < 1. Note that there is no contradiction between the conditions (15) and $\text{str}\left[Z\right]{< }\frac{1}{2}$ that an EW has to satisfy, since the relation $\text{str}\left[Z\right]{\leqslant}{\sum }_{i=1}^{N}\text{str}\left[{Z}_{i}\right]$ holds [16].

Nevertheless, the two formulations of the theorem above are equivalent. In order to show this we will use the results from reference [16] where it is proven that Tr[Zγ] ⩾ 1 for all separable CMs γ if and only if Z ⩾ 0 and ${\sum }_{j=1}^{k}\text{str}\left[{Z}_{j}\right]{\geqslant}\frac{1}{2}$. In addition, it is shown that Tr[Zγ] ⩾ 1 for all CMs γ if and only if Z ⩾ 0 and $\text{str}\left[Z\right]{\geqslant}\frac{1}{2}$. As Z ≱ 0 would contradict Tr[Zγ] ⩾ 1 for all separable CMs γ, it follows that if Tr[Zγ] < 1 for some CM γ then $\text{str}\left[Z\right]{< }\frac{1}{2}$. Conversely, if $\text{str}\left[Z\right]{< }\frac{1}{2}$ then there exists some CM γ such that Tr[Zγ] < 1.

The problem of finding an EW that most robustly detects entanglement in a given CM arises as a SDP (see reference [17] where the authors provide also numerical routines performing this task). Here we consider the situation when no information about the state is available, and we aim at constructing the EWs from given random measurements. For this purpose, the description of EWs given in the theorem above can serve as constraints in our optimization program.

However, the inequality (15) cannot be used in this form as a semidefinite constraint in an SDP, because its left-hand side cannot be regarded as a linear function since the symplectic eigenvalues of a matrix M ⩾ 0 are given by the eigenvalues of the matrix ${M}^{\frac{1}{2}}\left(\mathrm{i}{{\Omega}}_{N}\right){M}^{\frac{1}{2}}$ [5]. In the following, we propose a set of linear semidefinite constraints for EWs, which are stronger than conditions (15).

Proposition . For the EW Z of a k-partite entangled N-mode CM with ${\sum }_{j=1}^{k}{N}_{j}=N$, the inequalities (15) are satisfied if the following conditions are fulfilled:

$$\begin{align}\hfill & \qquad \qquad Z{\geqslant}0,\hfill \\ \hfill & {Z}_{j}+\mathrm{i}\frac{{x}_{j}}{{N}_{j}}{{\Omega}}_{{N}_{j}}{\geqslant}0,\quad {x}_{j}\in \mathbb{R},\quad j=1,\dots ,k-1,\hfill \\ \hfill & {Z}_{k}+\mathrm{i}\frac{1}{{N}_{k}}\left(\frac{1}{2}-\sum _{j=1}^{k-1}{x}_{j}\right){{\Omega}}_{{N}_{k}}{\geqslant}0.\hfill \end{align} \tag{ 16 }$$

Proof. 

• (a)  
  In the first part we prove the proposition for k = N. First, we prove this for N = 2, i.e. for a two-mode witness with the following block form:

  $$\begin{equation}Z=\left(\begin{matrix}\hfill {Z}_{1}\hfill & \hfill {Z}_{c}\hfill \\ \hfill {Z}_{c}^{\text{T}}\hfill & \hfill {Z}_{2}\hfill \end{matrix}\right),\end{equation} \tag{ 17 }$$

  where Z₁, Z₂ and Z_c are 2 × 2 matrices. Since Z is a positive semidefinite matrix, also the principal submatrices Z₁ and Z₂ are positive semidefinite. Let us assume the following inequality:

  $$\begin{equation}{Z}_{1}+\mathrm{i}x{{\Omega}}_{1}{\geqslant}0,\quad \text{where}\quad {{\Omega}}_{1}=\left(\begin{matrix}\hfill 0\hfill & \hfill 1\hfill \\ \hfill -1\hfill & \hfill 0\hfill \end{matrix}\right),\quad x\in \mathbb{R}.\end{equation} \tag{ 18 }$$

  By symplectic transformation S the positive matrix above can be brought to the Williamson normal form as follows ⁴ :

  $$\begin{equation}S\left({Z}_{1}+\mathrm{i}x{{\Omega}}_{1}\right){S}^{\mathrm{T}}={Z}_{1}^{w}+\mathrm{i}x{{\Omega}}_{1}=\left(\begin{matrix}\hfill {z}_{1}\hfill & \hfill x\mathrm{i}\hfill \\ \hfill -x\mathrm{i}\hfill & \hfill {z}_{1}\hfill \end{matrix}\right),\end{equation} \tag{ 19 }$$

  where ${Z}_{1}^{w}=\mathrm{diag}\left({z}_{1},{z}_{1}\right)$, with z₁ the positive symplectic eigenvalue of Z₁. The eigenvalues α of matrix (19) are determined from the equation:

  $$\begin{equation}{\left({z}_{1}-\alpha \right)}^{2}-{x}^{2}=\left({z}_{1}-\alpha -x\right)\left({z}_{1}-\alpha +x\right)=0,\end{equation} \tag{ 20 }$$

  and hence

  $$\begin{equation}{z}_{1}{\pm}x=\alpha {\geqslant}0.\end{equation} \tag{ 21 }$$

  Thus, the symplectic eigenvalue z₁ fulfills the inequality z₁ ⩾ ±x, or z₁ ⩾ |x|. A similar inequality can be formulated for the block matrix Z₂:

  $$\begin{equation}{Z}_{2}+\mathrm{i}\left(\frac{1}{2}-x\right){{\Omega}}_{1}{\geqslant}0,\end{equation} \tag{ 22 }$$

  from which we obtain the following condition for the symplectic eigenvalue z₂:

  $$\begin{equation}{z}_{2}{\geqslant}\left\vert \frac{1}{2}-x\right\vert .\end{equation} \tag{ 23 }$$

  Now, the sum of symplectic eigenvalues gives:

  $$\begin{equation}{z}_{1}+{z}_{2}{\geqslant}\vert x\vert +\left\vert \frac{1}{2}-x\right\vert {\geqslant}\left\vert x+\frac{1}{2}-x\right\vert =\frac{1}{2}.\end{equation} \tag{ 24 }$$

  The above inequality assures that the condition (15) is always fulfilled. The generalization to more modes is straightforward. For instance, consider a three-mode CM and we want an EW detecting three-partite entanglement. Then, according to the proposition, we need to impose constraints on the three block diagonal matrices of the witness, which amount to the following inequalities for the corresponding symplectic eigenvalues:

  $$\begin{align}\hfill & {z}_{1}{\geqslant}\vert {x}_{1}\vert ,\quad {x}_{1}\in \mathbb{R},\hfill \\ \hfill & {z}_{2}{\geqslant}\vert {x}_{2}\vert ,\quad {x}_{2}\in \mathbb{R},\hfill \\ \hfill & {z}_{3}{\geqslant}\left\vert \frac{1}{2}-{x}_{1}-{x}_{2}\right\vert .\hfill \end{align} \tag{ 25 }$$

  These inequalities imply the constraint (15).
• (b)  
  In the second part, we present the generalization of the proof for k-partite entanglement of N-mode CMs, with k < N. Consider, for simplicity, a three-mode state and the bipartition between the first and the other two modes. The witness Z is a 6 × 6 matrix where Z₁ is the 2 × 2 block diagonal matrix of Z acting on the first mode, and we denote by Z' the 4 × 4 block matrix acting on the other two modes. Then the corresponding constraints on the witness are:

  $$\begin{align}\hfill & \qquad \quad Z{\geqslant}0,\hfill \\ \hfill & {Z}_{1}+\mathrm{i}x{{\Omega}}_{1}{\geqslant}0,\quad x\in \mathbb{R},\hfill \\ \hfill & {Z}^{\prime }+\mathrm{i}\frac{1}{2}\left(\frac{1}{2}-x\right){{\Omega}}_{2}{\geqslant}0.\hfill \end{align} \tag{ 26 }$$

  If we denote by z₁ the symplectic eigenvalue of Z₁, and by z₁', z₂' the two symplectic eigenvalues of Z', then the conditions above are equivalent to:

  $$\begin{align}\hfill & {z}_{1}{\geqslant}\vert x\vert ,\quad x\in \mathbb{R},\hfill \\ \hfill & {z}_{1}^{\prime }{\geqslant}\frac{1}{2}\left\vert \frac{1}{2}-x\right\vert ,\hfill \\ \hfill & {z}_{2}^{\prime }{\geqslant}\frac{1}{2}\left\vert \frac{1}{2}-x\right\vert ,\hfill \end{align} \tag{ 27 }$$

  which satisfy the condition (15). Since this is a bipartite state, the lower bounds for the three symplectic eigenvalues depend on a single parameter x, while, according to the proposition, the detection of three-partite entanglement in a three-mode state would require two optimization parameters, x₁ and x₂ (see equation (25)). The generalization of the proof to N modes and k parties is straightforward. Note that the conditions in the proposition are stronger for k-partite entanglement (with k < N) than for genuine multipartite entanglement (i.e. k = N), where the optimization for every symplectic eigenvalue is done independently.

□

While the semidefinite inequalities proposed in the previous proposition present the advantage of being linear, the drawback of these constraints is that they are stronger than those required by the theorem characterizing the EWs based on second moments, and therefore some EWs will not satisfy conditions (16).

Here we will shortly present the physical set-up of the homodyne detection, that encodes the experimental settings measuring the variances of the state. Homodyne measurements are phase sensitive measurements which allow the detection of the moments of quadratures up to the second order [3, 4]. We denote by $\hat{k}$ and ${\hat{k}}^{{\dagger}}$ the mode operators of our state. A simple scheme for balanced homodyne measurements is composed of a balanced beam splitter superposing the signal mode to be measured $\hat{k}$ with a strong local oscillator field α_LO = |α_LO|e^iθ with phase θ, and two photon detectors, converting the electromagnetic modes into two output photon currents, i₁ and i₂. The actual quantity to be measured is the difference in the photon currents, given by:

$$\begin{equation}\delta i={i}_{1}-{i}_{2}=q\vert {\alpha }_{\text{LO}}\vert \enspace \langle {\hat{x}}_{\theta }\rangle ,\end{equation} \tag{ 28 }$$

with q being a constant, and ${\hat{x}}_{\theta }$ is the generalized quadrature operator of mode $\hat{k}$ defined as:

$$\begin{equation}{\hat{x}}_{\theta }=\frac{\mathrm{exp}\left(-\text{i}\theta \right)\hat{k}+\mathrm{exp}\left(\text{i}\theta \right){\hat{k}}^{{\dagger}}}{\sqrt{2}},\end{equation} \tag{ 29 }$$

which covers the whole continuum of quadratures for θ ∈ [0, π]. It was shown in reference [26] that in the strong local oscillator limit the homodyne detection performs projective measurements corresponding to the positive operator valued measure {|x_θ ⟩⟨x_θ |}, where |x_θ ⟩ is the eigenstate of the quadrature phase operator ${\hat{x}}_{\theta }$.

In two-mode homodyne detection, we rely on the detection scheme proposed in reference [27], where the two-mode states are characterized by a single homodyne detector. By denoting with $\hat{a}$ and $\hat{b}$ the initial modes to be detected, the mode $\hat{k}$ arriving at the detector can be expressed as [27]

$$\begin{equation}\hat{k}=\mathrm{exp}\left(\text{i}\varphi \right)\mathrm{cos}\phi \enspace \hat{a}+\mathrm{sin}\phi \enspace \hat{b},\end{equation} \tag{ 30 }$$

which corresponds to applying a phase shift of angle φ between the horizontal and vertical polarization components, a polarization rotator of angle ϕ, and a polarizing beam splitter (PBS) reflecting the vertically polarized component of the beam toward the detector [27]. Using repeated measurements of the quadratures for a set of identical states, the homodyne data are collected for which a probability distribution can be assigned with the variance given by:

$$\begin{equation}\langle {\hat{x}}_{\theta }^{2}\rangle -{\langle {\hat{x}}_{\theta }\rangle }^{2}=\mathrm{Tr}\left[P\gamma \right],\end{equation} \tag{ 31 }$$

where P is the matrix for the measurement of the quadrature variance of the mode $\hat{k}$:

$$\begin{equation}P=u{u}^{\text{T}},\quad {u}^{\text{T}}=\left(\begin{matrix}\hfill \mathrm{cos}\phi \mathrm{cos}\left(\theta -\varphi \right),\enspace \mathrm{cos}\phi \mathrm{sin}\left(\theta -\varphi \right),\enspace \mathrm{sin}\phi \mathrm{cos}\theta ,\enspace \mathrm{sin}\phi \mathrm{sin}\theta \hfill \end{matrix}\right).\end{equation} \tag{ 32 }$$

As P is a symmetric, real 4 × 4 matrix we can see that for 10 different combinations of angles θ, ϕ and φ the entire two-mode CM can be reconstructed (the number of unknown independent parameters in an N-mode CM is N(2N + 1)).

The extension of detection to N-mode CV states by a single homodyne detector can be achieved by applying the same two-mode combination scheme N − 1 times. For example, for the initial modes $\hat{a}$, $\hat{b}$ and $\hat{c}$, the generalized mode arriving at the detector is:

$$\begin{equation}\hat{k}=\mathrm{exp}\left(\text{i}{\varphi }_{1}\right)\mathrm{cos}\phi \enspace \hat{a}+\mathrm{exp}\left(\text{i}{\varphi }_{2}\right)\mathrm{sin}\phi \enspace \mathrm{cos}\psi \enspace \hat{b}+\mathrm{sin}\phi \enspace \mathrm{sin}\psi \enspace \hat{c},\end{equation} \tag{ 33 }$$

from where we can see that for ψ = 0 and φ₂ = 0 the two-mode case in equation (30) is obtained. We denote by P_j the matrix of the jth measurement.

4.1. Constructing witnesses

Random measurement directions in the case of two modes are given by random angles θ, ϕ, φ that are drawn from a uniform distribution in an interval:

$$\begin{equation}0{\leqslant}\theta {\leqslant}\pi ,\end{equation} \tag{ 34 }$$

$$\begin{equation}0{\leqslant}\phi {\leqslant}\pi ,\end{equation} \tag{ 35 }$$

$$\begin{equation}0{\leqslant}\varphi {< }2\pi .\end{equation} \tag{ 36 }$$

The problem of finding a witness operator Z, given the repeated independent measurements P_j on the CM, reduces to finding the coefficients c_j such that Z = ∑_j c_j P_j . Therefore, we apply the proposition in order to find the best witness for two-mode CMs, and propose the following SDP:

$$\begin{equation}\begin{aligned}\hfill \underset{x}{\mathrm{minimize}}\mathbf{\text{c}}\cdot \mathbf{\text{m}}\\ \hfill \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}Z& =\sum _{j}{c}_{j}{P}_{j}\hfill \\ \hfill Z& =\left(\begin{matrix}\hfill {Z}_{1}\hfill & \hfill {Z}_{c}\hfill \\ \hfill {Z}_{c}^{\text{T}}\hfill & \hfill {Z}_{2}\hfill \end{matrix}\right){\geqslant}0\hfill \\ \hfill & {Z}_{1}+\mathrm{i}x{{\Omega}}_{1}{\geqslant}0\hfill \\ \hfill & {Z}_{2}+\mathrm{i}\left(\frac{1}{2}-x\right){{\Omega}}_{1}{\geqslant}0,\hfill \end{aligned}\end{equation} \tag{ 37 }$$

where m = Tr( Pγ), with P being the vector of measurement matrices P_j . This SDP finds the matrix Z, given the experimentally obtained data, such that

$$\begin{equation}\mathbf{\text{c}}\cdot \mathbf{\text{m}}=\mathrm{Tr}\left[Z\gamma \right]\end{equation} \tag{ 38 }$$

takes its minimal value, while being an EW as defined in theorem above. If the obtained value in equation (38) is smaller than one, then the CV state with CM γ can be unambiguously identified as being entangled.

This SDP also allows for the identification of the minimal number of measurements that are required for entanglement assessment in arbitrary states. The number of measurements in a tomographically complete setting is given by N(2N + 1), where N is the number of modes. This is the maximal number of measurement settings required to detect entanglement. However, the set of EWs described in the proposition is more restrictive than the set of all EWs. The consequences will be discussed later.

The proposed SDP has the immediate advantage that it does not require any information about the state, except the number of modes N. We will now test the performance of this method by simulating its implementation on random two-mode entangled CV states, and on four-mode bipartite bound entangled states.

The entanglement of two-mode CMs with block structure given by:

$$\begin{equation}\gamma =\left(\begin{matrix}\hfill {\gamma }_{1}\hfill & \hfill {\varepsilon }_{1,2}\hfill \\ \hfill {\varepsilon }_{1,2}^{\text{T}}\hfill & \hfill {\gamma }_{2}\hfill \end{matrix}\right)\end{equation} \tag{ 39 }$$

is quantified by means of the logarithmic negativity [28]:

$$\begin{equation}E=\mathrm{max}\left\{0,-\frac{1}{2}\enspace {\mathrm{log}}_{2}\enspace f\right\},\end{equation} \tag{ 40 }$$

where

$$\begin{align}\hfill f& =\frac{1}{2}\left(\mathrm{det}\enspace {\gamma }_{1}+\mathrm{det}\enspace {\gamma }_{2}\right)-\mathrm{det}\enspace {\varepsilon }_{1,2}\hfill \\ \hfill & \quad -{\left({\left[\frac{1}{2}\left(\mathrm{det}\enspace {\gamma }_{1}+\mathrm{det}\enspace {\gamma }_{2}\right)-\mathrm{det}\enspace {\varepsilon }_{1,2}\right]}^{2}-\mathrm{det}\enspace \gamma \right)}^{1/2}.\hfill \end{align} \tag{ 41 }$$

An EW provides a lower bound for the logarithmic negativity measure when the PPT criterion of separability is necessary and sufficient [16]:

$$\begin{equation}E{\geqslant}{\mathrm{log}}_{2}\enspace \frac{1}{w},\end{equation} \tag{ 42 }$$

where w ∈ (0, 1) is the outcome of measuring an EW on CM γ: Tr[Zγ] = w. For two-mode CMs the logarithmic negativity corresponds to the minimal ⁵ EW Z_min giving the smallest possible value w_min.

In the following we investigate the efficiency of our method for detecting entanglement of arbitrary CV states, with respect to the minimal number of measurements required to accomplish this task. Thus, given an arbitrary unknown CM our algorithm first computes the variances of the generalized quadrature (31) for one random measurement direction in phase space and then carries out the SDP optimization to check if the state is entangled. If entanglement is not detected, additional random measurements are successively simulated and the optimization algorithm is executed each time until the entanglement is detected. At least two measurement settings are required in order to detect entanglement.

5.1. Detecting entanglement in squeezed vacuum states

The CMs of squeezed vacuum states have the form:

$$\begin{equation}\gamma =\left(\begin{matrix}\hfill \mathrm{cosh}\enspace 2r\hfill & \hfill 0\hfill & \hfill \mathrm{sinh}\enspace 2r\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathrm{cosh}\enspace 2r\hfill & \hfill 0\hfill & \hfill -\mathrm{sinh}\enspace 2r\hfill \\ \hfill \mathrm{sinh}\enspace 2r\hfill & \hfill 0\hfill & \hfill \mathrm{cosh}\enspace 2r\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -\mathrm{sinh}\enspace 2r\hfill & \hfill 0\hfill & \hfill \mathrm{cosh}\enspace 2r\hfill \end{matrix}\right),\end{equation} \tag{ 43 }$$

where r is the squeezing parameter. For such states the logarithmic negativity can be calculated using (40), obtaining a linear dependence on the squeezing parameter:

$$\begin{equation}E=2r\enspace {\mathrm{log}}_{2}\enspace e,\end{equation} \tag{ 44 }$$

where e is the Euler constant. Squeezed vacuum states are Gaussian states, which are naturally accessible in many experimental situations where spontaneous down conversion is involved, being also useful in many quantum optics applications [3, 4]. ⁶

In figure 1 we show the fraction of entanglement detection of squeezed vacuum states, using random EWs. Contrary to intuition, states with less entanglement are more easily detected, i.e. they require on average fewer measurements than states with higher entanglement. This is due to the fact that in this case the amount of entanglement is linked to the strength of quadrature squeezing. It is well known that it is difficult to measure high squeezing in CV states [29] (see also the explanation given in figure 2). The full tomography for two-mode CMs is reached by 10 independent measurements. The CM (43) of the squeezed vacuum state has some zero elements, and with this knowledge about the state one would need only 6 measurements to reconstruct the CM entirely. However, our method may require more than 6 measurements to assess entanglement since we assumed no information about the states, except the dimension of the CM. As a consequence of the stronger constraints imposed on the EWs in equation (16) our method requires, with very low probability (0.0094% in our example), more than 10 measurements, which correspond to full tomography.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 2 we show the variance of the generalized quadrature Tr[Pγ], see (31), for θ = 0, as a function of φ ∈ [−π, π] and ϕ ∈ [−π, π], for different values of the squeezing parameter, r = 0.2 (left) and r = 1 (right), of the squeezed vacuum state. The outcomes of the random measurements are represented by the points on this surface. The horizontal plane is given by Tr[Pγ] = 1, which holds for a separable vacuum state with r = 0. The areas below this plane, where Tr[Pγ] < 1, correspond to the region of parameters φ and ϕ for which entanglement is detected. We observe that the areas of the regions of entanglement detection are decreasing with increasing the squeezing. This corresponds to the fact that highly squeezed states occupy a smaller region in phase space in terms of the angles ϕ, φ. Thus, more random measurements are needed to detect the entanglement.

5.2. Detecting entanglement in random covariance matrices

Random CMs are produced as follows. Starting with a CM in diagonal form, with symplectic eigenvalues ν_i ⩾ 1 (i = 1, ..., N) randomly generated from a uniform distribution in a finite real interval [1, t], t > 1:

$$\begin{equation}{\gamma }_{\mathrm{t}\mathrm{h}}=\underset{i=1}{\overset{N}{\oplus }}\left(\begin{matrix}\hfill {\nu }_{i}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\nu }_{i}\hfill \end{matrix}\right),\end{equation} \tag{ 45 }$$

the general random CMs γ are created by applying random symplectic transformations $S\in Sp\left(2N,\mathbb{R}\right)$, as follows [5]:

$$\begin{equation}\gamma =S{\gamma }_{\mathrm{t}\mathrm{h}}{S}^{\text{T}}.\end{equation} \tag{ 46 }$$

The matrix (45) is the CM of thermal states, with the symplectic eigenvalue of every mode i related to the thermal photon number n_i as follows: ν_i = 2n_i + 1 [5]. Random symplectic matrices are generated using the Euler decomposition (7). First, unitary matrices X and Y in equation (9) are generated from the Haar distribution [24], and the symplectic orthogonal matrices K and L are formed as in equation (10). The one-mode squeezers defined in equation (8) are created by randomly choosing parameters r_i via a uniform distribution in a finite interval. For this purpose we implemented the Matlab code presented in reference [31].

In figure 3 we illustrate the efficiency of entanglement detection for general random two-mode CMs, created from thermal state CMs (45) with random symplectic eigenvalues ν_i ∈ [0, 5], by random symplectic transformations (7) with squeezing parameters r_i ∈ [0, 2]. Our method shows a slight improvement in the efficiency of entanglement detection for highly entangled states compared to less entangled states, and most of the time it does not require full tomography. However, because of the strength of our proposed linear constraints, we may need, with very low probability, more measurements than in the full tomography in order to detect entanglement in random two-mode states. The probability that entanglement is detected by 11–12 measurements in this case, is 0.05%.

Zoom In Zoom Out Reset image size

The evident difference in the efficiency of entanglement detection in random CMs compared to squeezed vacuum states may reside in the fact that highly squeezed states look classical in random measurement directions, which does not have to be the case for random states. In addition, squeezed vacuum states are a special class of states for which the logarithmic negativity has a linear dependence on the squeezing parameter alone (see equation (44)), while for a general two-mode CM the logarithmic negativity depends also on thermal photon number of the modes, and the simulation of entanglement detection shows a different behaviour.

In general, it is unlikely to draw randomly an entangled state with high logarithmic negativity, especially for states with a high number of modes. However, for the two-mode CMs, with the range of entanglement considered in figure 3, a substantial fraction of randomly generated CMs is entangled, which allowed us to perform the simulation.

5.3. Detecting bipartite bound entanglement

Since the proposed SDP algorithm can be easily generalized to multi-mode CV states, we provide an example of a four-mode CM with 12 independent parameters, mentioned in reference [32], which has bipartite bound entanglement:

$$\begin{equation}\gamma =\left(\begin{matrix}\hfill 2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 4\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 4\hfill \end{matrix}\right).\end{equation} \tag{ 47 }$$

The detection of bound entanglement by our method proves that the EWs defined in the theorem above, goes beyond the criteria which detect entanglement only in states with non-PPT. A general N-mode CM has N(2N + 1) independent variables, and for the four-mode CM in equation (47) by performing 36 measurements our algorithm provides the best estimate of entanglement, Tr[Z_min γ] = 0.8966, which is in agreement with the results of reference [17]. In figure 4 we depict the frequency of entanglement detection as a function of the number of random measurements composing the witness. The CM in equation (47) is of a rather simple form, however, the construction of the EW detecting bound entanglement requires 33 random measurements on average, since our SDP considers the number of modes of the state as the only available information.

Zoom In Zoom Out Reset image size

Until now we have considered only ideal measurements, where we used the exact variances ${m}_{i}=\mathrm{Tr}\left[{P}_{i}\gamma \right]={\left({\Delta}{\hat{x}}_{{\theta }_{i}}\right)}^{2}$ (see section 4) in order to construct the EW. In real experiments the accessible data are subject to statistical fluctuations. In the following we perform the statistical analysis for the case of Gaussian states, that is, we assume that the data obtained in homodyning, which represent the collection of outcomes ${X}_{ij}={\langle {\hat{x}}_{{\theta }_{i}}\rangle }_{j}$, (j = 1, ..., n_i ), from n_i repetitions of the measurement with the measurement direction given by θ_i , are governed by the normal probability distribution ${\mathcal{N}}_{i}\left({\mu }_{i},{m}_{i}\right)$ with the mean μ_i , and variance ${m}_{i}={\left({\Delta}{\hat{x}}_{{\theta }_{i}}\right)}^{2}$. Given the homodyne data from n_i measurements for a fixed measurement direction θ_i , the sample variance denoted as ${\bar{P}}_{i}$, which estimates the variance m_i , is given by ⁷ :

$$\begin{equation}{\bar{P}}_{i}=\frac{1}{{n}_{i}-1}\sum _{j=1}^{{n}_{i}}{\left({X}_{ij}-{\bar{X}}_{i}\right)}^{2},\end{equation} \tag{ 48 }$$

where ${\bar{X}}_{i}$ is the sample mean:

$$\begin{equation}{\bar{X}}_{i}=\frac{1}{n}\sum _{j=1}^{{n}_{i}}{X}_{ij}.\end{equation} \tag{ 49 }$$

In this case, the estimated value of our witness Tr[Zγ], denoted as $\bar{Z}$, is given by:

$$\begin{equation}\bar{Z}=\sum _{i}{c}_{i}{\bar{P}}_{i},\end{equation} \tag{ 50 }$$

where the index i is used to denote different measurement settings, and the coefficients c_i were introduced in equation (37). In the case when the data comes from a Gaussian probability distribution, the distribution of the sample variance follows the ${\chi }_{{n}_{i}-1}^{2}$ distribution [33]:

$$\begin{equation}\frac{{n}_{i}-1}{{m}_{i}}{\bar{P}}_{i}\sim {\chi }_{{n}_{i}-1}^{2},\end{equation} \tag{ 51 }$$

where ${\chi }_{{n}_{i}-1}^{2}$ is the chi-square distribution with n_i − 1 degrees of freedom, which by definition represents the distribution of sum of squares of n_i − 1 independent, standard normal random variables. The statistical error carried by ${\chi }_{{n}_{i}-1}^{2}$ is given by:

$$\begin{equation}{\Delta}{\chi }_{{n}_{i}-1}^{2}=\sqrt{\mathrm{Var}\left({\chi }_{{n}_{i}-1}^{2}\right)}=\sqrt{2\left({n}_{i}-1\right)},\end{equation} \tag{ 52 }$$

where $\mathrm{Var}\left({\chi }_{{n}_{i}-1}^{2}\right)=2\left({n}_{i}-1\right)$ is the variance of the chi-square distribution [33]. Using the error propagation formula the uncertainty of ${\bar{P}}_{i}$ satisfies:

$$\begin{equation}{\Delta}{\bar{P}}_{i}=\frac{\mathrm{d}{\bar{P}}_{i}}{\mathrm{d}{\chi }_{{n}_{i}-1}^{2}}{\Delta}{\chi }_{{n}_{i}-1}^{2}=\frac{{m}_{i}}{{n}_{i}-1}{\Delta}{\chi }_{{n}_{i}-1}^{2}={m}_{i}\sqrt{\frac{2}{{n}_{i}-1}}.\end{equation} \tag{ 53 }$$

Using again standard error propagation and considering that the number of measurement repetitions is equal for every measurement direction, i.e., n = n_i for every i, we obtain that the resulting error of $\bar{Z}$ defined in equation (50) has the following expression:

$$\begin{equation}{\Delta}\bar{Z}=\sqrt{\sum _{i}{\left(\frac{\mathrm{d}\bar{Z}}{\mathrm{d}{\bar{P}}_{i}}\right)}^{2}{\left({\Delta}{\bar{P}}_{i}\right)}^{2}}=\sqrt{\frac{2}{n-1}}\sqrt{\sum _{i}{c}_{i}^{2}{m}_{i}^{2}}.\end{equation} \tag{ 54 }$$

We stress the fact that, although by our method we can also detect entanglement in non-Gaussian states, this formula for the error of the value of the witness is valid only for Gaussian states. If the X_ij are not normally distributed, then the statistical analysis of EWs based on second moments will require also higher moments of the distribution.

In our method the coefficients c_i are derived from the variances m_i (see equation (37)), while equation (54) neglects the fact that they are not independent. To solve this difficulty one has to divide the homodyne data into two sets. First, one of them is used to derive the coefficients c_i , and then this witness is evaluated using the variances obtained from the other set of data [35]. In this way, the coefficients can be regarded as independent from the errors in the variances of the second set of data. With the quantity in equation (54) it is possible to decide whether it is better to perform additional repetitions of the measurements, or to add new measurements to detect entanglement. For example, consider the single detection of a low entanglement CM with Tr[Z_min γ] = 0.852. In figure 5 the 3σ-confidence of Tr[Zγ] is plotted as a function of the number n of measurement repetitions. It shows that a certification of Tr[Zγ] < 1 with 99.7%-confidence is possible for 6 measurements, which requires a high number of repetitions of the measurements. However, this number significantly decreases when adding another measurement setting.

Zoom In Zoom Out Reset image size

We have proposed a method to detect entanglement of unknown CV states, given only the dimension of their covariance matrices, using random homodyne measurements. Our method provides an alternative for performing full tomography. We characterize the EWs based on second moments using stronger semidefinite constraints than those presented in reference [17], and which account for obtaining a valid witness at all times. Therefore, a quantum state can be clearly considered entangled if it is detected by this criterion. As these constraints are linear, they can be implemented in an SDP. We studied the feasibility of this method in experimental situations, where the figure of merit is considered the number of measurements required to detect entanglement.

First, we tested the proposed algorithm for two-mode squeezed vacuum states, for which the logarithmic negativity linearly depends on the squeezing. We showed that the number of necessary random measurements is very likely to be smaller than for full tomography. We observed an increasing number of measurements required to detect highly entangled states, which is explained by the well-known difficulties in detecting high squeezing.

Our primary objective was to simulate the performance of this method for uniformly drawn random two-mode covariance matrices. Without adding any information about the states, we still found a reduction in the number of measurements needed to certify the presence of entanglement. The phenomenology of entanglement detection in random CV states is very similar to the case of decomposable witnesses for discrete systems [21]. Hence, a higher entanglement is easier to detect, but in our case this improvement is not as significant as in the discrete case. Only with very low probability our method needs more than a tomographically complete set of measurements in order to detect entanglement in random two-mode states.

Bound entangled CV states can also efficiently be detected by a random EW. Similarly to the previous cases entanglement is detected with less than a tomographically complete set of measurements.

As we provide a method of efficient entanglement detection in CV systems which does not require specific measurements, nor necessarily the full information about the CM, it is of high relevance for any experimental application which relies on the valuable resource of entanglement. The experimental scheme implementing our method for two-mode CV states consists of a phase shift in polarization basis, a rotator of polarization, a PBS and a homodyne detector, as e.g. presented in reference [27]. We also extended this scheme to multimode CV states. Our semidefinite optimisation program can easily be adapted to different experimental situations, like making specific measurements or including some additional information about the CM to be measured. Likewise, our method is also applicable to any experimental scheme measuring the quadrature variances, by directly inserting the results of the experiment into the optimisation program we provided. We investigated the statistical robustness of the method, and showed that it has a good robustness to statistical errors.

The authors thank Jochen Szangolies, Thomas Wagner and Matteo Paris for helpful discussions. Giulio Gianfelici has received funding from the German Federal Ministry of Education and Research (BMBF), within the Hardware-based Quantum Security (HQS) project.