Correlations for computation and computation for correlations

Abstract

Quantum correlations are central to the foundations of quantum physics and form the basis of quantum technologies. Here, our goal is to connect quantum correlations and computation: using quantum correlations as a resource for computation—and vice versa, using computation to test quantum correlations. We derive Bell-type inequalities that test the capacity of quantum states for computing Boolean functions within a specific model of computation and experimentally investigate them using 4-photon Greenberger–Horne–Zeilinger (GHZ) states. Furthermore, we show how the resource states can be used to specifically compute Boolean functions—which can be used to test and verify the non-classicality of the underlying quantum states. The connection between quantum correlation and computability shown here has applications in quantum technologies, and is important for networked computing being performed by measurements on distributed multipartite quantum states.

Introduction

Since the beginning of quantum theory, the puzzling and non-local aspects of the theory have been a major topic of research in theoretical^1,2,3,4 and experimental physics^5,6,7,8 with the demonstration of loophole-free Bell tests being a key achievement^9,10,11,12. Besides its fundamental nature, entanglement is one of the key ingredients of quantum technologies and forms the basis for quantum communication and quantum computing. In quantum communication, Bell inequalities and testing correlations have practical applications and ensure the security of protocols and devices^13,14. Quantum computing shows speed-ups in certain computational tasks and it is believed that it will have tremendous impact¹⁵. Although, an advantage of quantum computers over classical computers has been shown recently for the first time^16,17, current quantum devices are not yet at a stage where they can solve large-scale problems. However, beyond full-power quantum computing, achieving an advantage in some form of non-classical computation is highly desirable¹⁸. The main goal of this work is to demonstrate a quantum advantage in computing with simple quantum resources and to develop tools that quantify the usefulness of the resources (see Fig. 1).

Fig. 1: Duality.

Illustration of the dual link between quantum states and computing Boolean functions: We can derive Bell-like inequalities to test whether certain quantum states can be used as a resource for computation. Vice versa, we can use the computation of Boolean functions as a test of quantum correlations.

While the most common model of quantum computation is the circuit model, measurement-based quantum computing^15,19, is computationally equivalent. Here, one first generates a universal entangled quantum state and the computation is carried out by successive, adaptive measurements on that state—measurement results are processed by a classical control and determine the settings of future measurements²⁰. Crucially, the classical control only needs computation with XOR and NOT gates, called linear side-processing. In this setting, adaptivity of the measurements is crucial: removing it disables determinism and makes universal quantum computing impossible²¹.

However, it has been shown that non-adaptive measurements on entangled states are a resource for universal classical computation. For example, three-qubit GHZ states and linear side processing (XOR and NOT gates alone) are sufficient to implement (universal) NAND gates²². More generally, a linear side processing combined with non-adaptive measurements on entangled resources is sufficient to realize non-linear Boolean functions and thus allows universal classical computation^23,24. This model is also referred to as NMQC_⊕—non-adaptive measurement-based quantum computing with linear side-processing²³. Another motivation for studying this model is that it is experimentally challenging to maintain coherence, and in a photonic setting, store resources for long enough to allow for adaptive measurement. The setting of NMQC_⊕ gives a formal framework for studying resources for measurement-based quantum computation in current and near-term experiments, reducing the need for fully adaptive measurements. In the setting of NMQC_⊕, computing a Boolean function deterministically requires a number of qubits that scales exponentially with the length of the input bit string²⁵. However, in the case of probabilistic computation of Boolean functions there is an advantage using even small-scale quantum resources.

Here, we build on this advantage and show the computation of non-linear Boolean functions with quantum resources (see Fig. 1) in this specific setting of NMQC_⊕. We link the violation of certain Bell-like inequalities to the capacity of quantum states for being a resource for computation. We experimentally generate GHZ states, the optimal states for this task^26,27, and demonstrate the violation of different Bell inequalities that are related to computing certain non-trivial Boolean functions.

The appearance of Bell inequalities in this work underlines the deep connection between the measurement-based model of quantum computation in this specific context and non-locality, or non-classicality in general. Within this framework, the non-classical nature of resources can be studied from a computational perspective and be characterized by the computations that they perform. This work demonstrates this connection with an experimental demonstration of a quantum advantage.

Results

The setting

The basic model of computing by non-adaptive measurements on quantum resources (within the framework of NMQC_⊕) is shown in Fig. 2. Let x ∈ {0, 1}ⁿ be the input, we aim to compute the Boolean function f: {0, 1}ⁿ ↦ {0, 1} (upper layer). Assume, the input x is generated with probability p(x).

Fig. 2: NMQC and Boolean functions.

