Quantum steering with positive operator valued measures

TL;DR

This paper investigates quantum steering with POVMs, formulating steerability as a convex nesting problem, and provides numerical evidence that two-qubit Werner states are unsteerable with POVMs up to a certain mixing probability.

Contribution

It introduces a convex nesting framework for assessing steerability with POVMs and applies it to Werner states, strengthening theoretical understanding and providing numerical evidence.

Findings

Steerability with POVMs can be formulated as a convex nesting problem.

Werner states are unsteerable with POVMs up to a mixing probability of 1/2.

Numerical methods confirm unsteerability within a small error margin.

Abstract

We address the problem of quantum nonlocality with positive operator valued measures (POVM) in the context of Einstein-Podolsky-Rosen quantum steering. We show that, given a candidate for local hidden state (LHS) ensemble, the problem of determining the steerability of a bipartite quantum state of finite dimension with POVMs can be formulated as a nesting problem of two convex objects. One consequence of this is the strengthening of the theorem that justifies choosing the LHS ensemble based on symmetry of the bipartite state. As a more practical application, we study the classic problem of the steerability of two-qubit Werner states with POVMs. We show strong numerical evidence that these states are unsteerable with POVMs up to a mixing probability of $\frac{1}{2}$ within an accuracy of $10^{-3}$.

Equations

$${\mathrm{Tr}}_A[\rho (E_i\otimes {\mathbb{I}}_B)]=\int \mathrm{d}\omega (P)u(P)G_i(P)P,$$

$$K_i=\int \mathrm{d}\omega (P)u(P)G_i(P)P,$$

$$\int \mathrm{d}\omega (P)u(P)P={\mathbb{I}}_A^{\prime },$$

$$\underset{K\in {\mathcal{K}}^n(u)}{\max }⟨Z,K⟩=\int \mathrm{d}\omega (P)u(P)\underset{i}{\max }⟨Z_i,P⟩.$$

$$\int \mathrm{d}\omega (P)u(P)\underset{i}{\max }⟨Z_i,P⟩\ge \underset{E\in {\mathcal{M}}^n}{\max }\sum _{i=1}^n⟨Z_i,E_i^{\prime }⟩$$

$$\underset{Z\in {\mathcal{C}}^n}{\min }\left\{\int \mathrm{d}\omega (P)u(P)\underset{i}{\max }⟨Z_i,P⟩-\underset{E\in {\mathcal{M}}^n}{\max }\sum _{i=1}^n⟨Z_i,E_i^{\prime }⟩\right\}.$$

$$W_p=p|{\psi }^{-}⟩⟨{\psi }^{-}|+(1-p)\frac{{\mathbb{I}}_A}{2}\otimes \frac{{\mathbb{I}}_B}{2},$$

$$\underset{K\in {\mathcal{K}}^n(u)}{\max }⟨Z,K⟩=\underset{G}{\max }\int \mathrm{d}\omega (P)u(P)\sum _{i=1}^n{G}_i(P)⟨Z_i,P⟩,$$

$$I[\lambda (P)]$$

$$-\int \mathrm{d}\omega (P)\lambda (P)\left[\sum _{i=1}^n{G}_i(P)-1\right]\}$$

$$=\underset{G}{\max }\left\{\int \mathrm{d}\omega (P)\sum _{i=1}^n{G}_i(P)\left[u(P)⟨Z_i,P⟩-\lambda (P)\right]\right\}$$

$$+\int \mathrm{d}\omega (P)\lambda (P).$$

$$G_i^{\ast }(P)=\Theta [u(P)⟨Z_i,P⟩-\lambda (P)],$$

$$\sum _{i=1}^n\Theta [u(P)⟨Z_i,P⟩-\lambda (P)]=1.$$

$$\lambda (P)=u(P)\underset{i}{\max }⟨Z_i,P⟩.$$

$${\mathrm{Tr}}_A[\rho (E_i\otimes {\mathbb{I}}_B)]={\int }_{\Lambda }\mathrm{d}\nu (\lambda )G_i(\lambda )F(\lambda ),$$

$${\int }_{\Lambda }\mathrm{d}\nu (\lambda )\underset{i}{\max }⟨Z_i,F(\lambda )⟩\ge \underset{E\in {\mathcal{M}}^n}{\max }\sum _{i=1}^n⟨Z_i,E_i^{\prime }⟩$$

$${\int }_{\Lambda }\mathrm{d}\nu (\lambda )\underset{i}{\max }⟨Z_i,F(\lambda )⟩=\int \mathrm{d}\mu (P)\underset{i}{\max }⟨Z_i,P⟩,$$

$$\underset{E\in {\mathcal{M}}^n}{\max }\sum _{i=1}^n\mathrm{Tr}[\rho (E_i\otimes Z_i)]$$

$$=\underset{E\in {\mathcal{M}}^n}{\max }\sum _{i=1}^n\mathrm{Tr}[\rho U(g)E_i{U}^{†}(g)\otimes V(g)Z_i{V}^{†}(g)]$$

$$=\underset{E\in {\mathcal{M}}^n}{\max }\sum _{i=1}^n\mathrm{Tr}[\rho E_i\otimes V(g)Z_i{V}^{†}(g)],$$

$$\int \mathrm{d}\omega (P)u(P)\underset{i}{\max }⟨Z_i,P⟩\ge \underset{E\in {\mathcal{M}}^n}{\max }\sum _{i=1}^n⟨V(g)Z_i{V}^{†}(g),E_i^{\prime }⟩.$$

$$\int \mathrm{d}\omega (P)u(P)\underset{i}{\max }⟨V^{†}(g)Z_{i}V(g),P⟩\ge \underset{E\in {\mathcal{M}}^n}{\max }\sum _{i=1}^n⟨Z_i,E_i^{\prime }⟩$$

$$\int \mathrm{d}\omega (P)u(P)\underset{i}{\max }⟨V^{†}(g)Z_{i}V(g),P⟩$$

$$=\int \mathrm{d}\omega (P)u[V^{†}(g)PV(g)]\underset{i}{\max }⟨Z_i,P⟩,$$

$$E_i^{\prime }=\int \mathrm{d}\omega (P)u(P)G_i^E(P)P.$$

$$E_i^{\prime }=\int \mathrm{d}\omega (P)R_V(g)u(P)G_i^{U^{†}(g)EU(g)}[V^{†}(g)PV(g)]P.$$

$$E_i^{\prime }=\int \mathrm{d}\omega (P)P\int d\mu (g)R_V(g)u(P)G_i^{U^{†}(g)EU(g)}[V^{†}(g)PV(g)].$$

$$\int \mathrm{d}\mu (g)R_V(g)u(P)G_i^{U^{†}(g)EU(g)}[V^{†}(g)PV(g)]=$$

$${\bar{G}}_i(P)\int \mathrm{d}\mu (g)u[V^{†}(g)PV(g)].$$

$$E_i^{\prime }=\int \mathrm{d}\omega (P)u^{\ast }(P){\bar{G}}_i(P).$$

$$\left[\sum _{i=1}^n{\bar{G}}_i(P)\right]u^{\ast }(P)=u^{\ast }(P),$$

$$E_i^{\prime }$$

$$=\sum _{\alpha =0}^{K-1}\frac{1}{4\pi }\int \mathrm{d}S(P)G_i(\alpha ,P)c_{\alpha }P.$$

Full text

