Experimental robust self-testing of the state generated by a quantum network

Self-testing is a method of quantum state and measurement estimation that does not rely on assumptions about the inner working of the used devices. Its experimental realization has been limited to sources producing single quantum states so far. In this work, we experimentally implement two significant building blocks of a quantum network involving two independent sources, i.e. a parallel configuration in which two parties share two copies of a state, and a tripartite configuration where a central node shares two independent states with peripheral nodes. Then, by extending previous self-testing techniques we provide device-independent lower bounds on the fidelity between the generated states and an ideal state made by the tensor product of two maximally entangled two-qubit states. Given its scalability and versatility, this technique can find application in the certification of larger networks of different topologies, for quantum communication and cryptography tasks and randomness generation protocols.


INTRODUCTION
Within the last few years, a large number of quantum resource-based protocols have been designed, with a wide range of applications. However, it is crucial, and far from trivial, to discriminate the devices that are working correctly from those that are not. Indeed, two difficulties can emerge: on one hand, the task required by the user may be hard to verify (a notorious example being the Boson Sampling problem [1][2][3][4][5]) and, on the other, the devices may be affected by noise and imperfections that are unknown to the user. The latter case is especially relevant for tasks aimed to be secure against eventual adversaries, which could exploit such defects to obtain secret information or sabotage the operation of the devices. For instance, this is the case of private randomness generation or amplification and quantum key distribution protocols [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. Hence, the ability to certify that the device is operating properly, and possibly without relying on the knowledge of its internal working, is crucial for a wider application of quantum technologies. The approach where conclusions about the correctness of the device's operation are drawn only from input/output statistics, is known as device-independent (DI) [22] and typically relies on the quantum violation of Bell-like inequalities [23].
A key protocol in the device-independent scenario is that of self-testing [24]. There, a multipartite quantum state is subject to a number local measurements (a procedure called a Bell test), and the resulting statistics alone are enough to certify the specific form of the state and measurements. For instance, obtaining the maximum value of 2 √ 2 in a Clauser-Horne-Shimony-Holt (CHSH) Bell test [25] certifies that the state is equivalent to a two-qubit maximally entangled state. In recent years several self testing protocols have been proposed to certify different states and measurements [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41]. In this work we present experimental demonstrations of selftesting for two types of quantum network, each featuring two independent sources (see Fig. 1): (i) a network in which the sources are placed in a parallel configuration between two par-ties, and (ii) a network featuring three parties, where a central party shares a source with two peripheral parties. The design of our experimental setup follows the bipartite self-testing strategies recently proposed in [42], which we further adapt to the multipartite network (ii). To experimentally implement the network structures, we employ a flexible and versatile platform, introduced in [43], which allows one to easily change the quantum network topology. Precise lower bounds on the self-testing fidelity with the desired states are obtained from the experimental statistics via the SWAP method [31,41], a numerical tool based on semidefinite programming. Our results show that under realistic experimental conditions we can obtain non-trivial device-independent lower bounds on the fidelity between the actual state and ideal states. Moreover, the present techniques can in principle be extended to an arbitrary number of nodes and to an arbitrary target state, which makes them a promising tool for the certification of larger networks and in the implementation of quantum communication and cryptography tasks.

SELF-TESTING OF QUANTUM NETWORKS
In the device-independent scenario the measurement devices and sources are treated as black boxes, exchanging only classical communication with external users. Suppose the users (which we label as A, B, C, . . .) share some state ρ ABC··· that is unknown to them, and that they can prepare and measure the state in an independent, identically distributed (i.i.d.) manner. After the experiment is repeated many times the users can estimate the probabilities p(a, b, c, . . . |x, y, z, . . .) of obtaining measurement outcomes a, b, c, . . . , if measurements x, y, z, . . . are performed. According to quantum mechanics, such probabilities are given, through the Born rule as p(a, b, c, . . . |x, y, z, . . . Self-testing scenarios. The self-testing procedure consists in performing an experiment and analyse the produced data without assuming a particular implementation (a,b), i.e. by considering a black-box scenario in which we only have access to the conditional probability distributions of measurement results (labeled by a, b, and c) conditioned to measurement choices (labeled by x, y and z). Nothing is assumed on the shared quantum state and measurements. Then self-testing techniques are used to obtain the minimum fidelity between the real states produced in the experiment and the ideal situation shown in the right panel (c,d), where the sources produce perfect maximally entangled two-qubit states (e.g. |ψ − ) and each party applies the local Pauli measurements to each qubit, corresponding to a maximal violation of the CHSH inequality, e.g A0 = σz, A1 = σx, Here we perform self-testing analysis of the state produced in two geometries: c) a bipartite situation where we aim at certifying the presence of two copies of a maximally entangled state and d) a tripartite scenario in which a central party shares maximally entangled states with two peripheral parties.
where A x a , B y b , C z c . . . denote the local measurement operators. We say that the probabilities p(a, b, c, . . . |x, y, z, . . .) self-test the target state |ψ ABC··· if the observation of p(a, b, c, . . . |x, y, z, . . .) necessarily implies the existence of a local quantum channel Self-testing therefore certifies that the parties share the state |ψ ABC··· , in the sense that there exist local operations the parties could perform to extract the state from ρ. Note that this statement holds for any state ρ satisfying (1) (for some local measurements) and is thus a device-independent statement. As an example, it is known that any bipartite state ρ AB producing correlations that give the maximal quantum violation of the CHSH Bell inequality (with value 2 √ 2), there exists a channel such that ( In realistic scenarios, however, it is impossible to exactly meet the self-testing conditions (1), not only due to experimental noise, but also because the finite experiment time im-plies that one can only infer the probabilities up to a given confidence interval. For this reason, the self-testing statement has to be robust, i.e. to be able to say something about the underlying state even when the self-testing condition is only approximately met. To do this, we focus on lower bounding the fidelity between the extracted state and the target state given the experimental statistics. That is, one proves that for any state producing the experimental statistics there exists a local channel Ω such that F ≥ f (for some f ≤ 1), up to a given confidence level. Note f = 1 corresponds to the case of perfect self-testing given in (2). A useful method that we use for calculating such lower bounds is the SWAP method [31, 41], a numerical tool based on semidefinite programming and the Navascués-Pironio-Acín (NPA) hierarchy. In particular, this technique amounts to numerically swapping part of the generated state on a dummy register, in order to find a proper expression for the fidelity (3), as a function of the correlations terms p(a, b, c . . . |x, y, z . . . ). Then, a lower bound on this fidelity can be obtained through an SDP, over a superset of the quantum correlations set, mathematically defined by a level l of the NPA hierarchy. In order to get tighter bounds, further linear constraints can be added to the problem, e.g. the observed correlations p(a, b, c . . . |x, y, z . . . ). More details about this method can be found in the Supplementary Information.
In this work we report on the self-test of target states that correspond to two independent sources producing maximally entangled singlet states |ψ − . We focus on two scenarios, depicted in Fig. 1, which can be seen as two possible building blocks of a quantum network. The first scenario we consider features two parties (A and B), and the target state we self-test corresponds to preparing the two maximally entangled states in parallel (see Fig. 1c). That is, we self-test the state |Ψ 2 = |ψ − A1B1 ⊗ |ψ − A2B2 , where A i , B j denote local qubit Hilbert spaces of party A, B. The second structure we consider features three parties (A, B, C), and the target state corresponds to preparing the sources in an entanglement swapping network (see Fig. 1d). This configuration, although it has been investigated, both theoretically and experimentally, within the last years [44][45][46][47][48], was only very recently implemented exploiting truly independent sources [43,49] and closing the locality loophole [49]. Our target state in this scenario is therefore |Ψ 3 = |ψ − A1,B ⊗ |ψ − A2,C . We stress that, although the target states correspond to states produced by two independent sources, we do not assume this structure at the black-box level, that is, the parties can in principle share any multipartite state, but the self-testing statements ensure they have the desired product form.
Our self-testing protocol is inspired from the recent work [42], which presents a method for self-testing tensor products of copies of a state, while keeping the number of inputs constant. Such a method is desirable for self-testing quantum networks, since standard methods feature a number of inputs that grows exponentially with the number of copies, which will become a practical issue for larger networks. We consider scenarios in which all parties have two inputs; x, y = 0, 1 for scenario (a) and x, y, z = 0, 1 for scenario (b). The number of outputs is given by the local Hilbert space dimension of the target state: in scenario (a) both parties have four outputs, which we write as a = (a 1 , a 2 ) and b = (b 1 , b 2 ), where a i , b i take values ±1; in scenario (b) we have a = (a 1 , a 2 ) as before and b = ±1, c = ±1. The measurements are chosen so that in the ideal experiment, the marginal distributions maximally violate the CHSH Bell inequality. More precisely one has CHSH(p(a i , b i |x, y)) = 2 √ 2, i = 1, 2 for scenario (a) and CHSH(p(a 1 , b|x, y)) = CHSH(p(a 2 , c|x, z)) = 2 √ 2 for scenario (b). Following results from [42] such distributions are known to self-test the desired target states. The measurement strategy corresponds to the performing the standard CHSH measurements in parallel, which we elaborate on in the following section. More details about merging the SWAP method with the techniques from [42] can be found in the Supplementary Information.

EXPERIMENTAL APPARATUS
In the experimental implementation of the two scenarios of interest, we resort to the versatile photonic platform introduced in [43]. In particular, for the parallel self-testing case (Fig. 1c) we employ two separate laboratories equipped with independent quantum state sources (respectively, Λ 1 and Λ 2 ) and two measurement stations. Each measurement station is composed by a half-wave plate (HWP) and a polarizing beam splitter (PBS), which allows to perform polarization projective measurements of the form cos(4θ)σ z + sin(4θ)σ x , where σ x and σ z are Pauli operators, by simply rotating the HWP of the angle θ with respect to its optical axis. In the end, all the registered counts are sent to a central counter, which recognizes as coincidence events of distant detectors, the counts within a given time window. In our notation, the measurement stations in laboratory 1 represent Alice, while those in laboratory 2, Bob. Then, the two laboratories are connected through two ∼ 30 m long single-mode fibers, as shown in Fig. 2a. The source in laboratory 1 employs spontaneous parametric down-conversion (SPDC) of type II to generate a pair of polarization entangled photons of wavelength equal to λ = 785 nm, through a Beta-Barium Borate (BBO) crystal which is pumped, in pulsed regime, by a λ = 392.5nm U V laser beam. Instead, in laboratory 2, we have a periodically polarized Tytanil Phosphate, pumped in the continuouswave regime, which generates polarization entangled photon pairs at λ = 808 nm. Both sources generate a 2-qubit maximally entangled state, e.g. the singlet state |ψ − , where the computational basis (|0 and |1 ) is encoded in the horizontal and vertical photon polarization states (|H and |V ), hence . At this point, both laboratories send a photon to the other one and use one measurement station to perform projective measurements on the photon they kept and the second one to measure the photon they received. In de-tail, Alice's and Bob's operators will be those maximising the CHSH inequality violation, i.e., up to unitary transformations, , indicating with the superscript the source generating the subsystem.
In order to detect coincidence events between distant detectors, we designed a software to coordinate the counters located in the different laboratories. In particular, the software considered two coincidence time windows: one to detect the 2-fold coincidences generated by each source (set to 1.05 ns), w 1 , and another to reveal 4-fold coincidences, i.e. simultaneous 2-fold events occurring for both sources, w 2 . In other words, if a 2-fold event is registered for both sources, within w 2 , those are labelled as simultaneous and considered as a 4fold event. Then, the optimal value for w 2 is chosen through the corresponding CHSH values brought by the two sources, as depicted in Fig. 2c, and choosing the one giving the highest weighted average between the two values.
In the tripartite scenario ( Fig. 1d) we have 3 laboratories. Laboratories 1 and 3 (the peripherical nodes) are constituted by a quantum state source and one measurement station each, while laboratory 2 (the central node) has two measurement stations, as shown in Fig. 2b. The source in laboratory 1 sends one photon to Bob's measurement station and the other to Alice, while the source in laboratory 3 will send one photon to Alice the other one to Charlie. In this case, analogously to before, the measurement operators will be the following:

RESULTS
Our main goal is to experimentally obtain lower bounds on the fidelities between the produced states and the reference states, based on the statistics observed through the two apparatuses shown in Fig. 2a-b. However, we cannot simply apply the SWAP method using the raw statistics, because due to finite statistics the experimental frequencies do not correspond to physically allowed correlations (for instance, they violate the no-signalling conditions) and the optimization constraints imposed by the NPA would make the problem infeasible. Therefore, to overcome this problem, we proceed in the following way (see Supplementary Information for more details): 1. We use a regularization method in which we approximate the experimental frequencies f j by probability distributions belonging to the NPA set Q 4 [50], where j = (a, b, x, y) or j = (a, b, c, x, y, z) depending on the scenario. The solution of this method provides a no-signalling set of distributions P reg j that are guaranteed to be close to the set of quantum distributions. See Experimental apparatus and detection of four-fold coincidences. a) Parallel self-testing scenario implementation. In this scenario, the involved laboratories are two, both equipped with a quantum state source and two measurement stations. The measurement stations in Laboratory 1 represent Alice, while those in Laboratory 2, Bob. Both the sources generate a bipartite entangled state, send a subsystem to the other laboratory through a ∼ 30 m long single-mode fiber and keep the other one to measure it. The target state to be shared between the two parties involved in this network is the tensor product of two maximally entangled two-qubit states. b) Tripartite scenario implementation. In this scenario, the involved laboratories are three, two (1 and 3) equipped with quantum state source and a measurement station (representing, respectively, Bob and Charlie), and one, just with 2 measurement stations, constituting Alice. The source in laboratory 1, Λ1, sends one photon to Bob and one to Alice, while Λ2 sends one to Charlie and the other one to Alice, through a ∼ 30 m long single-mode fiber. c) In order to detect the couples of 2-qubit states generated by sources Λ1 and Λ2, we set a coincidence window within which two two-fold coincidence events must occur, to be recognized as a four-fold coincidence. The curves indicate the value of the CHSH of the states generated by the two sources (Λ1 blue curve, Λ2 orange curve), with the statistical uncertainty obtained considering the Poissonian distribution of the events (1 standard deviation), in terms of such a 4-fold coincidences window. The optimal time interval for this window is then chosen by evaluating the weighted mean of the two violations and choosing the highest one. In this case, it amounts to 1.033 µs.
Methods for details.
2. We then run the SWAP SDP using P reg j as inputs.
The solution of this SDP provides a linear functional d(P ) = j c * j P j (through its dual formulation) for which the value gives a lower bound on the self-testing fidelity.
3. We run another SWAP SDP in which we do not assume the actual value of the distributions, but we impose as a constraint the experimentally obtained value of the dual functional.
Using this method and a Monte Carlo simulation (assuming Poissonian statistics for the experiment) to calculate the uncertainties, we found that F (ρ AB , |Ψ 2 ) = 0.587 ± 0.053 for the parallel configurations and F (ρ ABC , |Ψ 3 ) = 0.863 ± 0.032 for the tripartite case. In detail, in Fig.3, we show the observable terms of the probability distribution given as solution by the SDP run in the third step, in comparison to the experimental frequencies.

