# Device-independent quantum key distribution from generalized CHSH inequalities

Pavel Sekatski, Jean-Daniel Bancal, Xavier Valcarce, Y.-Z Tan, Renato Renner, Nicolas Sangouard

► To cite this version:

Pavel Sekatski, Jean-Daniel Bancal, Xavier Valcarce, Y.-Z Tan, Renato Renner, et al., Device-independent quantum key distribution from generalized CHSH inequalities. Quantum, 2021, 5, pp.444. 10.22331/q-2021-04-26-444. cea-02955219

HAL Id: cea-02955219

https://cea.hal.science/cea-02955219

Submitted on 1 Oct 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Device-independent quantum key distribution from generalized CHSH inequalities

Pavel Sekatski<fnref n="1" />, Jean-Daniel Bancal<fnref n="2" />, Xavier Valcarce<fnref n="3" />, Ernest Y.-Z. Tan<fnref n="4" />, Renato Renner<fnref n="4" />, and Nicolas Sangouard<fnref n="1" /> <fnref n="3" />

September 21, 2020

Device-independent quantum key distribution aims at providing security guarantees even when using largely uncharacterised devices. In the simplest scenario, these guarantees are derived from the CHSH score, which is a simple linear combination of four correlation functions. We here derive a security proof from a generalisation of the CHSH score, which effectively takes into account the individual values of two correlation functions. We show that this additional information, which is anyway available in practice, allows one to get higher key rates than with the CHSH score. We discuss the potential advantage of this technique for realistic photonic implementations of device-independent quantum key distribution.

## 1 Introduction

The aim of quantum key distribution (QKD) is to give two parties — Alice & Bob — the possibility to generate a secret key when they share a quantum channel. For instance, in the implementation proposed by Ekert [8], the channel consists of a source producing entangled particles that are distributed to Alice & Bob. At each round, each of Alice & Bob measure one particle by choosing one out of several measurement settings. The claim that Alice's measurement results are secure, i.e. unknown to any third party — Eve — who may control the quantum channel, is guaranteed by inferring (from Alice and Bob's measurement results) that the source emits states close to pure bipartite entangled states. This ensures at the same time that Bob's results are correlated to Alice's ones if he chooses an appropriate measurement setting, i.e. Alice and Bob's measurement results can form a secret key.

Ekert suggested that the information about the key that may be available to an adversary can be quantified by choosing settings allowing Alice & Bob to

violate a Bell inequality. This idea was later progressively formalised and led to what is now called device-independent QKD (DIQKD). In its simplest version, DIQKD is implemented by letting Alice choose randomly between two measurement settings at each round, $A_x$ where $x \in \{0,1\}$, while Bob's measurement includes three possible settings, $B_y$ where $y \in \{0,1,2\}$. For settings $x, y \in \{0,1\}$, the results — which can possibly take many values — are post-processed locally and turned into binary values $A_x, B_y \in \{-1, +1\}$. After several iterations, Alice and Bob communicate classically to estimate the CHSH score

$$S = \langle A_0 \otimes (B_0 + B_1) + A_1 \otimes (B_0 - B_1) \rangle \quad (1)$$

where $\langle A_x \otimes B_y \rangle = p(A_x = B_y|x, y) - p(A_x \neq B_y|x, y)$ quantifies the correlation between the outcomes for measurement choices $x$ and $y$, respectively. The remaining measurement setting $y = 2$ is chosen to generate an outcome $B_2$ that minimises the uncertainty with respect to $A_0$. Alice then forms the raw key from the outcomes $A_0$ of the pairs that Bob measured with respect to $y = 2$.

We consider $n$ such rounds, over which the source produces a tripartite state $|\Psi_{ABE}\rangle$ shared between Alice, Bob and Eve. Ref. [2] showed that Eve's information is the same as in the case where the devices have no memory and behave identically and independently in each communication round of the protocol, up to corrections vanishing with $n$. In particular, we can write $|\Psi_{ABE}\rangle = |\psi\rangle_{ABE}^n$ where $|\psi\rangle_{ABE}$ is the tripartite state of a single round and consider the case where measurements are done successively on the state $|\psi\rangle_{ABE}$.

In the asymptotic limit of large $n$, the number of secret bits per round obtained after one-way error correction and privacy amplification (i.e. the key rate) is then given by [6]

$$r = H(A_0|E) - H(A_0|B_2), \quad (2)$$

<footnote n="1">¹Department of Physics, University of Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland</footnote>

<footnote n="2">²Department of Applied Physics, University of Geneva, Chemin de Pinchat 22, 1211 Geneva, Switzerland</footnote>

<footnote n="3">³Université Paris-Saclay, CEA, CNRS, Institut de physique théorique, 91191, Gif-sur-Yvette, France</footnote>

<footnote n="4">⁴Institute for Theoretical Physics, ETH Zürich, 8093 Zürich, Switzerland</footnote>
where $H$ is the von Neumann entropy. Ref. [1] showed that the conditional entropy $H(A_0|E)$ optimized over all states $\psi_{ABE}$ and measurements $A_x, B_y$ compatible with the observed CHSH score $S$ is lower bounded by

$$H(A_0|E) \geq 1 - h\left(\frac{1 + \sqrt{(S/2)^2 - 1}}{2}\right) \quad (3)$$

where $h$ denotes the binary entropy. This provides a lower bound on the key rate, as the conditional entropy $H(A_0|B_2)$ can be estimated directly from Alice and Bob measurement results associated to setting choices $A_0$ and $B_2$. Interestingly, this bound is obtained device-independent, i.e. without assumptions on the dimension of quantum states and the calibration of measurements. This is not the case for standard (non-device independent) QKD protocols which are not based on the violation of a Bell inequality and whose security guarantees rely on the assumption that the source and measurements carry out precisely the operations foreseen by the protocol. This assumption is hard to meet in practice and leads to vulnerabilities, as demonstrated by hacking experiments [10, 18, 35, 38]. The robustness of device-independent quantum key distribution against these attacks makes it appealing, and a race between several experimental groups is ongoing to report the first proof-of-principle distribution of a key with a fully device-independent security. Measurement-DIQKD, a precursor of DIQKD where device-independence only applies to the measurement devices, but not to those used for state preparation [3, 17], already admit a number of experimental implementations [5, 9, 15, 26, 32].

Let us note that the proof leading to the bound given in Eq. (3) only uses the knowledge of the CHSH score. This score is computed as a linear combination of the correlation functions $\langle A_x B_y \rangle$, but the additional information provided by considering these correlations individually – which is anyway available in practice – might help to facilitate a realisation of device-independent quantum key distribution. This motivation is at the core of this work.

Concretely, we consider the individual values of two terms appearing in the CHSH score, namely

$$\begin{aligned} X &= \langle A_0 \otimes (B_0 + B_1) \rangle, \\ Y &= \langle A_1 \otimes (B_0 - B_1) \rangle. \end{aligned} \quad (4)$$

The use of values of X and Y proved to be useful for device-independent state certification by improving the certified fidelity from the CHSH score [34]. It is also expected to be useful in DIQKD, as the knowledge of X and Y allows one to differentiate the contributions of the key generating measurement $A_0$ from the ones associated with $A_1$, from which no key is generally extracted (see [27] for a noticeable exception).

Finally, in implementations of DIQKD with non-unit detection efficiencies where no-detection events are attributed a fixed value $\pm 1$, no-detections on Bob's side can only contribute to one of these two correlation functions (X or Y). The main result of this work is to confirm the intuition that the use of individual values of X and Y improves the bounds on Eve's information derived from the CHSH score, and hence the key rate of DIQKD.

## 2 Formulation of the problem

**Generalization of the CHSH test**– We are interested in bounding Eve's conditional entropy $H(A_0|E)$ appearing in Eq. (2) as a direct function of observed quantities X and Y. Formally, this can be accomplished by considering all possible quantum models ($\psi_{ABE}, A_y, B_y$) that satisfy $X_{\text{model}} \ge X$ and $Y_{\text{model}} \ge Y$. However, it is clear that the set of points $(X,Y)$ for which Eve's conditional entropy is bounded from below by some constant is convex: two quantum models giving some amount of information to Eve can be joined into a new model on which Eve's conditional entropy is bounded by the weighted sum of entropy bounds associated to the individual models. It is thus equivalent to bound Eve's conditional entropy with linear constraints of the form

$$\frac{\cos(\Omega)}{2} X_{\text{model}} + \frac{\sin(\Omega)}{2} Y_{\text{model}} \ge \beta, \quad (5)$$

where $\beta$ is deduced from the observed quantities X and Y from the following formula

$$\beta = \frac{1}{2} (\cos(\Omega) X + \sin(\Omega) Y). \quad (6)$$

(Further in the text, we will use a compact notation for sine $S_\Omega = \sin(\Omega)$ and cosine $C_\Omega = \cos(\Omega)$.) Obviously, $\Omega = \pi/4$ reduces back to the CHSH constraint (up to normalization). Just like the CHSH score can be seen as the result of a test of the CHSH inequality, we can associate $\beta$ to the test of a Bell inequality – a generalisation of the CHSH inequality – that is characterized in Ref. [?]. This characterization is also done below for the sake of completeness.

**Reduction to qubits**– The score $\beta$ given in Eq. (6) is estimated when Alice and Bob choose the measurement $A_x, B_y, x, y = 0, 1$, which are observables with eigenvalues $\pm 1$. A simple lemma – the Jordan lemma [31] – tells us that such observables can be jointly block diagonalised with blocks of size $2 \times 2$, i.e.

$$A_x = \bigoplus_k A_{x,k} \qquad B_y = \bigoplus_{k'} B_{y,k'}, \quad (7)$$

where, without loss of generality, we can assume the restriction to each qubit block to be a real Pauli operator satisfying $A_{x,k}^2 = 1$ and $B_{y,k'}^2 = 1$. This means
that in each block labelled by $k$ and $k'$ respectively, the measurement is characterized by unit vectors $\mathbf{a}_x^k$, $\mathbf{b}_y^{k'}$ such that

$$A_{x,k} = \mathbf{a}_x^k \cdot \begin{pmatrix} \sigma_z \\ \sigma_x \end{pmatrix} \qquad B_{y,k'} = \mathbf{b}_y^{k'} \cdot \begin{pmatrix} \sigma_z \\ \sigma_x \end{pmatrix} \qquad (8)$$

where $\sigma_z$ and $\sigma_x$ are Pauli operators.

The state $|\psi\rangle_{ABE}$ can be enforced to take the form