Introduction. Ever since its first examination by Einstein, Podolsky, and Rosen (EPR) in 1935 Einstein et al. (1935), quantum nonlocality has been a puzzling phenomenon. In the EPR thought experiment, one observer, Alice, can perform a measurement on her half of an entangled pair to steer the other half (that belongs to a distant observer, Bob) to ensembles that conflict with the very intuition of classical locality Schrödinger (1936). This conflict was so profound that it prompted EPR to conclude that quantum theory was “incomplete” Einstein et al. (1935) and caused the longest debate in the history of quantum mechanics Bell (2004). A more “complete” theory would be supplemented by hidden variables; however, the seminal work of Bell Bell (1964) demonstrated that no such theory, when constrained by locality, is capable of explaining all quantum mechanical predictions for bipartite systems. Nowadays, quantum nonlocality is perceived as one of the hallmarks of quantum theory that sets it apart from classical notions and underlies numerous quantum information applications Brunner et al. (2014).

Bell’s work defined the first class of quantum nonlocality, now known as Bell nonlocality Brunner et al. (2014). Some 25 years later, Werner realised that Bell nonlocality and entanglement (nonseparability) were in fact two independent forms of quantum nonlocality Werner (1989). In 2007, Wiseman, Jones and Doherty Wiseman et al. (2007) recognised that the original idea of the EPR thought experiment is actually best captured by yet another form of quantum nonlocality – quantum steerability. Since then, quantum steerability has been successfully demonstrated experimentally in loophole-free tests Bennet et al. (2012); Smith et al. (2012); Wittmann et al. (2012). It has been employed in a range of practical quantum information tasks, including quantum cryptography Branciard et al. (2012), randomness certification Law et al. (2014); Passaro et al. (2015), and self-testing Šupić and Hoban (2016); Gheorghiu et al. (2017).

Among this surge of discoveries, a fundamental question remains: which bipartite states manifest quantum steerability? In fact, determining the steerability of a bipartite state when considering all possible measurements, i.e., positive operator valued measure (POVM) measurements, has been such a challenging task that it is unanswered for even the simplest case of two-qubit Werner states (oqp, 2017, Problem 39). The problem remains open in spite of many significant advances towards understanding quantum steering under particular subsets of POVMs, e.g., with projection valued measure (PVM) measurements Wiseman et al. (2007); Jevtic et al. (2015); Nguyen and Vu (2016a), with finite subsets of POVMs Cavalcanti et al. (2016); Hirsch et al. (2016), and with highly noisy POVMs or highly noisy states Quintino et al. (2015).

In this Letter, we are concerned with the problem of quantum steering with POVMs for bipartite systems of arbitrary (but finite) dimension. We demonstrate that for a given choice of local hidden state ensemble, the task of determining whether a quantum state is steerable can be considered as a nesting problem of two convex objects. As a consequence, we derive an inequality which allows a test of steerability for all measurements. Surprisingly, the inequality also reveals a fundamental aspect of quantum steering. Namely, in quantum steering, the choice of local hidden variable is no longer arbitrary as in Bell nonlocality, but can be limited to the set of Bob’s pure states. This in fact makes the study of quantum steering significantly simpler than its partner Bell nonlocality. In particular, one can strengthen the theorem (Lemma 1 of Ref. Wiseman et al. (2007)) which limits the choice of local hidden state ensemble based on the symmetry of the state. As the first application, we then apply the inequality to study the steerability of the two-qubit Werner states. Contrary to the fact that general POVMs provide an advantage over PVMs in many situations Vértesi and Bene (2010); Hirsch et al. (2013); Acín et al. (2016), we provide strong numerical evidence that POVMs and PVMs are in fact equivalent for steering two-qubit Werner states.

Quantum steerability. Suppose Alice and Bob share a bipartite quantum state $\rho$ over the finite-dimensional Hilbert space $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$. We use $\mathcal{A}^{H}$ (or $\mathcal{B}^{H}$) to denote the space of Hermitian operators over $\mathcal{H}_{A}$ (or $\mathcal{H}_{B}$). A POVM measurement with $n$ outcomes ($n$-POVM) $E$ implemented by Alice is an (ordered) collection of $n$ positive operators, $E=\{E_{i}\}_{i=1}^{n}$ with $E_{i}\in\mathcal{A}^{H}$, $E_{i}\geq 0$ and $\sum_{i=1}^{n}E_{i}=\mathbb{I}_{A}$, where $\mathbb{I}_{A}$ is the identity operator on $\mathcal{A}^{H}$. On performing the measurement, Alice steers Bob’s system to the steering ensemble $\{\operatorname{Tr}_{A}[\rho(E_{i}\otimes\mathbb{I}_{B})]\}_{i=1}^{n}$. However, despite the arbitrary choice of measurements, for certain bipartite states, the steering experiment can be locally simulated. More specifically, let $u$ be an ensemble (that is, a probability distribution) on the set of Bob’s pure states, denoted by $\mathcal{S}_{B}$. A state $\rho$ is then called $u$-unsteerable (always considered from Alice’s side) with respect to $n$-POVMs if, for any $n$-POVM $E$, Alice can find $n$ response functions $G_{i}$ (with $G_{i}(P)\geq 0$ and $\sum_{i=1}^{n}G_{i}(P)=1$ for all $P\in\mathcal{S}_{B}$) such that the steering ensemble can be simulated via a local hidden state model Wiseman et al. (2007),

[TABLE]

where the integral is taken over the Haar measure $\omega$ on Bob’s pure states $\mathcal{S}_{B}$. Equation (1) ensures that Bob, when performing state tomography conditioned on Alice’s outcomes, obtains the same result as if Alice were steering his system Wiseman et al. (2007). The ensemble $u$ is called a local hidden state (LHS) ensemble. In principle, the domain of the ensemble $u$ can be extended to mixed states. However, as a mixed state can be written as a convex combination of pure states, restricting the domain of $u$ to pure states causes no loss of generality. We say then that $\rho$ is unsteerable with $n$-POVMs if there exists $u$ such that $\rho$ is $u$-unsteerable with $n$-POVMs.

We note that this definition of quantum steering is slightly different from the original definition Wiseman et al. (2007). In the latter, the LHS ensemble is indexed by a local hidden variable. We will prove the two definitions are equivalent as parts of our results. Our seemingly minor simplification in fact has very important consequences, which will be discussed below.

The set of $n$-POVMs and its geometry. The key idea in our approach is that an $n$-POVM $E$ can be thought of as a point in the real vector space of composite operators $(\mathcal{A}^{H})^{\oplus n}=\oplus_{i=1}^{n}\mathcal{A}^{H}$. We therefore write $E=\oplus_{i=1}^{n}E_{i}$, which explicitly indicates that it is a composite operator in $(\mathcal{A}^{H})^{\oplus n}$ with components $E_{i}$ each bounded by $0\leq E_{i}\leq\mathbb{I}_{A}$. The space $(\mathcal{A}^{H})^{\oplus n}$ can be made Euclidean by defining an inner product $\left\langle X,Y\right\rangle=\sum_{i=1}^{n}\left\langle X_{i},Y_{i}\right\rangle$ for any composite operators $X$ and $Y$, where $\left\langle X_{i},Y_{i}\right\rangle$ denotes the Hilbert–Schmidt inner product of $\mathcal{A}^{H}$, $\left\langle X_{i},Y_{i}\right\rangle=\operatorname{Tr}(X^{\dagger}_{i}Y_{i})$. The set of $n$-POVMs is then a convex and compact subset of this space Chiribella et al. (2007), which we denote by $\mathcal{M}^{n}$. Since $\sum_{i=1}^{n}E_{i}=\mathbb{I}_{A}$, $\mathcal{M}^{n}$ in fact belongs to the linear manifold $\mathcal{P}^{n}=\{X|\sum_{i=1}^{n}X_{i}=\mathbb{I}_{A}\}$.