Device-Independent estimation of the uncertainty
The previous method for calculating experimental uncertainties are not fully device-independent as it assumes a Poissonian distribution for the measurement results. We now move a step forward in removing assumptions and quantify the confidence level on the fidelities bounds exploiting Hoeffding's inequality [51], which holds for independent variables that are in the range (0, 1).
For this second method, we relax the constraint that we obtain a given value for the SDP functional d(P ) and only assume bound to it given by d(P ) ≤ d exp + τ ( ) and d(P ) ≥ d exp − τ ( ). In this notation, 1 − constitutes the confidence level that the observed frequency f j are within a range t j ( ) from the real probability P j [51]. More specifically, such interval t j ( ) amounts to − ln( ) 2nj , where n j is the number of registered counts for configuration j. At this point, by the central limit theorem [52], the linear combination of frequencies d exp will be characterized by a Gaussian statistics, whose variance amounts to σ 2 = j c * 2 j V ar(f j ). Furthermore, given that Var(f j ) = 1 2nj , parameter τ ( ) was chosen as follows (see the Supplementary Information for full derivation): In the end, the confidence level that the true value of d(P ) is within a range of τ ( ) from d exp can be easily recovered from a standard normal table, considering that this interval amounts to − ln( ) standard deviations. In Fig. 4, we show the certifiable fidelities in the two studied scenarios, versus and indicate the corresponding number of standard deviations adding up to τ ( ) at the bottom of the bars. From such fidelities, we can extrapolate a lower bound on the Schmidt number of the state [53,54]. In particular, in the studied cases, a fidelity higher than k 4 implies a Schmidt number (SN) higher or equal to k + 1. In Fig. 4, we indicate the thresholds for SN = 2 (blue line), 3 (red line) and 4 (green line).

