The distribution and verification of quantum nonlocality across a network of users is essential for future quantum information science and technology applications. However, beyond simple point-to-point protocols, existing methods struggle with increasingly complex state preparation for a growing number of parties. Here, we experimentally demonstrate multiparty detection-loophole-free quantum steering, where one party simultaneously steers multiple spatially separate parties, using a multi-qubit state produced from a set of qubits of which only one pair is entangled. We derive loss-tolerant steering inequalities, enabling our experiment to close the detection loophole, and enabling us to show the scalability of this approach to rigorously verify quantum nonlocality across arbitrarily many parties. This provides practical tools for developing the future quantum internet.

Quantum nonlocality is a resource for secure communications and distributed information tasks.1–3 The latter include distributed quantum computing,4 quantum cryptography,5–7 randomness certification,8–10 quantum state teleportation,11,12 and long-range sensor nets such as extended baseline optical telescopes.13 

Nonlocal communication protocols prevent eavesdroppers or malicious parties from sabotaging or gaining information from sensitive communication and guarantee unconditional security.14–16 This is done by having networks of separate parties (usually two) measure correlations and violate a Bell inequality or steering inequality. Quantum (or Einstein–Podolsky–Rosen) steering17,18 is a form of nonlocality which distinguishes itself from Bell nonlocality in several ways.3,19 Notably, by only trusting a subset of parties, quantum steering is more robust to loss and noise on the untrusted channels.20–22 

Steering is often studied in point-to-point scenarios, but more than two spatially separate parties are required for most networking applications. Establishing multiparty quantum networks requires understanding new and more complex nonlocal and causal relationships beyond the well-studied two-party scenario described by Bell inequalities.23 For quantum steering in a multiparty scenario, many network topologies can arise, as any observer in the network could play the role of a trusted or untrusted party within a given steering task. This leads to novel nonlocality phenomena such as network quantum steering24 or collective steering,25–28 where subsets of parties jointly attempt to steer other subsets.

Increasing the number of parties in multiparty steering by directly extending traditional approaches is daunting, as it entails scaling up the quantity or complexity of entangled resources. Even a simple scenario in which one (untrusted) Alice steers N (trusted) Bobs, where Alice and each of the Bobs share an entangled pair, requires 2N total qubits.29 Alternative approaches that reduce this number rely on complex quantum states with entanglement depth larger than 2,30 such as N-qubit W, GHZ, cluster or graph states,31–34 or multi-mode squeezing.26 A different approach involves using a single entangled pair and sequential weak measurements.35–37 These works require all trusted parties to share a common quantum channel and also rely on post-selection, significantly increasing network complexity, and preclude nonlocal behavior from being observed simultaneously.

Moreover, real-world implementations require ruling out cheating strategies that untrusted parties may utilize, especially those that exploit known loopholes. The detection loophole is one such example, that is opened if the trusted party makes a fair sampling assumption on the statistics governing reported null results by the other party. Unlike with other loopholes, accounting for photonic loss requires a more sophisticated theoretical treatment. This is crucial, as failing to close the detection loophole entirely compromises the device-independent properties which are unique to the steering task.

In this work, we experimentally overcome the problems of growing entangled resource needs by utilizing a scalable and resource-efficient approach to multiparty steering. Following one proposal in Ref. 38, we build a multi-photon state from a single entangled photon pair plus other photons in a product state. We develop a theory for verifying multiparty steering with these resource-efficient states in the presence of loss, exploring a scenario where one party steers an arbitrary number N of spatially separated parties simultaneously. Our experiment demonstrates simultaneous quantum steering of N = 2 trusted parties with the detection loophole closed. Moreover, we predict steerability of up to N = 26 trusted parties with the noise of our currently produced states and prove loss-tolerant steerability for arbitrarily many parties with ideal detectors and only three measurement settings.

Consider a network involving N + 1 spatially distant parties, composed of one untrusted Alice and N trusted Bobs, as in Fig. 1. The parties in this network share a quantum state such that each observer has one out of N + 1 qubits. Through a quantum steering test, Alice attempts to convince each of the Bobs simultaneously that she shares nonlocality with them. The test result is evaluated based on measurement outcomes reported by the parties over repeated protocol runs.

Fig. 1.

Multiparty quantum steering scenario. N + 1 parties share a quantum state. (Untrusted) Alice receives measurement instructions x and reports outcomes a. Each (trusted) Bob makes Pauli projections from a tomographically complete set in each protocol run and can independently verify if Alice's measurements can steer their respective particle.

Fig. 1.

Multiparty quantum steering scenario. N + 1 parties share a quantum state. (Untrusted) Alice receives measurement instructions x and reports outcomes a. Each (trusted) Bob makes Pauli projections from a tomographically complete set in each protocol run and can independently verify if Alice's measurements can steer their respective particle.

Close modal

Since Alice is untrusted, we treat her measurement device as a black box, making no assumptions about how her outcomes are generated. Alice receives classical instructions labeled by x, specifying which out of a set of predetermined measurements to perform in each protocol run, and she broadcasts the corresponding outcomes a { + , , } to the other N parties. Here, represents the null outcome, corresponding to an event where Alice has received a measurement instruction but reports no outcome. The probability of reporting a non-null outcome per protocol run is called Alice's efficiency.

In every protocol run, each Bob (labeled Bn) can perform a projective measurement on his qubit (from a tomographically complete set). Here, no option of a null outcome is required because the trusted Bobs—who collectively decide what constitutes a run of the protocol—exclude instances where their measurement devices do not report an outcome. Over time, each Bob sorts his measurement statistics by Alice's announced results, thus creating ensembles of locally observed quantum states normalized by the probability of Alice observing a corresponding measurement outcome. These ensembles are commonly known as assemblages.19 Importantly, an assemblage contains all the information relevant to deciding the result of a quantum steering test.

For a given Alice–Bobs' bipartite state ρ, the n th Bob's assemblage is a collection of unnormalized quantum states σ a | x B n = Tr A [ ( E a | x I ) ρ A , B n ]. Here, ρ A , B n : = Tr ¬ A , B n [ ρ ] is the reduced system shared by Alice and Bob n, while { E a | x } a is the positive operator-valued measure (POVM)—over a—for each of Alice's settings x. Any bona fide assemblage must contain only positive semi-definite matrices that satisfy the no-signaling condition, x , x , a σ a | x B n = a σ a | x B n = : ρ B n. Here, we employ a convex optimization technique to reconstruct the assemblages ( Appendix C), which ensures these properties.

Alice has steered the n th Bob if every element of the assemblage cannot be written as a coarse-graining over local-hidden-states { ρ λ B n } λ detected by the n th Bob.17 That is, there does not exist a set of normalized states { ρ λ B n } λ, and probability distributions p λ and P λ ( a | x ) over λ and a, respectively, such that
a , x , σ a | x B n = λ p λ P λ ( a | x ) ρ λ B n .
(1)

If Alice's efficiency is below unity, a common approach is to postselect by ignoring the null outcomes and introducing the fair sampling assumption that Alice's reported measurement outcomes accurately represent the total statistical sample. However, as in any rigorous nonlocality test, such experimental assumptions open loopholes, allowing false nonlocality verification.

The loophole associated with the fair sampling assumption is the detection loophole; it allows Alice to cheat by not reporting some of her measurement outcomes to mimic a steerable assemblage. Alice's null outcome instances need to be taken into account to close this loophole and drop the need for the fair sampling assumption. This amounts to placing a lower bound on Alice's efficiency, which she must surpass to steer Bob via her non-null outcomes.

To close this loophole, our analysis will use inequalities of the following form, whose satisfaction certifies detection-loophole-free steering:
ε ( { σ a | x B n } ) < ε exp .
(2)
Here, ε exp is Alice's measured efficiency—the proportion of her non-null results—in the experiment, while ε is a cutoff efficiency. This is a function of the assemblage { σ a | x B n } a | x observed by the nth Bob when Alice announces any non-null outcome, i.e., a . The cutoff efficiency ε is the maximum efficiency that could allow a cheating Alice to exploit the fair sampling assumption. This cutoff efficiency can be computed through a single instance of a semi-definite program (SDP) (see  Appendix A).

Preparing an appropriate multiparty entangled quantum state is crucial to realizing a loss-tolerant multiparty steering demonstration. Here, we construct a practical quantum state that can demonstrate loophole-free steering from Alice to each of the Bobs independently, based on a single entangled pair of qubits and N – 1 pure ancillas. The ideal (lossless) theory of the type of entanglement we will create was recently introduced in Ref. 38, where it was referred to as semi-random pair entanglement. Here, we introduce a new theory for verifying steering in the presence of loss, which facilitates a robust demonstration of steering in this scenario with no detection loophole, as well as to prove loss-tolerant steerability for arbitrarily many parties.

The mixed global quantum state we prepare consists of N + 1 qubits, acting on a composite Hilbert space H A H B 1 H B N. Below, we omit system labels when there is no risk of confusion. Consider the family of two-qubit pure states | ψ α : = α | 11 + 1 α | 00 , that is entangled for α ( 0 , 1 ). From | ψ α , we construct the family of states between the N + 1 parties
ρ α , N : = 1 N n = 1 N V 1 , n ( | ψ α ψ α | | 0 0 | N 1 ) ,
(3)
where V m , n is a superoperator that swaps subsystems Bm and Bn: V m , n ( ρ B m ρ B n ) ρ B n ρ B m. As we show later, this state permits loophole-free quantum steering with an arbitrary number of parties and can be prepared using deterministic quantum gates.
For now, we concentrate on the case of N = 2, where we obtain the three-party state
ρ α , 2 : = 1 2 ( | ψ α ψ α | A , B 1 | 0 0 | B 2 + | ψ α ψ α | A , B 2 | 0 0 | B 1 )
(4)
distributed between the untrusted party, Alice, and the trusted parties, Bob 1 and Bob 2, as indicated by the system labels. The state represents a mixture of two cases in which one part of | ψ α remains at Alice's station, and the other is sent to either Bob 1 or Bob 2, with the remaining Bob receiving a single photon in the computational zero state. The cutoff efficiency ε ( { σ a | x B n } ) of the assemblage produced for Bob 1 and Bob 2, assuming Alice performs three measurements corresponding to Pauli operators, is illustrated as the purple solid line in Fig. 3. Quantum steering is in principle possible whenever ε ( { σ a | x B n } ) is below one. This occurs for α ( 0 , 2 / 3 ) and the minimal value of ε is achieved in the singular limit α 0, i.e., at vanishing entanglement of | ψ α (as measured, e.g., by concurrence).

Intuitively, this behavior is a result of saturating the steered states of | ψ α and the ancillas close to the boundary of space of quantum states. More specifically, when Alice measures the three Pauli observables, the steered states for each Bob [see Eqs. (E2)–(E4) of  Appendix E] are contained within a ball of radius 1 / N in the Bloch ball, which touches its surface once, at | 0 . Varying α changes the probabilities (and positions) with which the steered states appear inside this ball. For such an arrangement of steered states, the most difficult assemblages (for a given efficiency) to construct an LHS model for are those where the purest states appear with the highest probability, since the LHS ensemble must comprise strictly physical quantum states.