While the set of POVMs $\mathcal{M}^{n}$ is perhaps unfamiliar, it is similar to the classical probability simplexes. Figure 1 illustrates this similarity when $n=3$. The basis to construct a probability simplex is the probability range $[0,1]$. To obtain a probability simplex in $\mathbb{R}^{3}$, one aligns $3$ probability ranges $[0,1]$ along the $3$ axes and forms the triangle with vertices $(1,0,0)$, $(0,1,0)$, $(0,0,1)$. To construct the set of $3$-POVMs, the probability range $[0,1]$ is replaced by the set $0\leq X\leq\mathbb{I}_{A}$, which for qubits forms a double cone in $\mathcal{A}^{H}$ (illustrated in Figure 1) Nguyen and Vu (2016b). One then ‘aligns’ these $3$ sets $0\leq X\leq\mathbb{I}_{A}$ along the three orthogonal component spaces of $(\mathcal{A}^{H})^{\oplus 3}$. The set of $3$-POVMs is formed in between the points $(\mathbb{I}_{A},0,0)$, $(0,\mathbb{I}_{A},0)$, $(0,0,\mathbb{I}_{A})$. This analogy between a probability simplex and the set of $n$-POVMs applies in the same way for any $n$. There is, however, a crucial difference between classical probability simplexes and sets of POVMs: while $[0,1]$ is $1$-dimensional with $2$ extreme points [math] and $1$, the set $0\leq X\leq\mathbb{I}_{A}$ is generally high-dimensional and with more extreme points other than [math] and $\mathbb{I}_{A}$. As a result, the set of $n$-POVMs is also of high dimension and carries other extreme points apart from the special ones at the ‘corners’, which are of the form $\oplus_{i=1}^{n}\delta_{ik}\mathbb{I}_{A}$ with $k=1,2,\ldots n$.

The steering assemblage of $n$-POVMs. Now each POVM measurement $E$ performed on Alice’s side gives rise to a steering ensemble on Bob’s side, $\oplus_{i=1}^{n}\operatorname{Tr}_{A}[\rho(E_{i}\otimes\mathbb{I}_{B})]$. This is most easily implemented by the concept of a steering function $\rho^{A\to B}:\mathcal{A}^{H}\to\mathcal{B}^{H}$, which maps $X\in\mathcal{A}^{H}$ to $\operatorname{Tr}_{A}[\rho(X\otimes\mathbb{I}_{B})]\in\mathcal{B}^{H}$ Nguyen and Vu (2016b). This induces the map $(\rho^{A\to B})^{\oplus n}:(\mathcal{A}^{H})^{\oplus n}\to(\mathcal{B}^{H})^{\oplus n}$. For $X$ being an element or a subset of $\mathcal{A}^{H}$, we denote $X^{\prime}=\rho^{A\to B}(X)$. The same notation is used for composite vectors, namely, for $X$ being an element or a subset of $(\mathcal{A}^{H})^{\oplus n}$, $X^{\prime}=(\rho^{A\to B})^{\oplus n}(X)$.

Geometrically, the map $(\rho^{A\to B})^{\oplus n}$ maps a point in the set of POVMs $\mathcal{M}^{n}$ to a point in $(\mathcal{B}^{H})^{\oplus n}$. The set $(\mathcal{M}^{n})^{\prime}=(\rho^{A\to B})^{\oplus n}(\mathcal{M}^{n})$ is called the steering assemblage of $n$-POVMs. Being a linear image of $\mathcal{M}^{n}$, which is convex and compact Chiribella et al. (2007), $(\mathcal{M}^{n})^{\prime}$ is also convex and compact. Moreover, since $\mathcal{M}^{n}$ belongs to $\mathcal{P}^{n}$, $(\mathcal{M}^{n})^{\prime}$ belongs to $(\mathcal{P}^{n})^{\prime}$.

The capacity of an ensemble of Bob’s pure states. For an ensemble $u$ of Bob’s pure states $\mathcal{S}_{B}$, the $n$-capacity $\mathcal{K}^{n}(u)$ is the set of $n$-component ensembles that it can simulate. That is to say, the capacity $\mathcal{K}^{n}(u)$ is a subset of $(\mathcal{B}^{H})^{\oplus n}$ consisting of composite operators $K=\oplus_{i=1}^{n}K_{i}$, each component being given by

[TABLE]

with all possible choices of response functions $G_{i}$ that satisfy $G_{i}(P)\geq 0$, $\sum_{i=1}^{n}G_{i}(P)=1$. It is easy to show that the $n$-capacity $\mathcal{K}^{n}(u)$ is also a convex compact set, which has $n$ special extreme points of the form $\oplus_{i=1}^{n}\delta_{ik}\int\mathrm{d}\omega(P)u(P)P$ with $k=1,2,\ldots,n$.

Steerability as a nesting problem. With the above definitions, the following lemma is obvious.

We first consider the special extreme points of $\mathcal{M}^{n}$. It is easy to show that for their steering images to be in $\mathcal{K}^{n}(u)$, one has

[TABLE]

which is referred to as the minimal requirement for $u$ 111Since if $\int\mathrm{d}\omega(P)u(P)G_{i}(P)P=\mathbb{I}_{A}^{\prime}$ for some response function $G_{i}$, then by taking the trace of this equality, one obtains $\int\mathrm{d}\omega(P)u(P)G_{i}(P)=1$. Because $G_{i}(P)\geq 0$, this tells us that $G_{i}(P)=1$ almost everywhere with respect to the measure generated by $u$. This minimal requirement can also be obtained by summing both sides of equation (1) over $i$, bearing in mind the constraints $\sum_{i=1}^{n}E_{i}=\mathbb{I}_{A}$ and $\sum_{i=1}^{n}G_{i}(P)=1$. It can also be thought as the condition for $(\mathcal{M}^{1})^{\prime}\subseteq\mathcal{K}^{1}(u)$.

Once reformulated in terms of a nesting problem of convex objects (Lemma 1), one can apply nesting criteria to test steerability. The following lemma is such a nesting criterion based on a duality representation.

The idea behind this lemma is that if $\mathcal{X}$ contains $\mathcal{Y}$, then its projection onto any direction contains that of $\mathcal{Y}$ and vice versa (see Figure 2). A full proof is given in Appendix B.

To apply this lemma with $\mathcal{X}=\mathcal{K}^{n}(u)$, $\mathcal{Y}=(\mathcal{M}^{n})^{\prime}$, we need to solve two maximisation problems: $\max_{K\in\mathcal{K}^{n}(u)}\left\langle Z,K\right\rangle$ and $\max_{E\in\mathcal{M}^{n}}\left\langle Z,E^{\prime}\right\rangle$ for a given composite operator $Z$ in $(\mathcal{B}^{H})^{\oplus n}$. While the latter is a semidefinite program, the former is a linear constrained maximisation, which can be solved explicitly (see Appendix C for the details):

[TABLE]

From Lemma 1 and Lemma 2, the following theorem immediately follows.

[TABLE]

for all composite operators $Z=\oplus_{i=1}^{n}Z_{i}$ in $(\mathcal{B}^{H})^{\oplus n}$.