DISCUSSION
In this work, we experimentally implemented deviceindependent self-testing protocols in two scenario represent-ing two basic building blocks for quantum networks. In particular, we studied a parallel self-testing scenario, in which two parties share two copies of a bipartite state, and a tripartite one, in which two bipartite states are shared among two peripheral nodes and a central one [44][45][46][47][48][49]. In these two cases, we were able to obtain lower bounds on the fidelity of the generated states with the desired target states that demonstrate that both sources indeed produce entangled states. In detail, we found lower bounds on the Schmidt number of the generated states, certifying a SN ≥ 2 for the first case and a SN ≥ 3 for the second one.
We stress that we considered the scenario presented in Ref. [42] where the number of local measurement choices are kept constant independently of the number of quantum state copies that one aims at certifying. Because of this, in our implementation each party had to implement only two measurements basis each, although the self-testing protocol considered two copies of maximally entangled states. From an experimental perspective, this method provides an extra significant advantage, represented by the fact that this technique requires only separable measurements. In particular, this result shows that, with a simple platform and low resource expense, it is We show the lower bounds on the certifiable fidelities, in both of the studied scenarios, as a function of the probability that for all configurations j, fj / ∈ {Pj − tj( ), Pj + tj( )}, where fj are the observed frequencies, Pj the probabilities characterizing the probability distribution underlying the experiment and nj the registered counts corresponding to configuration j. Furthermore, by Hoeffding's inequality [51], tj( ) = − ln( ) 2n j . From those statistical uncertainties on the probabilities, we recovered the confidence level on the obtained lower bounds, which amounts to − ln( ) standard deviations of a Gaussian distribution (see Supplementary Information for a full derivation). In a), we report the case of parallel self-testing, while, in b), the tripartite one. The numbers on the bars indicate the number of standard deviations corresponding to the confidence level of such fidelities. Such confidence levels, then, can be found on a standard normal table. The dashed lines indicate the fidelity values that certify, respectively, that the state has a Schmidt number higher or equal to 2 (blue line), 3 (red line) and 4 (green line).
possible to obtain non-trivial device-independent estimates of the quality of quantum networks. We believe that the present tool will find future applications in quantum communication, in particular in cryptographic scenarios such as quantum key distribution and blind quantum computation.