When the family of states of Eq. (4) is modified slightly through the addition of noise, the singularity disappears, as shown by the purple dashed line in Fig. 3. Although counterintuitive, this finding that a small value of α is optimal aligns well with the fact that for the case of two parties, pure states with a small amount of entanglement exhibit detection-loophole-free nonlocality at a lower efficiency bound than a maximally entangled state in the cases of both Bell nonlocality39 and quantum steering40 tests.

We perform a photonic experiment, encoding qubits in the polarization degree of freedom of single photons. Our experimental setup is shown in Fig. 2. To prepare the state of Eq. (4), we first generate an entangled photon pair in the state | ψ α and a heralded single photon using two high-efficiency, telecom-wavelength photon pair sources based on group-velocity-matched spontaneous parametric down-conversion with a pulsed pump.41 One half of the entangled pair is sent to Alice, and the other half is probabilistically distributed to one of the trusted parties, Bob 1 or Bob 2, with the other Bob receiving the heralded single photon. Alice is instructed to perform one out of a set of three projective measurements on her photon. The Bobs each perform a quantum state tomographic measurement on their respective photon.

Fig. 2.

Experimental layout. Sources: the top photon source produces the entangled state | ψ α from Eq. (4) and distributes half to Alice and the other to the Swap stage. The bottom source produces the unentangled | 0 , sending one photon to the Swap stage and the other to the heralding detector. Swap: mirrors mounted on a linear translation stage displace the beam's trajectory with 0.5 probability, producing the state's mixture. While the swap operation occurs (mirror positions are shown in the figure), the light follows the dotted lines; otherwise, the mirrors are removed from the beam path, and the light will follow the dashed line. Path segments common in both cases are solid lines. Alice, Bob 1, and Bob 2: each party perform projective measurements with motorized wave plates and a polarizing beam splitter. Coincidences are measured with superconducting nanowire single-photon detectors (SNSPD). See  Appendix G for more details.

Fig. 2.

Experimental layout. Sources: the top photon source produces the entangled state | ψ α from Eq. (4) and distributes half to Alice and the other to the Swap stage. The bottom source produces the unentangled | 0 , sending one photon to the Swap stage and the other to the heralding detector. Swap: mirrors mounted on a linear translation stage displace the beam's trajectory with 0.5 probability, producing the state's mixture. While the swap operation occurs (mirror positions are shown in the figure), the light follows the dotted lines; otherwise, the mirrors are removed from the beam path, and the light will follow the dashed line. Path segments common in both cases are solid lines. Alice, Bob 1, and Bob 2: each party perform projective measurements with motorized wave plates and a polarizing beam splitter. Coincidences are measured with superconducting nanowire single-photon detectors (SNSPD). See  Appendix G for more details.

Close modal

Bob 1 and Bob 2 do not need to share a quantum channel but require the other party to indicate through a classical channel if they obtain an outcome. In the cases where both Bobs announce that they have obtained an outcome, the run goes ahead, regardless of whether Alice announces a measurement outcome or claims to have lost her photon.

To build up statistics, the protocol is repeated for each of Alice's measurement settings and various tomographic projections for the Bobs. Afterward, Bob 1 and Bob 2 reconstruct and analyze their assemblages, determine Alice's efficiency based on the proportion of runs she reported an outcome, and test the inequality of Eq. (2). Alice's efficiency includes all losses associated with her photon, from state preparation through to her detection efficiency. Further experimental details are provided in the Appendices.

The results of steering tests for different states are shown as the data points in Fig. 3. Multiparty steering is demonstrated when ε exp > ε Bob 1 and ε exp > ε Bob 2 , which is clearly observed at α values between 0.065 and 0.3. The most statistically significant steering occurs at α 0.1, where we measured a minimum efficiency 5.34 (5.35) standard deviations above the cutoff for Bob 1 (Bob 2). This is thanks to our high experimental heralding efficiencies above 0.69. Since the cutoff efficiencies are numerically determined from the experimentally reconstructed assemblages, the steering demonstrations are conclusive independent of how well the experimental states approximate the target states of Eq. (4).

Fig. 3.

Data demonstrating detection-loophole-free steering as a function of the quantum state parameter α, for N = 2. In the experiment, bipartite quantum steering takes place when ε exp > ε BobN . The green (blue) points illustrate the minimum efficiency required to demonstrate steering for the group Alice–Bob 1 (Alice–Bob 2). The orange squares correspond to Alice's minimum measured efficiency across her measurement settings, and the black dotted line is the mean across the different α values, with the gray-shaded region illustrating ± one standard deviation across all measured efficiencies. All data points include horizontal error bars, though frequently smaller than the data marker. The purple solid line depicts the theoretical lower bound for the cutoff efficiency, ε for the ideal target family of states from Eq. (4). The purple dashed line illustrates numerically determined cutoff efficiencies ε (see  Appendix A) for a more realistic family of states, consistent with experimental tomographic reconstructions of the target states before the multiparty steering test; these are well-modeled by the action of the depolarizing channel on each target, Δ η ( ρ ) = η ρ + ( 1 η ) I / 4 with η = 0.9931.

Fig. 3.

Data demonstrating detection-loophole-free steering as a function of the quantum state parameter α, for N = 2. In the experiment, bipartite quantum steering takes place when ε exp > ε BobN . The green (blue) points illustrate the minimum efficiency required to demonstrate steering for the group Alice–Bob 1 (Alice–Bob 2). The orange squares correspond to Alice's minimum measured efficiency across her measurement settings, and the black dotted line is the mean across the different α values, with the gray-shaded region illustrating ± one standard deviation across all measured efficiencies. All data points include horizontal error bars, though frequently smaller than the data marker. The purple solid line depicts the theoretical lower bound for the cutoff efficiency, ε for the ideal target family of states from Eq. (4). The purple dashed line illustrates numerically determined cutoff efficiencies ε (see  Appendix A) for a more realistic family of states, consistent with experimental tomographic reconstructions of the target states before the multiparty steering test; these are well-modeled by the action of the depolarizing channel on each target, Δ η ( ρ ) = η ρ + ( 1 η ) I / 4 with η = 0.9931.

Close modal
As previously mentioned, our steering scenario can be extended to N + 1 parties by considering the generalized state ρ α , N in Eq. (3). A method for creating the states is illustrated in Fig. 4, where a sequence of N − 1 deterministically implementable gates successively acts on one-half of the entangled pair and one of the ancilla qubits, which are initialized to | 0 . Each gate is a random-swap gate, a completely positive trace-preserving map between states acting on the Hilbert spaces of B1 and Bn, Φ n : = p n V 1 , n + ( 1 p n ) I 1 , n. Here, I n is the identity superoperator, and p n : = 1 / ( N n + 1 ) is the swapping probability for the gate. The ordered composition of these gates gives the desired output of the circuit,
ρ α , N Φ N 1 ° ° Φ n ° ° Φ 1 ( | ψ α ψ α | | 0 0 | N 1 ) .
(5)
Fig. 4.

Quantum circuit diagram for the construction of the global N + 1 qubit state Eq. (3). One half of the entangled pair | ψ α is sent to Alice, and the other is randomly sent to one of the N Bobs through a sequence of deterministically implementable random-swap gates.

Fig. 4.

Quantum circuit diagram for the construction of the global N + 1 qubit state Eq. (3). One half of the entangled pair | ψ α is sent to Alice, and the other is randomly sent to one of the N Bobs through a sequence of deterministically implementable random-swap gates.

Close modal

This state can be created in a resource-efficient way and with several appealing properties.

  • Number of qubits: Each party only requires one qubit, unlike the simple extension of two-party steering where an entangled pair of qubits gets distributed between Alice and each of the Bobs, which would involve 2N qubits.

  • 2-Producibility: Multiparty quantum states can be characterized using the concept of k-producibility.30 A pure quantum state | ψ is k-producible if it can be written as a composition of quantum states involving at most k parties. That is, | Ψ = | ψ 1 | ψ 2 | ψ 3 where each | ψ i is a state shared between k parties. Similarly, a mixed state is k-producible if it can be decomposed as a mixture of k-producible pure states. If a state is not (k − 1)-producible, but is k-producible, then its entanglement depth is k. The states from Eq. (3) are 2-producible and have entanglement depth 2, independently of the number of parties involved—unlike N + 1-party GHZ and W states which both have entanglement depth N + 1.

  • Number of entangled pairs: The state preparation only requires a single pair of entangled qubits, which is a stronger constraint on the required resources than 2-producibility alone.

  • Deterministic implementation of gates: The creation of the states involves gates that can be implemented on photons deterministically, in the sense that the gates do not require postselection or heralding, unlike controlled-NOT gates.

To understand how this class of states permits loophole-free multiparty steering, we can derive an exact expression for Alice's cutoff efficiency ε , when she measures three dichotomic Pauli observables. Its closed form, presented in Eq. (A3) of  Appendix A and derived in  Appendix E, is directly a function of N, and the degree of entanglement, α, present in the initial entangled pair.

From a theoretical point of view, this efficiency cutoff has various interesting properties. Detection-loophole-free steering could be observed when this expression is strictly below unity, which occurs for α ( 0 , 2 / ( N + 1 ) ). Moreover, we show in Sec. 3 of  Appendix E that the minimal value of ε is achieved in the singular limit α 0, the same limit in which the concurrence of | ψ α vanishes. This generalizes the earlier behavior evident from the N = 2 case (purple solid curve in Fig. 3). In this limit, the required detection efficiency is only ε = ( 1 + 2 / N ) 1, which is quite striking for a state produced from a single pair of entangled qubits and only three measurement settings.

Importantly, these bounds demonstrate that there always exists a value of Alice's efficiency, which permits arbitrarily many Bobs to be steered in a loophole-free way, using our construction. However, this is only true for an ideal scenario where noise is absent from the experiment. Therefore, we also investigate the scalability of steering for the states defined by Eq. (3), but with the addition of noise at the level we observed in our experiment. These states are well-modeled by a white noise model acting on the entangled state prior to swapping, Δ η ( | ψ α ψ α | ) = η | ψ α ψ α | + ( 1 η ) I / 4 with η = 0.9931. We implement the numerical program from Ref. 42 to certify steerability from Alice to each Bob under all qubit projective measurements. That is, in the limit of an infinite number of measurement settings. Denoting these noisy states by a tilde, for a fixed two-qubit state Δ η ( ρ ̃ α , 2 ), this algorithm computes a quantity called the critical radius R ( Δ η ( ρ ̃ α , 2 ) ). A state that achieves a value R < 1 implies it is steerable under all projective measurements made by Alice. Conversely, R 1 implies it cannot be steered under that set of measurements. Treating α and N as fixed parameters for each computation, we compute an upper (lower) bound on the critical radius using an outer (inner) approximation of the Bloch sphere by a polytope with 1514 vertices. The results of these simulations are illustrated in Fig. 5. The green (purple) line denotes where the upper (lower) bound on R equals unity for given α and N. The region shaded in gray denotes where these upper and lower bounds on R lie either side of 1, so steerability cannot decided, owing to the precision used in the simulation. Remarkably, we find that a judicious choice of α exists such that up to 26 Bobs could be steered simultaneously.