Inequality (5) is the main result in this Letter: it is valid for systems of arbitrary dimension and POVMs of arbitrary number of outcomes. We now discuss some of its important consequences; details of proofs and further discussions are given in Appendix D.

Had one started with the original definition of quantum steering with an “indexed LHS ensemble”, that is, an LHS ensemble indexed by some hidden variable, one would arrive at a similar inequality as (5). In that case, the integration is taken over the hidden variable instead (see Appendix D). However, one can rewrite it as an integral over the push-forward measure over Bob’s pure states $\mathcal{S}_{B}$ Bogachev (2007). This implies that it is only the push-forward measure on $\mathcal{S}_{B}$ that really determines the capacity of an LHS ensemble. Two indexed LHS ensembles generating the same measure on Bob’s pure states would have the same capacity. In other words, our definition of quantum steering where the local hidden variable is omitted is equivalent to the original definition of quantum steering (Corollary 1, Appendix D). Having eliminated the arbitrary choice of local hidden variable in quantum steering, the symmetry of the state directly has a stronger implication on the symmetry of LHS ensembles (see Theorem 2). In fact, it is this stronger implication of symmetry that actually renders many unsolvable problems in Bell nonlocality solvable in the context of quantum steering Wiseman et al. (2007).

More specifically, the state $\rho$ is said to have $(\mathcal{G},U,V)$-symmetry with $\mathcal{G}$ being a group and $U$ and $V$ being its two representations on $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$, respectively, if $\rho=U^{\dagger}(g)\otimes V^{\dagger}(g)\rho U(g)\otimes V(g)$ for all $g\in\mathcal{G}$. The action $V$ of $\mathcal{G}$ on Bob’s pure states $\mathcal{S}_{B}$ generates an action $R_{V}$ on the space of distributions on $\mathcal{S}_{B}$ defined by $[R_{V}(g)u](P)=u[V^{\dagger}(g)PV(g)]$. We then have a strengthened form of Lemma 1 of Ref. Wiseman et al. (2007) on the symmetry of LHS ensemble.

This theorem is applicable as well when measurements are restricted to PVMs. To understand the difference with Lemma 1 of Ref. Wiseman et al. (2007), we consider the example where $\mathcal{G}$ acts transitively on Bob’s pure states $\mathcal{S}_{B}$ (i.e., a single orbit covers all of $\mathcal{S}_{B}$). If $\rho$ is unsteerable, Lemma 1 of Ref. Wiseman et al. (2007) then states the existence of an indexed LHS ensemble, on which $\mathcal{G}$ acts covariantly on the indices and the states. However, due to the arbitrariness in choice of the local hidden variable, there exist in fact infinitely many different $\mathcal{G}$-covariant indexed LHS ensembles (see Appendix D). One then could not single out an unique choice of LHS ensemble. On the other hand, under the same conditions, Theorem 2 implies that the state is unsteerable with the unique uniform distribution on Bob’s pure states as an LHS ensemble.

Beyond revealing very general aspects of quantum steering, Theorem 1 can also be used to test steerability in practice. The most difficult part is to determine the existence of a LHS ensemble $u$. Even when measurements are limited to PVMs, the question is so far solved only for highly symmetric states, e.g., the Werner state Wiseman et al. (2007) and the two-qubit $T$-state (mixtures of Bell states) Jevtic et al. (2015); Nguyen and Vu (2016a). The implication of our approach on this problem will be discussed elsewhere. However, as we mentioned, even when $u$ is known, the problem of determining steerability with POVMs is still open (oqp, 2017, Problem 39). It is this latter problem that we are concerned with in the following. We will show that Theorem 1 can provide a strong numerical evidence for steering with POVMs with a given identified candidate for the LHS ensemble.

The gap function. We first note that for a LHS ensemble $u$ satisfying the minimal requirement (3), the inequality (5) is invariant with respect to the transformation $\oplus_{i=1}^{n}Z_{i}\to\frac{1}{\sqrt{D}}\oplus_{i=1}^{n}(Z_{i}-C)$, where $C=\frac{1}{n}\sum_{i=1}^{n}Z_{i}$ and $D=\sum_{i=1}^{n}\left\langle Z_{i}-C,Z_{i}-C\right\rangle$. We can therefore restrict $Z$ to the set of those satisfying $\sum_{i=1}^{n}Z_{i}=0$ and $\sum_{i=1}^{n}\left\langle Z_{i},Z_{i}\right\rangle=1$, denoted by $\mathcal{C}^{n}$. For clarity, we introduce the gap function $\Delta[(\mathcal{M}^{n})^{\prime},\mathcal{K}^{n}(u)]$, defined to be

[TABLE]

The gap function characterises the gap between the boundary of $(\mathcal{M}^{n})^{\prime}$ and that of $\mathcal{K}^{n}(u)$ from inside. The state $\rho$ is $u$-unsteerable with $n$-POVMs if and only if $\Delta[(\mathcal{M}^{n})^{\prime},\mathcal{K}^{n}(u)]\geq 0$.

Restriction to rank-$1$ POVMs. As all POVMs can be post-processed from those of rank-$1$, to test quantum steerability we can concentrate on the latter Barrett (2002). To this end, we define $\mathcal{\widetilde{\mathcal{M}}}^{n}=\{\oplus_{i=1}^{n}\alpha_{i}P_{i}\}$ where $P_{i}$ are (not necessarily independent) rank-$1$ projections and $0\leq\alpha_{i}\leq 1$ such that $\sum_{i=1}^{n}\alpha_{i}P_{i}=\mathbb{I}_{A}$. To test the steerability with $\mathcal{M}^{n}$, we therefore only need to calculate $\Delta[(\widetilde{\mathcal{M}}^{n})^{\prime},\mathcal{K}^{n}(u)]$.

Steerability of two-qubit Werner states with POVMs. Consider the two-qubit Werner state,

[TABLE]

which is a mixing between the singlet Bell state $\left|\psi^{-}\right\rangle=\frac{1}{\sqrt{2}}(\left|01\right\rangle-\left|10\right\rangle)$ and the maximally mixed state with mixing parameter $p$ ($0\leq p\leq 1$). When restricted to PVMs, by explicitly constructing the response functions, it was shown that the Werner state is unsteerable for the mixing probabilities $p\leq\frac{1}{2}$ Wiseman et al. (2007). By adapting Barrett’s model Barrett (2002) of local hidden variables, it was further possible to show that Werner states are unsteerable with POVMs when $p\leq\frac{5}{12}$ Quintino et al. (2015). We are to study the conjecture Werner (2014):

Although further analyses restricted to finite subsets of POVMs Cavalcanti et al. (2016); Hirsch et al. (2016) or POVMs with special symmetry Werner (2014) support unsteerability of the Werner state for $\frac{5}{12}\leq p\leq\frac{1}{2}$ Bavaresco et al. (2017); Werner (2014), there has not been a concrete evidence when one considers all POVMs. For $n=1$, the conjecture is obvious. For $n=2$, it has been proven by demonstrating the equivalence to steering with PVMs Nguyen and Vu (2016b). The proof for $n=3$ is also known Werner (2014). Finally, it is known that it is sufficient to consider the conjecture for $n=4$ D’Ariano et al. (2005); Werner (2014). Here, by computing the corresponding gap function for the Werner state, we provide strong numerical evidence for Conjecture 1 for $n=4$.