Experimental details
For the experimental setups of Fig. 2, the pump laser for source 1, with λ = 392.5 nm are produced by a second harmonic generation (SHG) process from a Ti:Sapphire mode locked laser with repetition rate of 76 MHz. Photon pairs entangled in the polarization degree of freedom are generated exploiting type-II SPDC in 2 mm-thick beta-barium borate (BBO) crystals. Source 2, instead, employs a continuous wave diode laser with wavelength of λ = 404 nm, which pumps a 20mm-thick periodically-poled KTP crystal inside a Sagnac interferometer, to generate photon pairs using a type-II degenerate SPDC process. The photons generated in all the sources are filtered in wavelength and spatial mode by using narrow band interference filters and single-mode fibers, respectively.

Coincidence counting
The photon detection events were collected and timed by a different time tagger device for each party, located in the corresponding laboratory [43]. For each 1 s of data acquisition the events were sent to a central server, along with a random clock signal shared between all the time-taggers, which was used to synchronize the timestamps of events relative to different devices. To filter out part of the noise the raw data was first pre-processed by keeping only double coincidence events for each photon source, using a narrow coincidence window of 1.05 ns. Then, the 4-fold coincidence events between the two sources were counted every time one of such double coincidence event was recorded for each source in a window of ∼ 1.033 µs.    In this section, we describe in further detail the self-testing approach we adopted in this work, firstly introduced in [1], which can be applied when experimental imperfections do not allow to reach the maximal quantum violation of a causal constraint, like the CHSH inequality. Indeed, in this case, the violation reveals the presence of non-classical correlations, but it does not single out the system producing it. However, through the techniques introduced in [1], it is possible to establish a lower bound on the fidelity between the target state and the unknown generated quantum state, resorting to the Navascués-Pironio-Acín (NPA) hierarchy [2,3].
Indeed, let us consider the simplest quantum scenario, with one source of bipartite entangled system and two parties (Alice and Bob) performing 2-output local measurements on a given subsystem, according to an input (x, y) ∈ (0, 1). The figure of merit for a self-testing of a two-qubit maximally entangled state is the CHSH inequality, but, if the maximal violation extent is not achieved, for instance due to experimental imperfections, the test is inconclusive.
Hence, we can resort to the following protocol: first of all, we use, as figure of merit, the fidelity with a target state |ψ target , defined as follows: However, without making assumptions on the dimension of the generated state, this fidelity is not properly defined, therefore, to have a fully device-independent protocol, we can resort to the SWAP operator [1]. As a first step, we consider an ancillary register, defined on an Hilbert space of the same dimension than that of the target, and we trust that each system is prepared in the dummy state, i.e. |0 . In our case we take dim(H target ) = 2. Then, we define the following local operator S = S AA ⊗ S BB , where S AA = U AA V AA and and and analogously for Bob (B). These operations are unitary if both O 0 and O 1 are unitary and Hermitian. Through this operator, we aim to swap part of the state ρ AB , which is seen as a black box, onto the ancillary register, and ρ swap will have the following form: , the entries of ρ swap , from the partial trace in Eq. (4), are given by linear combinations of correlation terms from the set c = (c I = T r(ρ AB I), ). Hence, we can solve the following semidefinite program (SDP): where Q l is a set which includes the one of quantum correlations and which corresponds to the l − th level of the NPA hierarchy [2,3]. In this way, by simply evaluating the CHSH inequality on the generated state, and putting it as a constraint in problem (5), we can lower bound the fidelity with the target state. To obtain higher bounds, we can add further constraints, for instance if the statistics corresponds to isotropic black boxes, we could add constraints of the following kind: A better lower bound can be obtained giving the full statistics to the SDP program, solving: given p(a, b|x, y) f = min ψ target | ρ swap |ψ target s.t. c ∈ Q n The bound on the certifiable fidelity, considering a generated state ρ AB = v|ψ ψ| + (1 − v)I/2 is plotted in Fig. 1 as a function of the visibility v (dotted curve). This method is general and can be extended to more parties and different topologies and the operators O 0 and O 1 must be chosen according to the target state, although general methods have been introduced, for the case in which one simply has some distributions p(a, b|x, y) and wishes to guess the state and measurements involved [4].
A. The SWAP method for parallel self-testing In the first part of our experiment, we use this protocol for a scenario involving two independent sources 1 and 2 generating singlets and two parties performing 4-output local measurements on the subsystem they get, according to an input (x, y) ∈ (0, 1). Both observers own dummy states with total dimension equal to the one of the target Hilbert space, which is dim(H target ) = 4, so Alice will have |00 A 1 A 2 , similarly to Bob. The fidelity to be bounded takes the form: where the target state is |ψ target = |ψ ⊗ |ψ and |ψ is defined as follows: which is maximally entangled and therefore equivalent to |ψ − up to local unitaries. This is chosen for simplicity of notation since this state reaches I CHSH = 2 √ 2 for the operators: to be inserted in the expression for the operators U and V defined in Eq. 2 and in Eq. 3, where σ x,z are Pauli operators. Alice and Bob then perform a swap between the dummy states on their ancillary register and the subsystem they receive, and the ρ swap density matrix takes the following form: where the SWAP operator is defined as S = S A1A 1 ⊗ S B1B 1 ⊗ S A2A 2 ⊗ S B2B 2 . Let us write the complete expression for the SWAP operator acted by Alice S AA = S A1A 1 ⊗ S A2A 2 , that is analogous for Bob: Since the operators U and V act on a subsystem of dimension 2, Alice and Bob can get 4 possible outcomes from their measurements, and this can be described by defining A AB Ax is the projector on the eigenspace corresponding to the eigenvalue (a = 0, 1, 2, 3) of operator A x . With these definitions, the expression of U and V operators has the following form: After finding the complete expression of ρ swap , the next step is to solve the semi-definite program defined in Eq. 6, where probabilities p(a 1 , b 1 |x, y) and p(a 2 , b 2 |x, y) are given. For this scenario, it is sufficient to choose Q n with n = 3, i.e. NPA hierarchy [2,3] at level 3, since the correlation term with highest degree appearing in the fidelity is c A0A1A0B0B1B0 = T r(ρ A1B1A2B2 A 0 A 1 A 0 B 0 B 1 B 0 ). The bound on the certifiable fidelity, considering a generated state ρ A1,2B1,2 = ρ A1B1 ⊗ ρ A2B2 , where ρ A1B1 = ρ A2B2 = v|ψ ψ| + (1 − v)I/d is plotted in Fig. 1 as a function of the visibility v (blue curve).
Supplementary Figure 1: Certifiable fidelities in the bilocality and parallel self-testing scenarios. The minimum fidelity certifiable with our protocol, respectively in the simple CHSH scenario (dashed curve) [1], in the bilocality scenario (red curve) and in the case of two parties sharing two singlets, i.e. parallel self-testing scenario, (blue curve) [5]. In the first case, the target state would be a 2-qubit maximally entangled state, |ψ , while in the other two cases, the target would be a product of two 2-qubit maximally entangled states, respectively |ψ A1B ⊗ |ψ A2C (bilocality scenario) and |ψ A1B1 ⊗ |ψ A2B2 (parallel self-testing scenario). The fidelities are plotted in terms of the visibilities of the generated states. In particular, we consider that, in the CHSH case, the generated state would have the following form: In the bilocality case, In the parallel self-testing case, analogously, The inset highlights the bilocality and parallel self-testing cases and the black dashed line indicates the trivial fidelity.