Fig. 5.

The limit to the number of Bobs N which could be steered, using noisy network states. This result is based on the two-qubit state used in the experiment, Eq. (4), with the presence of experimental noise, which is estimated as per the caption of Fig. 3. The simulation here assumes conditions that are, apart from the noise, ideal: unit detection efficiency and infinite measurement settings. The numerically uncertain region is shaded in gray.

Fig. 5.

The limit to the number of Bobs N which could be steered, using noisy network states. This result is based on the two-qubit state used in the experiment, Eq. (4), with the presence of experimental noise, which is estimated as per the caption of Fig. 3. The simulation here assumes conditions that are, apart from the noise, ideal: unit detection efficiency and infinite measurement settings. The numerically uncertain region is shaded in gray.

Close modal

In this work, we experimentally demonstrate a multiparty quantum nonlocality experiment with spatially separated parties in the discrete variable regime—the first such experiment to close the detection loophole. We achieve this by a novel, resource-efficient state preparation scheme that allows multiparty quantum steering. Remarkably, using just a single pair of entangled qubits, one party can, in theory, steer arbitrarily many other parties. Importantly, this approach is robust to the commonly overlooked effects of photon loss and noise: our steering inequality takes these effects into account and can be satisfied under demanding but nevertheless viable conditions. We apply our method to experimentally satisfy a detection-loophole-free quantum steering inequality in a network of three spatially separated parties by 5.35 standard deviations. We thereby demonstrate multiparty steering, where one untrusted party simultaneously steers two trusted parties. We show steering of up to N = 26 is possible with realistic experimental quantum state fidelities and ideal efficiencies. Unlike other methods, our approach is scalable as it does not require heralded or postselected gates to generate the steerable state, and neither does it require increasing entanglement for larger numbers of parties. Our protocol does not rely on sequential measurements, such as those used in continuous variable protocols. Thus, we do not require a quantum channel connecting our trusted parties. In the experiment, we focused on closing the detection loophole, which is crucial for real-world implementations.

Verification of quantum nonlocality is essential for the implementation of secure quantum networks. For the experimental design and data analysis, we used semi-definite programing to determine steering bounds and novel techniques to reconstruct quantum assemblages using maximum-likelihood estimation; these can also find wider applications in other quantum steering contexts.

Future directions include the implementation of a fully loophole-free protocol, implementing a larger network with more parties, and including a larger number of untrusted elements in various topologies. Our work demonstrates a realistic method for steering verification in a large-scale quantum network. We show how steering-based quantum networks of tens of users can be implemented. From a secure quantum communication application side, our protocol allows a trusted network to introduce an untrusted member, which may be useful in user authentication, such as banking, multi-factor authentication, and implementations of a quantum internet.

This work was supported by the Australian Research Council Centre of Excellence CE170100012 and the National Natural Science Foundation of China (No. 12288201). A.P., T.J.B., and Q.C.S. acknowledge support from the Australian Government Research Training Program (RTP). The authors acknowledge the support of the Griffith University eResearch Service and Specialized Platforms Team and the use of the High-Performance Computing Cluster “Gowonda” to complete this research.

The authors have no conflicts to disclose.