It is easy to see that the Werner state has $\operatorname{U}(2)$ symmetry as defined above. Moreover, since the action is transitive on Bob’s Bloch sphere, by Theorem 2, the candidate for the LHS ensemble can be limited to the uniform distribution, $u=1$. To simplify the notation, from now on we simply use $\Delta$ to denote the gap function without specifying the set of measurements and the capacity, which can be understood from the context. The computation of the gap function $\Delta$ generally requires global minimisation over $Z$, which is carried out by the standard (non-deterministic) simulated annealing algorithm (see Appendix E for details).

In Figure 3, we present values of $\Delta$ found for the mixing probability $p\approx\frac{1}{2}$. For comparison, we also present the numerical results of $\Delta$ when the measurements are limited to PVMs, which are in very good agreement with the analytical calculation (see Appendix F). For $p>\frac{1}{2}$, one observes that the gap function $\Delta$ for rank-$1$ $4$-POVMs is negative and coincides with the gap function for PVMs. For all $p\leq\frac{1}{2}-10^{-3}$, the best obtained values for $\Delta$ are less than $10^{-10}$ but persistently non-negative. The fact that the gap function tends to vanish for $p\leq\frac{1}{2}$ instead of attending finite positive values seems to be because of the high dimensionality of POVMs 222The effect of dimensionality can be understood directly from Figure 2. On the plane, the gap between the boundaries of the two objects is positive. However when imbedded in the 3D space, the gap will be zero in direction that is perpendicular to plane that contains the two objects. The presented results therefore support conjecture 1. Unfortunately, within this work, the ambiguity region of $\frac{1}{2}-10^{-3}\leq p\leq\frac{1}{2}$ cannot be resolved due to limits of numerical accuracy (see Appendix E).

Conclusion. Our work simplifies the definition of quantum steering, where we show that the local hidden variable indexing LHS ensemble can be omitted. As a direct consequence, a stronger theorem on the symmetry of the LHS ensemble is derived. We have thereby opened a general approach to studying quantum steerability with POVMs. Further works to strengthen the numerical evidence for the unsteerability of the Werner state at $p=\frac{1}{2}$ and testing steerability with POVMs of arbitrary two-qubit states are underway. Moreover, although the current illustrative applications are based on two qubits, our approach is not limited by the dimensionality of the systems. It is hoped that systematic tests for steering with POVMs (particularly for high dimensional systems) will give a complete answer to the fundamental question of the equivalence between POVMs and PVMs for steerability. Beyond quantum steering, we leave open the question of whether this approach can be extended to characterise Bell nonlocality with POVMs.

Appendix A The compactness of the capacity

In this appendix, we prove the compactness of the capacity $\mathcal{K}^{n}(u)$. On the Bloch sphere $\mathcal{S}_{B}$, the distribution $u$ defines a measure $\mu$. We assume that the response functions $G_{i}(P)$ are squared-integrable with respect to $\mu$. Consider the space $\mathcal{L}^{2}(\mathcal{S}_{B},\mu)$ of squared-integrable functions on $\mathcal{S}_{B}$ with respect to $\mu$, which is a Hilbert space (as usual, we ignore the difference on a zero-measured set) (Hu, 1967, Chapter V). We then construct the Hilbert space $[\mathcal{L}^{2}(\mathcal{S}_{B},\mu)]^{\oplus n}$ as usual. We define the subset of $\Omega=\{G=\oplus_{i=1}^{n}G_{i}\in[\mathcal{L}^{2}(\mathcal{S}_{B},\mu)]^{\oplus n}:0\leq G_{i}(P)\leq 1,\sum_{i=1}^{n}G_{i}(P)=1\}$, which is closed and convex, thus weakly closed (Conway, 1990, Chapter V, Corollary 1.5). Moreover, it is obviously bounded, thus weakly compact (Conway, 1990, Chapter V, Theorem 4.2).

Now consider the linear operator $T:[\mathcal{L}^{2}(\mathcal{S}_{B},\mu)]^{\oplus n}\to(\mathcal{B}^{H})^{\oplus n}$, defined by $T(G)=\oplus_{i=1}^{n}\int\mathrm{d}\mu(P)G_{i}(P)P$. It is obvious that $T$ is bounded, thus continuous and weakly continuous (Conway, 1990, Chapter VI, Theorem 1.1). It then follows directly that $\mathcal{K}^{n}(u)=T(\Omega)$ is compact.

Appendix B Nesting criterion by duality

In this Appendix, we provide the proof for the nesting criterion by duality.

Lemma 2** (Nesting criterion by duality).**

Let $\mathcal{X}$ be a convex compact subset of a finite-dimensional Euclidean space. Then a compact subset $\mathcal{Y}$ is contained in $\mathcal{X}$ if and only if $\max_{X\in\mathcal{X}}\left\langle Z,X\right\rangle\geq\max_{Y\in\mathcal{Y}}\left\langle Z,Y\right\rangle$ for all vectors $Z$ in the space.

Proof.

It is obvious that if $\mathcal{Y}\subseteq\mathcal{X}$ then $\max_{X\in\mathcal{X}}\left\langle Z,X\right\rangle\geq\max_{Y\in\mathcal{Y}}\left\langle Z,Y\right\rangle$ for all $Z$. Now suppose $\max_{X\in\mathcal{X}}\left\langle Z,X\right\rangle\geq\max_{Y\in\mathcal{Y}}\left\langle Z,Y\right\rangle$ for all $Z$ and $\mathcal{Y}\not\subseteq\mathcal{X}$. Because $\mathcal{Y}\not\subseteq\mathcal{X}$, there exists $A\in\mathcal{Y}$, $A\not\in\mathcal{X}$. Since $\mathcal{X}$ is a convex and compact set, by the separation theorem, $A$ is separated from $\mathcal{X}$ by a hyperplane, i.e., there exists a vector $Z$ such that $\left\langle Z,A\right\rangle>\max_{X\in\mathcal{X}}\left\langle Z,X\right\rangle$ Rockafellar (2015). It follows that $\max_{Y\in\mathcal{Y}}\left\langle Z,Y\right\rangle\geq\left\langle Z,A\right\rangle>\max_{X\in\mathcal{X}}\left\langle Z,X\right\rangle$, contradicting the assumption. ∎

Appendix C Solving the first optimisation problem

Here we provide details of the solution to the constrained maximisation problem (4). Using the definition of $\mathcal{K}^{n}(u)$, we have

[TABLE]

subject to the constraints $G_{i}(P)\geq 0$ and $\sum_{i=1}^{n}G_{i}(P)=1$. This is a linear maximisation problem with linear constraints, which can be solved easily by Lagrange’s multipliers. For every constraint $\sum_{i=1}^{n}G_{i}(P)=1$ for each $P\in\mathcal{S}_{A}$, we introduce a Lagrange’s multiplier $\lambda(P)$. This leads us to a modified unconstrained maximisation problem

[TABLE]

Note that the function under maximisation is a linear function of $G_{i}(P)$, which is bounded by $0\leq G_{i}(P)\leq 1$. Therefore $I[\lambda(P)]$ is saturated by

[TABLE]

where $\Theta$ is the Heaviside step function, $\Theta(x)=1$ if $x\geq 0$, and $\Theta(x)=0$ otherwise. One now needs to choose $\lambda(P)$ such that the constraint is satisfied,

[TABLE]

Consider some fixed $P$. The last equation means that $\lambda(P)$ must be such that out of $\{u(P)\left\langle Z_{i},P\right\rangle\}_{i=1}^{n}$, only one is larger than or equal to $\lambda(P)$. In other words, the suitable choice for $\lambda(P)$ is