B. The SWAP method for the tripartite scenario
Let us now consider the case in which the states generated by sources 1 and 2 are sent to three parties (Alice, Bob and Charlie) that perform local measurements according to an input (x, y, z) ∈ (0, 1). In detail, Bob and Charlie perform 2-output local measurements on their subsystem, respectively generated by source 1 and 2, while Alice performs 4-output local measurements on the subsystem generated by both sources. In such scenario, self-testing could be achieved for a state product of two-qubit maximally entangled states, either performing two separate CHSH tests or a bilocality test [6][7][8], when observing the maximum violation. However, if we violate the classical bound, but we do not observe the maximum violation extent or we want to self-test another state, these tests are not conclusive. Therefore, we wish to bound the fidelity between the generated state and the target state, i.e.: where the target is a state product of |ψ , defined in Eq. 8. In order to test the state, the participants own an ancillary register each, with a dummy state of dimension 2 for Alice, and 1 for Bob and Charlie. As in the previous case, the parties perform the swap, and the resulting ρ swap density matrix will take the form: The SWAP operator is defined as S = S A1A 1 ⊗ S BB ⊗ S A2A 2 ⊗ S CC , Alice's one having dimension 2 and thus defined as in Eq. 12 and Eq. 13, while Bob and Charlie's having dimension 1 and thus defined as in Eq. 2 and Eq. 3.
The SDP finds a lower bound starting from the probabilities p(a 1 , b|x, y) and p(a 2 , c|x, z), but the computational power required is higher in this case, since the correlation term with highest degree is c A0A1A0B0B1B0C0C1C0 = T r(ρ A1B1A2B2 A 0 A 1 A 0 B 0 B 1 B 0 C 0 C 1 C 0 ), and level 5 of the hierarchy would be necessary. The bound on the certifiable fidelity, considering a generated state ρ A1,2BC = ρ A1B ⊗ ρ A2C , where ρ A1B = ρ A2C = v|ψ ψ| + (1 − v)I/d is plotted in Fig. 1 as a function of the visibility v (red curve).