Alex Pepper: Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Software (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). Travis J. Baker: Conceptualization (supporting); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). Yuanlong Wang: Formal analysis (equal); Funding acquisition (supporting); Methodology (supporting); Software (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Qiu-Cheng Song: Formal analysis (supporting); Software (supporting); Writing – review & editing (supporting). Lynden K. Shalm: Resources (equal); Writing – review & editing (equal). Varun B. Verma: Resources (equal); Writing – review & editing (equal). Sae Woo Nam: Resources (equal); Writing – review & editing (equal). Nora Tischler: Investigation (supporting); Methodology (supporting); Project administration (supporting); Supervision (equal); Validation (equal); Writing – review & editing (equal). Sergei Slussarenko: Investigation (equal); Methodology (supporting); Project administration (supporting); Supervision (equal); Validation (equal); Writing – review & editing (equal). Howard M. Wiseman: Conceptualization (lead); Funding acquisition (supporting); Project administration (lead); Supervision (lead); Validation (equal); Writing – review & editing (equal). Geoff J. Pryde: Funding acquisition (lead); Project administration (equal); Resources (equal); Supervision (lead); Validation (supporting); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and data are available from corresponding author upon reasonable request.

Since imperfections such as noise are experimentally unavoidable and lead to deviations from the ideal behavior, we follow a multi-step process to adapt Alice's measurement directions and determine the correct cutoff efficiency ε . In the design stage prior to the steering test, we aim to prepare seven members of the ideal family of two-qubit states | Ψ α and perform quantum state tomography to characterize the level of noise. We find the average fidelity between the prepared and ideal states to be 0.987 ± 0.003. Based on this estimate, we perform a global differential evolution optimization routine to find Alice's measurement directions that minimize the ε we expect to obtain in the steering test. We then perform the steering test using these measurement settings, which provide the assemblages and Alice's efficiency ε exp. The assemblages are used to determine ε via a semi-definite program, and the outcomes of the steering test are decided based on Eq. (2), ε ( { σ a | x B n } ) < ε exp. The left side of this inequality is the maximum efficiency, i.e., average proportion of non-null results, that a dishonest Alice could have announced in an experiment preparing the assemblages { σ a | x B n }. This quantity is found by solving the SDP:
max { σ λ } λ ε , s . t . λ D ( a | x , λ ) σ λ = ε σ a | x B n , x , a , λ σ λ = ρ B n , σ λ 0 , λ .
(A1)
Here, { D ( a | x , λ ) } λ is the set of deterministic probability distributions mapping x to all outcomes (both null and non-null). This optimization problem can be interpreted as a functional that maps from the space of assemblages to real numbers in ( 0 , 1 ]. For a given input assemblage, exact values for ε can be found using efficient numerical solvers.
When the assemblage is prepared by Alice measuring on her part of ρ α , 2, the assemblage (and hence ε ) is a function of the measurements she chooses to make. Assuming these are dichotomic projective measurements { Π a | x } a , x on her qubit, we decompose them as Π a | x = ( I + a r x · σ ) / 2, with r x : = ( sin θ x cos ϕ x , sin θ x sin ϕ x , cos θ x ) and σ : = ( X , Y , Z ) is the vector of Pauli matrices. We then wish to find Alice's minimum detection efficiency by searching over these measurement settings. That is, for various values of α, we seek solutions to
min { Π a | x } a , x ε ( { Tr A [ ( Π a | x I ) ρ ̃ α , 2 ] } ) .
(A2)
The quantum state ρ ̃ α , 2 appearing in the objective function is experimentally determined from averaging over quantum state tomographies to account for noise. We perform this outer minimization over the measurement settings heuristically, using a global differential evolution optimization routine. The results of these searches are summarized in Table I.
Table I.

Alice's numerically found optimal measurement directions. Three measurements are used for each α in | Ψ α . Each direction is expressed as a pair of Euler angles (°) for the Bloch sphere: the angle from the azimuth ( θ n ) and the zenith ( ϕ n ).

α θ1 ϕ 1 θ2 ϕ 2 θ3 ϕ 3
0.015  57.2955  0.859 88  −3.3977  −88.0451  0.7153 
0.065  4.2398  84.1949  −48.6642  −88.781  48.0435 
0.100  39.2478  −60.7293  39.2834  57.8822  −41.0124 
0.185  −33.2033  −70.5399  −42.0172  64.8311  44.2747 
0.300  −1.1167  −89.9633  −2.5873  32.5804  57.2852 
0.400  −57.1865  −53.4757  −9.8684  36.7915  10.1887 
0.500  −13.0359  8.8413  57.2252  89.8365  0.50483 
α θ1 ϕ 1 θ2 ϕ 2 θ3 ϕ 3
0.015  57.2955  0.859 88  −3.3977  −88.0451  0.7153 
0.065  4.2398  84.1949  −48.6642  −88.781  48.0435 
0.100  39.2478  −60.7293  39.2834  57.8822  −41.0124 
0.185  −33.2033  −70.5399  −42.0172  64.8311  44.2747 
0.300  −1.1167  −89.9633  −2.5873  32.5804  57.2852 
0.400  −57.1865  −53.4757  −9.8684  36.7915  10.1887 
0.500  −13.0359  8.8413  57.2252  89.8365  0.50483 
Generally, when Alice attempts to steer N Bobs by performing three dichotomic Pauli measurements on her system, we are able to derive an exact expression for ε ( { σ a | x B n } ), when the noise is absent from the underlying network state. The derivation of this cutoff efficiency is given in  Appendix E. The main result is
ε ( { σ a | x B n } ) : = 2 N α ( N + 1 ) + 2 2 ( 1 α ) N α + α ( N 1 ) ( ( α + 4 ) N 5 α 4 2 ( 1 α ) N α ) 2 ( 2 2 ( 1 α ) N α + N + 2 α 2 α N ) .
(A3)
For N = 2 Bobs, this efficiency cutoff appears as the solid purple line in Fig. 3.

In Table II, we summarize the data appearing in Fig. 3 of the main text. Rows shaded in green and yellow indicate datasets that demonstrate detection-loophole-free steering of both Bobs.

Table II.

Data demonstrating detection-loophole-free steering. Green (yellow) backgrounds indicate steering by at least two (one) standard deviations. Experimental data with a red background do not show steering. εmean is the average value of ε exp across all measurements for a given α. Reported uncertainties correspond to one standard deviation of 200 Monte Carlo simulations adding Poisson noise to the raw data and repeating the MLE and steering test.

α ε exp ε mean ε Bob 1 ε Bob 2 ε ideal ( sim ) ε noisy ( sim )
0.015  0.7171(0.0100)  0.7276(0.0056)  0.7975(0.0993)  0.8023(0.0934)  0.5040  0.5978 
0.065  0.6899(0.0116)  0.7059(0.0062)  0.5858(0.0402)  0.5873(0.0513)  0.5171  0.5521 
0.100  0.7259(0.0108)  0.7390(0.0058)  0.5476(0.0333)  0.5394(0.0349)  0.5277  0.5561 
0.185  0.7282(0.0090)  0.7317(0.0049)  0.6021(0.0272)  0.6257(0.0258)  0.5553  0.5804 
0.300  0.7020(0.0086)  0.7160(0.0046)  0.6644(0.0182)  0.6591(0.0145)  0.6053  0.6303 
0.400  0.7206(0.0095)  0.7325(0.0051)  0.7393(0.0273)  0.7095(0.0186)  0.6646  0.6926 
0.500  0.6949(0.0099)  0.7079(0.0047)  0.7990(0.0227)  0.7780(0.0247)  0.7491  0.7807 
α ε exp ε mean ε Bob 1 ε Bob 2 ε ideal ( sim ) ε noisy ( sim )
0.015  0.7171(0.0100)  0.7276(0.0056)  0.7975(0.0993)  0.8023(0.0934)  0.5040  0.5978 
0.065  0.6899(0.0116)  0.7059(0.0062)  0.5858(0.0402)  0.5873(0.0513)  0.5171  0.5521 
0.100  0.7259(0.0108)  0.7390(0.0058)  0.5476(0.0333)  0.5394(0.0349)  0.5277  0.5561 
0.185  0.7282(0.0090)  0.7317(0.0049)  0.6021(0.0272)  0.6257(0.0258)  0.5553  0.5804 
0.300  0.7020(0.0086)  0.7160(0.0046)  0.6644(0.0182)  0.6591(0.0145)  0.6053  0.6303 
0.400  0.7206(0.0095)  0.7325(0.0051)  0.7393(0.0273)  0.7095(0.0186)  0.6646  0.6926 
0.500  0.6949(0.0099)  0.7079(0.0047)  0.7990(0.0227)  0.7780(0.0247)  0.7491  0.7807 
Here, we formulate our task as finding the most likely assemblages and Bobs' states to generate the given experimental data:
max { σ a | x } a , x , ρ B n a , b , x , y C ( b | y ) | a | x log Tr ( E b | y σ a | x ) , s . t . x , σ | x = ( 1 ε ) ρ B n , x , a σ a | x = ρ B n , Tr ρ B n = 1 , a , x , σ a | x 0 ,
(C1)
where C ( b | y ) | a | x is the experimental count of outcome b from POVM E b | y, and ε is calculated from x , y , b ; a C ( b | y ) | a | x / x , y , a , b C ( b | y ) | a | x before solving Problem (C1). Every quantity in (C1) except the indices is indexed by α and Bn, which we mostly omit to ease the burden of expression. We solve (C1) for each α and n independently, maximizing the logarithm of the likelihood of the unknowns generating the data in experiment. For fixed α and n, the maximum likelihood estimation (MLE) problem in (C1) is a conic optimization problem,43 which is more general than an SDP, and can be solved using standard software. Here, we employ the YALMIP toolbox44 and call the MOSEK solver45 under a MATLAB environment, and a selection of tomography results from experiment data is shown in Fig. 6.
Fig. 6.

The actual reconstructed assemblage tomography result versus the expected experimental assemblage for the largest steering demonstration α = 0.1. (a), (b), and (c) show the Alice–Bob 1 assemblage with the first group of POVM setting (see Table I) corresponding to Alice's three outcomes: +, −, and , respectively. Green (purple) bars correspond to the expected (reconstructed) assemblage. Imaginary components are not plotted as they are all <0.002.

Fig. 6.

The actual reconstructed assemblage tomography result versus the expected experimental assemblage for the largest steering demonstration α = 0.1. (a), (b), and (c) show the Alice–Bob 1 assemblage with the first group of POVM setting (see Table I) corresponding to Alice's three outcomes: +, −, and , respectively. Green (purple) bars correspond to the expected (reconstructed) assemblage. Imaginary components are not plotted as they are all <0.002.

Close modal

We then perform hypothesis testing to check the degree to which our tomography results are consistent with the experimental data (see  Appendix F). We obtain a test value of 572.0617 from our tomography results, smaller than the threshold of the critical region 595.1683 (corresponding to a significance level s α = 5 %), indicating an acceptance of the null hypothesis that the MLE result matches the true value of the assemblage.

The tomography for α = 0.015 was omitted from this testing as this data point has a large statistical uncertainty arising from diminishing photon counts as α decreases, which causes the hypothesis testing to fail.

The expected ensembles (Tables III and V) vs the MLE reconstructed ensemble (Tables IV and VI) tomography results for Alice–Bob 1 conditioned on Alice's first measurement. Here, these are shown for the two datasets α = 0.015, 0.1. Imaginary components are not displayed as they are all <0.002.

Table III.

Expected ensemble for α = 0.1.

Outcome | H H | | H V | | V H | | V V |
0.6017  0.0166  0.0166  0.0017 
−  0.1032  −0.0166  −0.0166  0.0183 
Null  0.2676  0.0076 
Outcome | H H | | H V | | V H | | V V |
0.6017  0.0166  0.0166  0.0017 
−  0.1032  −0.0166  −0.0166  0.0183 
Null  0.2676  0.0076 
Table IV.

MLE reconstruction for α = 0.1.

Outcome | H H | | H V | | V H | | V V |
0.6073  0.0121  0.0121  0.0027 
−  0.0911  −0.0370  −0.0370  0.0388 
Null  0.2454  −0.0088  −0.0088  0.0146 
Outcome | H H | | H V | | V H | | V V |
0.6073  0.0121  0.0121  0.0027 
−  0.0911  −0.0370  −0.0370  0.0388 
Null  0.2454  −0.0088  −0.0088  0.0146 
Table V.

Expected ensemble for α = 0.015.

Outcome | H H | | H V | | V H | | V V |
0.7104  −0.0313  −0.0313  0.0059 
−  0.0045  −0.0004  −0.0004  0.0078 
Null  0.2663  −0.0018  −0.0018  0.0051 
Outcome | H H | | H V | | V H | | V V |
0.7104  −0.0313  −0.0313  0.0059 
−  0.0045  −0.0004  −0.0004  0.0078 
Null  0.2663  −0.0018  −0.0018  0.0051 
Table VI.

MLE reconstruction for α = 0.015.

Outcome | H H | | H V | | V H | | V V |
0.7139  0.0025  0.0025  0.0026 
−  0.0055  0.0001  0.0001  0.0066 
Null  0.2680  0.0009  0.0009  0.0034 
Outcome | H H | | H V | | V H | | V V |
0.7139  0.0025  0.0025  0.0026 
−  0.0055  0.0001  0.0001  0.0066 
Null  0.2680  0.0009  0.0009  0.0034 
The reduced state between Alice and the nth Bob, ( n = 1 , , N) from Eq. (3) is
ρ α , N A B n : = Tr ¬ A B n [ ρ α , N ] = 1 N ( | ψ α ψ α | + ( N 1 ) ρ α | 0 0 | ) .
(E1)
Here, we show the bound in Eq. (A3) pertains to the assemblage that arises when Alice measures three projective measurements corresponding to the Pauli observables. By exchange symmetry between each of the trusted parties, each Bob is steered to ensembles within identical assemblages, and so we omit the labels Bn below. For ease of notation, below we refer to measurements corresponding to the Pauli observables Z : = | 0 0 | | 1 1 | , X : = | 0 1 | + | 1 0 | and Y : = i | 0 1 | + i | 1 0 | according to the ordered index x = 0, 1, 2, in that order. The assemblage detected by each of the N Bobs can, therefore, be defined as follows.
The ensembles conditioned for the first two measurements are defined by the positive operators
σ + 1 | 0 ( α , N ) = ( 1 α ) | 0 0 | ,
(E2)
σ 1 | 0 ( α , N ) = α 2 ( I + ( 1 2 N ) Z ) ,
(E3)
σ ± 1 | 1 ( α , N ) = 1 4 ( I ± 2 N α ( 1 α ) X + ( 1 2 α N ) Z ) ,
(E4)
and the third is given by a unitary rotation on the last of these, σ ± 1 | 2 = U z ( π / 2 ) σ ± | 1 U z ( π / 2 ), where U z ( θ ) denotes rotation about the z-axis through an angle θ. Note that
ρ B = 1 2 ( I + ( 1 2 α N ) Z ) .
(E5)
We begin by finding when each Bob's assemblage is steerable without consideration of Alice's reported non-detection events, in terms of the parameters α and N. The purpose is to determine when the assemblage demonstrates steering, so that appropriate ranges of these variables can be considered in the derivations to follow. For this purpose, we can examine the steering inequality derived in main text Ref. 46,
Tr [ X σ + 1 | 1 ] p + 1 | 1 1 2 ( p + 1 | 0 1 z + 1 | 0 2 + p 1 | 0 1 z 1 | 0 2 ) .
(E6)
The quantities appearing in this inequality can be calculated from the assemblage in Eqs. (E2)–(E4) directly via p a | x : = Tr [ σ a | x ] and z a | 0 : = Tr [ Z σ a | 0 ] / p a | 0. Therefore, we find Eq. (E6) is violated whenever
α ( 0 , 2 N + 1 ) .
(E7)
In this interval, I, we seek to derive an equation describing the behavior of the cutoff-efficiency.

To do this, we will find an analytic solution to the SDP in Eq. (A1) of  Appendix A. We proceed in Sec. 1 of  Appendix E by first simplifying this SDP into an equivalent form, by exploiting the symmetry present in the assemblages held by each Bob. Using the same symmetries, its dual program is also given. Then, we begin Sec. 2 of  Appendix E by conducting extensive numerical tests to guide an ansatz for the primal variables which achieve ε , as a function of α and N in I. We further discuss the equation in Sec. 3 of  Appendix E, providing tight lower bounds on it. In Sec. 4 of  Appendix E, we formulate an ansatz for the dual program. This section of the Appendix concludes by showing that these ansatzes achieve the same value in Sec. 5 of  Appendix E, and therefore are optimal.

1. Simplifications

Recall that the set of probability distributions { D ( a | x , λ ) } λ which encode the LHS constraints map x to all outcomes, both null and non-null. There are 3 3 = 27 such distributions, and as many operators { σ λ } λ over which the objective function must be maximized. To reduce the dimension of this search space, we make some observations.

The first simplification arises from the rotational symmetry present in the assemblage. Notice that the four assemblage members conditioned on settings 1 and 2 are all related by unitary rotations. For instance, consider the Lie group G generated by π Z / 4. Each group element (represented by k = 0 , , 3) corresponds to a unitary operator U k : = exp [ i π k Z / 4 ]. Let U k ( ) : = U k U k . Observe that the assemblage defined in Eqs. (E2)–(E4) is G-covariant, meaning that { σ a | x } a | x { U k ( σ a | x ) } a | x k. It is evident that the steered states conditioned on settings x = 1, 2 are related by
σ a | x = U k ( σ a | x ) ,
(E8)
for an appropriate rotation determined by k ( a , x , a , x ); the argument of this function will be omitted below for brevity. Moreover, the steered states conditioned on x = 0 are invariant under these rotations. This implies that the SDP constraint
λ D ( a | x , λ ) σ λ = ε σ a | x
(E9)
has interesting symmetry properties. Suppose there exists, for all x, and all a , an ensemble { σ λ } λ and ε for which this constraint is satisfied. From Eq. (E8), one can then always construct a G-covariant ensemble { U k ( σ λ ) / 4 } ( λ , k ), which can exactly reproduce the x = 0 steered states by the choice
ε σ a | 0 = k , λ D ( a | 0 , λ ) U k ( σ λ ) ,
(E10)
and the G-covariant steered states by
ε σ a | x = U k ( ε σ a | x )
(E11)
= U k ( λ D ( a | x , λ ) σ λ )
(E12)
= k , λ D ( f ( a , x , k ) | g ( x , k ) , λ ) U k ( σ λ ) ,
(E13)
for x = 1 , 2 and a = ± 1. The functions f and g satisfy this equation for the choice
a = f ( a , x , k ) : = { ( 1 ) h ( x , k ) a if a = ± 1 , x 0 , a else ,
(E14)
x = g ( x , k ) : = { ( x mod 2 ) + 1 , k odd , x 0 , x else .
(E15)
Here, we have defined the function h ( x , k ) : = ( ( k + x 1 ) mod 4 ) / 2 , which has binary outputs. Equations (E10) and (E13) show that all components of the assemblage (labeled by a and x ) can also be reproduced by the G-covariant ensemble { U k ( σ λ ) / 4 } ( λ , k ). This implies that, for the original problem, if we assume—without loss of generality—that the set { σ λ } λ is G-covariant, it is sufficient to confirm that, for any ε, the elements σ a | x for ( a , x ) t : = { ( + 1 , 0 ) , ( 1 , 0 ) , ( + 1 , 1 ) } are reproduced by it. The set of tuples t is minimal (but not unique) in this sense, in that it encodes a minimal set of constraints for the SDP, which are sufficient as a consequence of symmetry.
This symmetry can further be reflected in Alice's cheating strategies by partitioning them into equivalence classes (ECs). We say that two different strategies with labels λ , λ , are in the same equivalence class if
D ( a | x , λ ) D ( f ( a , x , k ) | g ( x , k ) , λ )
(E16)
for some k. Notice that these ECs have either K or 1 members. For brevity, we keep this multiplicity when referring to representatives of ECs with unit cardinality below, to simplify the writing of group averages over G. For our problem, the original 27 strategies can be partitioned into 9 ECs; these are characterized in Table VII. This means that, instead of labeling strategies by λ, they can be represented by the tuple (c, k), which uniquely specifies each strategy to be member k of the EC c. Moreover, instead of associating a different σ λ with each strategy, one can consider an operator σc for each EC, with an appropriate rotation indexed by k, σ λ U k ( σ c ), corresponding to each strategy (c, k). The optimization problem in Eq. (A1) can therefore be expressed in a form which encodes these symmetries:
max { σ c } c ε , s . t . c , k D ( f ( a , x , k ) | g ( x , k ) , c ) σ c , k = ε σ a | x , a , x t , c , k σ c , k = ρ B , σ c 0 , c .
(E17)
Table VII.

Canonical representatives σc for each of the nine equivalence classes of strategies, which comprise an LHS ensemble (if it exists). These are found by enumerating the original 3 3 = 27 strategies and grouping them according into ECs Eq. (E16). The final three columns indicate the outcomes announced for each of the three settings. The cardinality (number of strategies in each EC) is determined by having at least one non-null outcome announced for settings x = 1, 2.

EC representative Cardinality x = 0 x = 1 x = 2
σ0  +1  +1  +1 
σ1  +1  +1   
σ2  +1     
σ3  −1  +1  +1 
σ4  −1  +1   
σ5  −1     
σ6    +1  +1 
σ7    +1   
σ8       
EC representative Cardinality x = 0 x = 1 x = 2
σ0  +1  +1  +1 
σ1  +1  +1   
σ2  +1     
σ3  −1  +1  +1 
σ4  −1  +1   
σ5  −1     
σ6    +1  +1 
σ7    +1   
σ8       
The primal variables for the problem (E17) are the positive operators representative of each equivalence class, { σ c } c, and the real scalar ε. By introducing Hermitian operators F a | x , M , and Hc as dual variables corresponding to the three types of constraints, we can form the Lagrangian
L = ε + Tr [ a , x t F a | x ( c , k D ( f ( a , x , k ) | g ( x , k ) , c ) σ c , k ε σ a | x ) ] + Tr [ M ( ρ B c , k σ c , k ) ] + Tr [ c H c σ c ]
(E18)
= Tr [ M ρ B ] + ε ( 1 Tr a , x t F a | x ) + Tr [ c σ c ( k a , x t D ( f ( a , x , k ) | g ( x , k ) , c ) U k F a | x U k + H c k U k M U k ) ] .
(E19)
The dual program is then
min Tr M ρ B , s . t . a , x t Tr F a | x σ a | x = 1 , k a , x t D ( f ( a , x , k ) | g ( x , k ) , c ) U k F a | x U k + H c = M ¯ , c , H c 0 , c ,
(E20)
where M ¯ : = k U k M U k. Notice that all Hc matrices are slack variables. We explicitly include them in the dual formulation, because they will be important for deriving a certificate of optimality below.
2. Constructing the ansatz: Primal

The SDPs above can be straightforwardly solved, given fixed numeric values of α and N, using standard solvers. We will go one step further and derive a closed-form equation that explains the behavior of ε as a function of α and N, permitting us to take limits of the parameters analytically. To do this, we will be guided by numerical simulations to construct primal and dual sets of variables, which are both valid (feasible) and achieve the same values for the primal and dual objective functions, respectively. That is, they are solutions with zero duality gap and are thus optimal.

The main difficulty is that the primal problem takes into account all ECs of strategies for Alice, some of which are not used by her when the optimal value of ε is obtained. Here, and below, a denotes a primal or dual variable that is optimal. To construct an ansatz for primal variables, we first find the classes c for which we can expect σ c = 0. To this end, we examine all families of solutions to (E17) for which subsets of the operators σc are set to zero—or, equivalently, omitted from the problem formulation. To begin, we observe that σ8 represents the strategy for which null outcomes are always announced, so we set it to zero. For the remaining eight equivalence classes of strategies, there are 2 8 = 256 such combinations to consider. For each of these, we fix N = 2 and perform a sweep over 0 < α 2 / 3, solving the SDP for each value of α. These results are shown by the gray points in the left subfigure in Fig. 7. The solution to the original problem (E17) is reproduced exactly by the simulations for which σ 0 = σ 4 = σ 7 = 0, and Tr [ σ i ] > 0 for all others.

Fig. 7.

Left: Solutions ε for 28 families of SDP problems, defined by omitting all combinations of the matrices σ 0 , , σ 7 from the primal SDP formulation. For each such problem, we fix N = 2 and vary α to generate sets of gray points. The blue line corresponds to the solution for the sub-problem which encodes σ 1 = σ 4 = σ 7 = 0, which matches the solution giving the cutoff ε . Right: Cutoffs from Eq. (E39) for large numbers of Bobs, beyond the N = 2 curve (reproduced from Fig. 3 here in blue). This illustrates the remarkable observation that, for arbitrarily large N, there exists an interval of α such that ε is below one, facilitating loss-tolerant steering. Moreover, as shown in  Appendix E3, the smallest requirement on detector efficiency occurs in the limit α 0, where ε = ( 1 + 2 / N ) 1.

Fig. 7.

Left: Solutions ε for 28 families of SDP problems, defined by omitting all combinations of the matrices σ 0 , , σ 7 from the primal SDP formulation. For each such problem, we fix N = 2 and vary α to generate sets of gray points. The blue line corresponds to the solution for the sub-problem which encodes σ 1 = σ 4 = σ 7 = 0, which matches the solution giving the cutoff ε . Right: Cutoffs from Eq. (E39) for large numbers of Bobs, beyond the N = 2 curve (reproduced from Fig. 3 here in blue). This illustrates the remarkable observation that, for arbitrarily large N, there exists an interval of α such that ε is below one, facilitating loss-tolerant steering. Moreover, as shown in  Appendix E3, the smallest requirement on detector efficiency occurs in the limit α 0, where ε = ( 1 + 2 / N ) 1.

Close modal
Many of the constraints in the primal problem (E20) are redundant, being satisfied from the symmetry properties encoded in the problem. Defining group averages by U ¯ : = k U k, and removing the multiplicities on the ECs represented by σ 2 , σ 5 , and σ8, we can formulate the problem concisely:
max ε , s . t . U ¯ ( σ 0 ) + U ¯ ( σ 1 ) + σ 2 = ε σ + | 0 , U ¯ ( σ 3 ) + U ¯ ( σ 4 ) + σ 5 = ε σ | 0 , U ¯ ( σ 6 ) + U ¯ ( σ 7 ) + σ 8 = ( 1 ε ) ρ B , i = 0 , 3 , 6 σ i + U 1 ( σ i ) = ε σ + | 1 , c , σ c 0.
(E21)
Let σ c : = p c ( I + r ̂ c · σ ) / 2 be a rank-one decomposition of the matrix variables in this problem, for some Bloch vector r ̂ c : = ( sin θ c cos ϕ c , sin θ c sin ϕ c , cos θ c ). With this parametrization, the final matrix inequality constraint is guaranteed provided p c 0. It is evident from the left subfigure in Fig. 7 that the solution to the original SDP for the case of N = 2 Bobs occurs when we proceed by making the ansatz σ 1 = σ 4 = σ 7 = σ 8 = 0. Furthermore, we choose angles ϕ 3 = ϕ 6 = π / 4, and θ 5 = π. The latter is a natural choice because σ5 must be invariant under rotations generated by Z, and so must (when normalized) be located at the top or bottom of the Bloch sphere, meaning either θ 5 = 0 or θ 5 = π. We make the ansatz that this choice of primal variables will also be optimal for arbitrary N. The constraints of the problem (E21) become
4 p 0 + p 2 = ε ( 1 α ) ,
(E22)
4 p 3 + p 5 = ε α ,
(E23)
4 p 3 cos ( θ 3 ) p 5 = ( 1 2 N ) ε α ,
(E24)
4 p 6 = 1 ε ,
(E25)
4 p 6 cos ( θ 6 ) = ( 1 ε ) ( 1 2 α N ) ,
(E26)
2 p 3 sin ( θ 3 ) + 2 p 6 sin ( θ 6 ) = α ( 1 α ) ε N ,
(E27)
2 p 0 + 2 p 3 + 2 p 6 = ε 2 ,
(E28)
2 p 0 + 2 p 3 cos ( θ 3 ) + 2 p 6 cos ( θ 6 ) = ε 2 ( 1 2 α N ) .
(E29)
It is now straightforward to solve for the unknown variables by substitution. We have
p 6 = 1 ε 4 ,
(E30)
θ 6 = cos 1 ( 1 2 α N ) ,
(E31)
and by rearranging, p 5 = ε α 4 p 3 , p 2 = ε ( 1 α ) 4 p 0, and p 0 = 1 4 ( 2 ε 4 p 3 1 ). Upon substituting, only three constraints remain:
( N 1 ) ε α N = 2 p 3 ( 1 + cos ( θ 3 ) ) ,
(E32)
2 N p 3 sin ( θ 3 ) + ( 1 ε ) α ( N α ) = ε 2 α ( 1 α ) ,
(E33)
α ( 2 ε 1 ) N = 2 p 3 ( 1 cos ( θ 3 ) ) .
(E34)
From (E32), we find
p 3 = ( N 1 ) ε α 2 N ( 1 + cos ( θ 3 ) ) ,
(E35)
which is valid if θ 3 π, and so Eq. (E34) implies that
θ 3 = cos 1 ( N ε 3 ε + 1 N ε + ε 1 ) .
(E36)
The last remaining constraint is
α ( N 1 ) ε tan ( 1 2 cos 1 ( ( N 3 ) ε + 1 N ε + ε 1 ) ) + ( 1 ε ) α ( N α ) = ε 2 α ( 1 α ) .
(E37)
Using the identity tan ( 1 2 cos 1 ( x / y ) ) = ( y x ) / ( y + x ), this simplifies to
ε 2 α ( 1 α ) + ( ε 1 ) α ( N α ) = α 2 ε 1 ( N 1 ) ε ,
(E38)
which implicitly defines a solution for the remaining variable ε. This can be converted into a quadratic equation in ε, from which we can find a closed-form solution. This solution exactly matches the numerics and is given by ε ( α , N ) = e α , N, where
e α , N : = 2 N α ( N + 1 ) + 2 2 ( 1 α ) N α + α ( N 1 ) ( ( α + 4 ) N 5 α 4 2 ( 1 α ) N α ) 2 ( 2 2 ( 1 α ) N α + N + 2 α 2 α N ) .
(E39)
3. Tight lower bounds on e α , N

Here, we provide tight lower bounds on e α , N, for two reasons. The first is to justify the statement regarding the minimum efficiency being obtained in the singular limit α 0. Second, will require e α , N > 1 / 2 for the ansatz of the dual variables to be well-defined below.

Recall that we are interested in the interval I : = α ( 0 , 2 / ( N + 1 ) ). To obtain a lower bound on e α , N in I, which is tight for any integer N 2, we first observe that it is monotonically increasing. To this end, define
z 0 : = 2 N α ( N + 1 ) + 2 2 ( 1 α ) N α ,
(E40)
z 1 : = α ( N 1 ) ( ( α + 4 ) N 5 α 4 2 ( 1 α ) N α ) ,
(E41)
z 2 : = 2 ( 2 2 ( 1 α ) N α + N + 2 α 2 α N ) ,
(E42)
so that e α , N = ( z 0 + z 1 ) / z 2. We will show that both z 0 / z 2 and z 1 / z 2 are strictly increasing functions of alpha in I, for any integer N 2.
First, we observe that
α z 0 z 2 = ( N 1 ) ( 2 α ( N + 4 ) + 3 1 α N N α + 2 1 α N α + 2 ( 3 N + 2 ) ) 2 1 α N α ( α 2 α N + 2 2 2 α N α + N + 2 ) 2 .
(E43)
The denominator is positive in I. To see that the numerator is also positive, by omitting some strictly positive terms, we know that it is lower bounded by 2 ( 3 N + 2 ) 2 α ( N + 4 ) > 0. Hence, z 0 / z 2 is strictly increasing in I.
Now, for the second term, note that
z 2 α = 2 ( 2 2 α N α 2 N α 2 2 α 2 N 1 ) ,
(E44)
which implies that z 2 1 is strictly increasing (and positive). For the final term, note that
z 1 2 α = 2 ( N 1 ) [ ( α + 2 ) N + 2 α ( 3 4 α ) + 2 ( 3 α 2 ) N 5 1 α α N α 1 α N α ] .
(E45)
The term in square brackets is strictly positive for N 2 and 0 < α < 1, and which contains the interval of interest. This means that its square root, z1, is strictly positive and increasing too. We conclude that z 1 / z 2 must also be strictly increasing in I.
Now, the sum of any two strictly increasing functions is also strictly increasing. Therefore, e α , N = ( z 0 + z 1 ) / z 2 is strictly increasing in I. This implies that, for all values of N, min α I e α , N is achieved at the start of the interval. That is,
e α , N > lim α 0 + e α , N
(E46)
= 1 1 + 2 N .
(E47)
This also shows that e α , N > 1 / 2 for all valid N. Finally, we observe that this limit is singular, since for α = 0, the family entangled states in (3) of the main text are unentangled, therefore non-steerable, and so ε = 1.
4. Constructing the ansatz: Dual

Note that have not yet shown that this choice of primal variables is optimal. The set of primal variables achieving the objective value in Eq. (E39) above are only one feasible set of variables, i.e., they satisfy the constraints of the primal problem. However, we can prove optimality by analyzing the dual program and constructing a set of dual variables for which the dual objective function is equal to the primal objective function above. In other words, they will form a primal/dual pair with zero duality gap and are thus optimal.

For convenience, we summarize the primal variables forming the ansatz above:
σ 0 = α ( 1 e α , N ) + ( 2 α ) N e α , N N 4 N | 0 0 | ,
(E48)
σ 1 = 0 ,
(E49)
σ 2 = ( 1 e α , N ) ( N α ) N | 0 0 | ,
(E50)
σ 3 = α ( N 1 ) e α , N 2 N ( 1 + cos θ 3 ) · 1 2 ( I + r ̂ 3 · σ ) ,
(E51)
σ 4 = 0 ,
(E52)
σ 5 = α e α , N ( 1 + 2 ( 1 N ) N ( 1 + cos θ 3 ) ) | 1 1 | ,
(E53)
σ 6 = 1 e α , N 4 · 1 2 ( I + r ̂ 6 · σ ) ,
(E54)
σ 7 = 0 ,
(E55)
σ 8 = 0 ,
(E56)
where the Bloch vectors r ̂ i : = ( sin θ i cos ϕ i , sin θ i sin ϕ i , cos θ i ) are defined by
θ 3 = cos 1 ( ( N 3 ) e α , N + 1 ( N + 1 ) e α , N 1 ) ,
(E57)
ϕ 3 = π 4 ,
(E58)
θ 6 = cos 1 ( 1 2 α N ) ,
(E59)
ϕ 6 = π 4 .
(E60)
We first remove some unconstrained degrees of freedom in the dual program in (E20), by using the fact that the steered state σ + | z / Tr [ σ + | z ] is pure. Since pure states are extreme points in the space of two-qubit density matrices, this implies that the portion of the LHS ensemble which averages to it must also consist of the same pure state, σ i = p i | 0 0 | for i = 0, 1, 2. Using this fact, we can formulate (E20) into the equivalent dual program:
min Tr M ρ B ,
(E61)
s . t . Tr [ F + | 0 σ + | 0 + F | 1 σ | 1 + 4 F + | 1 σ + | 1 ] = 1 ,
(E62)
0 | M | 0 = 0 | ( F + | 0 + F + | 1 + U 1 ( F + | 1 ) ) | 0 + h 0 ,
(E63)
0 | M | 0 = 0 | ( F + | 0 + F + | 1 ) | 0 + h 1 ,
(E64)
0 | M | 0 = 0 | F + | 0 | 0 + h 2 ,
(E65)
M = F | 1 + F + | 1 + U 1 ( F + | 1 ) + H 3 ,
(E66)
M = F | 1 + F + | 1 + H 4 ,
(E67)
M = F | 1 + H 5 ,
(E68)
M = F + | 1 + U 1 ( F + | 1 ) + H 6 ,
(E69)
M = F + | 1 + H 7 ,
(E70)
M = H 8
(E71)
for c = 0 , 1 , 2 , h c 0
(E72)
for c = 3 , , 8 , H c 0.
(E73)
The variables in the dual program are the Hermitian matrices M , F + | 0 , F | 1 , F + | 1, three real scalars hi and the six positive-semidefinite matrices Hc. Motivated by numerical simulations, and examining the structure of the dual constraints, we formulate an ansatz now for the dual variables. We begin by assuming the first four of these matrices are of the form
M = ( μ 0 0 0 μ 1 ) ,
(E74)
F + | 0 = μ 0 | 0 0 | ,
(E75)
F | 1 = ( c 0 0 0 μ 1 ) ,
(E76)
F + | 1 = c 1 I + c 2 X c 1 Z
(E77)
for unknowns μ 0 , μ 1 , c 0 , c 1 , c 2. The structure of the Hc's will be guided by the so-called complementary slackness condition for optimality. These are necessary conditions that a primal/dual set of optimal variables must satisfy. Following Boyd and Vandenberghe43 (Sec. 5.5.2), for any primal and dual optimal set of variables that have zero duality gap, we know that
ε = Tr [ M ρ B ]
(E78)
ε + Tr [ a , x t F a | x ( ε σ a | x c , k D ( f ( a , x , k ) | g ( x , k ) , c ) σ c , k ) ] + Tr [ M ( ρ B c , k σ c , k ) ] + Tr [ λ H λ σ λ ]
(E79)
ε .
(E80)
The second line follows from the definition of L in (E18), and the third from the facts that primal feasibility holds for the primal variables, and Tr [ X Y ] 0 for any two positive semidefinite matrices X, Y. It can therefore be deduced that both inequalities in this chain hold with equality. This implies that the final term in (E79) vanishes,
λ Tr [ H λ σ λ ] = 0.
(E81)
Since each term in this sum is non-negative, all terms must vanish, and so for any primal/dual optimal set of variables, the ranges of H λ and σ λ must be pair-wise orthogonal. This implies we should express the dual variables H c in the form
H c : = h c ( I σ c Tr [ σ c ] ) ,
(E82)
when σ c is not zero, and where h λ are (undetermined) positive real scalars. In particular, from Eqs. (E51) and (E54), we define
H 3 = h 3 2 ( I sin θ 3 2 ( X Y ) + cos θ 3 Z ) ,
(E83)
H 6 = h 6 2 ( I sin θ 6 2 ( X Y ) + cos θ 6 Z ) .
(E84)
The other remaining Hc's are not important for deriving the dual optimal variables; only their existence as positive semidefinite matrices is important from the point of view of dual feasibility; we confirm this at the end.
The dual constraints which allow the unknowns to be determined are Eqs. (E66) and (E69):
F | 1 + F + | 1 + U 1 ( F + | 1 ) + H 3 = M ,
(E85)
F + | 1 + U 1 ( F + | 1 ) + H 6 = M .
(E86)
The latter of these, Eq. (E86), expressed in terms of the Pauli operators reads
( 2 ( 2 c 1 + h 6 ) ) I + ( 2 c 2 2 2 h 6 α ( N α ) N ) ( X Y ) ( 4 c 1 + 2 h 6 ( N 2 α ) N ) Z
(E87)
= ( μ 0 + μ 1 ) I + ( μ 0 μ 1 ) Z .
(E88)
Matching off-diagonal (in the Z-basis) terms, we find
h 6 = c 2 N 2 α ( N α ) ,
(E89)
from which matching the identity and Z terms implies
c 1 = 1 4 ( μ 0 + μ 1 2 c 2 N α ( N α ) ) ,
(E90)
c 2 = μ 0 α ( N α ) 2 α .
(E91)
Similarly, Eq. (E85) is
( c 0 + 4 c 1 + 2 h 3 + μ 1 ) I + ( 2 ( c 2 h 3 sin θ 3 2 ) ) ( X Y ) ( c 0 4 c 1 2 h 3 z 3 μ 1 ) Z
(E92)
= ( μ 0 + μ 1 ) I + ( μ 0 μ 1 ) Z .
(E93)
Once again matching terms, we find
h 3 = 2 c 2 sin θ 3 ,
(E94)
c 0 = μ 0 ( N 2 csc ( θ 3 ) α ( N α ) ) α μ 1 ,
(E95)
μ 1 = μ 0 ( N α + ( 1 + cos ( θ 3 ) ) csc ( θ 3 ) α ( N α ) ) α .
(E96)
Using the identity ( 1 + cos ( θ 3 ) ) csc ( θ 3 ) = 2 ( 1 + cos ( θ 3 ) ) / ( 1 cos ( θ 3 ) ), and the definition of θ3 in Eq. (E57), we can express the equation for μ1 in the form
μ 1 = μ 0 ( 2 ( N α ) 2 e α , N ( N 1 ) α ( N α ) 2 e α , N 1 ) 2 α .
(E97)
This only remaining degree of freedom is μ0, which we choose to satisfy the scalar constraint in Eq. (E62),
Tr [ F + | 0 σ + | 0 + F | 1 σ | 1 + 4 F + | 1 σ + | 1 ] = 1.
(E98)
For the choice of operators we have made, this reads
( 1 α ) μ 0 + 4 ( 1 2 c 1 ( 1 2 α N ) + ( 1 α ) α c 2 N + c 1 2 ) + α c 0 ( 1 1 N ) + α μ 1 N = 1 ,
(E99)
from which we deduce that
μ 0 = α e α , N N 2 e α , N ( α ( N α ) 2 α e α , N ( N 1 ) α ( N α ) 2 e α , N 1 + 2 ( 1 α ) α α ( N α ) ) + α e α , N ( N 1 ) α ( N α ) 2 e α , N 1 .
(E100)
This defines all relevant degrees of freedom in our construction for the dual variables. It remains to check (for dual feasibility) that the slack variables are valid for this choice; h c 0 for i = 0, 1, 2, and H c 0 for i = 3 , , 8. These matrices are defined through Eqs. (E63)–(E71). For the scalars, we have
h 0 = 0 | ( M F + | 0 F + | 1 U 1 ( F + | 1 ) ) | 0
(E101)
= 0 ,
(E102)
h 1 = 0 | ( M F + | 0 F + | 1 ) | 0
(E103)
= 0 ,
(E104)
h 2 = 0 | ( M F + | 0 ) | 0
(E105)
= 0.
(E106)
We now check the eigenvalues of the Hc matrices are all non-negative. The eigenvalues of
H 3 = h 3 2 ( I sin θ 3 2 ( X Y ) + cos θ 3 Z )
(E107)
are 0 and h3. From ((E94)) and ((E91)), we know that
h 3 = μ 0 α ( N α ) α sin θ 3 > 0 ,
(E108)
where the inequality follows from observing that μ 0 > 0, multiplies a term that is strictly positive in the interval I. The least eigenvalue of H 4 = M F | 1 F + | 1 is
1 2 ( μ 0 c 0 2 c 1 ( c 0 + 2 c 1 μ 0 ) 2 4 ( 2 c 0 c 1 2 c 1 μ 0 c 2 2 ) ) .
(E109)
This evaluates to zero; to see this, observe that
2 c 0 c 1 2 c 1 μ 0 c 2 2 =
(E110)
μ 0 2 ( N α ) [ α ( 1 2 e α , N ) e α , N ( N 1 ) + 2 e α , N 1 ( α ( N α ) e α , N ( 2 ( 1 α ) α ) + α ( N α ) ) ) ] 2 α 2 ( 2 e α , N 1 ) e α , N ( N 1 ) ,
(E111)
which evaluates to zero. This is because e α , N is a solution for ε in Eq. (E38), which upon rearranging and substituting causes the term in the square brackets to vanish. Now, H 5 = M F | 1 = diag ( μ 0 c 0 , 0 ), which has non-negative eigenvalues since
μ 0 c 0 = 2 e α , N 1 μ 0 α ( N α ) α e α , N ( N 1 ) ,
(E112)
which is strictly positive in I. Similarly, the eigenvalues of
H 6 = h 6 2 ( I sin θ 6 2 ( X Y ) + cos θ 6 Z )
(E113)
are 0 and h6, the latter of which is positive because h 6 = μ 0 N / ( 2 α ) > 0. Both eigenvalues of H7 are positive. To see this, the least eigenvalue of H 7 = M F + | 1 is
1 2 ( ( 2 c 1 μ 0 μ 1 ) 2 4 ( 2 c 1 μ 0 c 2 2 + μ 0 μ 1 ) 2 c 1 + μ 0 + μ 1 ) .
(E114)
To check that this eigenvalue is positive, it suffices to verify that
( 2 c 1 + μ 0 + μ 1 ) 2 > ( 2 c 1 μ 0 μ 1 ) 2 4 ( 2 c 1 μ 0 c 2 2 + μ 0 μ 1 ) ,
(E115)
which is equivalent to
μ 0 ( μ 1 2 c 1 ) > c 2 2 .
(E116)
Upon substitution, this becomes
α 2 e α , N ( 2 e α , N 1 ) μ 0 2 ( α 2 α N + 2 ( ( α 1 ) α ) α ( N α ) ) < 0 ,
(E117)
which is valid in I because e α , N > 1 / 2 , μ 0 > 0. Finally, since the eigenvalues of M = H8 are both strictly positive, we conclude that H8 is also strictly positive. Hence, all dual variables for the ansatz we have made satisfy the constraints of the dual program.
5. Certifying optimality
It remains to see that this choice implies a zero duality gap for the values achieved by the primal and dual objective functions. The duality gap for our ansatz is given by
Tr [ M ρ B ] e α , N = μ 0 ( 1 α N ) + α μ 1 N e α , N .
(E118)
Since we know that e α , N defines a solution for ε in Eq. (E38), we know, by rearranging, that it also satisfies
α e α , N ( N 1 ) 2 e α , N 1 = 2 ( 1 α ) α e α , N + ( e α , N 1 ) α ( N α ) 2 e α , N 1 .
(E119)
We can then simplify the eigenvalues of M into the forms
μ 1 = e α , N μ 0 ( α N + 2 1 α ( N α ) ) α ( 2 e α , N 1 )
(E120)
and
μ 0 = α e α , N ( 2 e α , N 1 ) N α 2 ( 1 3 e α , N ) + α ( 3 e α , N 1 ) N 2 e α , N ( 1 α ) α α ( N α ) .
(E121)
Upon substituting into Eq. (E118), it is straightforward to verify that
Tr [ M ρ B ] e α , N = 0 ,
(E122)
for the sets of primal and dual variables, we have defined above, implying that they form an optimal primal and dual pair. That is, we know that Eq. (E39) is a closed-form solution to the optimization problem and hence represents the cutoff efficiencies for each Bob's assemblage.

Denote θ the vector of all the unknown parameters to be estimated, and D the experiment data obtained. Suppose there are altogether R different measurement settings, accounting for different POVM groups and different initial states. The rth POVM group has a total measurement shots Sr for Or different outcomes, with the occurrence statistics D r : = ( D r 1 , D r 2 , , D r O r ). Hence, S r = i = 1 O r D r i.

Given D, one can use certain algorithm to obtain an estimation θ ̂. Then, the natural question is, to what degree is θ ̂ consistent with D? This is usually answered by performing hypothesis testing in statistics. More specifically, in our measurement-based scenario, we need to do a multinomial test47 as follows.

We would like to test the null hypothesis H 0 : θ = θ ̂. For the rth measurement group, this hypothesis will give a prediction of the normalized measurement probabilities f r = ( f r 1 , f r 2 , , f r O r ), where i = 1 O r f r i = 1. Then, the probability (likelihood) of observing D r under the null hypothesis is clearly
p ( D r | H 0 ) = S r ! i = 1 O r f r i D r i D r i ! .
(F1)
One then typically define an alternative hypothesis H1 where the predicted measurement probabilities, instead of being f r, should be ( p r 1 MLE , p r 2 MLE , , p r O r MLE ) where p r i MLE = D r i / S r is the maximum likelihood estimate from data. The probability of observing D r under the alternative hypothesis is then
p ( D r | H 1 ) = S r ! i = 1 O r ( p r i MLE ) D r i D r i ! .
(F2)

Now a likelihood ratio test can be done as LRT r : = p ( D r | H 0 ) / p ( D r | H 1 ). Assume H0 is true, then asymptotically ( S r ) the distribution of 2 ln ( LRT r ) = 2 i = 1 O r D r i ln ( f r i / p r i MLE ) converges to the χ 2 distribution with O r 1 degrees of freedom. Since measurements from different measurement settings are independent, we can employ the additivity of χ 2 distribution to assert that r = 1 R 2 ln ( LRT r ) = 2 r = 1 R i = 1 O r D r i ln ( f r i / p r i MLE ) should asymptotically follow the χ 2 distribution with r = 1 R ( O r 1 ) degrees of freedom, when H0 holds. One can thus select a significance level s α to decide whether H0 will be rejected, based on the calculated value of r = 1 R 2 ln ( LRT r ).

For our assemblage tomography part, r = 1 R ( O r 1 ) = r = 1 108 ( 6 1 ) = 540 and its significance level s α = 5 % corresponds to the critical region r = 1 R 2 ln ( LRT r ) 595.1683, while our tomography result gives r = 1 R 2 ln ( LRT r ) = 572.0617. Hence, the null hypothesis is accepted, and our tomography result matches the data.

The experiment, Fig. 2, uses two sources of single photon pairs based on the design of Ref. 41. One source provides a tunable entangled two-qubit state | Ψ α and the other the single qubit state | 0 (with the fourth photon detected to herald the photon that encodes the single qubit). We use a 775 nm wavelength Ti:sapph pump laser with 1 ps pulse length and 80 MHz repetition rate. The pump is split on a 50:50 beam splitter. The resulting beams are used to pump two 15 mm long periodically poled potassium titanyl phosphate crystals with approximately 100 mW pump power per crystal, each resulting in degenerate type-II spontaneous parametric down-conversion at telecom wavelength.

In each protocol run, the entangled state | Ψ α is shared between Alice and one other trusted party, and the third party receives the | 0 state. For each quantum state with a different value of α, approximately 7.5 × 10 4 runs are performed—a run being an event where both trusted parties and the heralding detector detect a photon. Our average data collection rate was 231 trials per second. Of those runs, in approximately half ( 0.53 ± 0.05 ) of the cases, Bob 1 receives part of the entangled state, and in the other cases, Bob 2 does. The swap operation, which determines whether Bob 1 or Bob 2 receives half of the entangled state, is implemented by a linear translation stage that moves mirrors into the beam paths, redirecting their spatial modes.

Each of the three parties can perform projective measurements using automated quarter-wave and half-wave plates, a polarizing beam splitter, and two superconducting nanowire single-photon detectors (SNSPD). Alice's detectors, the only detectors whose efficiency directly impacts the protocol's success, have an efficiency of 90 %.

The data for the experiment are recorded in 720-s batches. A batch consists of the linear translation stage remaining in one position while Alice and the Bobs sequentially perform their measurements in a predetermined order. We then average the detection efficiency of our SNSPDs by repeating these measurements with the wave plates rotated to swap which detector receives which outcome of the POVMs. Alice's measurement settings come from an optimization routine (see  Appendix A), and the Bobs perform tomographic measurements before moving on to the next combination of measurement settings. After each combination, the linear translation stage shifts (or does not shift) and measurement iterations begin again. The data files for each batch are summed to obtain the mixture, then organized into outcome groups between Alice and each of the Bobs, and we extract the heralding efficiencies ( ε a , x). We reconstruct the assemblages through conic optimization. From these assemblages, we use a SDP Eq. (A1) to calculate the cutoff efficiency for the trusted parties, ε Bob 1 and ε Bob 2 .

1.
S.
Slussarenko
and
G. J.
Pryde
,
Appl. Phys. Rev.
6
,
041303
(
2019
).
2.
M.
Nielsen
and
I.
Chuang
,
Quantum Computation and Quantum Information
(
Cambridge University Press
,
2000
).
3.
R.
Uola
,
A. C. S.
Costa
,
H. C.
Nguyen
, and
O.
Gühne
,
Rev. Mod. Phys.
92
,
015001
(
2020
).
4.
D.
Cuomo
,
M.
Caleffi
, and
A. S.
Cacciapuoti
,
IET Quantum Commun.
1
,
3
(
2020
).
5.
C.
Branciard
,
E. G.
Cavalcanti
,
S. P.
Walborn
,
V.
Scarani
, and
H. M.
Wiseman
,
Phys. Rev. A
85
,
010301
(
2012
).
6.
Y.
Xiang
,
I.
Kogias
,
G.
Adesso
, and
Q.
He
,
Phys. Rev. A
95
,
010101
(
2017
).
7.
D.
Mondal
,
C.
Datta
,
J.
Singh
, and
D.
Kaszlikowski
,
Phys. Rev. A
99
,
012312
(
2019
).
8.
Y. Z.
Law
,
L. P.
Thinh
,
J.-D.
Bancal
, and
V.
Scarani
,
J. Phys. A
47
,
424028
(
2014
).
9.
A.
Máttar
,
P.
Skrzypczyk
,
G. H.
Aguilar
,
R. V.
Nery
,
P. H.
Souto Ribeiro
,
S. P.
Walborn
, and
D.
Cavalcanti
,
Quantum Sci. Technol.
2
,
015011
(
2017
).
10.
D. J.
Joch
,
S.
Slussarenko
,
Y.
Wang
,
A.
Pepper
,
S.
Xie
,
B.-B.
Xu
,
I. R.
Berkman
,
S.
Rogge
, and
G. J.
Pryde
,
Phys. Rev. A
106
,
L050401
(
2022
).
11.
C. H.
Bennett
,
G.
Brassard
,
C.
Crépeau
,
R.
Jozsa
,
A.
Peres
, and
W. K.
Wootters
,
Phys. Rev. Lett.
70
,
1895
(
1993
).
12.
S. L. N.
Hermans
,
M.
Pompili
,
H. K. C.
Beukers
,
S.
Baier
,
J.
Borregaard
, and
R.
Hanson
,
Nature
605
,
663
(
2022
).
13.
D.
Gottesman
,
T.
Jennewein
, and
S.
Croke
,
Phys. Rev. Lett.
109
,
070503
(
2012
).
14.
B.
Hensen
,
H.
Bernien
,
A. E.
Dréau
,
A.
Reiserer
,
N.
Kalb
,
M. S.
Blok
,
J.
Ruitenberg
,
R. F. L.
Vermeulen
,
R. N.
Schouten
et al,
Nature
526
,
682
(
2015
).
15.
M.
Giustina
,
M. A. M.
Versteegh
,
S.
Wengerowsky
,
J.
Handsteiner
,
A.
Hochrainer
,
K.
Phelan
,
F.
Steinlechner
,
J.
Kofler
,
J.-Å.
Larsson
et al,
Phys. Rev. Lett.
115
,
250401
(
2015
).
16.
L. K.
Shalm
,
E.
Meyer-Scott
,
B. G.
Christensen
,
P.
Bierhorst
,
M. A.
Wayne
,
M. J.
Stevens
,
T.
Gerrits
,
S.
Glancy
,
D. R.
Hamel
et al,
Phys. Rev. Lett.
115
,
250402
(
2015
).
17.
H. M.
Wiseman
,
S. J.
Jones
, and
A. C.
Doherty
,
Phys. Rev. Lett.
98
,
140402
(
2007
).
18.
S. J.
Jones
,
H. M.
Wiseman
, and
A. C.
Doherty
,
Phys. Rev. A
76
,
052116
(
2007
).
19.
D.
Cavalcanti
and
P.
Skrzypczyk
,
Rep. Prog. Phys.
80
,
024001
(
2017
).
20.
A. J.
Bennet
,
D. A.
Evans
,
D. J.
Saunders
,
C.
Branciard
,
E. G.
Cavalcanti
,
H. M.
Wiseman
, and
G. J.
Pryde
,
Phys. Rev. X
2
,
031003
(
2012
).
21.
S.
Kocsis
,
M. J. W.
Hall
,
A. J.
Bennet
,
D. J.
Saunders
, and
G. J.
Pryde
,
Nat. Commun.
6
,
5886
(
2015
).
22.
M. M.
Weston
,
S.
Slussarenko
,
H. M.
Chrzanowski
,
S.
Wollmann
,
L. K.
Shalm
,
V. B.
Verma
,
M. S.
Allman
,
S. W.
Nam
, and
G. J.
Pryde
,
Sci. Adv.
4
,
e1701230
(
2018
).
23.
J. S.
Bell
,
Phys. Phys. Fiz.
1
,
195
(
1964
).
24.
B. D. M.
Jones
,
I.
Šupić
,
R.
Uola
,
N.
Brunner
, and
P.
Skrzypczyk
,
Phys. Rev. Lett.
127
,
170405
(
2021
).
25.
Q. Y.
He
and
M. D.
Reid
,
Phys. Rev. Lett.
111
,
250403
(
2013
).
26.
S.
Armstrong
,
M.
Wang
,
R. Y.
Teh
,
Q.
Gong
,
Q.
He
,
J.
Janousek
,
H.-A.
Bachor
,
M. D.
Reid
, and
P. K.
Lam
,
Nat. Phys.
11
,
167
(
2015
).
27.
I.
Kogias
,
Y.
Xiang
,
Q.
He
, and
G.
Adesso
,
Phys. Rev. A
95
,
012315
(
2017
).
28.
Y.
Cai
,
Y.
Xiang
,
Y.
Liu
,
Q.
He
, and
N.
Treps
,
Phys. Rev. Res.
2
,
032046
(
2020
).
29.
S. K.
Joshi
,
D.
Aktas
,
S.
Wengerowsky
,
M.
Lončarić
,
S. P.
Neumann
,
B.
Liu
,
T.
Scheidl
,
G. C.
Lorenzo
,
Ž.
Samec
et al,
Sci. Adv.
6
,
eaba0959
(
2020
).
30.
O.
Gühne
,
G.
Tóth
, and
H. J.
Briegel
,
New J. Phys.
7
,
229
(
2005
).
31.
D.
Cavalcanti
,
P.
Skrzypczyk
,
G. H.
Aguilar
,
R. V.
Nery
,
P. H.
Souto Ribeiro
, and
S. P.
Walborn
,
Nat. Commun.
6
,
7941
(
2015
).
32.
C.-M.
Li
,
K.
Chen
,
Y.-N.
Chen
,
Q.
Zhang
,
Y.-A.
Chen
, and
J.-W.
Pan
,
Phys. Rev. Lett.
115
,
010402
(
2015
).
33.
W.
McCutcheon
,
A.
Pappa
,
B. A.
Bell
,
A.
McMillan
,
A.
Chailloux
,
T.
Lawson
,
M.
Mafu
,
D.
Markham
,
E.
Diamanti
et al,
Nat. Commun.
7
,
13251
(
2016
).
34.
H.
Lu
,
C.-Y.
Huang
,
Z.-D.
Li
,
X.-F.
Yin
,
R.
Zhang
,
T.-L.
Liao
,
Y.-A.
Chen
,
C.-M.
Li
, and
J.-W.
Pan
,
Phys. Rev. Lett.
124
,
180503
(
2020
).
35.
F. J.
Curchod
,
M.
Johansson
,
R.
Augusiak
,
M. J.
Hoban
,
P.
Wittek
, and
A.
Acín
,
Phys. Rev. A
95
,
020102
(
2017
).
36.
M.
Choi
and
S.
Lee
,
Quantum Inf. Process.
20
,
47
(
2021
).
37.
J.
Zhu
,
M.-J.
Hu
,
C.-F.
Li
,
G.-C.
Guo
, and
Y.-S.
Zhang
,
Phys. Rev. A
105
,
032211
(
2022
).
38.
Q.-C.
Song
,
T. J.
Baker
, and
H. M.
Wiseman
,
Phys. Rev. A
108
,
012216
(
2023
).
39.
J. F.
Clauser
and
M. A.
Horne
,
Phys. Rev. D
10
,
526
(
1974
).
40.
G.
Vallone
,
Phys. Rev. A
87
,
020101
(
2013
).
41.
N.
Tischler
,
F.
Ghafari
,
T. J.
Baker
,
S.
Slussarenko
,
R. B.
Patel
,
M. M.
Weston
,
S.
Wollmann
,
L. K.
Shalm
,
V. B.
Verma
et al,
Phys. Rev. Lett.
121
,
100401
(
2018
).
42.
H. C.
Nguyen
,
H.-V.
Nguyen
, and
O.
Gühne
,
Phys. Rev. Lett.
122
,
240401
(
2019
).
43.
S. P.
Boyd
and
L.
Vandenberghe
,
Convex Optimization
(
Cambridge University Press
,
2004
).
44.
J.
Löfberg
, “
YALMIP: A toolbox for modeling and optimization in MATLAB
,” in
IEEE International Conference on Robotics and Automation
(
IEEE
,
2004
).
45.
MOSEK ApS
,
The MOSEK optimization toolbox for MATLAB manual. Version 10.0.46
(
MOSEK ApS
,
2022
).
46.
S. J.
Jones
and
H. M.
Wiseman
,
Phys. Rev. A
84
,
012110
(
2011
).
47.
Y. M. M.
Bishop
,
S. E.
Fienberg
, and
P. W.
Holland
,
Discrete Multivariate Analysis: Theory and Practice
(
Springer Science & Business Media
,
2007
).