[TABLE]

With this solution for $\lambda(P)$, substituting (10) to (8) one then obtains (4).

So far we actually ignored the case where for some $P$, $\max_{i}\left\langle Z_{i},P\right\rangle$ is attained by two indices, say, $i=i_{1}$ and $i=i_{2}$. In this case, one then has to slightly modify (10): at such a point $P$, while $G_{i}^{\ast}(P)=0$ for $i\neq i_{1},i_{2}$, $G_{i_{1}}^{\ast}(P)$ and $G_{i_{2}}^{\ast}(P)$ can take arbitrary values between [math] and $1$, provided that $G_{i_{1}}^{\ast}(P)+G_{i_{2}}^{\ast}(P)=1$. Similar modification is needed if $\max_{i}\left\langle Z_{i},P\right\rangle$ is attained by more indices. The maximal value (4) however remains the same.

Appendix D Corollaries of Theorem 1

D.1 The equivalence between our definition and the original definition of steering

The original definition of quantum steering Wiseman et al. (2007) goes as follows. Let $(\Lambda,\nu)$ be a probability measure space (we ignore the symbol which denotes the $\sigma$-algebra for the measure $\nu$). Let $F:\Lambda\to\mathcal{S}_{B}$ be a measurable function from the index space $(\Lambda,\nu)$ to the set of Bob’s pure states $\mathcal{S}_{B}$. A state $\rho$ is then called $(\Lambda,\nu)$-unsteerable (from Alice’s side) with respect to $n$-POVMs if, for any $n$-POVM $E$, Alice can find $n$ response functions $G_{i}$ (with $G_{i}(\lambda)\geq 0$ and $\sum_{i=1}^{n}G_{i}(\lambda)=1$ for all $\lambda\in\Lambda$) such that the steering ensemble can be simulated via a local hidden state model Wiseman et al. (2007),

[TABLE]

where the integral is taken over the index space $\Lambda$. In this case, we say $\rho$ admits an indexed LHS model.

If $\rho$ satisfies the definition of steering in the main text, it is cleared that it admits an indexed LHS model, where the index space is Bob’s pure states themselves. On the other hand, if $\rho$ admits an indexed LHS model, it is unclear that the existence of a response function on Bob’s pure states is guaranteed. This is particularly important if the indexing function $F$ is a many-to-one function. One of the strengths of our approach is that it provides a proof that a response function on Bob’s pure states does exist. We know of no other (constructive) proof at the moment.

Corollary 1**.**

Our definition of quantum steering is equivalent to the conventional definition of quantum steering where the LHS ensemble is indexed by a hidden variable.

Proof.

Suppose $\rho$ admits an indexed LHS ensemble, which is indexed by $(\Lambda,\nu)$. Then following the argument that leads to Theorem 1 in the main text, we arrive at the following statement:

Lemma**.**

A bipartite state $\rho$ is $(\Lambda,\nu)$-unsteerable if and only if

[TABLE]

for all composite operators $Z=\oplus_{i}^{n}Z_{i}$ in $(\mathcal{B}^{H})^{\oplus n}$.

Now denote by $\mu$ the push-forward measure generated by $F$ on the set of Bob’s pure states Bogachev (2007). Changing the variable in the integral, one has

[TABLE]

where the latter integral is taken over Bob’s pure states. Let $u$ be the distribution on Bob’s pure states generated by $\mu$ with respect the Haar measure $\omega$. This means the inequality (5) in the main text is satisfied for distribution $u$. According to Theorem 1, $\rho$ is $u$-unsteerable. ∎

D.2 The symmetry of LHS ensembles

For a given state $\rho$, we denote by $\Omega^{n}(\rho)$ the set of ensembles $u$ over Bob’s pure states such that $\rho$ is $u$-unsteerable with $n$-POVMs ($\Omega^{n}(\rho)$ is empty if $\rho$ is steerable). From inequality (5), it is easy to see that $\Omega^{n}(\rho)$ is convex (Corollary 2). Moreover, we show that the symmetry of $\rho$ implies the symmetry of $\Omega^{n}(\rho)$.

Corollary 2**.**

For a given state $\rho$, $\Omega^{n}(\rho)$ is convex.

Proof.

Suppose $u_{1}$ and $u_{2}$ are in $\Omega^{n}(\rho)$. That is to say, $u_{1}$ and $u_{2}$ satisfy inequality (5). It is then easy to check that inequality (5) is also satisfied for all convex combinations of $u_{1}$ and $u_{2}$. In other words, all convex combinations of $u_{1}$ and $u_{2}$ are in $\Omega^{n}(\rho)$. ∎

Let us recall from the main text that the state $\rho$ is said to have $(\mathcal{G},U,V)$-symmetry with $\mathcal{G}$ being a group and $U$ and $V$ being its two representations on $\mathcal{H}_{A}$ and $\mathcal{H}_{B}$, respectively, if $\rho=U^{\dagger}(g)\otimes V^{\dagger}(g)\rho U(g)\otimes V(g)$ for all $g\in\mathcal{G}$. The action $V$ of $\mathcal{G}$ on Bob’s pure states $\mathcal{S}_{B}$ generates an action $R_{V}$ on the space of distributions on $\mathcal{S}_{B}$ defined by $[R_{V}(g)u](P)=u[V^{\dagger}(g)PV(g)]$.

Corollary 3**.**

For a given state $\rho$ which is $(\mathcal{G},U,V)$-symmetric, then $\Omega^{n}(\rho)$ is $(\mathcal{G},V)$-symmetric, i.e., $\Omega^{n}(\rho)=R_{V}(g)\Omega^{n}(\rho)$ for all $g\in\mathcal{G}$.

Proof.

We need to show that if $\rho$ is $u$-unsteerable for some $u$, that is, if inequality (5) holds for $u$, then it also holds for $R_{V}(g)u$. Due to the symmetry of $\rho$, we have

[TABLE]

where the last equality is because $\mathcal{M}^{n}$ is symmetric under the action $U$ of $\mathcal{G}$. Then inequality (5) is equivalent to

[TABLE]

That this inequality holds for all $Z$ is then equivalent to

[TABLE]

for all $Z$. Now we manipulate the left-hand side,

[TABLE]

where the last inequality is a change of integration variable. By definition, $R_{V}(g)u(P)=u(V^{\dagger}(g)PV(g))$. Therefore, inequality (5) indeed holds for $R_{V}(g)u$. ∎

Theorem 2** (Symmetry of LHS ensemble).**

For a given state $\rho$ which is $(\mathcal{G},U,V)$-symmetric with a compact group $\mathcal{G}$, if $\rho$ is unsteerable with $n$-POVMs then it admits an LHS ensemble $u^{\ast}$ which is $(\mathcal{G},R_{V})$-invariant, i.e., $u^{\ast}=R_{V}(g)u^{\ast}$ for all $g$ in $\mathcal{G}$.

Proof 1.

Since $\rho$ is unsteerable, there exists $u\in\Omega^{n}(\rho)$. By Corollary 3, $R_{V}(g)u\in\Omega^{n}(\rho)$ for all $g\in\mathcal{G}$. Since $\Omega^{n}(\rho)$ is convex (Corollary 2), the average over the Haar measure of $\mathcal{G}$, i.e., $u^{\ast}=\int_{\mathcal{G}}\mathrm{d}\mu(g)R_{V}(g)u\in\Omega^{n}(\rho)$, also belongs to $\Omega^{n}(\rho)$. This averaged distribution $u^{\ast}$ is obviously invariant under the action of $\mathcal{G}$, i.e., $u^{\ast}=R_{V}(g)u^{\ast}$ for all $g$ in $\mathcal{G}$. ∎