Supplementary Note 2. NPA HIERARCHY
In this section we discuss the Navascués-Pironio-Acín hierarchy, first introduced in [2], which represents an efficient tool to solve optimization problems over the set Q of all quantum states and measurements of arbitrary dimension.
Let us consider the case of two parties, Alice and Bob, sharing a state ρ and performing local measurements that, without loss of generality, we consider to be projectors, indicated with E α,β respectively for Alice and Bob. Let us remind the properties of the projection operation: (i ) µ E µ = I, (ii ) E µ E ν = δ µ,ν for E µ and E ν belonging to the same measurement, and (iii ) [E α , E β ] = 0, namely projectors belonging respectively to Alice and Bob commute with each other. Now, let us assume that for a given probability distribution P αβ there exist a quantum state ρ and a set {E µ } of projection operations such that The implications coming from this assumption give necessary conditions that must hold for the set {E µ }. Indeed, let us consider the matrix Γ defined as follows: where S = {S 1 , S 2 , ..., S n } is the set made by all the operators of the form E α E β E α or α c α E α . It can be shown [2] that for Γ, together with Hermitianity, the following linearity conditions hold: and it is positive semidefinite: It follows that, if for some set S it is not possible to find a Γ matrix with the aforesaid properties, we must conclude that the probability distribution P αβ cannot be reproduced through local measurements on a quantum state. Additionally, it can be noted that if a set S can be written as linear combinations of operators in another set S , the conditions imposed by S are at least as constraining as the ones imposed by S, and that the set S m of all possible products of m projectors generates by linear combinations all the operators that are linear combinations of products of m projectors, with m ≤ m [2]. Hence, given a set of projectors {E µ }, the NPA method allows to build in a hierarchical way all the conditions, i.e. the sets S, to be satisfied in order to certify the quantumness of a given probability distribution. The hierarchy is based on constructing sets S n made up of products of the given operators until degree n: for instance, we can first consider the set containing just the projectors of the two parties, Alice and Bob, that we indicate by S 1 = {E α , E β }. If conditions (18), (19), (20) hold, we can go on and verify the constraints imposed by the set S 2 = {E µ E ν }, made up of all the products of the operators {E α , E β }, then iterating the process until some condition fails, or until the set S n is sufficiently large to represent all the conditions we need to check.
For a concrete example, let us consider Alice performing measurements X = 1, 2 and Bob performing measurements Y = 3, 4, where each measurement yields one of the two possible outcomes {a, b} ∈ {+1, −1}. Let us also define the correlation functions C XY = ab abP (a, b|X, Y ) and the marginal quantities C X = a aP (a, b|X, Y ) and C Y = b bP (a, b|X, Y ). The test of the NPA hierarchy corresponding to the first level, i.e. n = 1, is built on the set S 1 = {X 1 , X 2 , Y 1 , Y 2 }, and the Γ matrix will have the following form: where the lower symmetric part has been elided (remember that Γ is hermitian), and the entries u and v correspond to the correlation terms involving non-commuting measurements (both belonging to Alice or Bob). These entries are thus indeterminate and can be adjusted by use of a semi-definite program in order to get an SDP Γ matrix, if the other correlations {C X , C Y , C XY } are quantum [2].
A. Application to SDP characterization in parallel self-testing In the parallel self-testing scenario the function to be minimized, i.e. the fidelity defined in (7), is a linear combination of correlation terms containing products of 2 (one for Alice and one for Bob) up to 6 (three each) operators. To minimize such function on a quantum set, we need to impose these terms to be the entries of one positive semidefinite Γ matrix, and given that, as shown in the previous example, level 1 of the hierarchy considers constraints on correlation terms up to degree 2, in our case it is necessary to build the hierarchy up to level 3, at least. To do so, we exploit the Python library first introduced in [9], which allows to set both the objective function to be optimized and the level of the relaxation, and then to solve the SDP problem with MOSEK optimization software [10]. Now, considering that Alice and Bob own 8 projector operators each, indicated with Π a A/Bx , where a = {0, 1, 2, 3} are the possible eigenvalues, and x ∈ (0, 1), and given also property (i ) of the projection operation, yielding Π 3 , at level 1 the S 1 set contains the following 12 elements Given also that the number of the Γ entries grows with the level of the relaxation and with the number of operators involved, and decreases with the number of allowed substitutions, namely if entry , than the two are the same variable, due to property (ii ) of the projection operation, the total number of variables concerned in this optimization problem is found to be 22602.
However, it is possible to further reduce this quantity by noticing that our objective function is symmetric under exchange of parties A ↔ B, and that this transformation keeps also the matrix Γ unchanged. In particular, we added to the aforementioned substitutions those of the kind: where we took care to respect the commutation rules of Alice and Bob's operators when defining LHS and RHS of each substitutions.