$$|\psi\rangle_{ABE} = \bigoplus_{k,k'} \sqrt{p_{(k,k')}} |\psi\rangle_{ABE}^{(k,k')}, \qquad (9)$$

where $p_{(k,k')}$ is a probability distribution and $|\Psi\rangle_{ABE}^{(k,k')} \in \mathbb{C}_A^2 \otimes \mathbb{C}_B^2 \otimes \mathcal{H}_E^{(4)}$, see Refs [13, 23] for detailed discussions. Given models with such measurements and state, the quantity of interest can be expressed as

$$H(A_0|E) = \sum_{k,k'} p_{(k,k')} H_{(k,k')} (A_0|E) \qquad (10)$$

where $H_{(k,k')}(A_0|E)$ is Eve's conditional entropy for four-qubit models (including the two qubits from Eve's purification) involving real Pauli measurements. If the minimization of $H_{(k',k)}(A_0|E)$ over such models satisfying $\frac{\cos(\Omega)}{2}X_{\text{model}, k, k'} + \frac{\sin(\Omega)}{2}Y_{\text{model}, k, k'} \ge \beta$ provides a convex function of $\beta$, this function can be used directly as a lower bound on the quantity $H(A_0|E)$ through Eq. (10). If it is not convex, it can be convexified so as to apply to all possible mixtures of state and measurement, and thus again apply to Eq. (10). This convexity property allows us to reduce the general problem of finding the minimum of Eve's conditional entropy over all possible models to a minimization over four-qubit models with real Pauli measurements. We will come back to this convexification requirement later.

**Noisy preprocessing**– We consider a simple post-processing of the raw key, known as *noisy preprocessing* [14, 24, 25], which has been shown to be beneficial to reduce the requirement on the detection efficiency in photonic implementations of device-independent quantum key distribution [13]. Once the raw key is obtained, Alice is instructed to generate a new raw key $\hat{A}_0$ by flipping each bit of the initial raw key with a probability $p$. Such a pre-processing can be accounted by the model of the key generating measurement, by considering a new operator

$$\hat{A}_0 = (1-p)A_0 - pA_0 \qquad (11)$$

which is not a Pauli operator anymore. Note that we will often parametrise the amount of noise that Alice

adds with a parameter $q = (1 - 2p)^2$.

**Symmetrization**– In order to simplify the analysis, it is convenient to consider a symmetrization step in which both parties, Alice and Bob, flip the outcomes of the key generating measurements depending on a public random bit string. This guarantees that bits of the raw key are random, i.e. $H(A_0) = H(\hat{A}_0) = 1$. Importantly, one can show the equivalence of protocols with and without symmetrization, meaning that the symmetrization does not need to be implemented in practice, see [13] for a complete description of the symmetrization step in the presence of noisy preprocessing.

**Reduction to Bell diagonal states**– If the constraints appearing in the minimization problem do not depend on the marginal probabilities $p(A_x|x)$ and $p(B_y|y)$ of Alice and Bob respectively, the symmetrization step previously presented reduces the model of the state to a Bell-diagonal structure

$$|\psi\rangle_{ABE} = \sum_{i=1}^{4} \sqrt{L_i} |\Phi^i\rangle_{AB} |i\rangle_E \qquad (12)$$

where $|\Phi^i\rangle = \{\Phi^+, \Psi^-, \Phi^-, \Psi^+\}_{i=1}^4$, and without loss of generality, a partial ordering of the eigenvalues $L_1 \ge L_2$ and $L_3 \ge L_4$ can be imposed [23]. Note that the superscripts $k, k'$ are omitted in the tripartite state appearing in Eq. (12), i.e. $|\psi\rangle_{ABE} \rightarrow |\psi\rangle_{ABE}^{(k,k')}$. Until the end of this section and in the next section which is dedicated to the resolution of the optimization presented in Eq. (19), we remove the index $k, k'$ for making the notation simpler and ask the reader to keep in mind that we consider the restriction to four qubit models with real Pauli measurements in this two sections.

**Eve's conditional entropy**– Eve's conditional entropy can be expressed as

$$H(\hat{A}_0|E) = H(\hat{A}_0) - H(\rho_E) + \sum_{\hat{a}=\pm 1} p(\hat{a}) H(\rho_{E|\hat{a}}) \qquad (13)$$

where $\rho_E$ is the reduced state of Eve and $\hat{\rho}_{E|\hat{a}}$ corresponds to Eve's state conditioned on Alice's noisy key bit $\hat{A}_0$ being equal to $\hat{a}$, which occurs with probability $p(\hat{a})$. The equivalence of the protocol with the symmetrized one allows us to take $H(\hat{A}_0) = 1$ and $p(\hat{a}) = \frac{1}{2}$.

$H(\rho_E)$ is given by the entropy $H(\mathbf{L})$ of the probability vector $\mathbf{L} = (L_1, \dots, L_4)$ while the term $\rho_{E|\hat{a}=+1}^{(\alpha, \beta)}$ equals
$$ \begin{pmatrix} L_1 & 0 & C_\phi \sqrt{L_1 L_3 q} & S_\phi \sqrt{L_1 L_4 q} \\ 0 & L_2 & S_\phi \sqrt{L_2 L_3 q} & -C_\phi \sqrt{L_2 L_4 q} \\ C_\phi \sqrt{L_1 L_3 q} & S_\phi \sqrt{L_2 L_3 q} & L_3 & 0 \\ S_\phi \sqrt{L_1 L_4 q} & -C_\phi \sqrt{L_2 L_4 q} & 0 & L_4 \end{pmatrix} \quad (14) $$

where $\phi$ labels Alice measurement $A_0 = \cos(\phi)\sigma_z + \sin(\phi)\sigma_x$ (we use the notation $C_\phi = \cos(\phi)$ and $S_\phi = \sin(\phi)$). The two states $\rho_{E|\hat{a}=\pm 1}$ are related by a simple unitary transformation and therefore have the same entropy, see App.A.2 for details. The expressions of these entropic quantities provide an explicit way to compute $H(\hat{A}_0|E)$ as a function of the parameters $\mathbf{L}$ and $\phi$. Let us now turn our attention to the constraints.

**Quantum correlations in the (X,Y) plane** – As mentioned earlier, we are considering quantum models with the values of correlators X and Y given by

Eq. (4). Without loss of generality, we can assume X, Y $\ge$ 0, which can always be attained by relabelling the measurement outcomes of $A_1, B_0$ and $B_1$ (i.e. without touching the angle $\phi$).

In this positive quadrant of the plane, the local strategies are delimited by the CHSH inequality $X + Y \le 2$, i.e. the line connecting the deterministic strategies $(X, Y) = (2, 0)$ and $(X, Y) = (0, 2)$. This implies the following local bounds for the generalized CHSH tests

$$ \frac{1}{2} C_{\Omega} X + \frac{1}{2} S_{\Omega} Y \le B_{\Omega}^{L} = \max(C_{\Omega}, S_{\Omega}). \quad (15) $$

To identify the upper limit of the quantum set, we consider the expected values of the generalized CHSH operator

$$ B_{\Omega} = \left\langle \frac{C_{\Omega}}{2} A_0 \otimes (B_0 + B_1) + \frac{S_{\Omega}}{2} A_1 \otimes (B_0 - B_1) \right\rangle. \quad (16) $$

To find its maximum value, we use the qubit parametrization of measurements $A_y, B_y$ and parametrize the measurement angles on Bob's side as

$$ \begin{aligned} B_0 + B_1 &= (\mathbf{b}_0 + \mathbf{b}_1) \cdot \begin{pmatrix} \sigma_z \\ \sigma_x \end{pmatrix} = 2C_\theta \mathbf{c} \cdot \begin{pmatrix} \sigma_z \\ \sigma_x \end{pmatrix} \\ B_0 - B_1 &= (\mathbf{b}_0 - \mathbf{b}_1) \cdot \begin{pmatrix} \sigma_z \\ \sigma_x \end{pmatrix} = 2S_\theta \mathbf{c}_\perp \cdot \begin{pmatrix} \sigma_z \\ \sigma_x \end{pmatrix} \end{aligned} \quad (17) $$

with two arbitrary perpendicular unit vectors $\mathbf{c}$ and $\mathbf{c}_\perp$, and $\cos(2\theta) = \mathbf{b}_0 \cdot \mathbf{b}_1$. From the diagonalization of the operator on the right hand side of Eq. (16), one finds the quantum bound

$$ B_{\Omega}^{Q} = 1 \quad (18) $$

attained at $(X, Y) = (2\cos(\Omega), 2\sin(\Omega))$ by a maximally entangled two qubit state and measurement settings $\mathbf{a}_0 \cdot \mathbf{a}_1 = 0$ and $\mathbf{b}_0 \cdot \mathbf{b}_1 = \cos(2\Omega)$. It follows that Eve's information is constrained by the part of the quantum set lying between the line $X+Y=2$ and the circle $X^2 + Y^2 = 4$. At this point, we can already conclude that any quantum model with $(X, Y)$ lying on the circle satisfies $H(\hat{A}_0|E) = 1$ (except for the two points with $X+Y=2$), since the underlying state of Alice and Bob has to be pure. This is a straightforward improvement over the CHSH bound.

**Formulation of the problem to solve** – The reductions introduced so far invite us to first solve the

following optimization

$$ \begin{aligned} I(\beta; \Omega, q) &= \max_{\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1} & H(\mathbf{L}) - H(\rho_{E|\hat{a}=+1}) \\ & \text{s.t. } B_{\Omega}(\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1) \ge \beta \end{aligned} \quad (19) $$

and then consider directly the solution $I(\beta; \Omega, q)$ if it is concave in $\beta_\Omega$ or construct a concave function $\hat{I}(\beta; \Omega, q) \ge I(\beta\Omega, q)$ to bound Eve's uncertainty using

$$ H(\hat{A}_0|E) \ge 1 - \min_{\Omega} \hat{I} \left( \frac{C_{\Omega}X + S_{\Omega}Y}{2}; \Omega, q \right). \quad (20) $$

Note that from the symmetries of the goal function $H(\mathbf{L}) - H(\rho_{E|\hat{a}=+1})$ and the constraint $B_\Omega$, we can assume $\phi \in [0, \pi/4]$ for the key generating setting and $L_1 - L_2 \ge L_3 - L_4$ in addition to $L_1 \ge L_2$ and $L_3 \ge L_4$ for the state, see App.A.1 for the details. Further note that we will often use a parametrisation of the tripartite state given by the following 3 component vector

$$ \begin{pmatrix} T_z \\ T_x \\ T_p \end{pmatrix} = \begin{pmatrix} L_1 - L_2 + L_3 - L_4 \\ L_1 - L_2 - L_3 + L_4 \\ L_1 + L_2 - L_3 - L_4 \end{pmatrix}, \quad (21) $$
with $0 \le T_x \le T_z \le 1$ and $T_z + T_x - 1 \le T_p \le 1 - (T_z - T_x)$.

## 3 Bounding Eve's information with generalized CHSH tests

We are now ready to compute a bound on Eve's information as a function of the generalized CHSH score given in Eq. (6) by solving the optimisation problem given in Eq. (19). Among the parameters of the model in Eq. (19), the measurement setting $\mathbf{a}_1, \mathbf{b}_0$ and $\mathbf{b}_1$ only influence the constraint but not the goal function. Furthermore, it is shown in Ref. [13] that $H(\rho_{E|\hat{a}=+1})$ is a monotonic function in the key generating setting $\phi \in [0, \frac{\pi}{4}]$. We can thus decompose the maximization problem in two steps. First, for a fixed state $\mathbf{L}$, we find the lowest angle $\phi$ allowing to satisfy the constraint

$$ \phi_*(\mathbf{L}, \beta, \Omega) = \min \phi \\ \text{s.t.} \quad \max_{\mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1} \mathcal{B}_{\Omega}(\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1) \ge \beta. \quad (22) $$

Second, we fix $\phi = \phi_*(\mathbf{L}, \beta, \Omega)$ to the optimal value for Eve, and maximize her information with respect to the state, that is, we solve

$$ I(\beta; \Omega, q) = \max_{\mathbf{L}} H(\mathbf{L}) - H(\rho_{E|\hat{a}=+1}; \phi_*(\mathbf{L}, \beta, \Omega)). \quad (23) $$

We solve Eq. (22) in App. A.3. The expression of the optimal angle $\phi_*$ depends on whether the

parameter $\Omega$ exceeds $\frac{\pi}{4}$. We treat the two cases $\Omega \le \frac{\pi}{4}$ and $\Omega > \frac{\pi}{4}$ separately.

### 3.1 The simple case with $\Omega \le \frac{\pi}{4}$

For the Bell tests satisfying $\Omega \le \frac{\pi}{4}$, which include the CHSH test, the observed score $\beta$ does not constrain the key generating setting $\phi$ but only the state $\mathbf{L}$. As a result, there always exists a realization with the optimal angle $\phi_* = 0$ as long as the state is such that the Bell score can be attained, i.e. if

$$ C_{\Omega}^{2} T_{z}^{2} + S_{\Omega}^{2} T_{x}^{2} \geq \beta_{\Omega}^{2}. \quad (24) $$

The maximization of the entropy then becomes possible, as for $\phi = 0$ the conditional state $\rho_{E|\hat{a}=+1}$ is block diagonal and its entropy has a simple closed-form expression. Such a maximization has been done for the CHSH case ($\Omega = \frac{\pi}{4}$) in Ref. [13], but Ref. [37] pointed out that the analytical bounds on conditional entropies given in this reference assumes qubit attacks. The same bounds were derived using a different approach in [37], where it is also proved that these bounds are convex. The convexity results of [37], together with Jordan's lemma, imply that the obtained qubit bounds are in fact valid for any dimensions. The same convexity proof applies to the current situation with $\Omega \le \pi/4$. For the sake of completeness, we provide in App. A.5 an alternative proof of convexity, which directly applies to the present case and to [13]. We show in particular (see App. A.4)

$$ \begin{gathered} I(\beta; \Omega, q) = h_q(z) = h(z) - h(n_q(z)) \\ \text{with } n_q(z) = \frac{1 + \sqrt{1 - 4(1 - q)z(1 - z)}}{2} \\ \text{and } z = \frac{1}{2} \left( \frac{\sqrt{\beta^2 - C_{\Omega}^2}}{S_{\Omega}} + 1 \right), \end{gathered} \quad (25) $$

where $h$ is the binary entropy. The concavity of $I(\beta; \Omega, q)$ and hence the convexity of Eve's entropy $H(\hat{A}_0|E)$ follows from (see App. A.5)

$$ \frac{d^2}{d\beta^2} h_q(z(\beta)) = h_q''(z) (z'(\beta))^2 + h_q'(z) z''(\beta) \le 0. \quad (26) $$

Finally, it remains to determine the optimal inequality to use for a given point (X,Y). $h_q(z)$ being a monotonic function of $z$, we derive the expression

$$ z = \frac{1}{2} \left( \frac{\sqrt{\frac{C_{\Omega}X+S_{\Omega}Y}{2}}^2 - C_{\Omega}^2}{S_{\Omega}} + 1 \right), \quad (27) $$

with respect to $\Omega$ to find its local maximum and verify that this value satisfies $\Omega \le \frac{\pi}{4}$. We find that for

$$ \frac{4 - X^2}{XY} < 1: \quad z_{\text{opt}} = \frac{1}{2} \left( \frac{Y}{\sqrt{4 - X^2}} + 1 \right) \quad (28) $$
associated to the generalized CHSH test with

$$ \tan(\Omega) = \frac{4 - X^2}{XY}. \qquad (29) $$

For

$$ \frac{4 - X^2}{XY} \geq 1: \quad z_{\text{opt}} = \frac{1}{2} \left( 1 + \sqrt{\left(\frac{X+Y}{2}\right)^2 - 1} \right), \qquad (30) $$

which is the formula for the CHSH case $\Omega = \frac{\pi}{4}$ [13].

## 3.2 The complicated case with $\Omega > \frac{\pi}{4}$

### 3.2.1 Bounding Eve's information from a numerical optimization

For the remaining Bell tests, with $\Omega > \frac{\pi}{4}$, the situation is different. Here, the generalized CHSH score $\beta$ does not only constrain the state $\mathbf{L}$ but also the setting of the key generating measurement. In particular, we find that there are two different regions.

First, for the states falling in the range

$$ \mathcal{S}(\beta, \Omega) = \{\mathbf{L}|(C_{\Omega}^2 T_z^2 + S_{\Omega}^2 T_x^2) \le \beta^2 \le (C_{\Omega}^2 T_x^2 + S_{\Omega}^2 T_z^2)\}, \qquad (31) $$

the constraint $B_{\Omega} \ge \beta$ can be satisfied with the measurement angle $\phi \ge \phi_{*}(\mathbf{L}, \Omega, \beta)$, where $\cos^2(\phi_{*}) = c_{*}^2(\mathbf{L}, \Omega, \beta)$ is given by

$$ c_{*}^{2}(\mathbf{L}, \Omega, \beta) = \frac{(\beta^{2} - S_{\Omega}^{2} T_{x}^{2})(C_{\Omega}^{2} T_{x}^{2} + S_{\Omega}^{2} T_{z}^{2} - \beta^{2})}{C_{\Omega}^{2}(T_{z}^{2} - T_{x}^{2})(S_{\Omega}^{2} T_{z}^{2} + S_{\Omega}^{2} T_{x}^{2} - \beta^{2})}. \qquad (32) $$

The constraint on the angle only becomes trivial, $c_{*}^2 = 1$, on the boundary of the region $\mathcal{S}$, where $(C_{\Omega}^2 T_z^2 + S_{\Omega}^2 T_x^2) = \beta^2$.

Second, in the region where $(C_{\Omega}^2 T_z^2 + S_{\Omega}^2 T_x^2) > \beta^2$ the Bell score $\beta$ can also be attained with $\phi = 0$. However, there models provide less information to Eve as compared to these on the boundary $(C_{\Omega}^2 T_z^2 + S_{\Omega}^2 T_x^2) = \beta^2$, see the discussion at the end of App. A.4. So we can safely ignore this region.

To find the best strategy for Eve, it thus remains to solve

$$ I(\beta; \Omega, q) = \max_{\mathbf{L} \in \mathcal{S}(\Omega, \beta)} H(\mathbf{L}) - H(E|\hat{a} = +1; c_{*}^{2}(\mathbf{L}, \Omega, \beta)). \qquad (33) $$

This optimization only involves an analytic function of three parameters on a compact domain. It can be easily and time-efficiently solved by standard numerical methods. We now give an ansatz on the solution of the optimization given in Eq. (33), which allows one to speed up its numerical resolution even further.

### 3.2.2 Ansatz

First, we observe that the vector $\mathbf{L}$ saturating Eve's information only has two non-zero coefficients $L_1 = 1 - L_3$ and $L_2 = L_4 = 0$, or $T_z = 1$ and $T_x = T_p$. With this observation, the previous optimization problem becomes a scalar optimization, that is

$$ \tilde{I}_{\text{anz}}(\beta; \Omega, q) = \max_{\frac{\beta^2 - S_{\Omega}^2}{C_{\Omega}^2} \le T_x^2 \le \frac{\beta^2 - C_{\Omega}^2}{S_{\Omega}^2}} h\left(\frac{1 + T_x}{2}\right) - h\left(\frac{1 + \sqrt{T_x^2 + C_*^2 q(1 - T_x^2)}}{2}\right). \qquad (34) $$

Second, we see that the bound $\tilde{I}_{\text{anz}}(\beta; \Omega, q)$ is not concave for small $\beta$. However, we observe that its concave roof can be obtained by drawing a line which passes through the point

$$ \beta = B_{\Omega}^{L} = \sin(\Omega) \\ I_{L}(q) = I(B_{\Omega}^{L}; \Omega, q) = 1 - h\left(\frac{1 + \sqrt{q}}{2}\right) \qquad (35) $$
which is tangent to the curve $\tilde{I}_{\text{anz}}(\beta; \Omega, q)$. The value of the generalized CHSH score at the tangent point can be found by solving

$$ \beta_* = \operatorname{argmin}_{\beta} \frac{I_L(q) - \tilde{I}_{\text{anz}}(\beta; \Omega, q)}{\beta - \sin(\Omega)}. \qquad (36) $$

Labeling $I_*(\Omega, q) = \tilde{I}_{\text{anz}}(\beta_*; \Omega, q)$, this leads to

$$ \hat{I}_{\text{anz}}(\beta; \Omega, q) = \begin{cases} \frac{I_*(\Omega, q)(\beta - \sin(\Omega)) - I_L(q)(\beta - \beta_*)}{\beta_* - \sin(\Omega)} & \beta < \beta_* \\ \tilde{I}_{\text{anz}}(\beta; \Omega, q) & \beta \ge \beta_* \end{cases} \qquad (37) $$

At this stage, we further observe that the optimal value of $T_x^2$ for values of $\beta \ge \beta_*$ coincides with its maximum possible value $T_x^2 = \frac{\beta^2 - C_\Omega^2}{S_\Omega^2}$, which defines $I_{\text{anz}}(\beta; \Omega, q)$. In the Eqs. (36) and (37) we can now replace $\tilde{I}_{\text{anz}}(\beta; \Omega, q)$ with $I_{\text{anz}}(\beta; \Omega, q)$, which does not involve any nonlinear optimization.

While we believe this expression to be the true bound, we do not have a formal proof. Anyway, this conjectured expression helps to solve the optimization of interest.

### 3.2.3 Certified numerical solution

As mentioned before, the optimization in Eq. (33) can be easily solved by standard numerical methods. However, to provide a strict security guarantee to an actual implementation of DIQKD, such a numerical optimization would need to be done in a certified manner, with a formal proof that the obtained numbers lower bound Eve's conditional entropy on the whole domain. Below we present an algorithm which allows one to do such a certified optimisation based on the Lipshitz continuity of the goal function.

The algorithm is rather time-costly, but it only has to be used once the optimal experimental parameters are fixed through an ad hoc maximization of Eq. (33), cf. below.

Concretely, we present in this section an algorithm that approximates the set of possible strategies of Eve, delimited by the bound $\hat{I}(\beta, \Omega, q)$, from the outside. To avoid the issue of concavity posed by Eq. (19), we rewrite the problem in the dual form in which we look for the tangent lines

$$ \begin{gathered} f(t; \Omega, q) = \max_{\mathbf{L}, \phi} H(\rho_E) - H(\rho_E | \hat{\mathbf{A}}_0) + t \beta_{\max}(\mathbf{L}, \phi; \Omega) \\ \beta_{\max}(\mathbf{L}, \phi; \Omega) = \max_{\mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1} \mathcal{B}_{\Omega}(\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1), \end{gathered} \qquad (38) $$

to the curve $\hat{I}(\beta, \Omega, q)$ with different slopes $t$. In Eq. (38) we used the fact that it is only the Bell score that depends on the measurement setting $\mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1$, so it can be maximized straightforwardly to define $\beta_{\max}(\mathbf{L}, \phi; \Omega)$, see App. A.6 for its closed form expression.

Before giving the details on the way we solve this dual form, let us shortly discuss on how it shall be used. We consider an actual implementation of DIQKD with fixed values of $X_*, Y_*$ and $q_*$ and for which there is an optimal value $\Omega_*$ which saturates the minimum in Eq. (20)

$$ H(\hat{\mathbf{A}}_0 | E) \ge 1 - \hat{I}(\beta_*; \Omega_*, q_*), \qquad (39) $$

with $\beta_* = \frac{1}{2}(C_{\Omega_*} X_* + S_{\Omega_*} Y_*)$. Therefore, an optimal security guarantee for this particular implementation only requires the knowledge of the function $\hat{I}(\beta_*, \Omega_*, q_*)$ on a single point. The same lower bond on Eve's conditional entropy can be obtained from the value of the dual bound $f(t; \Omega_*, q_*)$ in Eq. (38) on a single point. Indeed, the concavity of $\hat{I}(\beta; \Omega, q)$

ensures that there exists a value $t_*$ for which the inequality

$$ f(t_*; \Omega_*, q_*) \ge \hat{I}(\beta; \Omega_*, q_*) + t_* \beta \qquad (40) $$

is saturated at $\beta = \beta_*$ where $f(t_*; \Omega_*, q_*) = \hat{I}(\beta; \Omega_*, q_*) + t_* \beta_*$. Using Eq. (38), we deduce that

$$ \begin{aligned} H(\hat{\mathbf{A}}_0 | E) &\ge 1 - f(t_*; \Omega_*, q_*) + t_* \beta_* \\ &= 1 - \hat{I}(\beta_*; \Omega_*, q_*). \end{aligned} \qquad (41) $$

Hence, for a fixed experimental implementation of the protocol, it is sufficient to certify a single value of the function $f(t_*; \Omega_*, q_*)$ in order to provide a strict and optimal security guarantee. Furthermore, the value $t_*$ is straightforward to find from the knowledge of $\Omega_*, \beta_*, q_*$ and the function $\hat{I}(\beta; \Omega, q)$.

We can now comment on the algorithm to provably upper bound the quantity in Eq. (38). While the algorithm used recently for the derivation of key rates in Ref. [28] could be easily adjusted to solve this problem, we propose a branch and bound approach.
Concretely, we first derive a parametrization of the set $(\mathbf{L}, \phi)$ for which the goal function in Eq. (38) is Lipshitz continuous with a constant that we compute. Then, we obtain an upper bound on its value on the whole domain by computing the value of the function on a grid of points. Finally, the algorithm subsequently refines the grid around the points where the value of the function is large in order to approach, step by step, the global maximum. Additional details are given in the next three sections.

**Lipshitz continuity of the entropy with respect to the angle** – A key ingredient in the implementation of the desired certified algorithm is that the von Neumann entropy $H(\rho)$ has a bounded gradient with respect to the following angle between the two states $\rho$ and $\sigma$ (see App. B.1)

$$A(\rho, \sigma) = \arccos(F(\rho, \sigma)), \quad (42)$$

where the fidelity is defined as $F(\rho, \sigma) = \text{tr}|\sqrt{\rho}\sqrt{\sigma}|$. This angle is a metric on the set of states [22]. More precisely, we show that for $n$-dimensional quantum states, the gradient of the entropy satisfies

$$\left|\frac{dH}{dA}\right| < \begin{cases} 4.023 & n=4 \\ 2 \log(n) & n \ge 5. \end{cases} \quad (43)$$

This contrasts with the trace distance, another metric on the set of quantum states, for which the entropy has an infinite slope around non-full-rank states.

**Bounding the gradient of the goal function**– To apply the continuity bound above to the goal function in Eq. (38), we use the following parametrization of the state

$$\sqrt{\mathbf{L}} = \begin{pmatrix} \cos(\alpha) \cos(\mu) \\ \cos(\alpha) \sin(\mu) \\ \sin(\alpha) \cos(\xi) \\ \sin(\alpha) \sin(\xi) \end{pmatrix}, \quad (44)$$

so that the quantum models are described by four angles

$$\mathbf{x} = (\alpha, \mu, \xi, \phi), \quad (45)$$

with $\mathbf{x} \in [0, \frac{\pi}{4}]^3 \times [0, \frac{\pi}{2}]$. Here, the condition $\mu, \xi \le \frac{\pi}{4}$ follow from $L_1 \ge L_2, L_3 \ge L_4$. The bound $\alpha \le \frac{\pi}{4}$ is a consequence of $L_1 + L_2 \ge L_3 + L_4$, which can be imposed on the states as an alternative to $L_1 - L_2 \ge L_3 - L_4$, see App. A.1.

To obtain a bound on the gradient of the goal function

$$G(\mathbf{x}) = H(\mathbf{L}) - H(\rho_E | \hat{a} = +1) + t \beta_{\max}(\mathbf{L}, \phi; \Omega) \quad (46)$$

we bound the gradient of each term independently, as described in App. B.2 and B.3. For the entropic

terms, we use Eq. (43), while the computation of the maximal gradient is quite straightforward for the CHSH score. Combining the three terms, we obtain a general bound

$$|\nabla_{\mathbf{x}} G| \le 12.7 + 7t \quad (47)$$

on the whole domain of $\mathbf{x}$.

With this global bound on the gradient, a certified maximization of the function $G(\mathbf{x})$ can be obtained with a slight modification of the code that we developed in Ref. [33], which is available on GitLab. A detailed description of this code can be found in the App. B.4.

## 4 Results

**Improved bound on Eve's conditional entropy**– In order to demonstrate the advantage of considering the pair of variables $(X, Y)$ when bounding Eve's conditional entropy, we compute the bound on $H(\hat{A}_0|E)$ for different values of $X, Y$ and $p$ and compare it to the bound obtained from the CHSH score [13]. Because of the two regimes identified earlier, we compute both the optimal bound assuming $\Omega \le \pi/4$ and the one assuming $\Omega > \pi/4$ for each value of $X$ and $Y$, and keep the best one. In the case $\Omega \le \pi/4$, the optimal choice of $\Omega$ is readily given by Eq. (29). In the other case, we optimize the bound over $\Omega \in (\pi/4, \pi/2]$. The difference between this optimal bound on $H(\hat{A}_0|E)$ given $X$ and $Y$ and Eq. (3) is shown in Fig. 1. It shows that our bound on $H(\hat{A}_0|E)$ is better than the one derived from the CHSH score, except along the line satisfying $X(X + Y) = 4$.

**Implication for a practical realization of DIQKD**– We now study the potential impact of our bound on practical realizations of DIQKD. In the limit of asymptotically many repetitions, an implementation is uniquely characterized by its key rate $r$, which is given in Eq. (2). In our case, this key rate is determined from three quantities: $X$, $Y$ and Bob's uncertainty about Alice's key generating bit as a function of the noisification parameter, given by $H(\hat{A}_0|B_2)$. Hence, in order to find the optimal design for an experimental implementation of DIQKD, we express the quantities

$$\text{setup} \approx (X, Y, H(\hat{A}_0|B_2)) \quad (48)$$

as a function of the model's parameters. Then, we maximize the key rate over these parameters:

$$\begin{aligned} \tilde{r} &= \max_{\text{setup}, \Omega, q} H(\hat{A}_0|E) - H(\hat{A}_0|B_2) \\ &= \max_{\text{setup}, \Omega, q} \left(1 - \hat{I}\left(\frac{C_\Omega X + S_\Omega Y}{2}; \Omega, q\right) - H(\hat{A}_0|B_2)\right) \end{aligned} \quad (49)$$
<figure><img src="image_2.png" /></figure>

Figure 1: Difference between the bound on Eve's conditional entropy $H(\hat{A}_0|E)$ computed as a function of X and Y, and as a function of the CHSH quantity $S = 2\sqrt{2}(X + Y)$, for $p = 0$. In presence of noisification (i.e. $p > 0$), the advantage follows a similar distribution, but is smaller in magnitude. The CHSH bound is only optimal along the curve $X(X+Y) = 4$. The advantage on the right-hand side of this curve is obtained with $\Omega < \pi/4$, and on the left-hand side with $\Omega > \pi/4$. The yellow (red) curve shows the trajectory in the X-Y plane which optimizes the key rate in an optical implementation of DIQKD for low (high) detection efficiencies. At the efficiency $\eta = 0.923$, it is better to switch from one curve to the other one (at the transition between full and dashed curves).

Solving this maximization gives a bound on the key rate $\bar{r}$, the values $(X_*, Y_*, H(\hat{A}_0|B_2))$ expected for a given implementation, as well as the optimal values of the parameters $\Omega_*$ and $q_*$ for this implementation.

## SPDC-based implementation of DIQKD

Photonic experiments using a source based on spontaneous parametric down conversion (SPDC) are one of the most promising setups for implementing DIQKD, as shown by recent experiments reporting on the violation of a Bell inequality without the fair sampling assumption [4, 11, 12, 16, 29, 30]. We consider such a setup in which an SPDC source is used to create and distribute polarization entanglement between distant parties who perform measurements as requested here in the proposed protocol. The main limitation in this setup is the overall detection efficiency, i.e. the possibility of losing a photon at any point between its creation at the source and its final detection. To reflect photon losses and non-unit detection efficiency, the transmission channel between the source and the parties is modeled as a

lossy channel with an overall transmission $\eta$. It is also important to include the statistics of an SPDC source which does not produce a two-qubit state, but a state that contains vacuum and multiple photon components. We invite the reader to look at Ref. [13] to get explicit expressions of the exact statistics created by this source as well as a description of tunable parameters.

When computing the key rate for an SPDC source with a security determined by the CHSH score, any values of X and Y with the same sum impose the same bound on Eve's conditional entropy $H(\hat{A}_0|E)$. In a Eberhard-like scenario where a significant fraction of the entangled particles can be lost before yielding their measurement result, it is advantageous for Alice to use two measurement settings with different overlap with her Schmidt basis in order to maximize the CHSH quantity [7]. It is then easier for Bob to guess the outcome of one of the two measurements (the one best aligned with the Schmidt basis). When the key rate is extracted from this measurement to minimize the cost of error correction (see Ref. [13] for more details), the value of the X quantity is then larger than Y, as shown in Fig. 1. But in this case, the (X,Y) values are very close to the line $X(X + Y) = 4$, for which there is no advantage. Therefore, we only expect a small improvement in the key rate here.

Given the values shown in Fig. 1, a better bound on Eve's information would be obtained if the same CHSH value was obtained with the contributions from X and Y being inverted, i.e. with $Y > X$. However, this requires Alice to define her key-generating measurement as being the one less aligned with the Schmidt basis, hence leading to an increase of the error correction cost. In this case, both conditional entropies $H(\hat{A}_0|E)$ and $H(\hat{A}_0|B_2)$ are larger, and we need to check which one increases the most in order to infer a possible gain on the key rate. As it turns out, the tradeoff between these two entropies increase depends on the detection efficiency.

Namely, there are two regimes, as shown in Fig. 2. When the detection efficiency is larger than ~0.923, relabelling the measurements in order to enter the region with $Y > X$ (where inequalities with $\Omega > \pi/4$ significantly improve the bound on $H(\hat{A}_0|E)$) is advantageous, because the increase in Bob's uncertainty $H(\hat{A}_0|B_2)$ is smaller. The red curve in Fig. 1 shows the corresponding trajectory in the X-Y plane. When $\eta < 0.923$, the cost of error correction become prohibitive compared to the potential increase in Eve's uncertainty, and it is better to stay in the region $X > Y$, as represented by the yellow curve in Fig. 1. There, a small increase of the key rate is still found because the correlations do not satisfy
<figure><img src="image_3.png" /><figcaption>Figure 2: Key rate achievable with a photonic setup as a function of the symmetric detection efficiency $\eta$. The dashed blue curve corresponds to the key rate achieved with noisification and a security based on the CHSH quantity alone; it is nonzero for efficiencies above ~0.828 [13]. The continuous red curve bases its security on the values of both X and Y (and also includes noisification). A significant increase of the key rate is possible for efficiencies above $\eta > 0.923$. At $\eta = 1$, the obtained key rate is 0.346 instead of 0.252.</figcaption></figure>

$X(X+Y) = 4$ exactly. However, this condition is only slightly violated, resulting in an increase in key rate smaller than ~$10^{-4}$, which is practically negligible. The critical detection efficiency then also remains at 82.8%, unchanged compared to a bound based on CHSH alone.

**Comparison of qubit vs SPDC bounds**— To illustrate the impact of the photon statistics of SPDC sources on DIQKD, we now consider a simpler model in which the state shared between Alice and Bob is a two-qubit state:

$$|\psi\rangle = \cos(\theta)|00\rangle + \sin(\theta)|11\rangle. \qquad (50)$$

We are not aware of physical setups allowing one to produce a state with $\theta = \pi/4$ but not allowing for a different value of $\theta$. Still, for the sake of the discussion, we distinguish between the cases where the state can be either constrained to be maximally entangled, i.e. $\theta = \pi/4$, or can have an arbitrary parameter $\theta \in [0, \pi/4]$.

In Fig. 3, we report the critical detection efficiencies corresponding to various security analyses applied on these implementations. Like for the SPDC model, no advantage on the critical detection efficiency is found when using arbitrary two-qubit systems. A small advantage is however present when

<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of the critical detection efficiencies for several setup models. Blue markers are for SPDC statistics, yellow ones for measurement on a maximally entangled two-qubit state, and red ones for measurement on an arbitrary two-qubit state. (a) is with security from CHSH following Eq. (3) and error correction based on the QBER [23], (b) Eq. (3) with error correction based on conditional entropy $H(A|B)$ [19], (c) security from CHSH with noisification and error correction based on $H(A|B)$ following [13], (d) security from X and Y with noisification and error correction based on $H(A|B)$.</figcaption></figure>

restricting to measurements on the singlet state. Still, this model is not optimal and it remains of course better to use partially entangled states. In fact, even partially entangled states produced by an SPDC source perform better.

In this respect, it is worth noticing here that the performance of an SPDC source is essentially comparable to that of an arbitrary two-qubit state once noisification is taken into consideration, i.e. the requirement on the detection efficiency is very similar. In the case with no noisification, the state produced by these physical sources do not have a better tolerance to losses than measurements on a maximally entangled state. Noisy preprocessing is thus a key ingredient for a first proof of principle implementation with an SPDC source [13].

# 5 Discussion

In this paper, we introduced a refinement of the usual CHSH-based analysis of DIQKD experiments: Instead of projecting the measurement statistics onto a single line giving the CHSH score $2\sqrt{2}(X + Y)$, we kept the information about the individual values of X and Y throughout the whole security analysis. We found that this refined analysis gives a more restrictive bound on the information available to the eavesdropper for almost all values X and Y in the quantum set.
When applying our results to photonic implemen-
tations of DIQKD with a SPDC source, we found
that the key rate is improved by a factor going up
to 37% for unit detection efficiency. On the other
hand, we could not find any improvement for the
critical detection efficiency as compared to the CHSH
protocol with noisy processing presented in Ref. [13].
However, we focused on a given photonic implemen-
tation, and the question of the most favorable optical
setup combining squeezing operations, displacement
operations, linear optical elements and photon count-
ing techniques is still open. Advanced techniques
using automated design of quantum experiments
based on reinforcement learning which already proved
to be useful to optimize the CHSH score [20] are
inspiring. Applying them to the proposed protocol in
order to reduce the required detection efficiency for
implementing DIQKD appears to be promising for
future work.

Finally, we would like to remark that the certi-
fied numerical techniques we proposed also open up
the possibility of bounding Eve's information reliably
when more correlators, or even the full measurement
statistics, are taken into account.

6 Note added

While writing this manuscript, we became aware of
another manuscript [37] reporting on similar results.

7 Acknowledgments

We thank Stefano Pironio for pointing out that the
convexity argument provided in [13] was not com-
plete, see the discussion in Sec. 3.1 and App. A.5. We
acknowledge funding by the Swiss National Science
Foundation (SNSF), through the Grants PP00P2-
179109 as well as via the National Center for Compe-
tence in Research for Quantum Science and Technol-
ogy (QSIT).

References

[1] Antonio Acín, Nicolas Gisin, and Lluis Masanes. From Bell's theorem to secure quantum key distribution. *Phys. Rev. Lett.*, 97:120405, Sep 2006. DOI: 10.1103/PhysRevLett.97.120405.

[2] Rotem Arnon-Friedman, Frédéric Dupuis, Omar Fawzi, Renato Renner, and Thomas Vidick. Practical device-independent quantum cryptography via entropy accumulation. *Nature Communications*, 9(1):459, 2018. DOI: 10.1038/s41467-017-02307-4.

[3] Samuel L. Braunstein and Stefano Pirandola. Side-channel-free quantum key distribu-

tion. *Phys. Rev. Lett.*, 108:130502, Mar 2012.
DOI: 10.1103/PhysRevLett.108.130502.

[4] BG Christensen, KT McCusker, JB Altepeter, B Calkins, Thomas Gerrits, Adriana E Lita, A Miller, Lynden K Shalm, Y Zhang, SW Nam, et al. Detection-loophole-free test of quantum nonlocality, and applications. *Physical review letters*, 111(13):130406, 2013. DOI: 10.1103/PhysRevLett.111.130406.

[5] L. Comandar, M. Lucamarini, and B. et al. Fröhlich. Quantum key distribution without detector vulnerabilities using optically seeded lasers. *Nature Photon*, 10:312–315, 2016. DOI: 10.1038/nphoton.2016.50.

[6] Igor Devetak and Andreas Winter. Distillation of secret key and entanglement from quantum states. *Proc. R. Soc. Lond. A*, 461:207–235, 2005. DOI: 10.1098/rspa.2004.1372.

[7] Philippe H. Eberhard. Background level and counter efficiencies required for a loophole-free einstein-podolsky-rosen experiment. *Phys. Rev. A*, 47:R747–R750, Feb 1993. DOI: 10.1103/PhysRevA.47.R747.

[8] Artur K. Ekert. Quantum cryptography based on Bell's theorem. *Phys. Rev. Lett.*, 67:661–663, Aug 1991. DOI: 10.1103/PhysRevLett.67.661.

[9] T. Ferreira da Silva, D. Vitoreti, G. B. Xavier, G. C. do Amaral, G. P. Temporão, and J. P. von der Weid. Proof-of-principle demonstration of measurement-device-independent quantum key distribution using polarization qubits. *Phys. Rev. A*, 88:052303, Nov 2013. DOI: 10.1103/PhysRevA.88.052303.

[10] Juan Carlos Garcia-Escartin, Shihan Sajeed, and Vadim Makarov. Attacking quantum key distribution by light injection via ventilation openings. preprint arXiv:1910.08152, 2019.

[11] Marissa Giustina, Alexandra Mech, Sven Ramelow, Bernhard Wittmann, Johannes Kofler, Jörn Beyer, Adriana Lita, Brice Calkins, Thomas Gerrits, Sae Woo Nam, et al. Bell violation using entangled photons without the fair-sampling assumption. *Nature*, 497(7448):227, 2013. DOI: 10.1038/nature12012.

[12] Marissa Giustina, Marijn AM Versteegh, Sören Wengerowsky, Johannes Handsteiner, Armin Hochrainer, Kevin Phelan, Fabian Steinlechner, Johannes Kofler, Jan-Åke Larsson, Carlos Abellán, et al. Significant-loophole-free test of Bell's theorem with entangled photons. *Physical review letters*, 115(25):250401, 2015. DOI: 10.1103/PhysRevLett.115.250401.

[13] M. Ho, P. Sekatski, E. Y.-Z. Tan, R. Renner, J.-D. Bancal, and N. Sangouard. Noisy preprocessing facilitates a photonic realization of device-independent quantum key distribution. *Phys. Rev. Lett.*, 124:230502, Jun 2020. DOI: 10.1103/PhysRevLett.124.230502.
[14] B. Kraus, N. Gisin, and R. Renner. Lower and upper bounds on the secret-key rate for quantum key distribution protocols using one-way classical communication. *Phys. Rev. Lett.*, 95:080501, Aug 2005. DOI: 10.1103/PhysRevLett.95.080501.

[15] Yang Liu, Teng-Yun Chen, Liu-Jun Wang, Hao Liang, Guo-Liang Shentu, Jian Wang, Ke Cui, Hua-Lei Yin, Nai-Le Liu, Li Li, Xiongfeng Ma, Jason S. Pelc, M. M. Fejer, Cheng-Zhi Peng, Qiang Zhang, and Jian-Wei Pan. Experimental measurement-device-independent quantum key distribution. *Phys. Rev. Lett.*, 111:130502, Sep 2013. DOI: 10.1103/PhysRevLett.111.130502.

[16] Yang Liu, Qi Zhao, Ming-Han Li, Jian-Yu Guan, Yanbao Zhang, Bing Bai, Weijun Zhang, Wen-Zhao Liu, Cheng Wu, Xiao Yuan, et al. Device-independent quantum random-number generation. *Nature*, 562(7728):548, 2018. DOI: 10.1038/s41586-018-0559-3.

[17] Hoi-Kwong Lo, Marcos Curty, and Bing Qi. Measurement-device-independent quantum key distribution. *Phys. Rev. Lett.*, 108:130503, Mar 2012. DOI: 10.1103/PhysRevLett.108.130503.

[18] Lars Lydersen, Carlos Wiechers, Christoffer Wittmann, Dominique Elser, Johannes Skaar, and Vadim Makarov. Hacking commercial quantum cryptography systems by tailored bright illumination. *Nature Photonics*, 4:686–689, 2010. DOI: 10.1038/nphoton.2010.214.

[19] Xiongfeng Ma and Norbert Lütkenhaus. Improved data post-processing in quantum key distribution and application to loss thresholds in device independent qkd. *Quantum Inf. Comput.*, 12:203–214, 2012. DOI: https://doi.org/10.26421/QIC12.3-4.

[20] Alexey A Melnikov, Pavel Sekatski, and Nicolas Sangouard. Setting up experimental Bell test with reinforcement learning. preprint arXiv:2005.01697, 2020.

[21] L. Mirsky. A trace inequality of John von Neumann. *Monatshefte für Mathematik*, 79(4):303–306, December 1975. DOI: 10.1007/bf01647331.

[22] Michael A. Nielsen and Isaac L. Chuang. *Quantum Computation and Quantum Information: 10th Anniversary Edition*. Cambridge University Press, USA, 10th edition, 2011. ISBN 1107002176.

[23] Stefano Pironio, Antonio Acin, Nicolas Brunner, Nicolas Gisin, Serge Massar, and Valerio Scarani. Device-independent quantum key distribution secure against collective attacks. *New Journal of Physics*, 11(4):045021, 2009. DOI: 10.1088/1367-2630/11/4/045021.

[24] Joseph M. Renes and Graeme Smith. Noisy processing and distillation of private quantum states. *Phys. Rev. Lett.*, 98:020502, Jan 2007. DOI: 10.1103/PhysRevLett.98.020502.

[25] Renato Renner, Nicolas Gisin, and Barbara Kraus. Information-theoretic security proof for quantum-key-distribution protocols. *Physical Review A*, 72(1):012332, 2005. DOI: 10.1103/PhysRevA.72.012332.

[26] A. Rubenok, J. A. Slater, P. Chan, I. Lucio-Martinez, and W. Tittel. Real-world two-photon interference and proof-of-principle quantum key distribution immune to detector attacks. *Phys. Rev. Lett.*, 111:130501, Sep 2013. DOI: 10.1103/PhysRevLett.111.130501.

[27] Rene Schwonnek, Koon Tong Goh, Ignatius W. Primaatmaja, Ernest Y.-Z. Tan, Ramona Wolf, Valerio Scarani, and Charles C.-W. Lim. Robust device-independent quantum key distribution. preprint arXiv:2005.02691, 2020.

[28] Rene Schwonnek, Koon Tong Goh, Ignatius W Primaatmaja, Ernest Y-Z Tan, Ramona Wolf, Valerio Scarani, and Charles C-W Lim. Robust device-independent quantum key distribution. preprint arXiv:2005.02691, 2020.

[29] Lynden K Shalm, Evan Meyer-Scott, Bradley G Christensen, Peter Bierhorst, Michael A Wayne, Martin J Stevens, Thomas Gerrits, Scott Glancy, Deny R Hamel, Michael S Allman, et al. Strong loophole-free test of local realism. *Physical review letters*, 115(25):250402, 2015. DOI: 10.1103/PhysRevLett.115.250402.

[30] Lijiong Shen, Jianwei Lee, Jean-Daniel Bancal, Alessandro Cerè, Antia Lamas-Linares, Adriana Lita, Thomas Gerrits, Sae Woo Nam, Valerio Scarani, Christian Kurtsiefer, et al. Randomness extraction from bell violation with continuous parametric down-conversion. *Physical review letters*, 121(15):150402, 2018. DOI: 10.1103/PhysRevLett.121.150402.

[31] Ivan Šupić and Joseph Bowles. Self-testing of quantum systems: a review. preprint arXiv:1904.10042, 2019.

[32] Zhiyuan Tang, Zhongfa Liao, Feihu Xu, Bing Qi, Li Qian, and Hoi-Kwong Lo. Experimental demonstration of polarization encoding measurement-device-independent quantum key distribution. *Phys. Rev. Lett.*, 112:190503, May 2014. DOI: 10.1103/PhysRevLett.112.190503.

[33] Xavier Valcarce, Pavel Sekatski, Davide Orsucci, Enky Oudot, Jean-Daniel Bancal, and Nicolas Sangouard. What is the minimum CHSH score certifying that a state resembles the singlet? *Quantum*, 4:246, March 2020. ISSN 2521-327X. DOI: 10.22331/q-2020-03-23-246.

[34] Xavier Valcarce, Julian Zivy, Nicolas Sangouard, and Pavel Sekatski. Self-testing two-qubit maximally entangled states with two binary measurements. To be submitted, 2020.

[35] Henning Weier, Harald Krauss, Markus Rau, Martin Fürst, Sebastian Nauerth, and Harald Weinfurter. Quantum eavesdropping without interception: an attack exploiting the dead time
of single-photon detectors. New Journal of Physics, 13(7):073024, 2011. DOI: 10.1088/1367-2630/13/7/073024.

[36] Erik Woodhead. Tight asymptotic key rate for the Bennett-Brassard 1984 protocol with local randomization and device imprecisions. Phys. Rev. A, 90:022306, Aug 2014. DOI: 10.1103/PhysRevA.90.022306.

[37] Erik Woodhead, Antonio Acín, and Stefano Piro-nio. Device-independent quantum key distribu-tion based on asymmetric CHSH inequalities. preprint arXiv:2007.16146, 2020.

[38] Yi Zhao, Chi-Hang Fred Fung, Bing Qi, Christine Chen, and Hoi-Kwong Lo. Quantum hacking: Experimental demonstration of time-shift attack against practical quantum-key-distribution systems. Phys. Rev. A, 78:042333, Oct 2008. DOI: 10.1103/PhysRevA.78.042333.

# A Analytical results

## A.1 Parametrization of two-qubit models

Following the logic of [23], we assume, without loss of generality that after the qubit reduction, the state shared by Alice, Bob and Eve is of the form

$$|\Psi\rangle_{ABE} = \sum_{i=1}^{4} \sqrt{L_i} |\Phi^i\rangle_{AB} |i\rangle_E \quad (51)$$

where we $|\Phi_i\rangle = \{\Phi^+, \Psi^-, \Phi^-, \Psi^+\}_{i=1}^4$, with nonneg-ative $L_1 \ge L_2$ and $L_3 \ge L_4$, and the measurements $A_x, B_y$ appearing in the constraint are in the x-z plane

$$A_x = \mathbf{a}_x^T \begin{pmatrix} \sigma_z \\ \sigma_x \end{pmatrix} \qquad B_y = \mathbf{b}_y^T \begin{pmatrix} \sigma_z \\ \sigma_x \end{pmatrix}. \quad (52)$$

The key generating setting is explicitly parametrized by an angle $\phi$

$$\mathbf{a}_0 = \begin{pmatrix} C_\phi \\ S_\phi \end{pmatrix}. \quad (53)$$

(As mentioned in the main text we use a compact notation $C_\phi = \cos(\phi)$ and $S_\phi = \sin(\phi)$ throughout the

paper.) One notes that the application of a unitary transformation $\sigma_Z$ on Alice's system is equivalent to changing the state and the measurements as

$$\begin{align} (L_1, L_2, L_3, L_4) &\rightarrow (L_3, L_4, L_1, L_2) \nonumber \\ \mathbf{a}_x &\rightarrow \sigma_z \mathbf{a}_x \tag{54} \\ \mathbf{b}_y &\rightarrow \mathbf{b}_y. \nonumber \end{align}$$

The two situations are completely equivalent for our purpose, so our parametrization of quantum models is actually redundant. To avoid this, we can introduce an order relation between the pairs of coefficients $(L_1, L_2)$ and $(L_3, L_4)$ (permuted by a basis transformation). In particular, we will impose

$$L_1 - L_2 \ge L_3 - L_4 \quad (55)$$

below, but will also use $L_1 + L_2 \ge L_3 + L_4$ for the certified numerical algorithm.

We also introduce a different parametrization of the probability simplex $\mathbf{L}$ with a vector $\mathbf{T} = (T_z, T_x, T_p)$ given by

$$\left\{ \begin{aligned} T_z &= (L_1 - L_2) + (L_3 - L_4) \\ T_x &= (L_1 - L_2) - (L_3 - L_4) \\ T_p &= L_1 + L_2 - L_3 - L_4. \end{aligned} \right. \quad (56)$$

The conditions $L_1 \ge L_2, L_3 \ge L_4$ and $L_1 - L_2 \ge L_3 - L_4$ enforce

$$\begin{gathered} 0 \le T_x \le T_z \le 1 \\ T_z + T_x - 1 \le T_p \le 1 - (T_z - T_x) \end{gathered} \quad (57)$$

## A.2 Entropies of $\rho_{E|\hat{a}=\pm 1}$

Having introduced a parametrization of quantum models, we now express the quantities of interest $H(\rho_E)$, $H(\rho_E|\hat{A}_0)$ and $B_\Omega$ as functions of the distribution $\mathbf{L}$ describing the state $|\Psi\rangle_{ABE}$ and the measurement settings $\mathbf{a}_x, \mathbf{b}_y$. The marginal state of Eve is straightforward to write down:

$$\rho_E = \begin{pmatrix} L_1 & & & \\ & L_2 & & \\ & & L_3 & \\ & & & L_4 \end{pmatrix}. \quad (58)$$

For the conditional states, we have

$$\rho_{E|\hat{a}=+1} = 2 \operatorname{tr}_{AB} \frac{1+\hat{A}_0}{2} \otimes \mathbb{1}_{BE} |\Psi\rangle\langle\Psi|_{ABE} = \begin{pmatrix} L_1 & 0 & C_\phi \sqrt{L_1 L_3 q} & S_\phi \sqrt{L_1 L_4 q} \\ 0 & L_2 & S_\phi \sqrt{L_2 L_3 q} & -C_\phi \sqrt{L_2 L_4 q} \\ C_\phi \sqrt{L_1 L_3 q} & S_\phi \sqrt{L_2 L_3 q} & L_3 & 0 \\ S_\phi \sqrt{L_1 L_4 q} & -C_\phi \sqrt{L_2 L_4 q} & 0 & L_4 \end{pmatrix}, \quad (59)$$

where the factor of two arises because the probability to observe the outcome $\hat{A}_0 = 1$ is simply $\frac{1}{2}$. The conditional state for the other outcome $\hat{A}_0 = -1$ can be easily obtained by noticing that the two outcomes are interchanged
by a mapping $\hat{A}_0 \to -\hat{A}_0$, i.e. with the inversion of the vector $(S_\phi) \to -(S_\phi)$. Consequently, we have

$$ \rho_{E|\hat{a}=-1} = \begin{pmatrix} 1 & & & \\ & 1 & & \\ & & -1 & \\ & & & -1 \end{pmatrix} \rho_{E|\hat{a}=1} \begin{pmatrix} 1 & & & \\ & 1 & & \\ & & -1 & \\ & & & -1 \end{pmatrix}. \qquad (60) $$

It follows that the entropies of both conditional states are equal and

$$ H(\rho_E|\hat{A}_0) = \frac{1}{2} H(\rho_{E|\hat{a}=+1}) + \frac{1}{2} H(\rho_{E|\hat{a}=-1}) = H(\rho_{E|\hat{a}=+1}). \qquad (61) $$

Analogously to Eq. (60), there are unitary transformations (with different positions of 1 and -1 on the diagonal) that correspond to the transformations $(S_\phi) \to (-S_\phi)$ and $(S_\phi) \to (S_\phi)$. This implies that $H(\rho_E|\hat{A}_0)$ and $H(\hat{A}_0|E)$ do not depend on these sign changes, i.e. they are only functions of $C_\phi^2$ and $S_\phi^2$.

### A.3 Optimal key generating setting $\mathbf{a}_0$

In this section, we give the minimal angle $\phi$ for which constraint $\mathcal{B}_\Omega \ge \beta$ can be fulfilled for a given $\mathbf{L}$. We first focus on the possible values of $\phi$. The expected Bell score is straightforward to compute:

$$ \begin{aligned} \mathcal{B}_\Omega &= \frac{1}{2} \mathrm{tr} \left( \left( (C_\Omega A_0 \otimes (B_0 + B_1) + S_\Omega A_1 \otimes (B_0 - B_1)) \otimes \mathbf{1}_E \right) |\Psi\rangle\langle\Psi|_{ABE} \right) \\ &= \frac{1}{2} C_\Omega \mathbf{a}_0^T \begin{pmatrix} T_z & \\ & T_x \end{pmatrix} (\mathbf{b}_0 + \mathbf{b}_1) + \frac{1}{2} S_\Omega \mathbf{a}_1^T \begin{pmatrix} T_z & \\ & T_x \end{pmatrix} (\mathbf{b}_0 - \mathbf{b}_1) \end{aligned} \qquad (62) $$

From this expression we notice that any of the following transformations of the key generating setting:

$$ \mathbf{a}_0 \to \begin{pmatrix} (-1)^{s_1} & \\ & (-1)^{s_2} \end{pmatrix} \mathbf{a}_0, \qquad (63) $$

with $s_1, s_2 \in \{0, 1\}$ can be compensated by applying the same transformation to the remaining settings $\mathbf{a}_1, \mathbf{b}_y$ to give the same Bell score $\mathcal{B}_\Omega$. Furthermore, we have seen that this transformation does not change Eve's conditional entropy $H(\hat{A}_0|E)$. Therefore we can always restrict $\mathbf{a}_0$ to the positive quadrant of the circle

$$ \phi \in [0, \frac{\pi}{2}] \qquad (64) $$

without loss of generality. We now express the Bell score with Bob's settings parametrized as in Eq. (17):

$$ \mathcal{B}_{\Omega} = \left\langle C_{\Omega} C_{\theta} \mathbf{a}_{0} \cdot \begin{pmatrix} Z \\ X \end{pmatrix} \otimes \mathbf{c} \cdot \begin{pmatrix} Z \\ X \end{pmatrix} + S_{\Omega} S_{\theta} \mathbf{a}_{1} \cdot \begin{pmatrix} Z \\ X \end{pmatrix} \otimes \mathbf{c}_{\perp} \cdot \begin{pmatrix} Z \\ X \end{pmatrix} \right\rangle. \qquad (65) $$

Computing the expected value of the operators on our Bell diagonal state gives

$$ \mathcal{B}_{\Omega} = C_{\Omega} C_{\theta} \mathbf{a}_{0}^{T} \begin{pmatrix} T_z & \\ & T_x \end{pmatrix} \mathbf{c} + S_{\Omega} S_{\theta} \mathbf{a}_{1}^{T} \begin{pmatrix} T_z & \\ & T_x \end{pmatrix} \mathbf{c}_{\perp}. \qquad (66) $$

We introduce an angle $\gamma$ to parametrize the vectors $\mathbf{c}$ and $\mathbf{c}_\perp$:

$$ \mathbf{c} = \begin{pmatrix} C_{\gamma} \\ S_{\gamma} \end{pmatrix}, \quad \mathbf{c}_{\perp} = \begin{pmatrix} -S_{\gamma} \\ C_{\gamma} \end{pmatrix}. \qquad (67) $$

The maximization with the second setting of Alice is straightforward, that is

$$ \begin{aligned} \max_{\mathbf{a}_1} S_{\theta} \mathbf{a}_1^T \begin{pmatrix} T_z & \\ & T_x \end{pmatrix} \mathbf{c}_{\perp} &= \| S_{\theta} \begin{pmatrix} T_z & \\ & T_x \end{pmatrix} \mathbf{c}_{\perp} \| \\ &= |S_{\theta}| \sqrt{T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2}. \end{aligned} $$

The Bell score (optimized with respect to $\mathbf{a}_1$) becomes

$$ \mathcal{B}_{\Omega} = C_{\Omega} C_{\theta} \begin{pmatrix} C_{\phi} \\ S_{\phi} \end{pmatrix}^{T} \begin{pmatrix} T_z & \\ & T_x \end{pmatrix} \begin{pmatrix} C_{\gamma} \\ S_{\gamma} \end{pmatrix} + S_{\Omega} |S_{\theta}| \sqrt{T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2}. $$

Since $T_z, T_x \ge 0$ and $\phi, \Omega \in [0, \frac{\pi}{2}]$, we can also restrict the angle $\theta$ and $\gamma$ to the interval $[0, \frac{\pi}{2}]$ without loss of generality, and drop the absolute value: $|S_\theta| = S_\theta$. The constraint $\mathcal{B}_\Omega \ge \beta$ takes the form

$$ C_{\Omega} \begin{pmatrix} C_{\phi} \\ S_{\phi} \end{pmatrix}^{T} \begin{pmatrix} T_z & \\ & T_x \end{pmatrix} \begin{pmatrix} C_{\gamma} \\ S_{\gamma} \end{pmatrix} \geq \frac{\beta - S_{\Omega} S_{\theta} \sqrt{T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2}}{C_{\theta}}. \qquad (68) $$
Recall that we wish to find the minimal $\phi$ for which this inequality can be fulfilled for at least one value of the free parameters $\theta$ and $\gamma$. We observe that if the right hand side (RHS) can become zero or negative by some choice of $\theta$ and $\gamma$, the constraint becomes trivial. Since $T_z \ge T_x$, this is possible for

$$ \beta^2 \le S_{\Omega}^2 T_z^2 \implies \phi_* = 0. \quad (69) $$

In the following we assume that this is not the case, i.e. $\beta^2 > S_{\Omega}^2 T_z^2$. The angle $\theta$ only appears on the right of the inequality, so to satisfy the inequality our best choice is to minimize the RHS with respect to $\theta$. The expression

$$ \frac{\beta - S_{\Omega}S_{\theta}\sqrt{T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2}}{C_{\theta}} \quad (70) $$

either has a local minimum that can be found by setting its derivative to zero, or there is no local minimum and the expression is minimal at the boundary $\theta = 0$ since it diverges for $\theta \to \frac{\pi}{2}$. Differentiating the expression with respect to $\theta$, we find that a local minimum does exist at

$$ S_{\theta} = \frac{S_{\Omega}\sqrt{T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2}}{\beta}, \quad (71) $$

(recall the assumption above). Plugging this value into Eq. (68) allows us to rewrite the constraint as

$$ C_{\Omega} \begin{pmatrix} C_{\phi} \\ S_{\phi} \end{pmatrix}^T \begin{pmatrix} T_z & 0 \\ 0 & T_x \end{pmatrix} \begin{pmatrix} C_{\gamma} \\ S_{\gamma} \end{pmatrix} \ge \sqrt{\beta^2 - S_{\Omega}^2(T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2)}, \quad (72) $$

which we rewrite as

$$ \begin{pmatrix} C_{\phi} \\ S_{\phi} \end{pmatrix}^T \mathbf{v}_{\gamma} \ge 1, \quad \text{with} \quad (73) $$

$$ \mathbf{v}_{\gamma} = \frac{C_{\Omega}}{\sqrt{\beta^2 - S_{\Omega}^2(T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2)}} \begin{pmatrix} T_z C_{\gamma} \\ T_x S_{\gamma} \end{pmatrix}. $$

Now it becomes simple to check if the constraint can be satisfied at all. The vector $\mathbf{v}_{\gamma}$ belongs to the positive quadrant of the plane with $\mathbf{v}_0 \parallel \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ and $\mathbf{v}_{\pi/2} \parallel \begin{pmatrix} 0 \\ 1 \end{pmatrix}$, hence the inequality $\mathbf{a}_0 \cdot \mathbf{v}_{\gamma} \ge 1$ can be satisfied if and only if the length of the vector $\mathbf{v}_{\gamma}$ reaches 1, i.e.

$$ |\mathbf{v}_{\gamma}|^2 = \frac{C_{\Omega}^2(T_z^2 C_{\gamma}^2 + T_x^2 S_{\gamma}^2)}{\beta^2 - S_{\Omega}^2(T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2)} \ge 1 \quad (74) $$

for some $\gamma$. We rewrite this inequality as

$$ \frac{1}{2} (T_z^2 + T_x^2 + (T_z^2 - T_x^2)C_{2\gamma}C_{2\Omega} - 2\beta^2) \ge 0. \quad (75) $$

Given that $T_z^2 \ge T_x^2$, the left hand side (LHS) is maximal for $\gamma = 0$ if $\Omega \le \frac{\pi}{4}$ and for $\gamma = \frac{\pi}{2}$ if $\Omega > \frac{\pi}{4}$. Hence, $B_{\Omega} \ge \beta$ can be fulfilled if and only if

$$ \begin{array}{ll} \text{for } \Omega \le \frac{\pi}{4}, & C_{\Omega}^2 T_z^2 + S_{\Omega}^2 T_x^2 \ge \beta^2 \\ \text{for } \Omega > \frac{\pi}{4}, & S_{\Omega}^2 T_z^2 + C_{\Omega}^2 T_x^2 \ge \beta^2. \end{array} \quad (76) $$

Consider the first case $\Omega \le \frac{\pi}{4}$. Setting $\gamma = 0$ one verifies that if the constraint can be fulfilled, then it can be fulfilled with $\phi = 0$:

$$ \begin{array}{l} \left( \begin{array}{c} 1 \\ 0 \end{array} \right)^T \mathbf{v}_0 = \frac{C_{\Omega} T_z}{\sqrt{\beta^2 - S_{\Omega}^2 T_x^2}} \ge 1 \quad (77) \\ \iff C_{\Omega}^2 T_z^2 + S_{\Omega}^2 T_x^2 \ge \beta^2. \end{array} $$

Hence, in this case the minimal possible angle is $\phi_* = 0$ as long as the Bell score can be attained as formalized by Eq. (76).

The other case $\Omega > \frac{\pi}{4}$ is less trivial. Assume that the constraint can be satisfied, $S_{\Omega}^2 T_z^2 + C_{\Omega}^2 T_x^2 \ge \beta^2$. We first check if $\phi = 0$ is a solution

$$ \left( \begin{pmatrix} 1 \\ 0 \end{pmatrix}^T \mathbf{v}_{\gamma} \right)^2 = \frac{C_{\Omega}^2 C_{\gamma}^2 T_z^2}{\beta^2 - S_{\Omega}^2(T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2)} \ge 1 \quad (78) \\ \iff C_{\gamma}^2(C_{\Omega}^2 T_z^2 + S_{\Omega}^2 T_x^2 - S_{\Omega}^2 T_z^2) + S_{\Omega}^2 T_z^2 \ge \beta^2 $$

As $S_{\Omega}^2 T_z^2 + C_{\Omega}^2 T_x^2 \ge \beta^2 \ge S_{\Omega}^2 T_z^2$, the LHS is maximal for $\gamma = 0$. Therefore $\phi_* = 0$ iff

$$ C_{\Omega}^2 T_z^2 + S_{\Omega}^2 T_x^2 \ge \beta^2. \quad (79) $$

We now assume

$$ S_{\Omega}^2 T_z^2 + C_{\Omega}^2 T_x^2 \ge \beta^2 > C_{\Omega}^2 T_z^2 + S_{\Omega}^2 T_x^2, \quad (80) $$

such that the constraint can be fulfilled but not with $\phi = 0$. To find the minimal $\phi$ that allows to do so we look on the dependence of the length of the vector $\mathbf{v}_{\gamma}$ on $\gamma$. We compute

$$ \frac{d}{d\gamma} |\mathbf{v}_{\gamma}|^2 = \frac{(T_z^2 - T_x^2)C_{\Omega}^2 S_{2\gamma}}{(\beta^2 - (T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2)S_{\Omega}^2)^2} ((T_z^2 + T_x^2)S_{\Omega}^2 - \beta^2), \quad (81) $$

here $(T_z^2 + T_x^2)S_{\Omega}^2 \ge S_{\Omega}^2 T_z^2 + C_{\Omega}^2 T_x^2 \ge \beta^2$. Thus the derivative is positive and the length of $\mathbf{v}_{\gamma}$ is increasing with $\gamma$.

We can now give a simple geometrical interpretation to our problem of finding the minimal $\phi$: for each value $\phi$ such that

$$ \begin{pmatrix} C_{\phi} \\ S_{\phi} \end{pmatrix}^T \mathbf{v}_{\gamma} \ge 1 \quad (82) $$
a line tangent to the unit circle at $(C_\phi/S_\phi)$ is also crossing the curve $v_\gamma$. So for the minimal value $\phi_*$ there is a line tangent to both $v_\gamma$ and the unit circle at $(C_{\phi_*}/S_{\phi_*})$. The equation of the line tangent to $v_\gamma$ reads $\ell(\lambda) = v_\gamma + \lambda v'_\gamma$. This line is tangent to the unit circle if and only if the equation

$$|v_{\gamma} + \lambda v'_{\gamma}|^2 = 1 \qquad (83)$$

With

$$v'_{\gamma} = \frac{C_{\Omega}}{(\beta^2 - S_{\Omega}^2(T_x^2 C_{\gamma}^2 + T_z^2 S_{\gamma}^2))^{3/2}} \begin{pmatrix} T_z S_{\gamma} (T_z^2 S_{\Omega}^2 - \beta^2) \\ -T_x C_{\gamma} (T_x^2 S_{\Omega}^2 - \beta^2) \end{pmatrix} \qquad (86)$$

lengthy but straightforward algebra gives

$$S_{\gamma}^{2} = - \frac{T_{x}^{2}(\beta^{2} - T_{x}^{2}S_{\Omega}^{2})(\beta^{2} - (T_{z}^{2} + T_{x}^{2})S_{\Omega}^{2})}{(T_{z}^{2} - T_{x}^{2})(\beta^{4} + S_{\Omega}^{2}(T_{x}^{2}T_{z}^{2} - 2\beta^{2}(T_{x}^{2} + T_{z}^{2}))) + S_{\Omega}^{4}(T_{x}^{4} + T_{z}^{4})} \qquad (87)$$

To find that the minimal angle $\phi_*$ note that for the tangent line $(C_{\phi_*}/S_{\phi_*}) \cdot v'_\gamma = 0$, and therefore $(-S_{\phi_*}/C_{\phi_*}) = \frac{v'_\gamma}{|v'_\gamma|}$. Plugging in the above equations we find

$$c_{*}^{2}(\mathbf{L}, \Omega, \beta) = \cos^{2}(\phi_{*}(\mathbf{L}, \Omega, \beta)) = \frac{(\beta^{2} - S_{\Omega}^{2}T_{x}^{2})(C_{\Omega}^{2}T_{x}^{2} + S_{\Omega}^{2}T_{z}^{2} - \beta^{2})}{C_{\Omega}^{2}(T_{z}^{2} - T_{x}^{2})(S_{\Omega}^{2}T_{z}^{2} + S_{\Omega}^{2}T_{x}^{2} - \beta^{2})} \qquad (88)$$

### A.4 Eve's maximum information for $\Omega \le \frac{\pi}{4}$

We here give details on the derivation of the formula (25) in the main text which corresponds to Eve's maximum information in the case $\Omega \le \frac{\pi}{4}$. For these inequalities we have seen that the constraint $C_\Omega^2 T_z^2 + S_\Omega T_x^2 \ge \beta^2$ can be fulfilled with $c_*^2(\mathbf{L}, \beta) = 1$. The choice of the measurement angle implies $C_\phi = 1$ implies a simple form for the state?

$$\rho_{E|\hat{a}=+1} = \begin{pmatrix} L_1 & 0 & \sqrt{L_1 L_3 q} & 0 \\ 0 & L_2 & 0 & -\sqrt{L_2 L_4 q} \\ \sqrt{L_1 L_3 q} & 0 & L_3 & 0 \\ 0 & -\sqrt{L_2 L_4 q} & 0 & L_4 \end{pmatrix}, \qquad (89)$$

with a closed form expression for its eigenvalues implying

$$H(\rho_{E|\hat{a}=1}) = H(\mathbf{p}) = \begin{pmatrix} \frac{1}{2} (L_1 + L_3 + \sqrt{4L_1L_3q + (L_1 - L_3)^2}) \\ \frac{1}{2} (L_1 + L_3 - \sqrt{4L_1L_3q + (L_1 - L_3)^2}) \\ \frac{1}{2} (L_2 + L_4 + \sqrt{4L_2L_4q + (L_2 - L_4)^2}) \\ \frac{1}{2} (L_2 + L_4 - \sqrt{4L_2L_4q + (L_2 - L_4)^2}) \end{pmatrix}. \qquad (90)$$

The constraint on the generalized CHSH score leads to the following constraint on the vector $\mathbf{L}$

$$\begin{pmatrix} T_Z^2 \\ T_X^2 \end{pmatrix} \cdot \begin{pmatrix} C_\Omega^2 \\ S_\Omega^2 \end{pmatrix} = 2(L_3 - L_4)(L_1 - L_2)C_{2\Omega} + (L_1 - L_2)^2 + (L_3 - L_4)^2 \ge \beta^2. \qquad (91)$$

Our goal is thus to find the components of the vector $\mathbf{L}$ maximizing $H(\mathbf{L}) - H(\mathbf{p})$ and satisfying the previous constraint. Inspired by Ref. [13], we first introduce the following parametrization

$$L_1 = P x \quad L_3 = P(1-x) \quad L_2 = (1-P)y \quad L_4 = (1-P)(1-y). \qquad (92)$$

The partial ordering of the $\mathbf{L}$ coefficients implies

$$(1-P)y \le P x \le (1-P)y + 2P - 1 \qquad (93)$$

has only one solution, where the derivative is with respect to $\gamma$. This is a quadratic equation

$$\lambda^2 |v'_{\gamma}|^2 + 2\lambda v'_{\gamma} \cdot v_{\gamma} + |v_{\gamma}|^2 - 1 = 0, \qquad (84)$$

which has a unique solution iff its determinant is zero

$$(v'_{\gamma} \cdot v_{\gamma})^2 = |v'_{\gamma}|^2 (|v_{\gamma}|^2 - 1). \qquad (85)$$
which requires $P \ge \frac{1}{2}$. The advantage of this parametrization comes from the fact that our figure of merit can be nicely rewritten as

$$
\begin{align}
H(\mathbf{L}) - H(\mathbf{p}) &= P h_q(x) + (1-P)h_q(y) \nonumber \\
h_q(z) &= h(z) - h(n_q(z)) \tag{94} \\
n_q(z) &= \frac{1 + \sqrt{1 - 4(1-q)z(1-z)}}{2} \nonumber
\end{align}
$$

where $h(z) = -z \log(z) - (1-z) \log(1-z)$ is the binary entropy function with the logarithm in base 2, while the constraint on the expected value of the generalized CHSH operator is given by

$$
S_{\Omega} (2P(x + y - 1) - 2y + 1)^2 + C_{\Omega} (1 - 2P)^2 \geq \beta^2 \quad (95)
$$

For a fixed $P$, the curve in the $(x, y)$-plan that cor-
responds to a constant value $\beta$ satisfies $P dx =
(1 - P) dy$. This remark allows one to maximize
$H(\mathbf{L}) - H(\mathbf{p})$ along this curve and find that it is op-
timal for Eve to set $x + y = 1$, see appendix C2 in
Ref. [13] for the detailed argument. The symmetry of
the function $h_q(x) = h_q(1 - x)$ allows us to write the
problem as

$$
\begin{align}
\max_{x,P} H(\mathbf{L}) - H(\mathbf{p}) &= h_q(x) \tag{96} \\
\text{s.t.} \ (2P - 1)^2 C_{\Omega}^2 + (1 - 2x)^2 S_{\Omega}^2 &\ge \beta^2. \nonumber
\end{align}
$$

As the goal function $h_q(x)$ does not depend on $P$, we
can set its value to $P=1$ because this is the value
maximizing the LHS of the constraint inequality and
allowing the largest possible interval for the remaining
variable $x$. We thus get

$$
\begin{align}
\max_{x,P} H(\mathbf{L}) - H(\mathbf{p}) &= h_q(x) \tag{97} \\
\text{s.t.} \ (1-2x)^2 &\ge \frac{\beta^2 - C_{\Omega}^2}{S_{\Omega}^2}. \nonumber
\end{align}
$$

Finally, as $h_q(x)$ is a monotonically decreasing func-
tion of $x$ (see Ref. [13]), it is optimal to set $x$ to the
least possible value compatible with the constraint.
This implies

$$
\begin{align}
I(\beta; \Omega, q) &= h_q(z) \nonumber \\
\text{with } z &= \frac{1}{2} \left( \frac{\sqrt{\beta^2 - C_{\Omega}^2}}{S_{\Omega}} + 1 \right). \tag{98}
\end{align}
$$

Let us now recall that the situation with $c_*^2(\mathbf{L}, \beta) = 1$ and $C_\Omega^2 T_z^2 + S_\Omega T_x^2 \ge \beta^2$ also occurs in the case $\Omega > \frac{\pi}{4}$. The above proof guarantees that it is optimal for Eve to use strategies where the inequality is saturated $C_\Omega^2 T_z^2 + S_\Omega T_x^2 = \beta^2$. Hence in the optimization problem for $\Omega > \frac{\pi}{4}$ we can ignore all the strategies with $C_\Omega^2 T_z^2 + S_\Omega T_x^2 > \beta^2$.

A.5 Concavity of $h_q \circ z(\beta)$

Recall that in order to use the bound $I(\beta; \Omega, q)$ derived for two-qubit strategies in the previous section as a universal bound (valid for strategies in any dimension), we have to show that this function is concave. In this case, for any mixture of qubit strategies (enforced by the Jordan’s lemma) with an average score $\bar{\beta} = \sum_i p_i \beta_i$, Eve’s information satisfies

$$
\bar{I}(\beta; \Omega, q) = \sum_{i} p_i I(\beta_i; \Omega, q) \leq I(\bar{\beta}; \Omega, q). \quad (99)
$$

The concavity of $I(\beta; \Omega, q) = h_q(z(\beta))$ follows from
the fact that its second derivative is negative

$$
\frac{d^2}{d\beta^2} h_q(z(\beta)) = h_q''(z) (z'(\beta))^2 + h_q'(z) z''(\beta) \le 0, \quad (100)
$$

which we are going to show below. In this section we
will use the natural algorithm instead of logarithm
in base 2. The function $h_q(z)$ takes a positive real
factor upon changing the base of the algorithm, so it
is irrelevant for its concavity. Note first that $z(\beta) \in$
$[\frac{1}{2}, 1]$ and

$$
\sqrt{\beta^2 - C_{\Omega}^2} = S_{\Omega}(2z - 1). \quad (101)
$$

Then consider the following identities

$$
\begin{align}
(z'(\beta))^2 &= \frac{\beta^2}{4S_{\Omega}^2(\beta^2 - C_{\Omega}^2)} = \frac{S_{\Omega}^2(2z-1)^2 + C_{\Omega}^2}{4S_{\Omega}^4(2z-1)^2} \tag{102} \\
z''(\beta) &= \frac{-C_{\Omega}^2}{2S_{\Omega}(\beta^2 - C_{\Omega}^2)^{3/2}} = \frac{-C_{\Omega}^2}{2S_{\Omega}^4(2z-1)^3}. \notag
\end{align}
$$

The identity (100) that we want to prove thus be-
comes

$$
\frac{h_q''(z)(2z-1)(S_{\Omega}^2(2z-1)^2 + C_{\Omega}^2) - 2h_q'(z)C_{\Omega}^2}{4S_{\Omega}^4(2z-1)^3} \le 0. \quad (103)
$$

Multiplying by a positive fraction $\frac{4S_{\Omega}^4(2z-1)^3}{C_{\Omega}^2}$ it can be straightforwardly simplified to the form

$$
h_q''(z)(2z-1)(T_{\Omega}^2(2z-1)^2+1)-2h_q'(z) \le 0. \quad (104)
$$

As $h_q''(z) \le 0$ was proven in [13, 36], we have the
inequality

$$
h_q''(z)(2z-1)^2 T_{\Omega}^2 \le 0. \quad (105)
$$
We use it to relax the inequality in Eq. (104) to

$$h_q''(z)\left(z-\frac{1}{2}\right)-h_q'(z)\le 0. \quad (106)$$

At the point $z = \frac{1}{2}$ the left hand side becomes zero, since $h'_{q}(\frac{1}{2}) = 0$ and $|h''_{q}(\frac{1}{2})| < \infty$, see below. From now on, we thus exclude the point $z = \frac{1}{2}$ and consider $z \in (\frac{1}{2}, 1]$. Now we can divide the whole expression by a strictly positive factor $(z - \frac{1}{2})^2$. We obtain

$$\frac{1}{2} \frac{h_q''(z)(z - \frac{1}{2}) - h_q'(z)(z - \frac{1}{2})'}{(z - \frac{1}{2})^2} \le 0, \quad (107)$$

or

$$\frac{1}{2} \left( \frac{h'_{q}(z)}{z - \frac{1}{2}} \right)' = \left( \frac{h'_{q}(z)}{2z - 1} \right)' \le 0. \quad (108)$$

In other words we want to show that the function

$$(*) \quad \frac{d}{dz} \left( \frac{h'_{q}(z)}{2z-1} \right) \le 0 \quad (109)$$

on the interval $z \in (\frac{1}{2}, 1]$. Let us now compute this function:

$$\frac{h'_{q}(z)}{2z-1} = \frac{h'(z)}{2z-1} - h'(n_{q}(z)) \frac{n'_{q}(z)}{2z-1}. \quad (110)$$

The last fraction can be simplified to

$$\begin{aligned} \frac{n'_{q}(z)}{2z-1} &= \frac{(2z-1)(1-q)}{(2z-1)\sqrt{1-4(1-q)z(1-z)}} \\ &= \frac{1-q}{2n_{q}(z)-1}. \end{aligned} \quad (111)$$

Therefore

$$\begin{aligned} \frac{h'_{q}(z)}{2z-1} &= \frac{h'(z)}{2z-1} - (1-q)\frac{h'(n_{q}(z))}{2n_{q}(z)-1} \\ &= g(z) - (1-q)g(n_{q}(z)), \end{aligned} \quad (112)$$

where

$$g(z) = \frac{h'(z)}{2z-1} = -\frac{\log\left(\frac{z}{1-z}\right)}{2z-1}. \quad (113)$$

To complete the proof we thus need to show that

$$(*) \quad g'(z) - (1-q)g'(n_q(z))n'_q(z) \le 0. \quad (114)$$

From Eq. (111), we have

$$n'_{q}(z) = (1-q) \frac{2z-1}{2n_{q}(z)-1}, \quad (115)$$

so the inequality to be shown can be rewritten as

$$\frac{g'(z)}{2z-1} - (1-q)^2 \frac{g'(n_q(z))}{2n_q(z)-1} \le 0. \quad (116)$$

For $q = 0$ we have $(1-q) = 1$ and $n_q(z) = z$ so the two terms are equal. The identity to be shown can then be expressed as

$$\begin{aligned} (*) \quad &f(z,0) - f(z,q) \le 0 && \text{with} \\ &f(z,q) = (1-q)^2 u(n_q(z)) \\ &u(z) = \frac{g'(z)}{2z-1} = \frac{\frac{1}{z-1} + \frac{1}{z} + 2\log\left(\frac{z}{1-z}\right)}{(2z-1)^3}, \end{aligned} \quad (117)$$

which holds for $q = 0$ trivially. To show that it holds for all $q$ it is sufficient to demonstrate that the function $f(z, q)$ is increasing with $q$, i.e.

$$(*) \quad \frac{d}{dq} f(z, q) \ge 0, \quad (118)$$

which we are going to show now. Using

$$\begin{aligned} \frac{d}{dq} n_q(z) &= \frac{(1-z)z}{\sqrt{1-4(1-q)(1-z)z}} \\ &= \frac{n_q(z)(1-n_q(z))}{(1-q)(2n_q(z)-1)} \end{aligned} \quad (119)$$

$$\begin{aligned} \frac{d}{dq} f(z,q) &= -2(1-q)u(n_q(z)) + (1-q)^2 u'(n_q(z)) \frac{d}{dq} n_q(z) \\ &= -2(1-q)u(n_q(z)) + (1-q)^2 u'(n_q(z)) \frac{n(1-n_q(z))}{2n_q(z)-1} \end{aligned} \quad (120)$$

for $n > \frac{1}{2}$, which is positive iff $q = 1$ or

$$(*) \quad u'(n) n(1-n) - 2u(n)(2n-1) \ge 0. \quad (121)$$

In the case $q < 1$, straightforward algebra allows to find a simple expression of the left hand side and rewrite the condition as

$$\begin{aligned} & \frac{6n^2 - 4n^3 - 1 + 4(n^4 - 2n^3 + 2n^2 - n) \log\left(\frac{n}{1-n}\right)}{(2n-1)^4 n(1-n)} \ge 0 \\ & \iff 6n^2 - 4n^3 - 1 + 4(n^4 - 2n^3 + 2n^2 - n) \log\left(\frac{n}{1-n}\right) \ge 0 \end{aligned} \quad (122)$$
Changing the variable to $x + 1 = \frac{n}{1-n}$ with $x \ge 0$ we express the above inequality as

$$ \frac{x(x^2 + 6x + 6)}{(x+2)^3} - 4 \frac{(x+1)(x^2 + 3x + 3)}{(x+2)^4} \log(x+1) \ge 0 \qquad (123) \\ \iff \log(1+x) \le \frac{x(x+2)(x^2 + 6x + 6)}{4(x+1)(x^2 + 3x + 3)}. $$

To prove this relation note that there is equality for $x = 0$, but the LHS increases slower than the RHS

$$ \begin{aligned} \log'(1+x) &\le \left( \frac{x(x+2)(x^2+6x+6)}{4(x+1)(x^2+3x+3)} \right)' &\iff \\ \frac{1}{1+x} &\le \frac{1}{4} \left( \frac{9(x+1)}{(x(x+3)+3)^2} + \frac{1}{(x+1)^2} \right. &\iff \\ &\qquad \left. + \frac{3}{x(x+3)+3} + 1 \right) &\iff \\ &\qquad \frac{x^4(x+2)^2}{4(x+1)^2(x(x+3)+3)^2} \ge 0, \end{aligned} \qquad (124) $$

which concludes the proof.

**Properties of $h'_q(\frac{1}{2})$ and $h''_q(\frac{1}{2})$** We start with the first derivative and want to show that $h'_q(\frac{1}{2}) = 1/2$. We have

$$ \begin{aligned} h'_q(z) &= (h(z) - h(n_q(z)))' \\ &= h'(z) - h'(n_q(z))n'_q(z) \end{aligned} \qquad (125) $$

The binary entropy hits a maximum at $z = 1/2$ so $h'(\frac{1}{2}) = 0$. For the second term, we use (115) to get

$$ \begin{aligned} h'(n_q(z))n'_q(z) &= -(1-q) \frac{2z-1}{2n_q(z)-1} \log \left( \frac{n_q(z)}{1-n_q(z)} \right), \qquad (126) \\ (2z-1) = 0 \text{ at } z = \frac{1}{2}, \text{ while} \end{aligned} $$

$$ \begin{aligned} & \left| \frac{1-q}{2n_q(z)-1} \log \left( \frac{n_q(z)}{1-n_q(z)} \right) \right|_{z=1/2} \\ & = \left| \frac{(1-q) \log \left( \frac{1+\sqrt{q}}{1-\sqrt{q}} \right)}{\sqrt{q}} \right| \end{aligned} \qquad (127) $$

Changing the variable to $x + 1 = \frac{1+\sqrt{q}}{1-\sqrt{q}}$ with $x \ge 0$ yields for the last expression

$$ \left| \frac{4(x+1) \log(x+1)}{x(x+2)} \right| < \infty. \qquad (128) $$

It is obviously bounded for $x > \epsilon$ with any $\epsilon$, and the fact that the limit $x \to 0$ exists can be seen by straightforward application of l'Hôpital's rule. Hence,

$$ h'_{q}(\frac{1}{2}) = 0. \qquad (129) $$

We also wish to show that the second derivative $h''_q(\frac{1}{2})$ is bounded. To do so we compute

$$ h''_{q}(\frac{1}{2}) = \frac{2(1-q) \log \left( \frac{1+\sqrt{q}}{1-\sqrt{q}} \right)}{\sqrt{q}} - 4. \qquad (130) $$

But as we have just shown that

$$ \left| \frac{(1-q) \log \left( \frac{1+\sqrt{q}}{1-\sqrt{q}} \right)}{\sqrt{q}} \right| < \infty, \qquad (131) $$

the desired result

$$ |h''_{q}(\frac{1}{2})| < \infty \qquad (132) $$

follows.

## A.6 Maximization of the generalized CHSH score with respect to auxiliary settings

In the goal function

$$ H(\rho_E) - H(\rho_E | \hat{a}_0(q)) + t \mathcal{B}_{\Omega}(\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1) \qquad (133) $$

that a priori appears in Eq. (38), it is only the Bell score which depends on the auxiliaty measurement settings $\mathbf{a}_1, \mathbf{b}_0$ and $\mathbf{b}_1$. As $t$ is always positive we can straightforwardly maximise the score with respect to these settings. We thus define

$$ \beta_{\max}(\mathbf{L}, \phi) = \max_{\mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1} \mathcal{B}_{\Omega}(\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1), \qquad (134) $$

which actually appears in Eq. (38).

Let us now compute this expression starting from Eq. (66), that we put in the form

$$ \mathcal{B}_{\Omega} = \begin{pmatrix} C_{\theta} \\ S_{\theta} \end{pmatrix} \cdot \begin{pmatrix} C_{\Omega} \mathbf{a}_{0}^{T} \begin{pmatrix} T_{z} C_{\gamma} \\ T_{x} S_{\gamma} \end{pmatrix} \\ S_{\Omega} \mathbf{a}_{1}^{T} \begin{pmatrix} -T_{z} S_{\gamma} \\ T_{x} C_{\gamma} \end{pmatrix} \end{pmatrix}. \qquad (135) $$
This form makes the maximization with respect to $\theta$ and $\mathbf{a}_1$ straightforward

$$ 
\begin{aligned}
& \max_{\theta, \mathbf{a}_1} \mathcal{B}_{\Omega} \\
& = \max_{\mathbf{a}_1} \sqrt{\left( C_{\Omega} \mathbf{a}_0^T \begin{pmatrix} T_z C_{\gamma} \\ T_x S_{\gamma} \end{pmatrix} \right)^2 + \left( S_{\Omega} \mathbf{a}_1^T \begin{pmatrix} -T_z S_{\gamma} \\ T_x C_{\gamma} \end{pmatrix} \right)^2} \\
& = \sqrt{C_{\Omega}^2 (C_{\phi} T_z C_{\gamma} + S_{\phi} T_x S_{\gamma})^2 + S_{\Omega}^2 (T_z^2 S_{\gamma}^2 + T_x^2 C_{\gamma}^2)}.
\end{aligned} 
\quad (136) $$

To find the maximum with respect to $\gamma$ or $\mathbf{c} = \begin{pmatrix} C_{\gamma} \\ S_{\gamma} \end{pmatrix}$ it is convenient to write the expression inside the square root as

$$ \left( \max_{\theta, \mathbf{a}_1} \mathcal{B}_{\Omega} \right)^2 = \mathbf{c}^T \begin{pmatrix} C_{\Omega}^2 C_{\phi}^2 T_z^2 + S_{\Omega}^2 T_x^2 & C_{\Omega}^2 C_{\phi} S_{\phi} T_z T_x \\ C_{\Omega}^2 C_{\phi} S_{\phi} T_z T_x & C_{\Omega}^2 S_{\phi}^2 T_x^2 + S_{\Omega}^2 T_z^2 \end{pmatrix} \mathbf{c}. \quad (137) $$

It is now obvious that the value is maximal if $\mathbf{c}$ is aligned with the eigenvector of the matrix which corresponds to the maximal eigenvalue. Therefore

$$ 
\begin{aligned}
\beta_{\max}(\mathbf{L}, \phi) &= \max_{\theta, \gamma, \mathbf{a}_1} \mathcal{B}_{\Omega}(\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1), \\
&= \sqrt{\mathrm{Eig}_+ \begin{pmatrix} C_{\Omega}^2 C_{\phi}^2 T_z^2 + S_{\Omega}^2 T_x^2 & C_{\Omega}^2 C_{\phi} S_{\phi} T_z T_x \\ C_{\Omega}^2 C_{\phi} S_{\phi} T_z T_x & C_{\Omega}^2 S_{\phi}^2 T_x^2 + S_{\Omega}^2 T_z^2 \end{pmatrix}} \\
&= \frac{1}{\sqrt{2}} \left( C_{\Omega}^2 (C_{\phi}^2 T_z^2 + S_{\phi}^2 T_x^2) + S_{\Omega}^2 (T_z^2 + T_x^2) + \sqrt{(C_{\Omega}^2 (C_{\phi}^2 T_z^2 - S_{\phi}^2 T_x^2) - S_{\Omega}^2 (T_z^2 - T_x^2))^2 + 4(C_{\Omega}^2 C_{\phi} S_{\phi} T_z T_x)^2} \right)^{1/2}.
\end{aligned} 
\quad (138) $$

# B Numerical tool

## B.1 Lipshitz continuity of von Neumann entropy

Consider two states $\rho$ and $\sigma$ on an $n$-dimensional Hilbert space, that are close in fidelity:

$$ F(\rho, \sigma) = \operatorname{tr}|\sqrt{\rho}\sqrt{\sigma}| = f. \quad (139) $$

Given the monotonicity of arccos in the range $[0, 1]$, this condition can be equivalently written in terms of the angle $A(\rho, \sigma) = \arccos(F(\rho, \sigma))$

$$ A(\rho, \sigma) = a = \arccos(f). \quad (140) $$

The angle is a metric on the space of density operators [22], in particular it satisfies the triangle inequality.

Next, note that the angle between two states is lower bounded by the angle between the ordered vectors made of their ordered eigenvalues

$$ a = A(\rho, \sigma) \ge A(\mathbf{p}, \mathbf{q}) = \arccos(\sqrt{\mathbf{p}} \cdot \sqrt{\mathbf{q}}), \quad (141) $$

with $\mathbf{p} = \mathrm{Eig}^{\downarrow}(\rho)$, $\mathbf{q} = \mathrm{Eig}^{\uparrow}(\sigma)$ such that $p_1 \ge p_2 \ge \dots$ This inequality follows from

$$ 
\begin{aligned}
F(\rho, \sigma) &= \operatorname{tr}|\sqrt{\rho}\sqrt{\sigma}| = \max_{U} \operatorname{tr}\sqrt{\rho}(\sqrt{\sigma}U) \\
&\le \sqrt{\mathbf{p}} \cdot \sqrt{\mathbf{q}},
\end{aligned} 
\quad (142) $$

where in second line we used the so-called von Neumann trace inequality [21]. This bound is useful because the entropies of the states match the entropies of the probability distributions

$$ H(\rho) = H(\mathbf{p}) \quad H(\sigma) = H(\mathbf{q}). \quad (143) $$

Let us now bound their difference

$$ \Delta H = |H(\rho) - H(\sigma)| = |H(\mathbf{p}) - H(\mathbf{q})|. \quad (144) $$

To do so, note that for any two unit vectors $\sqrt{\mathbf{p}}$ and $\sqrt{\mathbf{q}}$ on the $n$-sphere there exist a path $\gamma$ connecting the two and such that the integral along the path satisfies

$$ \int_{\sqrt{\mathbf{p}}}^{\sqrt{\mathbf{q}}} dA = A(\mathbf{p}, \mathbf{q}) \le a. \quad (145) $$

Let us bound the variation of the entropy along the path. To this end, we associate a probability distribution $\mathbf{r}$ to each vector $\mathbf{v}$ on the path $\gamma$, with $r_i = (v^{(i)})^2$ (note that the vectors along the curve remain in the positive part of the $n$-sphere $v^{(i)} \ge 0$). A step $dA$ from $\mathbf{v}$ along the path corresponds to some deformation of the vector given by

$$ \mathbf{v}_{dA} \rightarrow \mathbf{v} + \mathbf{v}_{\perp} dA, \quad (146) $$

with $\mathbf{v} \cdot \mathbf{v}_{\perp} = 0$. From

$$ H(\mathbf{r}) = - \sum_{i} (v^{(i)})^2 \log((v^{(i)})^2) \quad (147) $$
we compute the entropy variation for an infinitesimal increment of the angle

$$
\begin{align}
\left| \frac{dH}{dA} \right| &= 2 \left| \sum_i v^{(i)} v_{\perp}^{(i)} \left( \frac{1}{\log_e(2)} + \log((v^{(i)})^2) \right) \right| \nonumber \\
&= 2 \left| \sum_i v_{\perp}^{(i)} v^{(i)} \log((v^{(i)})^2) \right| \tag{148} \\
&= 2 \left| \mathbf{v}_{\perp} \cdot \mathbf{w} \right| \le 2 \left| \mathbf{w} \right| \nonumber
\end{align}
$$

where we defined a vector **w** as wᵢ = v⁽ᵢ⁾ log((v⁽ᵢ⁾)²) = √rᵢ log(rᵢ). Hence we obtain

$$
\left| \frac{dH}{dA} \right| \le 2 \sqrt{\sum_i r_i \log^2(r_i)}, \quad (149)
$$

and it remains to bound the expression on the RHS.
To do so, we will construct a concave upper bound on
the function

$$
c(r) = r \log^2(r) \quad (150)
$$

defined on $[0, 1]$. By computing the second derivative

$$
c''(r) = 2 \frac{1 + \log_e(r)}{r \log_e(2)^2} \quad (151)
$$

we see that the function is concave $c'' \le 0$ on the interval $r \in [0, \frac{1}{e}]$, and convex on the complement. To get a concave upper-bound we look for a line passing by $r=0$ and $c(1)=0$ and tangent to $c(r)$. In the $(r, c)$-plane the equation of a line tangent to $c(r)$ at $r$ is given by

$$
\ell_r(\lambda) = \begin{pmatrix} r \\ c(r) \end{pmatrix} + \begin{pmatrix} 1 \\ c'(r) \end{pmatrix} \lambda. \quad (152)
$$

It passes through the point $\ell_r(\lambda) = (1, 0)$ iff

$$
\begin{align}
c(r) + c'(r)(r-1) &= 0 \notag \\
\log(r)(-2r + \log(r) + 2) &= 0 \notag \\
\text{with solutions} \tag{153} \\
r_1 &= -\frac{1}{2}W_0(-\frac{2}{e^2}) \approx 0.203 \notag \\
r_2 &= 1 \notag
\end{align}
$$

$$
\begin{cases}
r_1 = -\frac{1}{2}W_0(-\frac{2}{e^2}) \approx 0.203 \\
r_2 = 1
\end{cases}
$$

where $W_0$ is the principal branch of the Lambert W-function. The second solution $r_2 = 1$ is irrelevant. From the first solution we can construct a concave upper bound

$$
c(r) \le \hat{c}(r) = \begin{cases} c(r) & r \le r_1 \\ c(r_1)^{\frac{1-r}{1-r_1}} & r > r_1. \end{cases} \quad (154)
$$

Note that $r_1 < \frac{1}{e}$.

Finally with the concave bound $\hat{c}$ it is easy to obtain

$$
\begin{align}
\sum_i r_i \log^2(r_i) &= \sum_i c(r_i) \nonumber \\
&\le \sum_i \hat{c}(r_i) = n \sum_i \frac{1}{n} \hat{c}(r_i) \tag{155} \\
&\le n \hat{c} \left( \sum_i \frac{r_i}{n} \right) = n \hat{c} \left( \frac{1}{n} \right). \nonumber
\end{align}
$$

So for entropy susceptibility we get

$$
\left| \frac{dH}{dA} \right| \le 2 \sqrt{n \hat{c} \left( \frac{1}{n} \right)} \quad (156)
$$

For $n \le 4$ one has $\frac{1}{n} > r_1$, while for $n \ge 5$ the contrary is true, $\frac{1}{n} < r_1$. Therefore we can have a more explicit expression

$$
\left| \frac{dH}{dA} \right| \le \begin{cases} 4\sqrt{r_1(1-r_1)}\sqrt{n-1} & n \le 4 \\ 2\log(n) & n \ge 5 \end{cases}, \quad (157)
$$

where we used the equation $-$log(r₁) = 2(1 − r₁) to
compute c(r₁). In particular we get for n = 4

$$
\hat{c}\left(\frac{1}{4}\right) = c(r_1) \frac{1-1/4}{1-r_1} < 1.011. \quad (158)
$$

and deduce

$$
\left| \frac{dH}{dA} \right| \le 2 \sqrt{4 \hat{c} \left( \frac{1}{4} \right)} < 4.023. \quad (159)
$$

To bound the global entropy difference, we simply integrate along the curve $\gamma$ to get

$$
\begin{align}
\Delta H &\le \int_{\sqrt{\hat{\mathbf{P}}}}^{\sqrt{\mathbf{q}}} |dH| \nonumber \\
&\Leftrightarrow \Delta H \le \int_{\sqrt{\hat{\mathbf{P}}}}^{\sqrt{\mathbf{q}}} \left| \frac{dH}{dA} \right| dA < 4.023 A(\rho, \sigma). \tag{160}
\end{align}
$$

As a side remark, we note that for $n \ge 5$, we get
$\frac{1}{n} < r_1$ and the bound becomes simpler

$$
\hat{c}\left(\frac{1}{n}\right)=\frac{\log ^{2}\left(n\right)}{n},
\quad(161)
$$

and

$$
\left| \frac{dH}{dA} \right| \le 2 \log(n). \quad (162)
$$

B.2 Continuity of the goal function

**The entropy term** – To apply the continuity bound previously described to our situation, let $\rho, \rho'$ be the states on $\hat{A}_0 E$ produced by measurements along angles $\phi, \phi'$ on the states $|\Psi\rangle_{ABE}$ with the weight $\mathbf{L}, \mathbf{L}'$ respectively. Our aim is to bound $|H(\hat{A}_0|E)_\rho - H(\hat{A}_0|E)_{\rho'}|$. We use

$$
\begin{align}
& |H(\hat{A}_0|E)_\rho - H(\hat{A}_0|E)_{\rho'}| \nonumber \\
& = |H(\rho_E) - H(E|\hat{A}_0)_\rho - H(\rho'_E) + H(E|\hat{A}_0)_{\rho'}| \tag{163} \\
& \le |H(\rho_E) - H(\rho'_E)| + |H(E|\hat{A}_0)_\rho - H(E|\hat{A}_0)_{\rho'}| \nonumber
\end{align}
$$

and bound the two last terms independently. For
the first term involving $H(\rho_E)_\rho$, things are straight-
forward. The states of Eve $\rho_E = \text{diag}(\mathbf{L})$ are 4-
dimensional and independent of $\phi$, so

$$
|H(\rho_E)_{\mathbf{L}} + H(\rho_E)_{\mathbf{L}'}| \le 4.023 \arccos(\sqrt{\mathbf{L}} \cdot \sqrt{\mathbf{L}'}) . \quad (164)
$$
For the other term, we note

$$
\begin{align}
& \left| H(E|\hat{A}_0)_\rho - H(E|\hat{A}_0)_{\rho'} \right| \tag{165} \\
&= \left| H(\hat{A}_0 E)_\rho - H(\hat{A}_0 E)_{\rho'} \right| \nonumber \\
&= \left| H(\rho) - H(\rho') \right| \nonumber
\end{align}
$$

as $H(\hat{A}_0) = 1$. Since $\hat{A}_0 E$ is 8-dimensional we get

$$
\begin{equation}
\begin{aligned}
& \left| H(E|\hat{A}_0)_\rho - H(E|\hat{A}_0)_{\rho'} \right| \\
& \le 2 \log(8) A(\rho, \rho') = 6 A(\rho, \rho').
\end{aligned}
\tag{166}
\end{equation}
$$

Given that $\rho = \rho(\mathbf{L}, \phi)$ and $\rho' = \rho(\mathbf{L}', \phi')$, we can use
the triangle inequality to write

$$
\begin{equation}
\begin{aligned}
& A(\rho, \rho') \\
& \le A(\rho(\mathbf{L}, \phi), \rho(\mathbf{L}', \phi)) + A(\rho(\mathbf{L}', \phi), \rho(\mathbf{L}, \phi')).
\end{aligned}
\tag{167}
\end{equation}
$$

For $A(\rho(\mathbf{L}, \phi), \rho(\mathbf{L}', \phi))$ we note that both states result from the action of the same CPTP map on two initial state $|\Psi(\mathbf{L})\rangle_{ABE}$ and $|\Psi(\mathbf{L}')\rangle_{ABE}$. By the data-processing inequality (inherited by the angle from the fidelity $F$) we have

$$
\begin{align}
& A(\rho(\mathbf{L}, \phi), \rho(\mathbf{L}', \phi)) \nonumber \\
& \le A(|\Psi(\mathbf{L})\rangle_{ABE}, |\Psi(\mathbf{L}')\rangle_{ABE}) \tag{168} \\
& = A(\mathbf{L}, \mathbf{L}') \nonumber \\
& = \arccos(\sqrt{\mathbf{L}} \cdot \sqrt{\mathbf{L}'}). \nonumber
\end{align}
$$

To bound the other term $A(\rho(\mathbf{L}', \phi), \rho(\mathbf{L}', \phi'))$ we note that if we apply a channel that performs a Z measurement followed by noisy preprocessing to the state $(e^{i\phi Y_A/2} \otimes \mathbf{1}_{BE}) |\Psi\rangle_{ABE}$, this produces exactly the state $\rho$ on $\hat{A}_0 E$; analogously, applying the same channel to the state $(e^{i\phi' Y_A/2} \otimes \mathbf{1}_{BE}) |\Psi\rangle_{ABE}$ produces the state $\rho'$. Therefore the data-processing inequality implies $F(\rho, \rho')$ is lower-bounded by

$$
\begin{align}
F &\left( \left( e^{i\phi Y_A/2} \otimes \mathbf{1}_{BE} \right) |\Psi\rangle, \left( e^{i\phi' Y_A/2} \otimes \mathbf{1}_{BE} \right) |\Psi\rangle \right) \nonumber \\
&= \left| \langle\Psi| \left( e^{i(\phi'-\phi)Y_A/2} \otimes \mathbf{1}_{BE} \right) |\Psi\rangle \right| \nonumber \\
&= \left| \cos \frac{\phi'-\phi}{2} + i \sin \frac{\phi'-\phi}{2} \left\langle \Psi \middle| (Y_A \otimes \mathbf{1}_{BE}) \middle| \Psi \right\rangle \right| \tag{169} \\
&\ge \left| \cos \frac{\phi'-\phi}{2} \right|, \quad \text{since } \langle\Psi| (Y_A \otimes \mathbf{1}_{BE}) |\Psi\rangle \in \mathbb{R}. \nonumber
\end{align}
$$

Putting these together, we conclude that (for $|\phi' - \phi| < \pi$)

$$
A(\rho(\mathbf{L}', \phi), \rho(\mathbf{L}', \phi')) \le \frac{|\phi' - \phi|}{2}. \quad (170)
$$

Combining everything together, we get

$$
\begin{align}
& \left| H(\hat{A}_0 | E)_\rho - H(\hat{A}_0 | E)_{\rho'} \right| \tag{171} \\
& \le 10.023 \arccos(\sqrt{\mathbf{L}} \cdot \sqrt{\mathbf{L}'}) + 3 |\phi - \phi'| \nonumber
\end{align}
$$

**The Bell score** – Next we wish to bound the increment $|d\beta_{\max}(\mathbf{L}, \mathbf{a}_0, \Omega)|$ for infinitesimal changes of the parameters ($d\mathbf{L}, d\phi$). It is actually straightforward to bound the gradient of the Bell score before the maximization with respect to $\mathbf{a}_0, \mathbf{b}_0, \mathbf{b}_1$, so we just need to be careful to apply this bound on $\beta_{\max}$.

First, note that $\mathcal{B}_{\Omega}(\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1)$ is a positive smooth infinitely differential function of all its parameters. For a fixed $\Omega$, let us group these parameters in two vectors $\mathbf{x} = (\mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1)$ and $\mathbf{y} = (T_z, T_x, \mathbf{a}_0)$. We then formally define

$$
\begin{align}
f(\mathbf{x}, \mathbf{y}) &= \mathcal{B}_{\Omega}(\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1) \nonumber \\
g(\mathbf{y}) &= \max_{\mathbf{x}} f(\mathbf{x}, \mathbf{y}) = \beta_{\max}(\mathbf{L}, \mathbf{a}_0, \Omega). \tag{172}
\end{align}
$$

We are interested in bounding $|dg(\mathbf{y})|$ as a function of $d\mathbf{y}$. Consider two values of the parameter, $\mathbf{y}_1$ and $\mathbf{y}_2$, and define

$$
\begin{align}
\bar{\mathbf{x}}_1 &= \operatorname*{argmax}_{\mathbf{x}} f(\mathbf{x}, \mathbf{y}_1) \nonumber \\
\bar{\mathbf{x}}_2 &= \operatorname*{argmax}_{\mathbf{x}} f(\mathbf{x}, \mathbf{y}_2), \tag{173}
\end{align}
$$

such that $g(\mathbf{y}_1) = f(\bar{\mathbf{x}}_1, \mathbf{y}_1)$ and $g(\mathbf{y}_2) = f(\bar{\mathbf{x}}_2, \mathbf{y}_2)$.
Without loss of generality we assume $g(\mathbf{y}_1) \ge g(\mathbf{y}_2)$
and consider the difference

$$
\begin{align}
|g(\mathbf{y}_1) - g(\mathbf{y}_2)| &= |f(\bar{\mathbf{x}}_1, \mathbf{y}_1) - f(\bar{\mathbf{x}}_2, \mathbf{y}_2)| \nonumber \\
&\le |f(\bar{\mathbf{x}}_1, \mathbf{y}_1) - f(\bar{\mathbf{x}}_1, \mathbf{y}_2)| \tag{174} \\
&\le \max_{\mathbf{x}} |f(\mathbf{x}, \mathbf{y}_1) - f(\mathbf{x}, \mathbf{y}_2)| \nonumber
\end{align}
$$

Taking the limit $\mathbf{y}_2 \to \mathbf{y}_1$ we get

$$
|dg(\mathbf{y})| \le \max_{\mathbf{x}} |\nabla_{\mathbf{y}} f(\mathbf{x}, \mathbf{y}) \cdot \mathbf{dy}|. \quad (175)
$$

Using the expression (66) for $\mathcal{B}_{\Omega}(\mathbf{L}, \phi, \mathbf{a}_1, \mathbf{b}_0, \mathbf{b}_1)$ we
then get

$$
|d\beta_{\max}(\mathbf{L}, \phi)| \le \begin{pmatrix} 1 \\ 1 \\ C_{\Omega} \end{pmatrix} \cdot \begin{pmatrix} |dT_z| \\ |dT_x| \\ |d\phi| \end{pmatrix}. \quad (176)
$$

B.3 Gradient of the goal function

We will now combine all the elements to upper bound
the gradient of the goal function

$$
G(\mathbf{L}, \phi; \Omega, q) = H(\rho_E) - H(E|\hat{a}_0(q)) + t\beta_{\max}(\mathbf{L}, \phi). \quad (177)
$$

We parametrize the vector **L** with the help of the angles as in Eq. (44), such that the “model” is described by four angles

$$
(\mathbf{L}, \phi) \approx \omega = (\alpha, \mu, \xi, \phi). \quad (178)
$$

We will bound the gradient of the goal function with
respect to this parametrization.

We start with the gradient of the angle. It satisfies

$$
dA = \begin{vmatrix} A(\mathbf{L}(\omega), \mathbf{L}(\omega + \mathbf{n}d\omega)) \\ |\nabla A(\omega)| \end{vmatrix} \quad (179)
$$
with a unit vector $\mathbf{n} = (n_\alpha, n_\mu, n_\xi, 0)$. For the fidelity one finds

$$F(\mathbf{L}(\omega), \mathbf{L}(\omega + \mathbf{n}d\omega)) = \sqrt{\mathbf{L}(\omega)} \cdot \sqrt{\mathbf{L}(\omega + \mathbf{d}\omega_i)} \quad (180) \\ \geq 1 - \frac{3}{4} d\omega^2 + O(d\omega^3)$$

for $d\omega_i = d\alpha, d\mu$ or $d\xi$, one easily finds

$$|\nabla A(\omega)| \leq \sqrt{\frac{3}{2}}. \quad (181)$$

We thus get from Eq. (171)

$$|d(H(\rho_E) - H(E|\hat{a}_0(q)))| \quad (182) \\ \leq 10.023 \sqrt{\frac{3}{2}} \sqrt{d\alpha^2 + d\mu^2 + d\xi^2} + 3d\phi$$

and

$$|\nabla(H(\rho_E) - H(E|\hat{a}_0(q)))| \quad (183) \\ \leq \sqrt{10.023^2 \frac{3}{2} + 9} < 12.7$$

For the gradient of the Bell score contribution $|\nabla\beta_{\max}(\mathbf{L}, \phi)|$ we note that

$$|dT_z|, |dT_x| \leq 2(|d\alpha| + |d\mu| + |d\xi|), \quad (184)$$

which implies from (176) that

$$|\nabla\beta_{\max}(\mathbf{L}, \phi)| \leq t\sqrt{3 \times 4^2 + C_{\Omega}^2} \leq 7t \quad (185)$$

Finally, for the goal function

$$|\nabla G| \leq 12.7 + 7t. \quad (186)$$

## B.4 Lipschitz function certified maximization on a compact space

**General approach** – In this section we explain how we obtain a certified maximum of a Lipschitz function $f: \mathbb{R}^n \to \mathbb{R}$ with a Lipschitz constant $\Lambda$, i.e. a bounded gradient $|\nabla f| \le \Lambda$, on a closed domain

$$\mathcal{D} = \{x = \begin{pmatrix} x_1 \\ \vdots \\ x_n \end{pmatrix} \in \mathbb{R}^n \mid x_i \in [\tau_i, \mu_i], \forall i\}. \quad (187)$$

We start by meshing $\mathcal{D}$ into an hypercube grid graph $G(s)$ with element of size $s$. If possible, we take $s$ so that $s$ evenly divides $\mu_i - \tau_i$ for all $i$. We label the center of each hypercube $\vec{c}$ so that each first neighbor in a given direction is separated by $s$, e.g. for $\vec{c}$ and its neighbor $\vec{c}'$ in the direction $+\vec{e}_1$

$$\vec{c}' = \begin{pmatrix} c'_1 \\ c'_2 \\ \vdots \\ c'_n \end{pmatrix} = \vec{c} + s.\vec{e}_1 = \begin{pmatrix} c_1 + s \\ c'_2 \\ \vdots \\ c'_n \end{pmatrix}. \quad (188)$$

A hypercube of center $\vec{c}$, $h(\vec{c})$, element of $G$, is thus defined as

$$h(\vec{c}) = \{x \in \mathbb{R}^n | c_i - \frac{s}{2} \le x_i \le c_i + \frac{s}{2}, \forall i \in \{1, \dots, n\} \}. \quad (189)$$

Given the Lipschitz constant $\Lambda$, the function $f$ in a given hypercube $h(\vec{c})$ is upper bounded by

$$f(x) \leq f(\vec{c}) + \frac{\sqrt{n}\Lambda}{2}s, \quad \forall x \in h(\vec{c}). \quad (190)$$

Trivially, an upper bound on the maximum of $f$ over $\mathcal{D}$ is given by,

$$\max_{x \in \mathcal{D}} f(x) \leq \max_{\vec{c} \in G(s)} f(\vec{c}) + \frac{\sqrt{n}\Lambda}{2}s. \quad (191)$$

The maximum of $f$ over $\mathcal{D}$ can thus be obtained by taking the smallest possible step, i.e.

$$\max_{x \in \mathcal{D}} f(x) = \lim_{s \to 0} \left( \max_{\vec{c} \in G(s)} \left( f(\vec{c}) + \frac{\sqrt{n}\Lambda}{2}s \right) \right). \quad (192)$$

**Numerical realisation** – While Eq. (192) gives the optimal maximum on $f$, it is however impracticable using numerical resources, due to the obvious need of discretization. A first naive approach is to relax the problem by setting a lower bound on $s$. This is, with the lower bound $s \ge \varepsilon$,

$$\max_{x \in \mathcal{D}} f(x) \leq \lim_{s \to \varepsilon} \left( \max_{\vec{c} \in G(s)} \left( f(\vec{c}) + \frac{\sqrt{n}\Lambda}{2}s \right) \right). \quad (193)$$

Resource wise, this method is cumbersome. Indeed, setting $\varepsilon$ to a small value compared to the domain space, i.e. $\varepsilon \ll \min_i(\tau_i - \mu_i)$, will result in a high number of hypercubes to explore. Furthermore, for a given step $s$, the number of hypercubes scales exponentially with the dimension $n$.

**Speed-up using a guess on the maximum** – Providing a guess on the maximum of $f$ may significantly speed up the previous method. Such a maximum can be found using an optimization algorithm without guarantees of optimality (BFGS, CMA, and others...). Denote such a maximum as $\nu$.

We start by setting a not-too-small step, $s_0$, so that the number of hypercubes composing the grid graph, $G_0(s_0)$, is reasonable. For each hypercube $h_0(\vec{c})$ we start by computing the potential maximum,

$$\xi^0(\vec{c}) = f(\vec{c}) + \frac{\sqrt{n}\Lambda}{2}s_0. \quad (194)$$

We then compare this value to $\nu$. In the case where $\nu > \xi^0(\vec{c})$, we pass to the next hypercube since the
guessed maximum is higher than the potential maximum value of $f$ in $h_0(\vec{c})$. Otherwise, we mesh $h_0(\vec{c})$ into a new hypercube grid graph $G_1(s_1)$ of element of size $s_1 = s_0/2$. For all of the new generated hypercubes, $h_1(\vec{c}) \in G_1(s_1)$, we compute the potential maximum $\xi^1(\vec{c})$ using Eq. (195). This method is applied recursively until either all new hypercubes of graph $G_m^k$ satisfy $\nu > h_m(\vec{c})$, or we reach the minimum step $s_m \le \varepsilon$. An upper bound of $f$ is thus given by the maximum $\xi^m(\vec{c})$

$$ \max(\xi^m(\vec{c})) = f(\vec{c}) + \frac{\sqrt{n\Lambda}}{2}s_m \quad (195) $$

where $s_m = s_0/2^m$.