a Concept of our setting to compute a Boolean function f(x) as described in the main text. The input x is transformed into a bit string s which determines the measurement settings M_j(s_j) on the physical resource. The outcomes of these measurements m_j determine the results of the computation: ${z:= \bigoplus }_{j = 1}^{l}{m}_{j}$ is generated by the parity of the outcomes from the quantum measurements. b Truth table for the three functions considered in this work with input string x = (x₁, x₂, x₃) and the bit string s = (s₁, s₂, s₃, s₄) (see main text for details).

First, the input x is processed by a linear side processor only, as in the middle layer. Here, the input bitstring x is transformed into another bitstring s ∈ {0, 1}^l with l ≥ n and

$${s}_{j} = \mathop {\bigoplus}\limits_{k = 1}^{n}{a}_{jk}{x}_{k}$$

(1)

for each j^th bit of s, where a_jk ∈ {0, 1} and ⨁ is summation modulo 2. The values a_jk can be seen as elements in an l-by-n binary matrix A (see Fig. 2a) and we can write s = (Ax)_⊕.

The j^th bit s_j now determines the settings for the measurement M_j(s_j) on the j^th qubit of the physical resource (bottom layer). For each measurement M_j(s_j), we obtain a measurement outcome m_j ∈ {0, 1}, associated with the eigenvalues ${(-1)}^{{m}_{j}}$. All outcomes m_j are collected in an outcome bitstring m ∈ {0, 1}^l. Note that the number of bits in the input x is distinct from the number of parties l in the physical resource. For example, in this work, we will focus on the case of n = 3 inputs and l = 4 parties. The reason we focus on this setting is because in the case of n = 2, the set of non-linear functions one can consider is limited to the NAND gate (up to additive terms modulo 2). The set of non-linear functions for n = 3 is less limited, and as we will demonstrate, allowing for a larger resource state (in this case l = 4) boosts the probability for quantum resources to compute a non-linear function.

We now ask ourselves how and when ${z:= \bigoplus }_{j = 1}^{l}{m}_{j}$, the parity of all outcomes, is equal to f(x), the designated Boolean function. To answer this question, as shown in Ref. ²⁵, one can determine the success probability for obtaining z = f(x) to be

$$p(z=f(x))=\frac{1}{2}(1+\beta ).$$

(2)

with

$$\beta =\sum _{x}p(x){(-1)}^{f(x)}E(x)$$

(3)

being a weighted sum of expectation values E(x) ≔ p(z = 0∣x) − p(z = 1∣x). Therefore, if β = 1, then E(x) = (−1)^f(x) for all x, and the function f(x) can be computed deterministically.

From Eq. (3), we obtain a Bell-like inequality, where the upper limit is determined by the physical resource:

$$\beta \le \left\{\begin{array}{ll}c\,,&\,\text{for classical resources}\,\\ q\,,&\,\text{for quantum resources}\,\end{array}\right..$$

(4)

Classical resources within this specific setting could simply be arbitrary measurements on an n-partite separable quantum state, where the statistics are convex mixtures of local probabilities. Alternatively, we can assume a local hidden variable model^22,25, or a non-contextual hidden variable model²³. These definitions of a classical resource are motivated by the assumption that there is no communication between the resources within this setting, and operations are local or that local measurements on one qubit commute with local measurements on another, respectively. The crucial point here is that all of these definitions of classical resources give rise to the same experimental predictions. Equivalently, we can assume the classical outcomes m_j are solely determined by the choice s_j and shared random variables between the parties. If we call the (finite) set of possible shared random variables Λ, which are distributed according to a probability distribution p(λ) where λ ∈ Λ, then the probability of observing the string of outcomes m giving the choices s will be

$$p(m| s)=\sum _{\lambda \in {{\Lambda }}}p(\lambda )\mathop{\prod }\limits_{j=1}^{l}p({m}_{j}| {s}_{j},\lambda ).$$

(5)

Given this notion of classical resource, it has been shown that the only functions f(x) that can be computed deterministically by classical resources are linear Boolean functions^23,25. In this way, classical resources have the same computational power as the linear side processor. As mentioned, quantum resources can have an advantage when the function f(x) is non-linear.

One can show that the maximal quantum bound q is achieved by GHZ states $\left|{\text{GHZ}}^{(l)}\right\rangle =({\left|0\right\rangle }^{\otimes l}+{\left|1\right\rangle }^{\otimes l})/\sqrt{2}$ and measurement of observables in the X-Y plane of the Bloch sphere ${M}_{j}({s}_{j})=\cos ({s}_{j}{\phi }_{j})X+\sin ({s}_{j}{\phi }_{j})Y$ for appropriately chosen angles²⁶. More details on the derivation of the equations above is given in the SI.

Inequalities for computation

We consider several functions, which are listed in Fig. 2b. These functions were chosen because they are examples of functions of different degree, i.e. h₃(x) is a degree-two function, where OR₃(x) and OR₃(x) ⊕ x₁x₃ are both degree-three functions, thus allowing for different forms of non-linear function. We choose one function and show how to derive the corresponding inequality in detail; and we list the results for the other functions.

Let us start with the function:

$${h}_{3}(x)={x}_{1}({x}_{2}\oplus {x}_{3}\oplus 1)\oplus {x}_{2}({x}_{3}\oplus 1)\oplus {x}_{3},$$

(6)

which leads us to the truth table shown in Fig. 2b and which is closely related to the pairwise AND function in Ref. ²⁵. First, the input bit string x is transformed by a linear side processor into measurement instructions s

$$\begin{array}{r}{s}_{1}={x}_{1},{s}_{2}={x}_{2},{s}_{3}={x}_{3},\,\text{and}\,{s}_{4}={x}_{1}\oplus {x}_{2}\oplus {x}_{3}.\end{array}$$

(7)

We will use this pre-processing for all examples of 3-bit Boolean functions in this work.

Now, our aim is to derive an inequality that tells us whether a certain physical resource is suitable for computing h₃(x). In order to do this, we make use of Eq. (3) and choose the uniform distribution p(x) = 1/8:

$${\beta }_{{\text{h}}_{3(\text{x})}}=\frac{1}{8}\left(\sum _{{x}_{1}={x}_{2}={x}_{3}}E(x)-\sum _{x\setminus ({x}_{1}={x}_{2}={x}_{3})}E(x)\right)\le \left\{\begin{array}{l}c\\ q\end{array}\right.\,,$$

(8)

with (−1)^f(x) according to the truth table in Fig. 2b. The maximal value of c can be obtained by maximizing Eq. (8) with E(x) = E(x₁)E(x₂)E(x₃) and enforcing that E(x_i) = ±1 for i = 1, 2, 3.

For the quantum case measurements are made on a four-qubit GHZ state. We have that for s_j equal to 0 or 1 the corresponding observables are given by the Pauli operators X or Y respectively, meaning that e.g. (s₁, s₂, s₃, s₄) = (0, 0, 1, 1) corresponds to a measurement of XXYY and we obtain the inequality shown in Fig. 4a. We obtain the classical and quantum bounds:

$${\beta }_{{\text{h}}_{3(\text{x})}}\le 1/2\,\,\text{vs.}\,\,1\,(c\,\,\text{vs.}\,\,q).$$

(9)

This means that if a physical resource violates the classical bound of this inequality, it is better suited for computing the function h₃(x) than classical resources, meaning it has a higher success probability to obtain the correct result. Quantum resources can deterministically compute this function if they have at least four qubits; for three qubits or less l = n = 3, the bound q is equal to $1/\sqrt{2}$²⁵.

Another function we consider is the three-bit OR function

$${\text{OR}}_{3}(x)={x}_{1}\vee {x}_{2}\vee {x}_{3}$$

(10)

which is only 0 for x₁ = x₂ = x₃ = 0 and 1 otherwise. With the distribution p(x = (0, 0, 0)) = 3/10 and p(x ≠ (0, 0, 0)) = 1/10, and measurement bases X/Y as above, we obtain:

$${\beta }_{{\text{OR}}_{3(\text{x})}}\le 4/10\,\,\text{vs.}\,\,8/10\,(c\,\,\text{vs.}\,\,q)\,,$$

(11)

where the value for q has been calculated according to Fig. 4a.

A similar example is the function

$${\text{OR}}_{3}(x)\oplus {x}_{1}{x}_{3},$$

(12)

for which we obtain

$${\beta }_{{\text{OR}}_{3(\text{x})}\oplus {\text{x}}_{1}{\text{x}}_{3}}\le 9/16\,\,\text{vs.}\,\,14/16\,(c\,\,\text{vs.}\,\,q),$$

(13)

with a distribution p(x) ∈ {1/16, 3/16} (see Fig. 4a), and again, measurement observables X and Y.

Finally, we aim at computing the two-bit AND function

$${\text{NAND}}_{2}(x)={x}_{1}{x}_{2}\oplus 1.$$

(14)

Choosing s₁ = x₁, s₂ = x₂, s₃ = x₁ ⊕ x₂ ⊕ 1 and s₄ = 1 and a uniform distribution p(x), we obtain the bounds

$${\beta }_{{\text{NAND}}_{2}(\text{x})}\le 1/2\,\,\text{vs.}\,\,1\,(c\,\,\text{vs.}\,\,q).$$

(15)

This computation is equivalent to the computation of a NAND using a three-qubit GHZ state is shown in Ref. ²². All these inequalities show that a quantum resource can violate the classical bounds for all Boolean functions considered here (see also Fig. 4). This means that the probability to compute the correct result is higher than with classical resources according to Eq. (2). For details on the derivations, see SI.

Computation for testing correlations

In the previous section, we used Bell-like inequalities to test whether certain physical resources are suitable for computing certain Boolean functions. Now, we would like to use computation to probe the non-classicality of the resource state. In other words, we perform computation in our model (see Fig. 1) and if we obtain the correct result with a certain probability, given by the inequalities above, we know our resource has to be non-classical in a particular, formal way.

Using Eq. (2) we can convert the classical and quantum bounds above into success probabilities

$${h}_{3}(x):\ \ \ 0.750\,\ \ \ \,\text{vs.}\,\,1.000$$

(16)

$${\text{OR}}_{3}(x):\ \ \ 0.700\,\ \ \ \,\text{vs.}\,\,0.900$$

(17)

$${\text{OR}}_{3}(x)\oplus {x}_{1}{x}_{3}:\ \ \ 0.813\,\ \ \ \,\text{vs.}\,\,0.938$$

(18)

$${\text{NAND}}_{2}(x):\ \ \ 0.750\,\ \ \ \,\text{vs.}\,\,1.000.$$

(19)

Here, the first value in each row indicates the maximum probability to obtain the correct results when the function is computed using classical resources, the second value indicates the probability for computing with quantum resources.

Experiment

For exploring relation between computation and Bell inequalities experimentally, we generate four-photon GHZ states using an all-optical setup that is shown and described in Fig. 3. The state we obtain in our experiment is

$$\left|\,\text{GHZ}\,^{\prime} \right\rangle =\left(\left|H,V,V,H\right\rangle -\left|V,H,H,V\right\rangle \right)/\sqrt{2}$$

(20)

with $\left|H\right\rangle \widehat{=}\left|0\right\rangle$ and $\left|V\right\rangle \widehat{=}\left|1\right\rangle$ denoting horizontal and vertical polarization. Note that the state $\left|\,\text{GHZ}\,^{\prime} \right\rangle$ is related to state $\left|{\text{GHZ}}^{(4)}\right\rangle$ by local unitary transformations, e.g. $\left|{\text{GHZ}}^{(4)}\right\rangle ={\mathbb{1}}XXZ\left|\,\text{GHZ}\,^{\prime} \right\rangle$. We verify the state obtained in the experiment through quantum state tomography²⁸. The reconstructed density matrix ρ_exp shows a fidelity $F=\left\langle \,\text{GHZ}^{\prime} | {\rho }_{\text{exp}}| \text{GHZ}\,^{\prime} \right\rangle$ of F = 0.82 ± 0.01 (see SI).

Fig. 3: Experimental setup.

A fs-pulsed Ti:sapphire laser at 780 nm is first frequency doubled and then passes through two non-linear β-barium-borate (BBO) crystals, each of which produces spatially separated single photons by type-II spontaneous parametric down-conversion (SPDC) in states: $\left|{\psi }^{-}\right\rangle =\left(\left|H,V\right\rangle -\left|V,H\right\rangle \right)/\sqrt{2}$. Half-wave plates (HWP) and additional BBO crystals compensate walk-off effects and allow to adjust the relative phase. The photons in modes 2 and 3 of the two states $\left|{\psi }_{12}^{-}\right\rangle \left|{\psi }_{34}^{-}\right\rangle$ are sent to a polarizing beam splitter (PBS). Upon post-selection to one photon in each of the output port, we obtain the state $\left|\,\text{GHZ}\,^{\prime} \right\rangle =\left(\left|H,V,V,H\right\rangle -\left|V,H,H,V\right\rangle \right)/\sqrt{2}$.

The values of β^exp we obtain for the individual Boolean functions are listed in Fig. 4b, together with the classical and quantum bounds. All values are clearly above the classical limit by more than at least 17 standard deviations.

Fig. 4: Results of computation.

a. Overview of the tested multipartite Bell inequalities with their corresponding Boolean functions. b–c. Classical, quantum, and experimentally obtained bounds for computing Boolean functions shown in the panels below. For both the bounds in Fig. 4b. and the probabilities Fig. 4c. the grey area of the bars indicate the regions completely obtainable with linear operations on classical resources. The coloured regions with the measurement points (error bars indicate the standard deviation) on top demonstrate that these limits have been surpassed in each case by ${\beta }_{{\text{h}}_{3}(\text{x})}:29\sigma$, ${\beta }_{{\text{OR}}_{3}(\text{x})}:19\sigma$, ${\beta }_{{\text{OR}}_{3}(\text{x})\oplus {\text{x}}_{1}{\text{x}}_{3}}:17\sigma$ and, ${\beta }_{{\text{NAND}}_{2}(\text{x})}:21\sigma$; showing that the correct computation of the non-linear Boolean function is more probable with quantum resources. The underlying white bars indicate the optimal quantum limits for 4 qubit measurements, while the grey dashed lines mark the quantum bound for three qubits, hence we see that with increased entangled resources the limit moves away form the classical bound unless it is already deterministic.

If we in turn use the GHZ state generated in the experiment to perform computation, we can quantify the probability to get the correct output. The corresponding probabilities are shown in Fig. 4c. This confirms that, for all probabilities, we lie above the classical values, which verifies that our physical resource must be quantum. In addition the grey dashed line in Fig. 4 highlights the limits of the computation on 3 qubits.

The discrepancy to the quantum bounds arises due to experimental imperfections. First of all, our resource state is not perfect. The quality of the 4-photon entanglement is limited by purity of the generated two-photon entangled states (two-photon fidelity F ≅ 0.96) as well as by the interference of the photons in modes 2 and 3 where we measured a Hong-Ou-Mandel dip visibility of V = 0.80 ± 0.02. In addition, imperfections in the polarization states, polarization drifts, and higher-order emissions (about 8% of the fourfold coincidences) reduce the quality of the generated GHZ state.

Discussion

In this work, we link a deeply fundamental question—the violation of a Bell inequality—to computing classical functions within a specific computational setting. We investigate this connection from two angles: verifying correlations through Bell tests quantifies the ability of a certain physical resource for computation. Furthermore, doing computation can be used as a tool to test non-classicality itself. We demonstrate this connection in a quantum optics experiment and show that already a four-qubit quantum state can provide an advantage.

The beauty of this connection between classicality and linear Boolean functions is that no matter how large the classical resource, its computational power does not change. However, as we increase the number of qubits in a quantum resource the computational power increases. Furthermore, as the number of qubits can be increased, it is possible to consider other families of non-linear Boolean functions as the input string to the computation could get larger. An interesting family is that of the Bent functions, which in the study of Boolean functions, are those that are, in some sense, the furthest away from the linear functions and exist only for an even number of bitstring inputs (cf.²³).

An interesting question is to study further types of classical resources in our setting. We could, for example, allow communication between measurement sites, as in study of multipartite non-locality^29,30. The computed Boolean functions have various implications on the amount of non-classicality of the resources. For one thing, these extra resources could enable the computation of non-linear Boolean functions, but perhaps not all functions. The amount of non-linearity needed could be a measure of how non-classical quantum resources can be. Given these or other additional powers, the central technical question is how the success probabilities of the enhanced classical resources compare to quantum resources. Some initial results in this direction were obtained in Ref. ³¹.

The relation between non-classicality and computing investigated here is related to the connection of Bell inequalities and quantum games³². It is also related to work on contextuality and the use of single-qubit operations for classical computation^33,34,35,36.

Even, if no fully fledged quantum computer is available, our work demonstrates the advantages of quantum resources for computation. In particular, our work has implications for quantum networks. Although, our approach here has been computational and not cryptographic, the quantum advantage in our work can be applied to a cryptographic setting if the shared resource state is distributed among agents in a network. For example, our methods could be directly used to self-test GHZ states³⁷ and generate randomness³⁸, both in a device-independent manner. Furthermore our quantum advantages for computation can be turned into an advantage for communication complexity³⁹. Thus, our work is a further example how the power of modern quantum technologies lies in fundamental quantum physics.

Methods

In the following sections we give some more details on the derivation of the different multipartite Bell inequalities presented in the main manuscript and show a more detailed data analysis. For a more general discussion of the theory we refer to Ref. ²⁵. First we outline what we mean by classical resources.

Definition of classical resources

A classical resource means that the probabilities (correlations) p(m∣s) of getting outcomes m: (m₁, . . . , m_l) given measurement choices s can be written as

$$p(m| s)=\mathop {\sum}\limits_{\lambda }p(\lambda )\mathop{\prod }\limits_{j=1}^{l}{\delta }_{{m}_{j},\mu ({s}_{j},\lambda )},$$

(21)

where {λ} is a set of shared random variables with probability distribution {p(λ)}_λ and μ(s_j, λ) ∈ {0, 1} is a map from the measurement choice s_j and λ to a bit-value. Equation 21 is arrived at from Eq. 5 by replacing the probabilities $p(m_j|s_j,\lambda)$ with indicator functions $\delta_{m_j,\mu(s_{j},\lambda)}$; this can be done by appealing to convexity of the probabilities in Eq. 5.

As already established, the measurement settings s are linear functions of the input bit-string x, and the outcome bit-value z is a linear function of the outcomes m. Thus from the statistics p(m|s) we obtain the probability p(z|x): the probability of getting bit-value z given input x. Furthermore given the form of the probabilities p(m|s) in Eq. 21, it can be derived that p(z|x) will only be a mixture of delta functions. Since the string s and bit-value z result from linear computation on x and m respectively, it can be seen that p(z∣x) will only be a mixture of delta functions δ_z,g(x), where g(x) is a linear Boolean function in x. When we wish to find the optimal classical bound of the inequalites below, by convexity, we only need to consider these delta functions $\delta_{z,g(x)}$. More formally, the set of classical correlations is a convex set and the optimal value of an inequality will be given by extreme points of the set. These extreme points are the deterministic correlations.

Such a classical model as above can be motivated in many ways: since, in the quantum case, there is no communication between the resources, and operations are local, a local hidden variable model is the natural classical analogue²⁵; since a local measurement on one qubit commutes with a local measurement on another, this motivates a non-contextual hidden variable model²³. Furthermore, if we associate a classical resource state with a separable quantum state, then the statistics will be convex mixtures of local probabilities.

Derivation of the inequalities

In the following, we show in detail how to derive the Bell inequalities for the Boolean functions f: xⁿ ↦ x given in the main manuscript:

$$\beta =\sum _{x}\beta (x)E(x)=\sum _{x}p(x){(-1)}^{f(x)}E(x)\le \left\{\begin{array}{ll}c\,,&\,\text{classic}\,\\ q\,,&\,\text{quantum}\,\end{array}\right.\,.$$

(22)

The function f(x) = h₃(x): The first example is the function h₃(x) = x₁(x₂ ⊕ x₃ ⊕ 1) ⊕ x₂(x₃ ⊕ 1) ⊕ x₃ and all input strings x are uniformly distributed p(x) = 2⁻³. The function satisfies f(0, 0, 0) = f(1, 1, 1) = 0, else it will yield the result 1 as shown in the table of Fig. 2 of the main section. Accordingly, from Eq. (22) we get the relation

$$\begin{array}{l}\beta =\frac{1}{8}\left(E(0,0,0)-E(0,0,1)-E(0,1,0)-E(0,1,1)\right.\\ - \, \left.E(1,0,0)-E(1,0,1)-E(1,1,0)+E(1,1,1)\right).\end{array}$$

(23)

For computing the largest classical bound, we maximize the sum in Eq. (23) over all possible values with E(x₁, x₂, x₃) = E(x₁)E(x₂)E(x₃) and the expectation values confined to E(x_i) = (±1) for i ∈ {1, 2, 3}. In Eq. (23) we obtain a maximal classical value of $\frac{1}{2}$.

The value for the quantum bound depends on the number of sites l≥n. In our measurements we have l = 4 and the measurement choices are encoded by the following linear map

$$\left(\begin{array}{l}{s}_{1}\\ {s}_{2}\\ {s}_{3}\\ {s}_{4}\end{array}\right)=\left(\begin{array}{lll}1&0&0\\ 0&1&0\\ 0&0&1\\ 1&1&1\end{array}\right){\left(\begin{array}{l}{x}_{1}\\ {x}_{2}\\ {x}_{3}\end{array}\right)}_{\oplus }\,.$$

(24)

Therefore, choosing for s = 0 observable X and for s = 1 observable Y the inequality Eq. (23) becomes

$$\frac{1}{8}\langle XXXX-XXYY-XYXY-XYYX-YXXY-YXYX-YYXX+YYYY\rangle \le 1\,.$$

(25)

The maximum in Eq. (25) is obtained exactly by $\left|{\text{GHZ}}^{(4)}\right\rangle =(\left|0,0,0,0\right\rangle +\left|1,1,1,1\right\rangle )/\sqrt{2}$, naturally the values of all other states will be within the region bounded by this GHZ state. Note that this value of 1 for the quantum bound is the maximum allowed algebraically.

The function f(x) = OR₃(x): The second example is the OR function. At this point it is also worthwhile to rewrite the function in algebraic normal form (ANF) so it can obviously be seen as non-linear:

$${\text{OR}}_{3}(x)={x}_{1}{x}_{2}{x}_{3}\oplus {x}_{1}{x}_{2}\oplus {x}_{1}{x}_{3}\oplus {x}_{2}{x}_{3}\oplus {x}_{1}\oplus {x}_{2}\oplus {x}_{3}.$$

(26)

The distribution for this function is chosen as $p(x)=\frac{1}{10}$, except for $p(0,0,0)=\frac{3}{10}$. We thus obtain the inequality

$$\beta =\frac{3}{10}E(0,0,0)-\frac{1}{10}\sum _{x\ne (0,0,0)}E(x).$$

(27)

In the same way as above we can compute the classical bound to be $\frac{4}{10}$ and the quantum bound

$$\frac{3}{10}\langle XXXX\rangle -\frac{1}{10}\langle XXYY+XYXY+XYYX+YXXY+YXYX+YYXX+YYYY\rangle \le \frac{8}{10}\,.$$

(28)

This bound can be readily confirmed to be the maximum allowed for all possible quantum resources (both states and measurements) using the methods described by Werner, Wolf, Żukowski and Brukner^26,27.

The function f(x) = OR₃(x) ⊕ x₁x₃: The third example is f(x) = OR₃(x) ⊕ x₁x₃ and the distribution $p(0,0,0)=p(0,0,1)=p(1,0,1)=p(1,1,1)=\frac{1}{16}$ and $p(0,1,0)=p(0,1,1)=p(1,0,0)=p(1,1,0)=\frac{3}{16}$. It should be noted that this function is still clearly non-linear after converting it into ANF using the identity described above. The correct signs can be read off from the truth table and we get

$$\begin{array}{l}\beta =\frac{1}{16}\left(E(0,0,0)-E(0,0,1)+E(1,0,1)+E(1,1,1)\right)\\- \, \frac{3}{16}\left(E(0,1,0)+E(0,1,1)+E(1,0,0)+E(1,1,0)\right)\end{array}$$

(29)

with a classical bound of $\frac{9}{16}$

$$\frac{1}{16}\langle XXXX-XXYY+YXYX+YYYY\rangle -\frac{3}{16}\langle XYXY+XYYX+YXXY+YYXX\rangle \le \frac{14}{16}\,.$$

(30)

Again, this bound can be readily confirmed to be the maximum for all quantum resources (both states and measurements) using the methods described by Werner, Wolf, Żukowski and Brukner^26,27.

The function f(x) = NAND₂(x): In the work of Anders and Browne²², it was demonstrated that a three-qubit GHZ state can be used to compute the NAND function of two bits, which we can write as NAND₂(x) = x₁x₂ ⊕ 1. Since this is just a function on two bits, things will be somewhat simplified. The distribution over these two bits is 2⁻² for all values of x ≔ (x₁, x₂). The inequality can be written as:

$$\beta =\frac{1}{4}\left(-E(0,0)-E(0,1)-E(1,0)+E(1,1)\right)\,.$$

(31)

The classical bound for this inequality is 1/2. In our setting for l = 4 parties we can also compute the function NAND₂(x) with the following linear map to generate the four inputs to the parties:

$$\left(\begin{array}{l}{s}_{1}\\ {s}_{2}\\ {s}_{3}\\ {s}_{4}\end{array}\right)=\left(\begin{array}{ll}1&0\\ 0&1\\ 1&1\\ 0&0\end{array}\right){\left(\begin{array}{c}{x}_{1}\\ {x}_{2}\end{array}\right)}_{\oplus }.$$

(32)

Now we can modify the setup so that the third and fourth parties measure Y for s_j = 0, and X for s_j = 1. Equivalently, we could have applied a NOT to the values of s_j. Thus the inequality in Eq. (31) can be rewritten as:

$$\frac{1}{4}\langle -XXYY-XYXY-YXXY+YYYY\rangle \le 1$$

(33)

The bound on the right-hand-side is attained with the GHZ state $(\left|0,0,0,0\right\rangle +\left|1,1,1,1\right\rangle )/\sqrt{2}$, as above; this bound is also the maximum allowed algebraically.

Bounds on three-qubit entanglement: Here we derive bounds on the inequalities above when the four parties are limited to sharing tripartite quantum resources. In this case three of the four parties can share a quantum state and the fourth party’s outcome is local, i.e. determined by the measurement choice s_j and shared randomness between the parties. Therefore, for each inequality, we can choose the party whose outcome m_j will be a deterministic function ${(-1)}^{{r}_{j}{s}_{j}}$, where r_j ∈ {0, 1}. Then we can optimize over all measurements and states for the other three parties with the same pre-processing as outlined above for each function.

To summarize, for all the functions above, and pre-processing outlined above, we can limit one of the four parties to having a deterministic outcome in the inequalities. In addition to this, we can adapt the numerical techniques in Refs. ^26,27, to find the optimal violation for the remaining three parties. The optimal violation for the remaining three parties will be attained by a three-qubit GHZ state. This optimization can be done for each choice of the party that does not share quantum resources with the other three.

For the function f(x) = x₁(x₂ ⊕ x₃ ⊕ 1) ⊕ x₂(x₃ ⊕ 1) ⊕ x₃, the bound on β for the inequality Eq. (23) for tripartite quantum resources is $\beta \le 1/\sqrt{2}$. Notably this bound relates to the maximal quantum violation of the Svetlichny inequality⁴⁰, as discussed in Ref. ²⁵. For the function OR₃(x), the bound on β for Eq. (27) for tripartite quantum resources is β ≤ 2/3. For inequality Eq. (29), the bound β ≤ 0.70235 holds. For the function NAND₂(x), the function can be computed deterministically by three qubits, as shown in Ref. ²².

Further experimental details and data analysis

For the experiment as illustrated in Fig. 3a 4 W fs-pulsed Ti:sapphire laser at 780 nm is first frequency doubled and then directed onto a non-linear β-barium-borate (BBO) crystal which produces spatially separated single photons by type-II spontaneous parametric down-conversion (SPDC). The two BBO crystals produce a polarization entangled state

$$|{\psi }^{(\varphi )}\rangle =\left(\left|H,V\right\rangle +{e}^{i\varphi }\left|V,H\right\rangle \right)/\sqrt{2}\,.$$

(34)

After the entangled photon pairs are emitted in conic sections the photons in each spatial mode pass a combination of half-wave plate (HWP) and another BBO-crystal which compensate previously induced walk-off effects and allow to adjust the relative angle φ in Eq. (34). By fixing φ = π, a product of two Bell-states $\left|{\psi }_{12}^{-}\right\rangle \left|{\psi }_{34}^{-}\right\rangle$ is obtained. Interference of modes 2 and 3 via a polarizing beam splitter (PBS) becomes

$${\text{PBS}}_{23}\left|{\psi }_{12}^{-}\right\rangle \left|{\psi }_{34}^{-}\right\rangle =\frac{1}{\sqrt{2}}\left(\left|\,\text{GHZ}\,^{\prime} \right\rangle +\left|\chi \right\rangle \right)\,.$$

(35)

The final state is thus a superposition of $\left|\,\text{GHZ}\,^{\prime} \right\rangle =\left(\left|H,V,V,H\right\rangle -\left|V,H,H,V\right\rangle \right)/\sqrt{2}$, a four-photon GHZ state and $\left|\chi \right\rangle =i\left(\left|H,HV,0,H\right\rangle +\left|V,0,HV,V\right\rangle \right)/\sqrt{2}$ . Any 4-fold coincidence measured at the 8 avalanche photodiodes (APDs) is thus only obtained by the GHZ state.

To check the quality of the experimental state ρ_exp we conducted a regular quantum state tomography by recording the 3⁴ combinations of the expected values Tr(O₁O₂O₃O₄ρ_exp), O_i = {X_i, Y_i, Z_i} being the Pauli operators. To calculate a state fidelity we first get an estimation of the distribution ρ_exp from the measured data by maximizing a likelihood function and then compare to the optimal state $F=\left\langle \,\text{GHZ}^{\prime} | {\rho }_{\text{exp}}| \text{GHZ}\,^{\prime} \right\rangle$²⁸. The tomography data was obtained at 500 mW pump power with a maximal 4-fold coincidence rate of ca. 1 Hz. We obtained a value of F = 0.824 ± 0.014; the corresponding density matrix is shown in Fig. 5. The error has been estimated from a Monte Carlo simulation, i.e. random sampling from a Poisson distribution and iterating the fidelity calculation 100 times. The expectation values E(x) that were recorded for the violation of the respective inequalities are exhibited in Fig. 5e.

Fig. 5: Tomography and expectation value.

Panels a–d. show the measured real and imaginary part of the quantum state ρ_exp calculated by a maximum likelihood estimation from the measured 4-fold coincidence counts as well as their ideal cases. e. Calculated expected values obtained by the probability measurements of 4² = 16 combinations to measure H or V polarized photons. The data was measured at 100 mW pump power.