For completeness, we also provide an alternative proof of this theorem without the use of inequality (5). This proof is more similar to the original proof in Ref. Wiseman et al. (2007); to get the stronger statement, one has to apply the so-called mean value theorem for integrals Bogachev (2007), though.

Proof 2.

Suppose $\rho$ is $u$-unsteerable with $n$-POVMs, then for a POVM $E$, there exists response function $G$ such that

[TABLE]

Here, to track the dependence of the response function $G$ on the measurement, we introduce the superscript $E$ for $G$. Now due to the symmetry of the state, we also have

[TABLE]

Since the right-hand-side is independent of $g$, we can take the average over $g$ with respect to the Haar measure $\mu$ of $\mathcal{G}$,

[TABLE]

According to the mean value theorem Bogachev (2007), there exists a function $\bar{G}_{i}(P)$ with $0\leq\bar{G}_{i}(P)\leq 1$ such that

[TABLE]

Let $u^{\ast}(P)=\int\mathrm{d}\mu(g)u[V^{\dagger}(g)PV(g)]$, which is obviously $R_{V}-$covariant. Then

[TABLE]

To see that $\bar{G}$ satisfies the normalisation, we sum (23) over $i$:

[TABLE]

which means $\sum_{i=1}^{n}\bar{G}_{i}(P)=1$ almost everywhere with measure generated by $u^{\ast}(P)$. Therefore $\bar{G}$ is a proper response function for measurement $E$ with LHS ensemble $u^{\ast}$. Thus $\rho$ is also $u^{\ast}$-unsteerable. ∎

One can easily check that both proofs work equally well when the measurements are restricted to PVMs. To better understand the relation of Theorem 2 with Lemma 1 of Ref. Wiseman et al. (2007), we come back to the example of the two-qubit Werner state $W_{p}$ in the main text. According to Lemma 1 of Ref. Wiseman et al. (2007), if $W_{p}$ is unsteerable, there exists a covariant indexed LHS ensemble for which $W_{p}$ is unsteerable, and then it is deduced that this singles out the uniform distribution over the hidden variables as the “optimal” LHS ensemble $u^{\ast}$. In fact, there exist infinitely many indexed ensembles that are covariant under the $\mathrm{U}(2)$ action. For example, consider and index space $\Lambda=Z_{K}\times\mathcal{S}_{B}$ with $Z_{K}=\{0,1,2,...,K-1\}$, i.e., we have a composite hidden variable $\lambda=(\alpha,P)$ with $\alpha=0,1,2,...,K-1$ and $P\in\mathcal{S}_{B}$. The measure $\nu$ on $\Lambda$ is generated by the distribution $u$ on $\Lambda$, defined by $u(\alpha,P)=c_{\alpha}$ independent of $P$ with $c_{\alpha}\geq 0$, $\sum_{\alpha=0}^{K-1}c_{\alpha}=1$.

The action of $g\in\mathrm{U}(2)$ on $\Lambda$ can be defined as $g\lambda=g(\alpha,P)=(\alpha,gPg^{\dagger})$. This action is not transitive on $\Lambda$; it has $K$ orbits indexed by $\alpha$. As a result, despite the fact that $u$ is $\mathcal{G}$-invariant, $u(\lambda)=u(g\lambda)$, it is not uniform over $\Lambda$ if $c_{\alpha}$ are distinct numbers.

Now consider the indexed LHS ensemble given by the indexing function $F:\Lambda\to\mathcal{S}_{B}$, $(\alpha,P)\mapsto F(\alpha,P)=P$ (which is many-to-one). The LHS ensemble is apparently covariant, that is, $F(\lambda)=g^{\dagger}F(g\lambda)g$ or even $u({\lambda})F(\lambda)=u(g\lambda)g^{\dagger}F(g\lambda)g$.

Being unsteerable with respect to this indexed LHS ensemble implies that for any $n$-POVM $E$, there exist response functions $G_{i}(\alpha,P)$ such that

[TABLE]

where $S$ is surface measure on Bob’s Bloch sphere. On the face of it, this does not imply that one can choose the uniform distribution on Bob’s Bloch sphere to be the LHS ensemble. The latter requires that there exists response function $\bar{G}_{i}(P)$ such that

[TABLE]

The existence of $\bar{G}_{i}(P)$ only follows upon applying the mean value theorem as in Proof 2, which states that there exist $\bar{G}(P)$ such that

[TABLE]

More complicated examples can be easily constructed by replacing $Z_{K}$ with a more complicated measurable space. Corollary 1 then implies that all these different constructions are actually equivalent when one concerns with simulating steering assemblages, since they generate the same uniform distribution on the Bloch sphere. Thus one sees that Lemma 1 of Ref. Wiseman et al. (2007), when augmented with Corollary 1, can also identify the uniform distribution as the optimal choice for LHS ensemble as stated directly in Theorem 2.

Appendix E Simulated annealing and computation of the gap function

Simulated annealing is a standard heuristic algorithm to solve a generic global optimisation problem Press et al. (2007). In our case, we wish to compute

[TABLE]

with

[TABLE]

where $S$ is the sureface measure of the Bloch sphere (which is different from the Haar measure by a factor $\frac{1}{4\pi}$). The simulated annealing algorithm goes as follows. One first regards $F(Z,E)$ as an energy function of a system in the state space $(Z,E)$. Simulated annealing couples this system to an effective heat bath, whose temperature is then lowered slowly, so that configurations with decreasing energy are explored. At each temperature the system follows stochastic dynamics leading to equilibrium with the heat bath. The system is cooled down slowly to sufficiently small temperature $T_{f}$. It is known that if the temperature schedule is sufficiently slow then the system converges to a global minimum of $F(Z,E)$ Granville et al. (1994). However, the required cooling schedule is too slow that it is not useful in practice and an alternative cooling schedule is used. Here we use an exponential cooling scheme, i.e., in each step the temperature is cooled down by a factor $f$. The system can in principle become stuck in a local minimum at $T_{f}$. It is then necessary to repeat the cooling procedure multiple times.

Coordinisation of variables. Note that for two qubit systems, $\mathcal{A}=\mathcal{B}=M(\mathbb{C},2)$, where $M(\mathbb{C},2)$ is the algebra of $2\times 2$ complex matrices. We use the Pauli basis $\{\sigma_{i}\}_{i=0}^{3}=\{\mathbb{I},\sigma_{x},\sigma_{y},\sigma_{z}\}$ to coordinate the real subspace $M^{H}(\mathbb{C},2)$. Each operator $X\in M^{H}(\mathbb{C},2)$ is therefore characterised by $4$ (real) coordinates $x_{i}$,

[TABLE]

The boundary of the positive cone of $M^{H}(\mathbb{C},2)$ is given by $x_{0}^{2}-x_{1}^{2}-x_{2}^{2}-x_{3}^{2}=0$ with $x_{0}\geq 0$, consisting of vectors of the form $\alpha\begin{pmatrix}1\\ \boldsymbol{n}\end{pmatrix}$ with $\alpha\geq 0$.