B. Application to SDP characterization in the tripartite scenario
Let us now consider the scenario in which Alice shares one singlet with Bob and another with Charlie. Here, the fidelity to be bounded is the one in (14), which contains correlation terms of at least 3 (one for each party), and at most 9 (three for each party) operators. In this case, the level 1 of the hierarchy would be built over the }, containing 10 projection operators, then at the minimum required level, that is 5, the problem would be too computationally requiring for a normal computer. To simplify the problem, it is possible to stop at level 3 of the hierarchy, though considering all the monomials with degree higher than 6 appearing in the objective function as extramonomials. In this way, not all the possible product combinations of the operators are taken into account as variables, for correlation terms of degree > 6, but only those given as extramonomials. The number of variables concerning the optimization is, considering also the standard substitutions due to the projection properties, is 10115. Besides this, it is possible to further reduce the number of the Γ entries by using, again, the symmetries of the objective function. This time, the symmetry holds for the exchange of the peripheral parties, i.e. Bob and Charlie, B ↔ C , and swapping the second and third output of the central party, i.e. 01 ↔ 10. By considering this symmetry when setting the list of the substitutions, the number of SDP variables decreases to 7670, while the time for the optimization decreases from 136.28s to 78.16s.

Supplementary Note 3. FREQUENCY REGULARIZATION AND DEFINITION OF THE DUAL EQUALITY CONSTRAINT
Navascués, Pironio and Acín hierarchy is just one of the several examples of theoretical tools taking into account the full quantum distribution for device-independent characterization, which all share the common assumption of non-signaling. However, due to finite statistics, raw distributions obtained by the experimental frequencies generally do not satisfy this condition, to the extent that it is almost always impossible to use raw data within these methods. In [11], a general tool called the device-independent least-square method is introduced, which can be used to estimate the non-signaling probability distribution, belonging to a superset of the quantum correlation set, Q l , which is the closest to the experimental frequencies, in terms of a norm-2 distance (|| f − P || 2 ). In particular, we indicate with f the experimental frequencies and P is the probability distribution satisfying the no-signalling condition. Let us briefly remind the non-signaling condition: a distribution of probability P (a, b|x, y) is non-signaling if the following two relations are simultaneously satisfied: P (a|x, y) ≡ b P (a, b|x, y) = P (a|x, y ), ∀a, x, y, y P (b|x, y) ≡ a P (a, b|x, y) = P (b|x , y), ∀b, y, x, x In detail, this technique is shown in [11] to be equivalent to performing a projection P Π ( f ) of f onto an affine subspace N of R 16 (16 is the dimension of the frequency vector in a simple CHSH 2-parties scenario), which contains only P s satisfying the non-signaling condition, followed by the minimization of the 2-norm distance between P Π and Q.
Formally, the method amounts to finding the unique minimizer of the least-square problem: Where, when the quantum set Q is approximated by a superset relaxation Q l that admits a semi-definite programming characterization, through the NPA hierarchy. Problem (24) can be cast as an SDP, although the norm-2 distance is not a linear function, in the following formulation, using the characterization of positive semidefinite matrices via their Shur's complement: where I is the identity matrix having the same dimension of the column vector f , and Γ is defined as in Eq.(21), depending on the chosen NPA hierarchy level. Indeed, by Theorem 7.7.7 of [12], given that ( and hence s || f − P || 2 . This problem only involves an objective function and matrix constraints that are linear in the SDP variables s, P and in the variables of the Γ matrix, and is shown [11] to be unique and totally equivalent to the one defined is Eq. (24).
At this point, the obtained probability distribution P is injected in the SDP program as linear constraints on the moments of the Γ matrix, as described in Eq. (6) of the Supplementary information or Eq. (10) of main text, where the p are the regularised probabilities. Now, consider that given C ∈ M n , A i ∈ M n , i = 1, 2, ..., m, and b ∈ R m , the semidefinite programming problem is to find a matrix X ∈ M n for the following optimization problem: inf C · X subject to A i · X = b i | i = 1, ..., m | X 0 (28) which individuates the primal solution to the problem, i.e. searches for the minimum solution from above. The corresponding dual problem can be written as: Here, y ∈ R m and S ∈ M n constitute the solution of the dual problem, namely the best approximation of the result from below, which coincides with the primal, due to strong duality. Note that the objective function of the dual problem, b T y, is, in our case, a linear combination of the regularized frequencies m i c i P i . At this point, we took the dual solution of the SDP, namely, the vector of coefficients y, and evaluated the linear combination D exp = b T exp y by putting as b exp the observed frequencies. In the end, this new linear combination of the matrix moment is given as equality constraint to a new SDP which finally gives our device-independent estimation of the minimal certifiable fidelity, as follows: This method allows to extract an experimental lower bound on the fidelity, directly from the experimental data and not from the regularised one, avoiding a possible overestimation of it.