The composite operator $Z=\oplus_{i=1}^{4}Z_{i}$ and $E=\oplus_{i=1}^{4}E_{i}$ are thought of as $4\times 4$ matrices, in which each column is the coordination of $Z_{i}$ and $E_{i}$ respectively. From now on, we will use $Z$ and $E$ to denote these matrices, and $Z_{i}$, $E_{i}$ to denote the $i^{\text{th}}$ column.

To implement the constraint $\mathcal{C}^{4}$ to $Z$, we write $Z=XR$, where

[TABLE]

and $X$ satisfies $\operatorname{Tr}(X^{T}X)=2$, $X_{4}=(0,0,0,0)^{T}$.

To implement the constraint $\mathcal{N}^{4}$ on $E$, we note that every component $E_{i}$ is on the boundary of the positive cone, $E_{i}=\alpha_{i}\begin{pmatrix}1\\ \boldsymbol{n}_{i}\end{pmatrix}$, where $\alpha_{i}\geq 0$ and $\boldsymbol{n}_{i}$ is a unit vector. The constraint $\sum_{i=1}^{4}E_{i}=\mathbb{I}$ can be considered as a constraint for $\{\alpha_{i}\}_{i=1}^{4}$. In fact, if $\begin{pmatrix}1\\ \boldsymbol{n}_{i}\end{pmatrix}$ are independent, $\{\alpha_{i}\}_{i=1}^{4}$ are uniquely determined. Note that the set of $E$ where $\begin{pmatrix}1\\ \boldsymbol{n}_{i}\end{pmatrix}$ are dependent are zero-measured in $\mathcal{N}^{4}$. In practice, we therefore do not need to worry about the case that $\begin{pmatrix}1\\ \boldsymbol{n}_{i}\end{pmatrix}$ are dependent if the linear solver is relatively stable. Here we use the Householder linear solver, provided by Eigen 3 Guennebaud et al. (2010). Further, we use a common technique Press et al. (2007) to take care of the constraint $0\leq\alpha_{i}\leq 1$ by assigning infinite values to the energy if the solution $\alpha_{i}$ are outside $[0,1]$. This allows us to describe the (extended) set $\mathcal{N}^{4}$ by $4$ vectors $\{\boldsymbol{n}_{i}\}_{i=1}^{4}$.

Annealing. At temperature $T$, the stochastic dynamics of the system is simulated by the Metropolis algorithm Press et al. (2007): at every time step, the system tries an elementary step, which will be accepted with probability $\min\{1,e^{-\Delta F/T}\}$, where $\Delta F$ is the change in energy due to the trial step. In each elementary trial step, either $Z$ or $E$ is updated with equal probability. If $Z$ is updated, we choose randomly two elements of $X$, say $X_{ij},X_{kl}$, where $j,l<4$, and perform a rotation $Q(\theta)\in SO$(2) on the vector $(X_{ij},X_{kl})^{T}$ by a random angle $\theta$ normally distributed with mean [math] and standard deviation $2\pi\sqrt{T}$. The components of the vector $Q(\theta)(X_{ij},X_{kl})^{T}$ replace the $ij$ and $kl$ elements of $X$. Note that the constraints $\operatorname{Tr}(X^{T}X)=2$ and $X_{4}=(0,0,0,0)^{T}$ are respected in the new $X$. Then $Z$ is updated as $Z=XR$. If $E$ is updated, we choose one of the vectors $\boldsymbol{n}_{i}$ randomly, and rotate it around one of the $3$ axes $x$, $y$, $z$ by a random angle normally distributed with mean [math] and standard deviation $2\pi\sqrt{T}$. At each temperature, the number of the simulated steps are at least $100$ times the degree of freedom.

Cooling schedule. After annealing the system at temperature $T$, the temperature is decreased by a factor $f$; here $f=0.95$. This is known as exponential temperature scheduling Press et al. (2007). The initial temperature $T_{i}$ is chosen as the maximal value minus the minimal value of the energy function sampled at $1000$ times the degree of freedom points. The algorithm is stopped at temperature $T_{f}=10^{-9}$. Lowering the final temperature does not significantly improve the results.

Replicas. We repeat the cooling procedure $M=512$ times. 64 such typical cooling trajectories are presented in Figure 4. As shown in Figure 5, all $512$ replicas produce a very similar minimum energy, suggesting that there is no major local minimum in the energy landscape. This provides confidence that the system indeed converges close to a global minimum.

Numerical accuracy. As seen in Figure 4, the gap can be overestimated by some order of $10^{-4}$. This makes it difficult to study very small gap values when $p$ is around $\frac{1}{2}\pm 10^{-4}$. Other cooling schedules Press et al. (2007), parallel tempering Swendsen and Wang (1986) or more subtle global optimisation techniques Arora et al. (1995) can be considered to increase the accuracy. However, the numerical accuracy is limited by another critical factor: the accuracy of the spherical integral in (30). Here, we used Lebedev’s quadrature with $5810$ points to compute spherical integrals. This effectively replaces the optimal uniform LHS ensemble by a suboptimal discrete distribution at $5810$ quadrature points. Accordingly, the expected transition probability $p$ is shifted by some value of $2\times 10^{-4}$ to the left of $\frac{1}{2}$, smearing out the accuracy of the simulated annealing optimisation as seen in Figure 6. Similar problems occur for $T$-states at a resolution of $\left|\epsilon\right|\approx 10^{-3}$ around the surface of unsteerable states.

Appendix F The gap function for steering the Werner states with PVMs

The numerical calculation of the gap function for the Werner states with PVMs is carried out similarly to Appendix E. The analytical calculation of the gap function for steering a Werner state with PVMs is rather straightforward. We are to calculate

[TABLE]

where $(P_{1},P_{2})$ forms a projective measurement, i.e., $P_{1}$, $P_{2}$ are orthogonal projections such that $P_{1}+P_{2}=\mathbb{I}$.

Since $(Z_{1},Z_{2})\in\mathcal{C}^{2}$ implies that $Z_{1}+Z_{2}=0$, we can set $Z_{1}=X$ and $Z_{2}=-X$. Moreover, because of the $\operatorname{U}(2)$ symmetry of the problem, we can suppose $X=\lambda_{0}\left|0\right\rangle\left\langle 0\right|+\lambda_{1}\left|1\right\rangle\left\langle 1\right|$ with $\lambda_{0}\geq\lambda_{1}$. Because $\left\langle Z_{1},Z_{1}\right\rangle+\left\langle Z_{2},Z_{2}\right\rangle=1$, we have $\left\langle X,X\right\rangle=\frac{1}{2}$, or $\lambda_{0}^{2}+\lambda_{1}^{2}=\frac{1}{2}$.

If we write the projections in Pauli coordinates as $P=\begin{pmatrix}1\\ \boldsymbol{n}\end{pmatrix}$, then $\left\langle Z_{1},P\right\rangle\geq\left\langle Z_{2},P\right\rangle$ is equivalent $n_{z}\geq-\frac{a}{b}$, with $a=\lambda_{0}+\lambda_{1}$ and $b=\lambda_{0}-\lambda_{1}$. Therefore

[TABLE]

which evaluates to $\frac{1}{4b}$.

On the other hand

[TABLE]

which evaluates to $\frac{pb}{2}$.

Thus we have $\Delta=\min_{b}\left\{\frac{1}{4b}-\frac{pb}{2}\right\}$. Note that $\lambda_{0}^{2}+\lambda_{1}^{2}=\frac{1}{2}$ implies that $a^{2}+b^{2}=1$, and so $b\leq 1$. Hence we obtain the gap function for PVMs as $\Delta=\frac{1}{4}-\frac{p}{2}$.