Supplementary Note 4. BOUNDING THE CONFIDENCE LEVEL ON THE FIDELITY BY USE OF HOEFFDING'S INEQUALITY
Hoeffding's inequality [13] represents a very useful tool when there is need to estimate uncertainties on experimental probabilities, taking into account finite statistics. It asserts that given n independent random variable X 1 , ..., X n , defined in [0, 1], with the following mean, expected value and variance: S = (X 1 + X 2 + ... + X n ) X = S/n µ = EX = ES/n σ 2 = Var(S)/n Then, the following inequality holds, for 0 ≤ t ≤ 1 − µ: where t depends on the arbitrary confidence interval assigned to the variable y = X − µ according to: t( ) = − ln( )/(2n). From this bound, it is possible to extract information about the effects of finite statistics on the estimation of the dual constraint, described in Supplementary Note 3. Let us consider, in fact, the worst-case scenario of Hoeffding's inequality (31), i.e. = e −2nt 2 . The RHS of this equality corresponds to one minus the cumulative distribution of the variable y, namely: For the Fundamental Theorem of Calculus, the probability distribution function P (y) for y > 0 can be obtained by differentiating one minus (32), thus yielding to: which is a well-defined, positive and normalized distribution of probability, when considering y > 0. For the case y < 0, Hoeffding's inequality still holds symmetrically [13], leading to the following an upper bound, for t > 0: This means that we can consider the two-sided variant of Hoeffding's bound, relative to the absolute value of the variable y, i.e.: P r( X − µ ≥ t) = 2e −2nt 2 Supplementary Figure 2: Hoeffding's distribution function. The black curve represents the normalised probability distribution of the variable y = X − µ > 0, P (y) = 2n |y| e −2ny 2 , while the red curve represents 1 minus the cumulative probability of such a distribution, i.e. the probability that y lies outside the range [0, t], given by Hoeffding's inequality [13]. Since, for the case y < 0, Hoeffding's inequality holds simmetrically and, hence, P (y) = P (−y), so, the probability that y lies outside the range [−t, t] amounts to P (|y| ≥ t, t > 0) = −t −∞ −2nye −2ny 2 dy + ∞ t 2nye −2ny 2 dy = e −2nt 2 . The variance of the black function is V ar(y) = 1/(2n). We used the variance of this distribution to compute the statistical uncertainty on the constraint c T (y) of problem (30), by asserting that it lies inside the interval d exp − τ ( ) ≤ c T y ≤ d exp + τ ( ) with a probability that can be recovered from a standardized normal table, considering that such interval amounts to − ln( ) standard deviations.
by defining normalized probability distribution over all real values of y, which is P (y) = 2n |y| e −2ny 2 and whose variance is defined as V ar(y) = 1 2n . For the central limit theorem, a linear combination N k=1 a k y k of N such independent, finite variance variables has a Gaussian distribution characterised by the following variance: Since the experimental value of the dual d exp is, actually, a linear combination of the experimental frequencies, its probability distribution can be approximated by a Gaussian, and we can estimate that the experimental dual lies within the following interval: where having defined σ 2 = N k=1 a 2 k Var(y k ), with a confidence level that can be computed from a standardized normal table considering that such interval amounts to − ln( ) standard deviations.