# Device-independent lower bounds on the conditional von Neumann entropy

Peter Brown¹,³, Hamza Fawzi², and Omar Fawzi¹

¹Univ Lyon, ENS Lyon, UCBL, CNRS, Inria, LIP, F-69342, Lyon Cedex 07, France

²DAMTP, University of Cambridge, United Kingdom

³Télécom Paris, LTCI, Institut Polytechnique de Paris, 91120 Palaiseau, France

The rates of several device-independent (DI) protocols, including quantum key-distribution (QKD) and randomness expansion (RE), can be computed via an optimization of the conditional von Neumann entropy over a particular class of quantum states. In this work we introduce a numerical method to compute lower bounds on such rates. We derive a sequence of optimization problems that converge to the conditional von Neumann entropy of systems defined on general separable Hilbert spaces. Using the Navascués-Pironio-Acín hierarchy we can then relax these problems to semidefinite programs, giving a computationally tractable method to compute lower bounds on the rates of DI protocols. Applying our method to compute the rates of DI-RE and DI-QKD protocols we find substantial improvements over all previous numerical techniques, demonstrating significantly higher rates for both DI-RE and DI-QKD. In particular, for DI-QKD we show a minimal detection efficiency threshold which is within the realm of current capabilities. Moreover, we demonstrate that our method is capable of converging rapidly by recovering all known tight analytical bounds up to several decimal places. Finally, we note that our method is compatible with the entropy accumulation theorem and can thus be used to compute rates of finite round protocols and subsequently prove their security.

## 1 Introduction

Quantum cryptography enables certain cryptographic tasks to be performed securely with their security guaranteed by the physical laws of nature as opposed to assumptions of computational hardness that are used in more conventional cryptography. However, standard quantum cryptography protocols, for instance BB84 [1], still require a strong level of trust in the hardware used and its implementation. Faulty hardware or malicious attacks on the implementation can compromise the security of these protocols, rendering them useless [2]. Whilst better security analyses and improved hardware checks can make it more difficult for an adversarial party to eavesdrop on the protocol, it is always possible that there are unknown side-channel attacks that remain exploitable.

Fortunately, quantum theory also offers a way to remove strong assumptions on the hardware and implementation by instead running a so-called device-independent protocol. Device-independent protocols offer the pinnacle of security guarantees. By relying on minimal assumptions, they remain secure even when the devices used within the protocol are completely untrusted. The central idea behind many device-independent protocols, including randomness expansion (RE) and quantum key distribution (QKD), is that there are certain correlations between multiple separate systems that (i) could only have been produced by entangled quantum systems and (ii) are intrinsically random. Intuitively, these nonlocal correlations then act as a certificate that guarantees the systems produced randomness [3]. Furthermore, nonlocal correlations have additional applications such as certifying the dimension of a system [4], reducing communication complexity [5] and in

<footnote>Peter Brown: peter.brown@telecom-paris.fr</footnote>
certain cases nonlocal correlations can even certify the exact systems used (up to some unavoidable symmetries) [6].

Beginning with the works of [7, 8], significant effort has been placed into developing new device-independent protocols and subsequently proving their security. For tasks like DI-RE and DI-QKD, security is now well understood, with tools like the entropy accumulation theorem (EAT) [9, 10] and quantum probability estimation (QPE) [11] enabling relatively simple proofs of security against all-powerful quantum adversaries [12, 13]. In both cases a security proof can be readily established once one has developed a quantitative relationship between the nonlocal correlations observed in the protocol and the quantity of randomness generated by the devices used. In particular, one is required to understand the *rate* of a protocol, i.e. the number of random bits or key bits that are generated per round. The asymptotic rates of both DI-RE and DI-QKD protocols are quantified in terms of the conditional von Neumann entropy. Moreover, using the EAT, one can also bound the total randomness produced in a finite round protocol in terms of the conditional von Neumann entropy.<fnref n="1" />

A central problem in device-independent cryptography then remains of how one actually computes this conditional von Neumann entropy? More specifically, one is required to compute (or lower bound) the minimum conditional von Neumann entropy produced by the devices in a single round of the protocol conditioned on the adversary's side information. This computation should be device-independent in the sense that there are no restrictions on what systems were used except for some constraints on the expected correlations produced by the parties' devices, e.g., a Bell-inequality violation. For the purpose of the introduction, we use the following informal notation $inf_ρ H[ρ]$ for this quantity of interest, where the infimum is understood to be over all states that are compatible with the observed correlations, and $H$ refers to the conditional von Neumann entropy. We refer to Section 2.2 for a more formal definition. At first glance this optimization would appear very challenging: the optimization itself is non-convex and furthermore we place no restrictions on the dimensions of the quantum systems used. Nevertheless, in spite of these apparent difficulties, significant progress has still been made towards solving such an optimization problem. For example, when the devices in the protocol have only binary inputs and outputs, it is often possible to exploit Jordan's lemma [14] in order to restrict the optimization to qubit systems. Using this technique, the work of [15] was able to give an exact analytical solution to the optimization problem for the one-sided conditional von Neumann entropy when the devices were constrained by a violation of the CHSH inequality [16]. More recently, other works have also managed to obtain analytical lower bounds on devices satisfying MABK inequality violations [17], tight bounds for the Holz inequality [18] and tight bounds for the asymmetric CHSH inequality [19]. These analytical results are difficult to obtain and also rely on the reduction to qubit systems, meaning that it is not possible to perform such an analysis when the inputs and outputs of the devices are not binary. As such, general numerical techniques were also developed that could handle more complicated protocols and general linear constraints on the expected correlations.

For example, in [20, 21] it was shown that the analogous optimization for the min-entropy [22] can be straightforwardly lower bounded using the Navascués-Pironio-Acín (NPA) hierarchy [23, 24]. As the min-entropy is never larger than the von Neumann entropy the results of these computations give lower bounds on the rates. However, one major drawback of this method is that the resulting lower bounds are in general quite loose. More recently, two other proposals for general numerical approaches have been developed. In [25], the authors give an explicit method to lower bound the conditional von Neumann entropy via a non-commutative polynomial optimization problem that can then be lower bounded using the NPA hierarchy. In [26], the authors introduced new entropies that are all lower bounds on the conditional von Neumann entropy and which can, like the min-entropy, be optimized using the NPA hierarchy. Both works improve upon the min-entropy method but neither has been shown to give tight bounds on the actual rate of a protocol and in general there appears to be significant room for improvement. As such, the question remains as to whether one can give a computationally tractable method to compute tight lower bounds on the rates of protocols.

<footnote n="1">¹For the QPE one must compute a related quantity, namely construct a quantum estimation factor but we leave such investigations to future work.</footnote>
## 1.1 Contributions of this work

In this work, we introduce a new method to solve this problem. To achieve this, we develop a family of variational expressions indexed by a positive integer $m$ that approximate the conditional von Neumann entropy to arbitrary accuracy as $m$ grows. For any $m$, the variational expression is determined by a noncommutative polynomial $P_m$ and for a state given by the density operator $\rho$ takes the form

$$H_m[\rho] := \inf_{Z_1, \dots, Z_n} \mathrm{Tr}[\rho P_m(Z_1, \dots, Z_n, Z_1^*, \dots, Z_n^*)], \quad (1)$$

where $Z_1, \dots, Z_n$ are bounded operators. For any choice of $m$, this expression is a lower bound on the von Neumann entropy i.e., $H[\rho] \ge H_m[\rho]$, and in particular can be used to lower bound the rates of DI protocols. Such variational expressions are closely related to Kosaki's expressions for the quantum relative entropy [27, 28, 29]. However, Kosaki's variational expression has a number of variables ($n$ in (1)) that is infinite and thus it is not suitable for computations. The variational expressions we obtain can be seen as well-chosen finite approximations of Kosaki's formula. Note that we need to be careful in the choice of approximation as we need to obtain expressions which always give lower bounds on the conditional von Neumann entropy.

The point of obtaining expressions of the form (1) is that when further taking the infimum over all states $\rho$ that are compatible with the observed correlations, it is possible to relax such an optimization problem to a semidefinite program (SDP) using the NPA hierarchy to yield a computationally tractable lower bound. Moreover, we show that for any fixed $m$, the resulting noncommutative polynomial optimization problem can be made to satisfy the property that all the variables have a fixed upper bound on their operator norm. As a consequence, we obtain that the sequence of SDPs given by the NPA hierarchy converges in the limit to $\inf_{\rho} H_m[\rho]$ where $\rho$ is compatible with the observed statistics. We also remark that unlike previous numerical techniques that were developed, our analysis applies when the systems are also infinite dimensional.

We then apply this to compute the rates of both DI-RE and DI-QKD protocols. Compared to the other numerical techniques we show significant improvements on the calculated rates. We give improved bounds on the DI randomness certifiable with qubit systems which could be used to yield more efficient experiments for DI randomness [30, 31]. In addition we also give new lower bounds on the minimal detection efficiency required to perform DI-QKD with qubit systems. This gives a promising approach to conduct DI-QKD experiments with current technologies. We also find that in practice our method can converge quickly as we demonstrate that we can recover (up to several decimal places) tight analytical bounds from [15], [18] and [19]. We also remark that the computations we ran in several of our examples were also vastly more efficient than the two numerical approaches of [25, 26]. Finally, we also explain how our technique can be used directly with the entropy accumulation theorem in order to compute non-asymptotic rates of protocols and subsequently prove their security.

The article is structured as follows. In Section 2 we begin by stating the main results relevant to device-independent cryptography and their application. Note that in this section we will restrict to a special case of our more general result that is sufficient for those interested solely in applications. Later in Section 3 we will state the general technical result along with its proof.

## 2 Results

Before we begin let us establish some notation. Let $H$ be a separable Hilbert space, we denote the set of bounded operators on $H$ to itself by $B(H)$. A state on $H$ is a positive semidefinite, trace-class operator $\rho$ such that $\mathrm{Tr}[\rho] = 1$. We denote the set of states on $H$ by $D(H)$. Given two positive semidefinite operators $\rho, \sigma$ on $H$ we write $\rho \ll \sigma$ if $\mathrm{ker} \sigma \subseteq \mathrm{ker} \rho$ where $\mathrm{ker} X := \{|v\rangle \in H : X|v\rangle = 0\}$. Throughout this work $\mathbb{R}_+$ will denote the set of nonnegative real numbers, $\mathbb{N}$ will denote the set of strictly positive integers, $\ln$ will denote the natural logarithm and $\log_2$ will denote the logarithm base two. In this section, for simplicity of presentation, we assume that $H$ is finite-dimensional. See Section 3 for the setting where $H$ is infinite-dimensional.
device. In a device-independent protocol such as DI-QKD or DI-RE the statistics of many rounds will be analyzed in order to determine whether or not randomness is being produced from the devices or whether a secret key can be distilled.

We assume that the devices are constrained by quantum theory and that they act in the following way. Let $Q_A$, $Q_B$ and $Q_E$ be the three separable Hilbert spaces of Alice's device, Bob's device and Eve's device respectively. At the beginning of each round a tripartite state $\rho_{Q_A Q_B Q_E}$ on $Q_A \otimes Q_B \otimes Q_E$ is shared between the three systems. In response to an input, Alice and Bob's devices will measure some preselected POVMs $\{\{M_a|x\}_a\}_x$, $\{\{N_b|y\}_b\}_y$ on the parts of the state that they received and return the measurement outcome. The state and measurements are unknown to Alice and Bob but may be known by Eve. Overall, the joint conditional distribution of their outputs may be described via the Born rule as

$$p(a, b|x, y) = \mathrm{Tr} \left[ \rho_{Q_A Q_B Q_E} (M_{a|x} \otimes N_{b|y} \otimes I) \right]. \qquad (6)$$

Immediately after the round, the joint state between the classical information recorded by the honest parties and Eve's quantum system may be described by the cq-state

$$\rho_{ABXYQ_E} = \sum_{abxy} p(x, y) |abxy\rangle\langle abxy| \otimes \rho_{Q_E}^{abxy}, \qquad (7)$$

where

$$\rho_{Q_E}^{abxy} = \mathrm{Tr}_{Q_A Q_B} [\rho_{Q_A Q_B Q_E} (M_{a|x} \otimes N_{b|y} \otimes I)] . \qquad (8)$$

In spot-checking protocols, secret key and randomness are only extracted from rounds with a particular input [33]. As such, we henceforth consider some distinguished inputs $(x^*, y^*)$ and the post-measurement state $\rho_{ABQ_E} = \sum_{ab} |ab\rangle\langle ab| \otimes \rho_{Q_E}^{abxy}$. The exact choices of $(x^*, y^*)$ will be made explicit when we compute rates for given protocols. For a spot-checking DI-RE protocol, the asymptotic rate is characterized by

$$H(AB|X = x^*, Y = y^*, Q_E) \qquad (9)$$

if the randomness is extracted from both Alice and Bob's outputs and

$$H(A|X = x^*, Q_E) \qquad (10)$$

if the randomness is only extracted from the outputs of Alice's device. For a spot-checking DI-QKD protocol with one-way error correction, the asymptotic rate is given in terms of the Devetak-Winter bound [34]

$$H(A|X = x^*, Q_E) - H(A|B, X = x^*, Y = y^*). \qquad (11)$$

To obtain the asymptotic<fnref n="4" /> rate of the device-independent protocols we must consider a minimization of the above quantities over all possible quantum systems that could have produced the statistics we observed on expectation. More formally, let $\mathcal{A}, \mathcal{B}, \mathcal{X}, \mathcal{Y}$ denote the finite sets from which the random variables $A, B, X, Y$ take their values. Then given a finite set $C$ let $C : \mathcal{A}\mathcal{B}\mathcal{X}\mathcal{Y} \to C$ be some function which will act as a statistical test on $ABXY$. Finally, we consider a probability distribution $q(c)$ on $C$. We then say a distribution $p(a, b, x, y)$ on $ABXY$ is compatible with the pair $(C, q)$ if

$$\sum_{(a,b,x,y): C(a,b,x,y)=c} p(a, b, x, y) = q(c). \qquad (12)$$

In other words, when we apply our statistical test $C$ to the random variables $ABXY$ we obtain a new random variable whose distribution is $q$. For example, let $\mathcal{A} = \mathcal{B} = \mathcal{X} = \mathcal{Y} = C = \{0, 1\}$,

$$C(a, b, x, y) = \begin{cases} 1 & \text{if } a \oplus b = xy \\ 0 & \text{otherwise} \end{cases} \qquad (13)$$

<footnote n="4">⁴For protocols with a finite number of rounds a more careful analysis of the rate is required. However general techniques like the entropy accumulation theorem [9, 10] allow us to lift the asymptotic rate to the rate of a protocol with a finite number of rounds.</footnote>
and let $q(1) = \omega$. Then for $p(x, y) = 1/4$ the pair $(C, q)$ imposes the constraint that the distribution $p(a, b, x, y)$ should achieve an expected score of $\omega$ in the CHSH game.

More generally we say that a tuple $(Q_A Q_B Q_E, \rho, \{M_{a|x}\}, \{N_{b|y}\})$ is compatible with the constraints $(C, q)$ if the probability distribution $p(a, b, x, y) = p(x, y) \text{Tr}[\rho(M_{a|x} \otimes N_{b|y} \otimes I)]$ is compatible with $(C, q)$. We refer to such a tuple as a *strategy*. For each strategy we can also associate a post-measurement state via (7). Then for a fixed $(C, q)$ the *device-independent* rate of randomness produced by Alice's device on input $X = x^*$, if the devices satisfy the constraints imposed by $(C, q)$, is given by

$$ \inf H(A|X = x^*, Q_E) \qquad (14) $$

where the infimum is taken over all strategies compatible with $(C, q)$ and the conditional von Neumann entropy is evaluated on the post-measurement state induced by the strategy. The rates for the randomness produced by both devices and the device-independent Devetak-Winter rate can be defined analogously by replacing the objective function with the appropriate quantity. Note that by considering appropriate dilations we can restrict the optimization to strategies wherein the measurements are projective and the state $\rho_{Q_A Q_B Q_E}$ is pure.

The following lemma demonstrates how to use the upper bounds on the relative entropy from Theorem 2.1 in order to lower bound this infimum by a converging sequence of optimizations that can be subsequently lower bounded using the NPA hierarchy.

**Lemma 2.3.** Let $m \in \mathbb{N}$ and let $t_1, \dots, t_m$ and $w_1, \dots, w_m$ be the nodes and weights of an m-point Gauss-Radau quadrature on $[0, 1]$ with an endpoint $t_m = 1$, as specified in Theorem 3.8. Let $\rho_{Q_A Q_B Q_E}$ be the initial quantum state shared between the devices of Alice, Bob and Eve and let $\{M_{a|x^*}\}$ denote the measurements operators performed by Alice's device in response to the input $X = x^*$. Furthermore for $i = 1, \dots, m-1$ let $\alpha_i = \frac{3}{2} \max\{\frac{1}{t_i}, \frac{1}{1-t_i}\}$. Then $H(A|X = x^*, Q_E)$ is never smaller than

$$ c_m + \sum_{i=1}^{m-1} \frac{w_i}{t_i \ln 2} \sum_a \inf_{Z_a \in B(Q_E)} \text{Tr} \left[ \rho_{Q_A Q_E} \left( M_{a|x^*} \otimes (Z_a + Z_a^* + (1-t_i)Z_a^* Z_a) + t_i (I_{Q_A} \otimes Z_a Z_a^*) \right) \right] \\ \text{s.t.} \quad \|Z_a\| \le \alpha_i \qquad (15) \\ \text{where } c_m = \sum_{i=1}^{m-1} \frac{w_i}{t_i \ln 2}. \text{ Moreover these lower bounds converge to } H(A|X = x^*, Q_E) \text{ as } m \to \infty. $$

*Proof.* Let $\rho_{AQ_E} = \sum_a |a\rangle\langle a| \otimes \rho_{QE}(a, x^*)$ be the cq-state after Alice has performed her measurement corresponding to the input $x^*$, i.e., $\rho_{QE}(a, x^*) = \text{Tr}_{QAQ_B} [\rho_{QAQ_BQ_E}(M_{a|x^*} \otimes I \otimes I)]$. Then for any $m \in \mathbb{N}$, using the relation $H(A|X = x^*, Q_E) = -D(\rho_{AQ_E} \| I_A \otimes \rho_{QE})$ and Theorem 2.1 we have that $H(A|X = x^*, Q_E)$ is never smaller than

$$ c_m + \sum_{i=1}^{m-1} \frac{w_i}{t_i \ln 2} \inf_{Z \in B(H)} \text{Tr} [\rho_{AQ_E}(Z + Z^* + (1-t_i)Z^*Z)] + t_i \text{Tr} [(I_A \otimes \rho_{QE})ZZ^*] \qquad (16) \\ \text{s.t.} \quad \|Z\| \le \frac{3}{2} \max \left\{ \frac{1}{t_i}, \frac{\lambda}{1-t_i} \right\} $$

where $c_m = -\frac{\lambda}{m^2 \ln 2} + \sum_{i=1}^m \frac{w_i}{t_i \ln 2}$ and $\lambda$ is some real number such that $\rho_{AQ_E} \le \lambda I_A \otimes \rho_{QE}$. Now as system A is finite dimensional, we can write the operator $Z = \sum_{ab} |a\rangle\langle b| \otimes Z_{(a,b)}$ for some operators $Z_{(a,b)} \in B(Q_E)$. Then for the first term we have

$$ \begin{aligned} \text{Tr} [\rho_{AQ_E}(Z + Z^*)] &= \sum_a \text{Tr} [\rho_{QE}(a, x^*)(Z_{(a,a)} + Z_{(a,a)}^*)] \\ &= \sum_a \text{Tr} [\text{Tr}_{QA} [\rho_{QAQ_E}(M_{a|x^*} \otimes I_{Q_E})] (Z_{(a,a)} + Z_{(a,a)}^*)] \qquad (17) \\ &= \sum_a \text{Tr} [\rho_{QAQ_E}(M_{a|x^*} \otimes (Z_{(a,a)} + Z_{(a,a)}^*))] \end{aligned} $$
where on the first line we traced out the A system, on the second line we substituted in the definition of $\rho_{AQ_E}(a, x^*)$ and on the third line we used the identity $\text{Tr}_A[X_{AB}(I \otimes Y_B)] = \text{Tr}_A[X_{AB}]Y_B$. Repeating this for the second term we find

$$ \begin{aligned} \text{Tr} [\rho_{AQ_E} Z^* Z] &= \sum_{ab} \text{Tr} [\rho_{Q_E}(a, x^*) Z^*_{(b,a)} Z_{(b,a)}] \\ &\geq \sum_{a} \text{Tr} [\rho_{Q_E}(a, x^*) Z^*_{(a,a)} Z_{(a,a)}] \qquad (18) \\ &= \sum_{a} \text{Tr} [\rho_{Q_A Q_E}(M_{a|x^*} \otimes Z^*_{(a,a)} Z_{(a,a)})] \end{aligned} $$

where on the second line we noted that $\sum_b Z^*_{(b,a)} Z_{(b,a)} \geq Z^*_{(a,a)} Z_{(a,a)}$. Finally we get a similar relation for the final term

$$ \text{Tr}[(I_A \otimes \rho_{Q_E})ZZ^*] \geq \sum_a \text{Tr} [\rho_{Q_A Q_E}(I_{Q_A} \otimes Z_{(a,a)} Z^*_{(a,a)})]. \qquad (19) $$

Inserting these three rewritings into the lower bound on $H(A|X = x^*, Q_E)$ and relabelling $Z_{(a,a)}$ to $Z_a$ we recover the objective function stated in the lemma. Note that the above rewritings and the fact that we are minimizing implies that we need only consider operators Z that are block diagonal in the sense that $Z = \sum_a |a\rangle\langle a| \otimes Z_{(a,a)}$.

As $\rho_{AQ_E}$ is a cq-state we have $\rho_{AQ_E} \leq I_A \otimes \rho_{Q_E}$ and we can set $\lambda = 1$. This recovers the form of $c_m$ stated in the lemma (noting that $w_m = \frac{1}{m^2}$ and $t_m = 1$). As it is sufficient to take $Z = \sum_a |a\rangle\langle a| \otimes Z_a$, we must have $\alpha_i = \frac{3}{2} \max\{\frac{1}{t_i}, \frac{\lambda}{1-t_i}\} \geq \|\sum_a |a\rangle\langle a| \otimes Z_a\| = \max_a \|Z_a\|$. Finally the convergence statement follows immediately from the convergence proven in Theorem 3.11. $\square$

**Remark 2.4 (Adapting to other entropies).** The above lemma only describes a bound for $H(A|X = x^*, Q_E)$. However the proof can be easily adapted to the case of the global entropy $H(AB|X = x^*, Y = y^*, Q_E)$ or for non-fixed inputs, e.g., $H(A|XQ_E)$. For example, the global entropy $H(AB|X = x^*, Y = y^*, Q_E)$ can be lower bounded by replacing the inner summation in (15) with

$$ \sum_{ab} \inf_{Z_{ab} \in B(Q_E)} \text{Tr} \left[ \rho_{Q_A Q_B Q_E} \left( M_{a|x^*} \otimes N_{b|y^*} \otimes (Z_{ab} + Z^*_{ab} + (1-t_i)Z^*_{ab} Z_{ab}) + t_i (I_{Q_A} \otimes Z_{ab} Z^*_{ab}) \right) \right]. \qquad (20) $$

Similarly one could also adapt the proof to bound $H(A|XQ_E)$, allowing the entropy to averaged over the inputs X. For this one would replace the inner summation with

$$ \sum_{ax} \inf_{Z_{ax} \in B(Q_E)} p(x) \text{Tr} \left[ \rho_{Q_A Q_E} \left( M_{a|x} \otimes (Z_{ax} + Z^*_{ax} + (1-t_i)Z^*_{ax} Z_{ax}) + t_i (I_{Q_A} \otimes Z_{ax} Z^*_{ax}) \right) \right]. \qquad (21) $$

The above lemma and remark provide a converging sequence of lower bounds on the conditional von Neumann entropy. In order to turn these into lower bounds on the rate of a device-independent protocol we must also include the optimizations of all states, measurements and Hilbert spaces, subject to any constraints on the devices' joint probability distribution that we wish to impose. Suppose for $1 \leq j \leq r$ for some $r \in \mathbb{N}$ we impose on Alice and Bob's devices a collection of constraints

$$ \sum_{abxy} c_{abxyj} p(a, b|x, y) \geq v_j \qquad (22) $$

where $c_{abxyj}, v_j \in \mathbb{R}$. Then using Lemma 2.3 we can compute a lower bound on $H(A|X = x^*, Q_E)$ for all possible devices that satisfy the above constraints by solving the following optimization
problem:

$$
\begin{align}
c_m + \inf \quad & \sum_{i=1}^{m-1} \frac{w_i}{t_i \ln 2} \sum_a \operatorname{Tr} \left[ \rho_{Q_A Q_B Q_E} \left( M_{a|x^*} \otimes I_{Q_B} \otimes (Z_{a,i} + Z_{a,i}^* + (1-t_i)Z_{a,i}^* Z_{a,i}) + t_i (I_{Q_A Q_B} \otimes Z_{a,i} Z_{a,i}^*) \right) \right] \nonumber \\
\text{s.t.} \quad & \sum_{abxy} c_{abxyj} \operatorname{Tr} \left[ \rho_{Q_A Q_B Q_E} (M_{a|x} \otimes N_{b|y} \otimes I) \right] \ge v_j && \text{for all } 1 \le j \le r \nonumber \\
& \sum_a M_{a|x} = \sum_b N_{b|y} = I && \text{for all } x, y \nonumber \\
& M_{a|x} \ge 0 && \text{for all } a, x \nonumber \\
& N_{b|y} \ge 0 && \text{for all } b, y \nonumber \\
& \|Z_{a,i}\| \le \alpha_i && \text{for all } a, i = 1, \dots, m-1 \nonumber \\
& M_{a|x} \in B(Q_A), \quad N_{b|y} \in B(Q_B), \quad Z_{a,i} \in B(Q_E) && \text{for all } a, b, x, y, i \nonumber \\
& \rho_{Q_A Q_B Q_E} \in D(Q_A Q_B Q_E) && \tag{23}
\end{align}
$$

where the infimum is taken over the all collections $(Q_A Q_B Q_E, \rho_{Q_A Q_B Q_E}, \{M_{a|x}\}, \{N_{b|y}\}, \{Z_{a,i}\})$
satisfying the constraints of the problem. Note that when rearranging the objective function we
have used the fact that the inner summation in (15) commutes with the infimum over the $Z_a$
operators. We can then further pull the infimum outside the outer summation by reparametrizing
the variables as $Z_{a,i}$ for each $i$ in the outer sum. In order to compute a lower bound on (23)
we employ the NPA hierarchy to relax this problem to an SDP. To do this we first drop the
tensor product structure and instead include commutation relations on the relevant variables of
the problem. In doing so we end up with the following noncommutative polynomial optimization
problem which gives a lower bound on (23)

$$
\begin{align*}
c_m + \inf \quad & \sum_{i=1}^{m-1} \frac{w_i}{t_i \ln 2} \sum_a \left\langle \psi | M_{a|x^*} (Z_{a,i} + Z_{a,i}^* + (1-t_i)Z_{a,i}^* Z_{a,i}) + t_i Z_{a,i} Z_{a,i}^* | \psi \right\rangle \\
\text{s.t.} \quad & \sum_{abxy} c_{abxyj} \left\langle \psi | M_{a|x} N_{b|y} | \psi \right\rangle \ge v_j && \text{for all } 1 \le j \le r \\
& \sum_a M_{a|x} = \sum_b N_{b|y} = I && \text{for all } x, y \\
& M_{a|x} \ge 0 && \text{for all } a, x \\
& N_{b|y} \ge 0 && \text{for all } b, y \\
& Z_{a,i}^* Z_{a,i} \le \alpha_i^2 && \text{for all } a, i = 1, \dots, m-1 \\
& Z_{a,i} Z_{a,i}^* \le \alpha_i^2 && \text{for all } a, i = 1, \dots, m-1 \\
& [M_{a|x}, N_{b|y}] = [M_{a|x}, Z_{b,i}^{(*)}] = [N_{b|y}, Z_{a,i}^{(*)}] = 0 && \text{for all } a, b, x, y, i \\
& M_{a|x}, N_{b|y}, Z_{a,i} \in B(H) && \text{for all } a, b, x, y, i \tag{24}
\end{align*}
$$

where we have recalled that it is sufficient to consider pure states. Note that in both of the above
optimizations we can also include projective measurement constraints without loss of generality.
Using the NPA hierarchy we can then relax this optimization to a sequence of SDPs that give us a
converging sequence of lower bounds on the optimal value and in turn a lower bound on the rate
of the protocol.

**Remark 2.5** (Commuting operator versus tensor product strategies). It is immediate that (24) is never larger than (23) and thus our subsequent relaxations of (24) will always give lower bounds on the rates of protocols as we defined them previously using the tensor product framework. However, due to recent work [35] it may not be the case that (23) and (24) are equal.

Let $R_m$ be the optimal value in (24) and $r_{m,k}$ be the optimal value of its $k^{\text{th}}$ level NPA relaxation. As we have explicit bounds on the operator norms of the variables we know that our NPA relaxations of (24) will converge to the optimal value of (24), i.e., $R_m = \lim_{k \to \infty} r_{m,k}$. We believe that $\sup_m R_m$ will correspond to infimum of the conditional von Neumann entropy for commuting
operator strategies. However proving this would require one to formally define commuting operator strategies (analogous to the tensor product strategies introduced earlier) and checking that (24) can be derived in a similar fashion to how (23) was derived. We leave this to future work.

**Remark 2.6** (Faster lower bounds). There are several ways to speed up the SDP relaxations of (24). We will also note in the caption of each figure which speedups were used.

1. Often including the operator inequalities $Z_{a,i}^*Z_{a,i} \le \alpha_i^2$ and $Z_{a,i}^*Z_{a,i} \le \alpha_i^2$ does not improve the lower bound. Hence we always removed them when performing the computations.

2. The choice of monomial indexing set for the moment matrices can greatly affect the accuracy and speed of the computations. We found that the local level 1 set, i.e., monomials of the form $ABZ$ where $A \in \{I\} \cup \{M_{a|x}\}_{a,x}$, $B \in \{I\} \cup \{N_{b|y}\}_{b,y}$ and $Z \in \{I\} \cup \{Z_{c,i}, Z_{c,i}^*\}_{c,i}$ performed particularly well.

3. It is also possible to commute the outer summation and the infimum, i.e., compute

$$ \sum_{i=1}^{m-1} \frac{w_i}{t_i \ln 2} \inf \sum_a \langle \psi | M_{a|x^*} (Z_{a,i} + Z_{a,i}^* + (1-t_i)Z_{a,i}^* Z_{a,i}) + t_i Z_{a,i} Z_{a,i}^* | \psi \rangle. \quad (25) $$

This can only decrease the lower bound on the entropy and hence it is sufficient for the purpose of lower bounding the rate of a protocol. The main advantage of doing so is that it drastically reduces the number of variables in the NPA hierarchy relaxations. Rather than running a single SDP with a $Z_{a,i}$ variable for each value of $a$ and $i$, we can instead run $m$ much smaller SDPs with only a $Z_{a,i}$ variable for each $a$. This significantly reduces the runtime of the SDPs and results in the runtime scaling linearly with the number of nodes in the Gauss-Radau quadrature. However we did notice in certain cases that the lower bounds computed in this manner were not converging to tight lower bounds. In such cases we did not include this speedup.

4. The moment matrix of the NPA relaxation can without loss of generality be taken to be a real symmetric matrix.

## 2.3 Numerical results

We will now apply our method of computing rates to several DI-RE and DI-QKD scenarios and compare our technique with known analytical results [15, 18, 19] and other numerical techniques [25, 26, 36]. We will only concern ourselves with the asymptotic rates in this work. However, as noted earlier, our technique can be combined with the entropy accumulation theorem in a relatively straightforward manner, similar to [37], and thus one could also use it to compute finite round rates for protocols. We discuss this further in Appendix E. The semidefinite relaxations were generated using the python package NCPOL2SDPA [38] and the resulting SDPs were solved using MOSEK [39]. As NCPOL2SDPA is no longer maintained, we used a maintained fork of the original package [40]. We also provide example python scripts that implement some of the computations, these can be found at the github repository [41].

### 2.3.1 Randomness from the CHSH game

To begin we consider the simplest possible setting of bounding $H(A|X = 0, Q_E)$ when Alice and Bob's devices are constrained to achieve some minimal score in the CHSH game. In Figure 1 we demonstrate how our bounds improve as we increase the number of nodes $m$ in the Gauss-Radau quadrature. We compare this with a known tight analytical bound in this setting [15]. In the figure we see that for $m = 8$ our numerical technique effectively recovers the known tight analytical bound. As far as we are aware this is the first numerical technique to do so without resorting to algebraic simplifications afforded by Jordan's lemma. Importantly, this also demonstrates that in certain settings our technique can converge very quickly in the NPA hierarchy and in the size of the Gauss-Radau quadrature. Furthermore, as we are able to run our computations at low levels of the NPA hierarchy, the computations are also relatively fast for this setting with each SDP taking less than a second to run. We include additional plots demonstrating the convergence of our technique for
<figure><img src="image_1.png" /><figcaption>Figure 1: Recovering the local CHSH bound. Comparison of lower bounds on $H(A|X=0, Q_E)$ for quantum devices that constrained to achieve some minimal CHSH score. Numerical bounds were computed using speedups (1), (3) and (4) from Remark 2.6 at a relaxation level $2+ABZ+AZZ$. A single SDP took less than one second to run.</figcaption></figure>

other known tight analytical bounds [18, 19] in Appendix D. These include multipartite scenarios useful for bounding the rates of DI conference key agreement protocols [42].

For randomness expansion protocols it is beneficial to consider the randomness generated from both devices. Therefore, in Figure 2 we plot lower bounds on $H(AB|X=0, Y=0, Q_E)$. We compare this again with the analytical bound on $H(A|X=0, Q_E)$ and we also compare it with the numerical technique of [26] and a known upper bound from [36]. As is clear from the figure, the global entropy can be significantly larger than the local entropy leading to much better rates for randomness expansion protocols based on the CHSH game. In comparison with the numerical technique of [26], we see that our new method vastly outperforms it. In the plot we also compare with some numerical upper bounds from [36] which used Jordan's lemma to reduce the problem to qubits and then optimized over explicit qubit strategies. As we see from the plot our rate curve almost coincides with this upper bound and we expect that by increasing $m$ further the gap would be further reduced. Note that we could also have used speedup (3) from Remark 2.6 in this case which would change the runtime from hours to minutes. However, we found that when using speedup (3) we were unable to recover the bound from [36] and so we elected to not use it.

### 2.3.2 Randomness from the full distribution

By increasing our knowledge about the conditional distribution with which Alice and Bob's devices operate we further constrain the possible devices that could produce the statistics we observe in the protocol. In turn we can hope that this leads to larger bounds on the randomness produced by the devices. At the extreme end of this, we would have knowledge of the entire distribution characterizing the devices. In the following we assume that we have access to the complete distribution and we impose these constraints in the SDP.

To make the plots more experimentally relevant we also assume the devices are affected by inefficient detection events. That is, with some probability $\eta \in [0,1]$ the devices will operate normally and with probability $1-\eta$ the devices will fail and deterministically output 0. The conditional
<figure><img src="image_2.png" /><figcaption>Figure 2: Global randomness from the CHSH game. Comparison of lower bounds on $H(AB|X = 0, Y = 0, Q_E)$ for quantum devices that constrained to achieve some minimal CHSH score. Our technique was computed using speedups (1) and (2) from Remark 2.6. For $m = 8$ a single data point can take hours to run. We also compare with the iterated mean entropy $H_{(4/3)}^*(AB|X = 0, Y = 0, Q_E)$ from [26] and an upper bound from [36].</figcaption></figure>

distribution for such devices takes the form

$$p(a, b|x, y) = \begin{cases} \eta^2 q(a, b|x, y) + \eta(1-\eta)(q(a|x) + q(b|y)) + (1-\eta)^2 & \text{if } a = 0 = b \\ \eta^2 q(a, b|x, y) + \eta(1-\eta)q(a|x) & \text{if } a \neq 0 = b \\ \eta^2 q(a, b|x, y) + \eta(1-\eta)q(b|y) & \text{if } a = 0 \neq b \\ \eta^2 q(a, b|x, y) & \text{otherwise} \end{cases} \quad (26)$$

where $q(a, b|x, y)$ is the conditional distribution of the devices when $\eta = 1$.

In Figure 3 we compare lower bounds on the entropy $H(AB|X = 0, Y = 0, Q_E)$ when Alice and Bob's devices have inefficient detectors. For the curve representing our technique we used a Gauss-Radau decomposition with $m = 8$ nodes and at each data point we selected a two-qubit distribution in order to maximize the entropy produced by the devices. We compare our technique with the iterated mean entropy $H_{(2)}^*(AB|X = 0, Y = 0, Q_E)$ from [26] and the TSGPL method from [25] both of which also constrained the devices by their full distribution. Compared with the other methods, our technique is everywhere larger and on the whole the difference is substantial. We also note that our technique is again significantly faster than the other numerical techniques presented here.

### 2.3.3 Better bounds on DI-QKD key rates

Thus far we have only concerned ourselves with bounds on local and global entropies. However, we can use these bounds to compute the asymptotic rates of some DI-QKD protocols. The asymptotic rate of a DI-QKD protocol with one-way error correction is given by the Devetak-Winter bound [34]

$$H(A|X = x^*, Q_E) - H(A|B, X = x^*, Y = y^*), \quad (27)$$

where again we are assuming a spot checking protocol and we consider the rate as the spot checking probability tends to 0 and the number of rounds tends to infinity. The second term in the rate can be directly estimated from the statistics in the protocol and the first term can be lower bounded using our technique.
<figure><img src="image_3.png" /><figcaption>Figure 3: Global randomness in the 2222-scenario. Comparison of lower bounds on $H(AB|X = 0, Y = 0, Q_E)$ for quantum devices that are constrained by a full distribution. The numerical bounds for our technique were computed using speedups (1) and (3) from Remark 2.6 at relaxation level $2 + ABZ$ including all monomials present in the objective function. A single SDP takes less than a minute to run at this level. We also compare with the iterated mean entropy $H_{(2)}^*(AB|X = 0, Y = 0, Q_E)$ from [26] and the numerical technique from [25] which we refer to as the TSGPL bound.</figcaption></figure>

We consider DI-QKD protocols in the 2322-scenario and we take $(x^*, y^*) = (0, 2)$, i.e., Bob's third input acts as his key generating input. In the same setup as the previous figure we consider devices with inefficient detectors and we constrain them by their full distribution.<fnref n="5" /> However, for QKD we allow Bob to record a device failure with an additional symbol $\perp$ when he receives his key generating input $y^* = 2$. By doing this, he collects more detailed information about his device's behaviour and in turn this allows him to reduce the size of $H(A|B, X = x^*, Y = y^*)$ slightly. We refer the reader to [43] for a more detailed discussion of this post-processing of no-click events.

To further boost the rates we also include preprocessing of the raw key [44]. Loosely, we allow for Alice to add additional noise to the outputs of her devices. In certain circumstances, this can increase the value of $H(A|X = x^*, Q_E)$ more than it increases the value of $H(A|B, X = x^*, Y = y^*)$ and therefore increasing the rate overall. More specifically, after Alice and Bob have collected their raw key (the outputs of their devices when $(X, Y) = (0, 2)$) Alice will independently flip each of her key-bits with some fixed probability $q \in [0, 1/2]$. For example, the post-measurement state for a single key-generating round is transformed after this preprocessing to

$$ 
\begin{align}
\rho_{ABQ_E} &= \sum_{ab} ((1-q)|a\rangle\langle a| + q|a \oplus 1\rangle\langle a \oplus 1|) \otimes |b\rangle\langle b| \otimes \rho_{Q_E}(a, b, x^*, y^*) \nonumber \\
&= \sum_{ab} |a\rangle\langle a| \otimes |b\rangle\langle b| \otimes ((1-q)\rho_{Q_E}(a, b, x^*, y^*) + q\rho_{Q_E}(a \oplus 1, b, x^*, y^*)) \tag{28} \\
&= \sum_{ab} |a\rangle\langle a| \otimes |b\rangle\langle b| \otimes \mathrm{Tr}_{Q_A Q_B} [\rho_{Q_A Q_B Q_E} (((1-q)M_{a|x^*} + qM_{a \oplus 1|x^*}) \otimes N_{b|y^*} \otimes I)]. \nonumber
\end{align}
$$

Thus this preprocessing can be seen as equivalent to Alice transforming her measurement from $\{M_0, M_1\}$ to $\{(1-q)M_0 + qM_1, (1-q)M_1 + qM_0\}$ on key generating rounds. It follows that we can then model this preprocessing in our numerical computations by modifying Alice's measurement operators in the objective function of (23) appropriately.

<footnote n="5">⁵Actually, we only constrain the devices by their distribution for inputs $(X, Y) \in \{0, 1\}^2$.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 4: **QKD 2322 rates.** Comparison of lower bounds on the key rate $H(A|X = 0, Q_E) - H(A|B, X = 0, Y = 2)$. The curves related to our technique were computed using speedups (1) and (3) from Remark 2.6 at a relaxation level 2 + ABZ and a single SDP at this relaxation level takes a few seconds to run. We also compare with the iterated mean entropy $H^}_{(8/7)}$ from [26] and with the analytic bound for a CHSH based DI-QKD protocol with a preprocessing step from [44, 19].</figcaption></figure>

In Figure 4 we plot the key rates achievable for devices with inefficient detectors. We compute our technique using a Gauss-Radau quadrature of $m = 12$ nodes. We compare our technique to bounds on the rate given by the iterated mean entropy $H^}_{(8/7)}$ from [26] (we do not include preprocessing for this curve as it did not give improvements) and an analytical bound from [44, 19] which is for a protocol based on the CHSH game and includes the preprocessing step. For both of the curves that incorporate the preprocessing, at each data point we optimized the probability $q$ with which Alice performs the her bitflip in order to maximize the obtained rates. We see that for the entirety of the plot the curve computed using our technique outperforms the other rate curves. In particular this shows the advantage that one can gain by changing the knowledge collected about the devices from a Bell-inequality violation to the full distribution. In the inset plot we zoom in on the region $[0.8, 0.88] \times [10^{-5}, 0.1]$. Looking at the inset plot we can inspect where the various rate curves vanish, in other words the minimal detection efficiency required to perform DI-QKD according to that curve. For the curve computed using the iterated mean entropy $H^}_{(8/7)}$ the minimal detection efficiency is around 0.845. The red curve based on the analytic bound actually vanishes around 0.826 [19]. However, we see that by using our method (purple curve) we are able to substantially reduce the threshold detection efficiency down to under 0.8. This threshold is now within a regime that is experimentally achievable. This warrants a more thorough analysis, using a more realistic model and computing finite round rates, in order to ascertain the performance capabilities of a current DIQKD experiment. We leave such an analysis to future work.

## 3 Methods

The objective of this section is to derive the variational upper bounds on the relatively entropy that were used in the previous section (see Theorem 3.11 below). However, first we need to define the quantum relative entropy and more generally quasi-relative entropies in the framework of von Neumann algebras and establish some of their properties (Section 3.1). We will then also describe rational approximations of the logarithm function (Section 3.2), which are a crucial ingredient in the derivation.
## 3.1 Quasi-relative entropies

Quasi-relative entropies can be defined in general for two positive semidefinite linear functionals on a von Neumann algebra $\mathcal{A}$. A linear functional $\rho: \mathcal{A} \to \mathbb{C}$ is said to be positive semidefinite if $\rho(a^*a) \ge 0$ for all $a \in \mathcal{A}$. We refer to a positive semidefinite linear functional satisfying $\rho(I) = 1$ as a state. We will mostly focus on the setting where $\mathcal{A} = B(H)$ for some separable Hilbert space $H$ and positive semidefinite linear functionals $\rho$ defined by some trace-class operator $\tilde{\rho}$ by $\rho(a) = \text{Tr}[\tilde{\rho}a]$ for all $a \in B(H)$. To simplify notation, we will often use the same symbol $\rho$ for both the positive semidefinite linear functional $\rho$ and the trace-class positive semidefinite operator $\tilde{\rho}$.

**Functional calculus of quadratic forms** We review briefly the functional calculus of sesquilinear (or quadratic) forms introduced by Pusz and Woronowicz in [45] and further applied in [46]. Let $\alpha$ and $\beta$ be two positive semidefinite quadratic forms on a complex vector space $U$. A function $\alpha: U \to \mathbb{R}$ is called a quadratic form if it satisfies $\alpha(\lambda u) = |\lambda|^2\alpha(u)$ for any $\lambda \in \mathbb{C}$ and $u \in U$ as well as the parallelogram identity $\alpha(u + v) + \alpha(u - v) = 2\alpha(u) + 2\alpha(v)$. We say that it is positive semidefinite if $\alpha(u) \ge 0$ for all $u \in U$. The theory of Pusz and Woronowicz allows one to define a new quadratic form $F(\alpha, \beta)$ for any function $F: \mathbb{R}_+^2 \to \mathbb{R}$ that is measurable (with respect to the $\sigma$-algebra of Borel subsets of $\mathbb{R}_+^2$), positive homogeneous (i.e., $F(\lambda x, \lambda y) = \lambda F(x, y)$ for all $\lambda \ge 0$) and locally bounded from below (i.e., bounded from below on all compact sets). For example if $F(x, y) = \sqrt{xy}$, this theory allows us to define the geometric mean of two positive semidefinite quadratic forms.

The main idea of the theory is that one can always represent a pair of positive semidefinite quadratic forms by two positive semidefinite, commuting operators. We call a representation of $(\alpha, \beta)$ a tuple $(H, A, B, h)$ where $H$ is a Hilbert space, $A, B$ are positive semidefinite and commuting bounded operators on $H$, and $h: U \to H$ is a linear map onto a dense subset of $H$ such that $\alpha(u) = \langle h(u), Ah(u) \rangle_H$, $\beta(u) = \langle h(u), Bh(u) \rangle_H$ for all $u \in U$. As shown in [45, Theorem 1.1], such a representation always exists. Then one defines the quadratic form $F(\alpha, \beta)$ on $U$ by

$$F(\alpha, \beta)(u) = \langle h(u), F(A, B)h(u) \rangle_H = \int_{\mathbb{R}_+^2} F(x, y) \langle h(u), dE(x, y)h(u) \rangle, \quad (29)$$

where $E$ is the joint spectral measure of the commuting operators $A$ and $B$; see [47, Chapter 5] for the definition of the joint spectral measure and [46] for the justification of the existence of this integral. The main result of the theory [45, Theorem 1.2] is that the quadratic form $F(\alpha, \beta)$ defined above is independent of the choice of representation $(H, A, B, h)$ of the pair $(\alpha, \beta)$.

**Defining quasi-relative entropies using the functional calculus of quadratic forms** Following the work of [46], we define the quasi-relative entropy of two positive semidefinite linear functionals $\rho, \sigma$ on a von Neumann algebra $\mathcal{A}$. Given $\rho$ and $\sigma$, we define two positive semidefinite quadratic forms on $\mathcal{A}$ (note that $\mathcal{A}$ is also a complex vector space) by

$$L_{\sigma}(a) = \sigma(aa^*) \quad \text{and} \quad R_{\rho}(a) = \rho(a^*a) \quad \text{for } a \in \mathcal{A}. \qquad (30)$$

Then using the functional calculus of Pusz-Woronowicz, one can define $F(L_{\sigma}, R_{\rho})$, as a quadratic form on $\mathcal{A}$. This leads to the following definition for the $F$-quasi-relative entropy.

**Definition 3.1** ($F$-quasi-relative entropy). Let $\rho$ and $\sigma$ be two positive semidefinite linear functionals on the von Neumann algebra $\mathcal{A}$. Then for any measurable, positive homogeneous and locally bounded from below function $F: \mathbb{R}_+^2 \to \mathbb{R} \cup \{+\infty\}$ the $F$-quasi-relative entropy between $\rho$ and $\sigma$ is

$$D_F(\rho \| \sigma) := F(L_{\sigma}, R_{\rho})(I) \qquad (31)$$

where

$$F(L_{\sigma}, R_{\rho})(I) := \langle h(I), F(A, B)h(I) \rangle_{K} \qquad (32)$$

with $(K, A, B, h)$ being a representation of the two positive semidefinite quadratic forms $(L_{\sigma}, R_{\rho})$ and this expression is defined by the integral in (29). We denote by $\nu_{\rho, \sigma}$ the positive measure
$\nu_{\rho,\sigma} = \langle h(I), Eh(I) \rangle$ on $\mathbb{R}_+^2$, where $E$ is the joint spectral measure of $A$ and $B$. Note that $\nu_{\rho,\sigma}$ only depends on $\rho$ and $\sigma$ and not on the function $F$. With this notation we have

$$D_F(\rho\|\sigma) = \int_{\mathbb{R}_+^2} F(x, y)d\nu_{\rho,\sigma}(x, y), \qquad (33)$$

with $\int_{\mathbb{R}_+^2} x d\nu_{\rho,\sigma}(x, y) = \sigma(I)$ and $\int_{\mathbb{R}_+^2} y d\nu_{\rho,\sigma}(x, y) = \rho(I)$.

We remark that for the applications we consider, it is important to allow the function $F$ to be infinite on some points. As described in [46], the functional calculus can readily be extended to this setting. In this setup, the quadratic form $F(L_\sigma, R_\rho)$ can take the value $+\infty$ on some points and thus, as expected, the corresponding $F$-quasi-relative entropy can be infinite.

**Definition 3.2** (Quantum relative entropy, also called Umegaki divergence [48]). Let $\rho$ and $\sigma$ be two positive semidefinite linear functionals on the von Neumann algebra $\mathcal{A}$. Then the relative entropy between $\rho$ and $\sigma$, written $D(\rho\|\sigma)$ is defined as the $F$-quasi-relative entropy with $F(x, y) = y \log_2(y/x)$:

$$D(\rho\|\sigma) := D_F(\rho\|\sigma). \qquad (34)$$

Note that when $\mathcal{A} = B(H)$ and $\sigma$ is the trace, we obtain the von Neumann entropy $H(\rho) = -D(\rho\|\sigma)$. The conditional von Neumann entropy can also be obtained in a similar way from the quantum relative entropy, this is described in Section 2.2.

**Definition 3.3** ($\alpha$-quasi-relative entropy). Let $\rho$ and $\sigma$ be two positive semidefinite linear functionals on $\mathcal{A}$ and let $\alpha \in (0, 1) \cup (1, \infty)$. Then the $\alpha$-quasi-relative entropy written $Q_\alpha(\rho\|\sigma)$ is defined as $F$-quasi-relative entropy with $F(x, y) = y^\alpha x^{1-\alpha}$.

**Remark 3.4** (Finite-dimensional case and the relative modular operator). It is instructive to consider the setting where $H$ is a $d$-dimension Hilbert space. In this case $B(H)$ is the algebra of $d \times d$ matrices with elements in $\mathbb{C}$. Then the two positive semidefinite linear functionals $\rho$ and $\sigma$ can be represented by positive semidefinite operators $\tilde{\rho}$ and $\tilde{\sigma}$, i.e., $\rho(a) = \text{Tr}[\tilde{\rho}a]$ and $\sigma(a) = \text{Tr}[\tilde{\sigma}a]$ for all $a \in B(H)$.

Let $\mathcal{L}_\sigma$ and $\mathcal{R}_\rho$ be left and right-multiplication operators by $\tilde{\sigma}$ and $\tilde{\rho}$, that is $\mathcal{L}_\sigma, \mathcal{R}_\rho : B(H) \to B(H)$ where $\mathcal{L}_\sigma(a) = \tilde{\sigma}a$ and $\mathcal{R}_\rho(a) = a\tilde{\rho}$. Note that $\mathcal{L}_\sigma$ and $\mathcal{R}_\rho$ are two commuting operators on $B(H)$ that are self-adjoint with respect to the Hilbert-Schmidt inner product on $B(H)$, $\langle a, b \rangle_{HS} = \text{Tr}[a^*b]$. These operators realize the quadratic forms $L_\sigma$ and $R_\rho$ defined in (30), in the sense that $L_\sigma(a) = \langle a, \mathcal{L}_\sigma a \rangle_{HS}$ and $R_\rho(a) = \langle a, \mathcal{R}_\rho a \rangle_{HS}$. Moreover, it is clear from their actions that $[\mathcal{L}_\sigma, \mathcal{R}_\rho] = 0$. In particular, our representation is $(B(H), \mathcal{L}_\sigma, \mathcal{R}_\rho, \text{id})$ where $\text{id} : B(H) \to B(H)$ is the identity map. Thus, for an $F$ satisfying the conditions of Definition 3.1, the $F$-quasi-relative entropy between $\rho$ and $\sigma$ is given by

$$F(L_\sigma, R_\rho)(I) = \langle I, F(\mathcal{L}_\sigma, \mathcal{R}_\rho)I \rangle_{HS}. \qquad (35)$$

This is precisely the definition used in [49]. Indeed, if we introduce the relative modular operator $\Delta = \mathcal{L}_\sigma \mathcal{R}_\rho^{-1}$ and we define $f(x) := F(x, 1)$ then for positive homogeneous $F$ we have $F(x, y) = yf(x/y)$ and the above is equal to

$$\langle I, \mathcal{R}_\rho f(\mathcal{L}_\sigma \mathcal{R}_\rho^{-1}) I \rangle_{HS} = \langle \mathcal{R}_\rho^{1/2} I, f(\Delta) \mathcal{R}_\rho^{1/2} I \rangle_{HS} = \langle \tilde{\rho}^{1/2}, f(\Delta) \tilde{\rho}^{1/2} \rangle_{HS} \qquad (36)$$

where we used the fact that $\mathcal{R}_\rho^{1/2}$ is self-adjoint with respect to the Hilbert-Schmidt inner product, and that $\mathcal{R}_\rho^{1/2}a = a\tilde{\rho}^{1/2}$. Note that when $\rho$ is not invertible the operator $\Delta$ multiplies on the right by the generalized inverse of $\tilde{\rho}$. It can be verified that in the finite-dimensional case, the quasi-relative entropy is given by (see e.g., [50, Equation 15])

$$D_F(\rho\|\sigma) = \sum_{j,k} F(q_k, p_j) |\langle \phi_k | \psi_j \rangle|^2 \qquad (37)$$

where $\tilde{\rho} = \sum_j p_j |\psi_j\rangle\langle\psi_j|$ and $\tilde{\sigma} = \sum_k q_k |\phi_k\rangle\langle\phi_k|$ are spectral decompositions of the density matrices $\tilde{\rho}$ and $\tilde{\sigma}$ respectively. Finally, choosing $F(x, y) = y \log_2(y/x)$ we recover the usual expression for
the quantum relative entropy

$$D(\rho||\sigma) = \mathrm{Tr}[\rho(\log_2 \rho - \log_2 \sigma)] \qquad (38)$$

and for the $\alpha$-quasi-relative entropies, $F(x, y) = x^{1-\alpha}y^{\alpha}$, we find $Q_{\alpha}(\rho||\sigma) = \mathrm{Tr}[\rho^{\alpha}\sigma^{1-\alpha}]$ which is the quantity within the logarithm of the Petz-Rényi divergences [51].

In order to obtain our variational expression for the quantum relative entropy, we will use a variational representation of $D_{F_t}$ when $F_t(x, y) = y \frac{x-y}{t(x-y)+y}$ for $t \in (0, 1]$. The relevance of $F_t$ to the relative entropy is that if we let $F(x, y) = y \log_2(y/x)$, then $F(x, y) = -\frac{1}{\ln 2} \int_0^1 F_t(x, y)dt$.

**Proposition 3.5.** Let $F_t(x, y) = y \frac{x-y}{t(x-y)+y}$ for $t \in (0, 1]$, let $\rho$ and $\sigma$ be positive semidefinite linear functionals on a von Neumann algebra $\mathcal{A}$. Then

$$D_{F_t}(\rho||\sigma) = \frac{1}{t} \inf_{a \in \mathcal{A}} \{\rho(I) + \rho(a + a^*) + (1-t)\rho(a^*a) + t\sigma(aa^*)\} . \qquad (39)$$

Furthermore, if $t < 1$, and $\rho$ and $\sigma$ are two trace-class positive semidefinite operators on the separable Hilbert space $H$ satisfying $\rho \le \lambda\sigma$ for some $\lambda \in \mathbb{R}_+$, then the infimum in (39) is achieved at an element with norm bounded by $\alpha := \max(\frac{3}{2(1-t)}, \frac{3\lambda}{2t})$, i.e., we can write

$$D_{F_t}(\rho||\sigma) = \frac{1}{t} \inf_{\substack{Z \in B(H) \\ \|Z\| \le \alpha}} \{\mathrm{Tr}[\rho] + \mathrm{Tr}[\rho(Z + Z^*)] + (1-t)\mathrm{Tr}[\rho Z^*Z] + t\mathrm{Tr}[\sigma ZZ^*]\} . \qquad (40)$$

**Remark 3.6.** Note that a constant bound on the norm of $Z$ in (40) can only be guaranteed when $t < 1$. Indeed, observe that for that $t = 1$ we have $F_1(x, y) = y - y^2/x$ so that $D_{F_1}(\rho||\sigma) = \rho(I) - Q_2(\rho||\sigma)$, and the variational expression becomes:

$$D_{F_1}(\rho||\sigma) = \inf_{Z \in B(H)} \{\mathrm{Tr}[\rho] + \mathrm{Tr}[\rho(Z + Z^*)] + \mathrm{Tr}[\sigma ZZ^*]\} . \qquad (41)$$

Assuming $\rho$ and $\sigma$ are finite-dimensional operators the infimum can be shown to be attained at $Z = -\sigma^{-1}\rho$ with an optimal value equal to $\mathrm{Tr}[\rho] - \mathrm{Tr}[\rho^2\sigma^{-1}]$. One can see that the operator norm of $\sigma^{-1}\rho$ can be made arbitrarily large even if we assume that $\rho \le \lambda\sigma$ (e.g., it suffices to take $\rho$ to be a pure state).

*Proof.* Proceeding in a similar way as [45, Lemma], if $A, B$ are two commuting Hermitian operators on a Hilbert space $H$, for any $y, z$ with $y + z = x$, a simple calculation gives

$$\langle x, F_t(A, B)x \rangle = \frac{1}{t} \langle x, \frac{t(A-B)B}{t(A-B)+B}x \rangle = \frac{1}{t} (\langle y, t(A-B)y \rangle + \langle z, Bz \rangle - \langle u, (t(A-B)+B)u \rangle) , \qquad (42)$$

where $u = (t(A - B) + B)^{-1}Bx - y$. As $t(A - B) + B$ is a positive semidefinite operator and we may choose $y$ so that $u = 0$, it follows that

$$\langle x, F_t(A, B)x \rangle = \frac{1}{t} \min_{y+z=x} \langle y, t(A-B)y \rangle + \langle z, Bz \rangle . \qquad (43)$$

Now we consider the positive semidefinite quadratic forms $L_{\sigma}$, $R_{\rho}$ for the positive semidefinite functionals $\sigma$ and $\rho$ as defined in (30) and we take a corresponding representation $(H, A, B, h)$. We get

$$F_t(L_{\sigma}, R_{\rho})(I) = \frac{1}{t} \min_{y+z=h(I)} \langle y, t(A-B)y \rangle + \langle z, Bz \rangle . \qquad (44)$$

Since $h: \mathcal{A} \to H$ is onto a dense subset of $H$, we can replace the minimum with an infimum to $y, z \in H$ such that $y = h(a)$ and $z = h(b)$. Then using the fact that $\langle h(a), Ah(a) \rangle = \sigma(aa^*)$ and $\langle h(a), Bh(a) \rangle = \rho(a^*a)$ we get

$$F_t(L_{\sigma}, R_{\rho})(I) = \frac{1}{t} \inf_{h(a)+h(b)=h(I)} \{t\sigma(aa^*) - t\rho(a^*a) + \rho(b^*b)\} . \qquad (45)$$
Note that if we have $a_1$ and $a_2$ satisfying $h(a_1) = h(a_2)$, then $\sigma(a_1a_1^*) = \langle h(a_1), Ah(a_1) \rangle = \langle h(a_2), Ah(a_2) \rangle = \sigma(a_2a_2^*)$ and similarly for $\rho$. As such we can replace the condition in the infimum by simply $a+b=I$, which gives

$$ 
\begin{aligned}
D_{F_t}(\rho\|\sigma) &= \frac{1}{t} \inf_{a+b=I} \{t\sigma(aa^*) - t\rho(a^*a) + \rho(b^*b)\} \\
&= \frac{1}{t} \inf_{a \in \mathcal{A}} \{\rho(I) - \rho(a + a^*) + (1-t)\rho(a^*a) + t\sigma(aa^*)\} \quad (46) \\
&= \frac{1}{t} \inf_{a \in \mathcal{A}} \{\rho(I) + \rho(a + a^*) + (1-t)\rho(a^*a) + t\sigma(aa^*)\}.
\end{aligned}
$$

where on the final line we made the substitution $a \mapsto -a$. This proves (39).

To prove the second part we first derive sufficient conditions for $a_0$ to achieve the infimum. Call $\phi(a)$ the expression to be infimized in (39). Note that $\phi$ is convex since $\rho$ and $\sigma$ are positive semidefinite (note for example that the restriction of $\phi$ to any line $a + sb$ is a convex quadratic in $s \in \mathbb{R}$). If $a_0$ achieves the infimum of $\phi$ it must be that $\frac{d}{ds}\phi(a_0 + sb)|_{s=0} = 0$ for any $b \in \mathcal{A}$. This equation gives

$$ \rho(b + b^*) + (1 - t)(\rho(a_0^*b + b^*a_0)) + t\sigma(a_0b^* + ba_0^*) = 0 \quad \forall b \in \mathcal{A}. \quad (47) $$

Since $\phi$ is convex this condition is also sufficient for optimality.

Now assume that $\rho$ and $\sigma$ are trace-class positive semidefinite operators on the separable Hilbert space $H$ so that $\rho(a) = \text{Tr}[\rho a]$ and $\sigma(a) = \text{Tr}[\sigma a]$. Equation (47) says that a necessary and sufficient condition for $A_0 \in B(H)$ to achieve the infimum is that

$$ \text{Tr}[\rho(B + B^*)] + (1 - t)\text{Tr}[\rho A_0^* B + \rho B^* A_0] + t\text{Tr}[\sigma A_0 B^* + \sigma B A_0^*] = 0 \quad \forall B \in B(H). \quad (48) $$

Letting $M = \rho + (1-t)A_0\rho + t\sigma A_0$, this is the same as $\text{Tr}[MB^*] + \text{Tr}[M^*B] = 0$ for all $B \in B(H)$ which is equivalent to $M = 0$. For convenience, we let $Z = -A_0$ and arrive at the following operator Sylvester equation

$$ (1 - t)Z\rho + t\sigma Z = \rho. \quad (49) $$

The existence of a bounded solution to this equation is guaranteed by the following lemma.

**Lemma 3.7.** Let $t \in (0, 1)$ and let $\rho$ and $\sigma$ be two trace-class positive semidefinite operators on the separable Hilbert space $H$ satisfying $\rho \le \lambda\sigma$ for some $\lambda \in \mathbb{R}_+$. Then the operator equation

$$ (1 - t)Z\rho + t\sigma Z = \rho \quad (50) $$

has a solution $Z \in B(H)$ which satisfies $\|Z\| \le \alpha := \max\{\frac{3}{2(1-t)}, \frac{3\lambda}{2t}\}$.

*Proof.* See Appendix B. $\square$

This completes the proof of Proposition 3.5. $\square$

## 3.2 Rational lower bounds on the logarithm

The quantum relative entropy is defined as the F-quasi relative entropy for $F(x, y) = y \log_2 y - y \log_2 x$. In this section we define a sequence of rational upper bounds on $F$ that are expressed as a finite sum of the functions $F_t$ for $t \in (0, 1]$.

The natural logarithm function (written $\ln$) has the following integral representation

$$ \ln(x) = \int_{0}^{1} f(t, x)dt \quad (51) $$

where

$$ f(t, x) = \frac{x-1}{t(x-1)+1}. \quad (52) $$
To get a rational approximation of ln, we can discretize the integral (51) with nodes $t_1, \dots, t_m \in [0, 1]$ and weights $w_1, \dots, w_m > 0$ to get a function

$$r(x) = \sum_{i=1}^{m} w_i f(t_i, x). \qquad (53)$$

For example, in [52] the authors used Gaussian quadrature to choose the weights and nodes of $r(x)$. The resulting function agreed with the first $2m+1$ derivatives (i.e., derivatives $0, \dots, 2m$) of the logarithm function. However, this approximation to the logarithm is unsuitable for the current work as it is not a global lower bound on ln.<fnref n="6" />

It turns out though that if we use Gauss-Radau quadrature then we get rational functions that are global bounds on ln. Gauss-Radau quadrature [53, p.103] is a variant of Gaussian quadrature where one of the endpoints of the integral interval is required to be a node of the quadrature formula.

**Theorem 3.8 (Gauss-Radau quadrature).** For any integer $m \ge 1$, there exist nodes $t_1, \dots, t_{m-1} \in (0, 1)$ and weights $w_1, \dots, w_m > 0$ such that the quadrature formula

$$\int_{0}^{1} g(t)dt \approx \sum_{i=1}^{m-1} w_i g(t_i) + w_m g(1) \qquad (54)$$

holds with equality for all polynomials $g$ of degree up to $2m-2$. Furthermore $w_m = 1/m^2$.

We note that such nodes and weights can be expressed in terms of properties of Legendre polynomials [53, p.103] and can be numerically computed efficiently [54]. If we use this quadrature formula in Theorem 3.8 for the integral representation of the logarithm (51) we get the following rational approximation of ln:

$$r_m(x) := \sum_{i=1}^{m-1} w_i f(t_i, x) + w_m f(1, x) = \sum_{i=1}^{m-1} w_i \frac{x-1}{t_i(x-1)+1} + w_m \frac{x-1}{x}. \qquad (55)$$

The next proposition shows that the rational functions $r_m$ are lower bounds on $\ln(x)$ that monotonically converge to $\ln(x)$.

**Proposition 3.9.** For any positive integer $m$, the function $r_m$ satisfies

1. $\ln(x) - r_m(x) = O((x-1)^{2m})$ for $x \to 1$,

2. $r_m(x) \to \ln(x)$ as $m \to \infty$ for any $x > 0$,

3. For any $x > 0$, $r_1(x) = 1 - 1/x \le r_m(x) \le r_{m+1}(x) \le \ln(x)$

*Proof.* See Appendix A. □

**Remark 3.10.** • One can obtain quantitative bounds on the approximation $\ln(x) - r_m(x)$ showing that the convergence is uniform on any compact segment $[\epsilon, \epsilon^{-1}] \subset (0, \infty)$ with a convergence rate $\approx (1 - \sqrt{\epsilon})^{2m}$. This can be proved using similar tools as [52, Proposition 6]. Furthermore, one can also show relative approximation error bounds that hold for all $x > 0$, such as

$$0 \le \ln(x) - r_m(x) \le \frac{1}{m^2} \frac{(x-1)^2}{x} \quad \forall x > 0. \qquad (56)$$

See [55, Proposition 2.2] and [32] for more details.

• Since $\ln$ satisfies $\ln(x^{-1}) = -\ln(x)$ it is natural to ask what is $r_m(x^{-1})$? Note that the rational functions $f(t, x)$ satisfy $f(t, x^{-1}) = -f(1 - t, x)$. As such one can verify that $r_m$ satisfies

$$r_m(x^{-1}) = -\bar{r}_m(x) \qquad (57)$$

<footnote n="6">⁶In fact, it is neither a lower bound nor an upper bound.</footnote>
where $\bar{r}_m(x)$ is the rational function obtained by applying an $m$-point Gauss-Radau quadra- ture on (51) with an endpoint at 0 (instead of 1). The resulting function $\bar{r}_m(x)$ is an upper bound on $\ln(x)$ and also satisfies $\bar{r}_m(x) - \ln(x) = O((x-1)^{2m})$ for $x \to 1$. Noting that $f(0, x) = x - 1$ we see that $\bar{r}_m(x)$ is a $(m, m-1)$ rational function$^7$, and as such it is the corresponding Padé approximant of $\ln$ at $x = 1$. Many properties of $\bar{r}_m(x)$ are discussed in detail in [56].

* One can obtain similar rational approximations for Petz divergences. The function $x \mapsto 1 - x^{-\alpha}$ is operator monotone for $\alpha \in (0, 1)$, and has the following integral representation

$$1 - x^{-\alpha} = \int_{0}^{1} \frac{x-1}{1+t(x-1)} d\mu_{\alpha}(t) \qquad (58)$$

where $d\mu_\alpha(t) = \frac{\sin(\alpha\pi)}{\pi} t^\alpha (1-t)^{-\alpha}$. To obtain rational approximations for $1 - x^{-\alpha}$, we can also apply the Gauss-Radau quadrature for the integral (58), this time with respect to the measure $d\mu_\alpha(t)$ (see e.g., [57]). This yields nodes $\{t_{\alpha,i}\}$ and weights $\{w_{\alpha,i}\}$ for $i = 1, \dots, m$, with $t_{\alpha,m} = 1$ that can be used to define the rational function

$$r_{\alpha,m}(x) = \sum_{i=1}^{m} w_{\alpha,i} f(t_{\alpha,i}, x). \qquad (59)$$

In [32] it is shown that the functions $r_{\alpha,m}(x)$ satisfy $r_{\alpha,m}(x) \le r_{\alpha,m+1}(x) \le 1 - x^{-\alpha}$ for all $x > 0$. Furthermore, if one uses the Gauss-Radau quadrature with an endpoint fixed at $t = 0$ (instead of $t = 1$), then the resulting rational functions are upper bounds on $1 - x^{-\alpha}$. We refer to [32] for further properties about these rational functions, including bounds on the approximation error as a function of $m$.

### 3.3 Variational expressions approximating the quantum relative entropy

We now state and prove the main result of this work, a sequence of converging variational upper bounds on the relative entropy.

**Theorem 3.11.** Let $\rho, \sigma$ be two positive semidefinite linear functionals on a von Neumann algebra $\mathcal{A}$ such that $Q_2(\rho\|\sigma) < \infty$ (see Definition 3.3). Then for any $m \in \mathbb{N}$ and the choice of $t_1, \dots, t_m \in (0, 1]$ and $w_1, \dots, w_m > 0$ as in Theorem 3.8, we have

$$D(\rho\|\sigma) \le -\sum_{i=1}^{m} \frac{w_i}{t_i \ln 2} \inf_{a \in \mathcal{A}} \{\rho(I) + \rho(a + a^*) + (1-t_i)\rho(a^*a) + t_i\sigma(aa^*)\} . \qquad (60)$$

Moreover, the RHS converges to $D(\rho\|\sigma)$ as $m \to \infty$.

In the special case where $\mathcal{A} = B(H)$ for a separable Hilbert space $H$ and $\rho$ and $\sigma$ are defined via trace-class operators on $H$ (also denoted by $\rho$ and $\sigma$) satisfying $\rho \le \lambda\sigma$ for some $\lambda \in \mathbb{R}_+$, we can give an explicit bound on the norm of the operators appearing in the optimization:

$$D(\rho\|\sigma) \le -c_m - \sum_{i=1}^{m-1} \frac{w_i}{t_i \ln 2} \inf_{Z \in B(H), \|Z\| \le \alpha_i} \{\mathrm{Tr}[\rho(Z + Z^*)] + (1-t_i)\mathrm{Tr}[\rho Z^* Z] + t_i \mathrm{Tr}[\sigma Z Z^*]\}, \qquad (61)$$

where $c_m = \mathrm{Tr}[\rho] \left( \sum_{i=1}^{m} \frac{w_i}{t_i \ln 2} - \frac{\lambda}{m^2 \ln 2} \right)$ and $\alpha_i = \frac{3}{2} \max\{\frac{1}{t_i}, \frac{\lambda}{1-t_i}\}$. Moreover, the RHS converges to $D(\rho\|\sigma)$ as $m \to \infty$.

**Remark 3.12.** This result can be seen as computationally tractable approximations of the variational expression for the quantum relative entropy due to Kosaki [27], a variant of which is also

<footnote>$^7$I.e., it can be written as the ratio of two polynomials $F/G$ where $\deg F = m$ and $\deg G = m - 1$</footnote>
present in the work of Donald [28]; we refer to the book [29] for more information on these expressions. For example in Donald [28, Section 4] it is shown that

$$D(\rho\|\sigma) = - \inf_{\substack{0<\delta<\epsilon<1 \\ a(t),b(t) \in \mathcal{A} \\ a(t)+b(t)=I}} \int_{\delta}^{\epsilon} \rho \left( \frac{1}{t} a(t) a(t)^* - b(t) b(t)^* \right) + \sigma (b(t)^* b(t)) dt, \quad (62)$$

where the functions $a(t), b(t) : (\delta, \epsilon) \to \mathcal{A}$ are piecewise constant. Importantly however, each one of our approximations is guaranteed to be an upper bound on the relative entropy.

**Proof.** The relative entropy $D$ is defined via the quasi-relative entropy $D_F$ for the function $F(x, y) = y \log_2 y - y \log_2 x$. Using the property of $\nu_{\rho,\sigma}$ (33), we have

$$D(\rho\|\sigma) = -\frac{1}{\ln 2} \int_{\mathbb{R}_+^2} y \ln(x/y) d\nu_{\rho,\sigma}(x, y). \quad (63)$$

Using Proposition 3.9, we can lower bound $\ln(x/y) \ge r_m(x/y)$ (defined in (55) and nodes and weights in (3.8)) obtaining

$$\begin{aligned} D(\rho\|\sigma) &\le -\sum_{i=1}^{m} \frac{w_i}{\ln 2} \int_{\mathbb{R}_+^2} y \frac{x/y-1}{t_i(x/y-1)+1} d\nu_{\rho,\sigma}(x, y) \\ &= -\sum_{i=1}^{m} \frac{w_i}{\ln 2} D_{F_{t_i}}(\rho\|\sigma) \end{aligned} \quad (64)$$

where on the second line we have defined the quasi-relative entropy $D_{F_{t_i}}(\rho\|\sigma)$ with $F_{t_i}(x, y) := y \frac{x-y}{t_i(x-y)+y}$. Now applying Proposition 3.5 to each $D_{F_{t_i}}$ we get

$$D(\rho\|\sigma) \le -\sum_{i=1}^{m} \frac{w_i}{t_i \ln 2} \inf_{a \in \mathcal{A}} \{\rho(I) + \rho(a + a^*) + (1 - t_i)\rho(a^*a) + t_i\sigma(aa^*)\}. \quad (65)$$

This concludes the proof of the upper bounds.

We now turn to the proof of convergence. Let $R_m(x, y) := \sum_{i=1}^m w_i F_{t_i}(x, y)$. Then we have already shown that $D(\rho\|\sigma) \le \sum_i -\frac{w_i}{\ln 2} D_{F_{t_i}}(\rho\|\sigma) = \frac{1}{\ln 2} D_{-R_m}(\rho\|\sigma)$. Furthermore, we have

$$\frac{1}{\ln 2} D_{-R_m}(\rho\|\sigma) - D(\rho\|\sigma) = \frac{1}{\ln 2} \int_{\mathbb{R}_+^2} y (\ln(x/y) - r_m(x/y)) d\nu_{\rho,\sigma}(x, y). \quad (66)$$

For any fixed $(x, y) \in \mathbb{R}_+^2$ we know from Proposition 3.9 that $y(\ln(x/y) - r_m(x/y))$ converges monotonically to 0 as $m \to \infty$. Thus, by the monotone convergence theorem we get that $\frac{1}{\ln 2} D_{-R_m}(\rho\|\sigma) - D(\rho\|\sigma) \to 0$ as desired.

To obtain the second statement with explicit bounds on the operators, we use the second part of Proposition 3.5. It suffices to use the variational expression in (40) for the terms $i \in \{1, 2, \dots, m-1\}$ instead of the general form. We get:

$$D(\rho\|\sigma) \le -\sum_{i=1}^{m-1} \frac{w_i}{t_i \ln 2} \inf_{\substack{Z \in B(H) \\ \|Z\| \le \alpha_i}} \{\mathrm{Tr}[\rho] + \mathrm{Tr}[\rho(Z + Z^*)] + (1-t_i)\mathrm{Tr}[\rho Z^* Z] + t_i \mathrm{Tr}[\sigma Z Z^*]\} - \frac{w_m}{\ln 2} D_{F_1}(\rho\|\sigma). \quad (67)$$

It now only remains to find a lower bound on $D_{F_1}(\rho\|\sigma)$ using the property $\rho \le \lambda\sigma$. Note that

$$\begin{aligned} D_{F_1}(\rho\|\sigma) &= \inf_{Z \in B(H)} \{\mathrm{Tr}[\rho] + \mathrm{Tr}[\rho(Z + Z^*)] + \mathrm{Tr}[\sigma Z Z^*]\} \\ &\ge \inf_{Z \in B(H)} \left\{\mathrm{Tr}[\rho] + \mathrm{Tr}[\rho(Z + Z^*)] + \frac{1}{\lambda}\mathrm{Tr}[\rho Z Z^*]\right\} \\ &= \inf_{Z \in B(H)} \left\{\mathrm{Tr}[\rho] - \lambda\mathrm{Tr}[\rho] + \mathrm{Tr}\left[\rho(\sqrt{\lambda}I + \frac{1}{\sqrt{\lambda}}Z)(\sqrt{\lambda}I + \frac{1}{\sqrt{\lambda}}Z)^*\right]\right\} \\ &= \mathrm{Tr}[\rho] - \lambda\mathrm{Tr}[\rho], \end{aligned} \quad (68)$$

which together with the fact that $w_m = \frac{1}{m^2}$ leads to the desired bound. $\square$
# 4 Conclusion

In this work we derived a converging sequence of upper bounds on the relative entropy between two positive semidefinite linear functionals on a von Neumann algebra. We then demonstrated how to use this sequence of upper bounds to derive a sequence of lower bounds on the conditional von Neumann entropy. The resulting optimization problems could then be relaxed to a convergent hierarchy of semidefinite programs using the NPA hierarchy. Overall this gives a computationally tractable method to compute lower bounds on the rates of DI-RE and DI-QKD protocols.

We applied our method to compute lower bounds on the asymptotic rates of various device-independent protocols. We compared the rates derived with our technique to other numerical techniques [25, 26, 36] and to known tight analytical bounds [15, 18, 19]. We found everywhere substantial improvements over the previous numerical techniques and we also demonstrated that our technique could recover all known tight analytical bounds, making it the first general numerical technique to do so and showing that the method can converge quickly. Compared with the previous numerical techniques, not only does our technique derive higher rates but it can also be much faster due to the fact that the noncommutative polynomial optimization that we derive is of low degree and can be ran without additional operator inequalities without affecting the rates in a substantial way. Furthermore, our derivation of this optimization problem was done in the general setting of bounded operators on a separable Hilbert space, allowing it to compute the rates of device-independent protocols even when the systems are infinite dimensional. This is in contrast to the previous general numerical techniques of [25, 26] that assumed finite dimensional Hilbert spaces in their analyses.

Computing key rates for DI-QKD, we also found significant improvements on the minimal detection efficiency required to generate secret key using a pair of entangled qubits. In particular we found a minimal detection efficiency lower than 0.8 which is now well within the regime of current device-independent experiments [30, 31, 58, 59, 60]. With further effort we can hope to increase the rates of these experiments further and we leave a thorough investigation of this question to future work.

We also demonstrated that min-tradeoff functions, necessary for the entropy accumulation theorem [9, 10], could be derived directly from our numerical computations. Therefore our technique can also be used to compute the rates of actual finite round protocols and to subsequently prove their security. For example, this then opens the possibility of using the global randomness bounds (see Figure 2) for the CHSH game to improve the rates of the recent DI-RE experiments [30, 31].

There remain a few open questions from this work. It would be interesting to investigate how we could make our computations more efficient. In particular, when exploring more complex scenarios the complexity of the SDPs can quickly grow (see the difference in runtimes between Fig. 1 and Fig. 2). To combat this one could search for other variational forms that converge faster to the von Neumann entropy. Otherwise, what are the best monomial sets to include in our relaxations or can we exploit symmetries to reduce the size of the SDPs [61]? When computing the NPA relaxations we made various simplifications to the problem (see Remark 2.6 and the captions of the various figures), in what settings can these simplifications be made without losing tightness for the rate curves? In a different direction, one could also look at applying our bounds on the relative entropy to bound other entropic quantities.

# Acknowledgements

We thank Roger Colbeck and Federico Grasselli for supplying data and information that facilitated some of the comparisons in the plots. We also thank Mateus Araújo for feedback on an earlier version of the manuscript, we thank Jean-Daniel Bancal for useful discussions that led to some improved numerics and we thank Glaucia Murta for pointing out an error in a previous version. This work is funded by the European Research Council (ERC Grant AlgoQIP, Agreement No. 851716).
References

[1] Charles H. Bennett and Gilles Brassard. "Quantum cryptography: Public key distribution and coin tossing". In Proceedings of IEEE International Conference on Computers, Systems, and Signal Processing. Pages 175-179. IEEENew York (1984). arXiv:2003.06557.

[2] 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). arXiv:1008.4593.

[3] Lluís Masanes, Antonio Acín, and Nicolas Gisin. "General properties of nonsignaling theories". Physical Review A 73, 012112 (2006). arXiv:quant-ph/0508016.

[4] Nicolas Brunner, Stefano Pironio, Antonio Acín, Nicolas Gisin, André Allan Méthot, and Valerio Scarani. "Testing the dimension of Hilbert spaces". Physical review letters 100, 210503 (2008). arXiv:0802.0760.

[5] Harry Buhrman, Richard Cleve, Serge Massar, and Ronald De Wolf. "Nonlocality and communication complexity". Reviews of modern physics 82, 665 (2010). arXiv:0907.3584.

[6] Ivan Šupić and Joseph Bowles. "Self-testing of quantum systems: a review". Quantum 4, 337 (2020). arXiv:1904.10042.

[7] Dominic Mayers and Andrew Yao. "Quantum cryptography with imperfect apparatus". In Proceedings of the 39th Annual Symposium on Foundations of Computer Science (FOCS-98). Pages 503-509. Los Alamitos, CA, USA (1998). IEEE Computer Society. arXiv:quant-ph/9809039.

[8] Artur K. Ekert. "Quantum cryptography based on Bell's theorem". Physical Review Letters 67, 661-663 (1991).

[9] Frédéric Dupuis, Omar Fawzi, and Renato Renner. "Entropy accumulation". Communications in Mathematical Physics 379, 1-47 (2020). arXiv:1607.01796.

[10] Frédéric Dupuis and Omar Fawzi. "Entropy accumulation with improved second-order term". IEEE Transactions on information theory 65, 7596-7612 (2019). arXiv:1805.11652.

[11] Emanuel Knill, Yanbao Zhang, and Honghao Fu. "Quantum probability estimation for randomness with quantum side information" (2018). arXiv:1806.04553.

[12] 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, 459 (2018).

[13] Rotem Arnon-Friedman, Renato Renner, and Thomas Vidick. "Simple and tight device-independent security proofs". SIAM Journal on Computing 48, 181-225 (2019). arXiv:1607.01797.

[14] C. Jordan. "Essai sur la géométrie à n dimensions". Bulletin de la S. M. F. 3, 103-174 (1875).

[15] Stefano Pironio, Antonio Acín, Nicolas Brunner, Nicolas Gisin, Serge Massar, and Valerio Scarani. "Device-independent quantum key distribution secure against collective attacks". New Journal of Physics 11, 045021 (2009). arXiv:0903.4460.

[16] J. F. Clauser, M. A. Horne, A. Shimony, and R. A. Holt. "Proposed experiment to test local hidden-variable theories". Physical Review Letters 23, 880-884 (1969).

[17] Federico Grasselli, Gláucia Murta, Hermann Kampermann, and Dagmar Bruß. "Entropy bounds for multiparty device-independent cryptography". PRX Quantum 2, 010308 (2021). arXiv:2004.14263.

[18] Federico Grasselli, Gláucia Murta, Hermann Kampermann, and Dagmar Bruß. "Boosting device-independent cryptography with tripartite nonlocality". Quantum 7, 980 (2023). arXiv:2209.12828.
[19] Erik Woodhead, Antonio Acín, and Stefano Pironio. "Device-independent quantum key distribution with asymmetric CHSH inequalities". *Quantum* 5, 443 (2021). arXiv:2007.16146.

[20] Jean-Daniel Bancal, Lana Sheridan, and Valerio Scarani. "More randomness from the same data". *New Journal of Physics* 16, 033011 (2014). arXiv:1309.3894.

[21] Olmo Nieto-Silleras, Stefano Pironio, and Jonathan Silman. "Using complete measurement statistics for optimal device-independent randomness evaluation". *New Journal of Physics* 16, 013035 (2014). arXiv:1309.3930.

[22] Robert König, Renato Renner, and Christian Schaffner. "The operational meaning of min- and max-entropy". *IEEE Transactions on Information Theory* 55, 4337-4347 (2009). arXiv:0807.1338.

[23] Miguel Navascués, Stefano Pironio, and Antonio Acín. "Bounding the set of quantum correlations". *Physical Review Letters* 98, 010401 (2007). arXiv:quant-ph/0607119.

[24] Stefano Pironio, Miguel Navascués, and Antonio Acín. "Convergent relaxations of polynomial optimization problems with noncommuting variables". *SIAM Journal on Optimization* 20, 2157-2180 (2010). arXiv:0903.4368.

[25] Ernest Y-Z Tan, René Schwonnek, Koon Tong Goh, Ignatius William Primaatmaja, and Charles C-W Lim. "Computing secure key rates for quantum cryptography with untrusted devices". *npj Quantum Information* 7, 158 (2021). arXiv:1908.11372.

[26] Peter Brown, Hamza Fawzi, and Omar Fawzi. "Computing conditional entropies for quantum correlations". *Nature communications* 12, 1-12 (2021). arXiv:2007.12575.

[27] Hideki Kosaki. "Relative entropy of states: a variational expression". *Journal of Operator Theory* Pages 335-348 (1986). url: https://www.jstor.org/stable/24714805.

[28] Matthew J Donald. "On the relative entropy". *Communications in Mathematical Physics* 105, 13-34 (1986).

[29] Masanori Ohya and Dénes Petz. "Quantum entropy and its use". Springer Science & Business Media. (2004).

[30] Wen-Zhao Liu, Ming-Han Li, Sammy Ragy, Si-Ran Zhao, Bing Bai, Yang Liu, Peter J Brown, Jun Zhang, Roger Colbeck, Jingyun Fan, et al. "Device-independent randomness expansion against quantum side information". *Nature Physics* Pages 1-4 (2021). arXiv:1912.11159.

[31] Lynden K Shalm, Yanbao Zhang, Joshua C Bienfang, Collin Schlager, Martin J Stevens, Michael D Mazurek, Carlos Abellán, Waldimar Amaya, Morgan W Mitchell, Mohammad A Alhejji, et al. "Device-independent randomness expansion with entangled photons". *Nature Physics* Pages 1-5 (2021). arXiv:1912.11158.

[32] Oisín Faust and Hamza Fawzi. "Rational approximations of operator monotone and operator convex functions" (2023). arXiv:2305.12405.

[33] Carl A. Miller and Yaoyun Shi. "Universal security for randomness expansion from the spot-checking protocol". *Siam Journal of Computing* 46, 1304-1335 (2017). arXiv:1411.6608.

[34] Igor Devetak and Andreas Winter. "Distillation of secret key and entanglement from quantum states". *Proceedings of the Royal Society A: Mathematical, Physical and engineering sciences* 461, 207-235 (2005). arXiv:quant-ph/0306078.

[35] Zhengfeng Ji, Anand Natarajan, Thomas Vidick, John Wright, and Henry Yuen. "MIP*=RE". *Communications of the ACM* 64, 131-138 (2021). arXiv:2001.04383.

[36] Rutvij Bhavsar, Sammy Ragy, and Roger Colbeck. "Improved device-independent randomness expansion rates using two sided randomness". *New Journal of Physics* 25, 093035 (2023). arXiv:2103.07504.
[37] Peter J Brown, Sammy Ragy, and Roger Colbeck. "A framework for quantum-secure device-independent randomness expansion". *IEEE Transactions on Information Theory* **66**, 2964–2987 (2019). arXiv:1810.13346.

[38] Peter Wittek. "Algorithm 950: Npol2sdpa—sparse semidefinite programming relaxations for polynomial optimization problems of noncommuting variables". *ACM Transactions on Mathematical Software (TOMS)* **41**, 1–12 (2015). arXiv:1308.6029.

[39] MOSEK ApS. "The MOSEK optimization suite 10.1.11". (2023). url: https://docs.mosek.com/latest/intro/index.html.

[40] Peter Brown. "Maintained fork of npol2sdpa". https://github.com/peterjbrown519/npol2sdpa.

[41] Peter Brown. "Example scripts for computing rates of device-independent protocols". Available at https://github.com/peterjbrown519/DI-rates (2021).

[42] Gláucia Murta, Federico Grasselli, Hermann Kampermann, and Dagmar Bruß. "Quantum conference key agreement: A review". *Advanced Quantum Technologies* **3**, 2000025 (2020). arXiv:2003.10186.

[43] Xiongfeng Ma and Norbert Lütkenhaus. "Improved data post-processing in quantum key distribution and application to loss thresholds in device independent qkd". *Quantum Information & Computation* **12**, 203–214 (2012). arXiv:1109.1203.

[44] M Ho, P Sekatski, EY-Z Tan, R Renner, J-D Bancal, and N Sangouard. "Noisy preprocessing facilitates a photonic realization of device-independent quantum key distribution". *Physical Review Letters* **124**, 230502 (2020). arXiv:2005.13015.

[45] Wiesław Pusz and Stanisław L Woronowicz. "Functional calculus for sesquilinear forms and the purification map". *Reports on Mathematical Physics* **8**, 159–170 (1975).

[46] Wiesław Pusz and Stanisław L Woronowicz. "Form convex functions and the WYDL and other inequalities". *Letters in Mathematical Physics* **2**, 505–512 (1978).

[47] Konrad Schmüdgen. "Unbounded self-adjoint operators on Hilbert space". *Volume 265*. Springer Science & Business Media. (2012).

[48] Hisaharu Umegaki. "Conditional expectation in an operator algebra, iv (entropy and information)". In Kodai Mathematical Seminar Reports. *Volume 14*, pages 59–85. Department of Mathematics, Tokyo Institute of Technology (1962).

[49] Edward G Effros. "A matrix convexity approach to some celebrated quantum inequalities". *Proceedings of the National Academy of Sciences* **106**, 1006–1008 (2009). arXiv:0802.1234.

[50] Naresh Sharma. "Equality conditions for quantum quasi-entropies under monotonicity and joint-convexity". In 2014 Twentieth National Conference on Communications (NCC). *Pages 1–6*. (2014).

[51] Dénes Petz. "Quasi-entropies for finite quantum systems". *Reports on Mathematical Physics* **23**, 57–65 (1986).

[52] Hamza Fawzi, James Saunderson, and Pablo A. Parrilo. "Semidefinite approximations of the matrix logarithm". *Foundations of Computational Mathematics* **19**, 259—296 (2019). arXiv:1705.00812.

[53] Philip J Davis and Philip Rabinowitz. "Methods of numerical integration". *Academic Press*. (1984). 2nd edition.

[54] Gene H Golub. "Some modified matrix eigenvalue problems". *SIAM review* **15**, 318–334 (1973).

[55] Hamza Fawzi and Omar Fawzi. "Semidefinite programming lower bounds on the squashed entanglement" (2022). arXiv:2203.03394.
[56] Flemming Topsøe. "Some bounds for the logarithmic function". Inequality theory and applications 4 (2007).

[57] Walter Gautschi. "Gauss-Radau formulae for Jacobi and Laguerre weight functions". Mathematics and Computers in Simulation 54, 403-412 (2000).

[58] Wei Zhang, Tim van Leent, Kai Redeker, Robert Garthoff, René Schwonnek, Florian Fertig, Sebastian Eppelt, Wenjamin Rosenfeld, Valerio Scarani, Charles C-W Lim, et al. "A device-independent quantum key distribution system for distant users". Nature 607, 687-691 (2022). arXiv:2110.00575.

[59] DP Nadlinger, P Drmota, BC Nichol, G Araneda, D Main, R Srinivas, DM Lucas, CJ Ballance, K Ivanov, EY-Z Tan, et al. "Experimental quantum key distribution certified by Bell's theorem". Nature 607, 682-686 (2022). arXiv:2109.14600.

[60] Wen-Zhao Liu, Yu-Zhe Zhang, Yi-Zheng Zhen, Ming-Han Li, Yang Liu, Jingyun Fan, Feihu Xu, Qiang Zhang, and Jian-Wei Pan. "Toward a photonic demonstration of device-independent quantum key distribution". Physical Review Letters 129, 050502 (2022). arXiv:2110.01480.

[61] Denis Rosset. "Symdpoly: symmetry-adapted moment relaxations for noncommutative polynomial optimization" (2018). arXiv:1808.09598.

[62] Gert K Pedersen. "On the operator equation $HT+TH=2K$". Indiana University Mathematics Journal 25, 1029-1033 (1976). url: https://www.jstor.org/stable/24891356.

[63] Timo Holz, Hermann Kampermann, and Dagmar Bruß. "Genuine multipartite Bell inequality for device-independent conference key agreement". Physical Review Research 2, 023251 (2020). arXiv:1910.11360.

[64] Renato Renner. "Security of quantum key distribution". PhD thesis. Swiss Federal Institute of Technology, Zurich. (2005). arXiv:quant-ph/0512258.

[65] Stephen Boyd and Lieven Vandenberghe. "Convex optimization". Cambridge university press. (2004).

A Proof of Proposition 3.9 on the rational functions $r_m$

1. We start by proving that $\ln(x) - r_m(x) = O((x-1)^{2m})$. We can write the Taylor expansion of $f(t, x)$ (defined in (52)) about $x = 1$ as

$$f(t, x) = (x - 1) \sum_{k=0}^{\infty} (-1)^k (x - 1)^k t^k, \qquad (69)$$

which is valid for $|t(x-1)| < 1$ and similarly for $\ln(x) = (x-1) + \sum_{k=1}^{\infty} (-1)^k \frac{(x-1)^{k+1}}{k+1}$ which is valid for $|x-1| < 1$. Let $\nu_m = \sum_{i=1}^{m} w_i \delta_{t_i}$ be the discrete measure coming from the Gauss-Radau quadrature, then we have for all $|x-1| < 1$

$$\begin{aligned}
\ln(x) - r_m(x) &= (x-1) \sum_{k=0}^{\infty} (-1)^k (x-1)^k \left( \int_0^1 t^k dt - \int_0^1 t^k d\nu_m(t) \right) \\
&= (x-1) \sum_{k=2m-1}^{\infty} (-1)^k (x-1)^k \left( \int_0^1 t^k dt - \int_0^1 t^k d\nu_m(t) \right) \\
&= (x-1) \sum_{k=0}^{\infty} (-1)^{k+2m-1} (x-1)^{k+2m-1} \left( \int_0^1 t^{k+2m-1} dt - \int_0^1 t^{k+2m-1} d\nu_m(t) \right) \\
&= (x-1)^{2m} \sum_{k=0}^{\infty} (-1)^{k-1} (x-1)^k \left( \int_0^1 t^{k+2m-1} dt - \int_0^1 t^{k+2m-1} d\nu_m(t) \right)
\end{aligned} \qquad (70)$$
where on the second line we used the fact that the Gauss-Radau formula (54) is exact for all polynomials of degree up to $2m-2$ means that the terms $k=0, \dots, 2m-2$ in the sum are equal to 0. Thus this gives $\ln(x) - r_m(x) = O((x-1)^{2m})$.

2. For any fixed $x > 0$ the function $t \mapsto f(t, x)$ is continuous on $[0, 1]$ and so can be approximated arbitrarily closely over $[0, 1]$ by polynomials. The claim then follows from the fact that the $m$-point quadrature rule is exact for all polynomials of degree $2m-2$.

3. This item can be shown using the fact mentioned in Remark 3.10 that $r_m(x) = -\bar{r}_m(x^{-1})$ where $\bar{r}_m$ is the $(m, m-1)$ Padé approximant of $\ln$ at $x = 1$. Indeed it was shown in [56, Equation (27)] that the functions $\bar{r}_m$ satisfy $\bar{r}_m(x) - \bar{r}_{m+1}(x) = \frac{(x-1)^{2m}}{S(x)}$ where $S(x)$ is a polynomial positive on $(0, \infty)$. This shows that $\bar{r}_m(x) \ge \bar{r}_{m+1}(x)$ and as such $r_m(x) \le r_{m+1}(x)$. Since $\lim_{m \to \infty} r_m(x) = \ln(x)$, this also shows $r_m(x) \le \ln(x)$.

## B Proof of Lemma 3.7 concerning the Sylvester equation (50)

We use a result of [62] to show the existence of a bounded operator $Z$ satisfying (50). Given the positive semidefinite operators $\rho$ and $\sigma$ on $H$, we construct the operators on the direct sum $H \oplus H: X = \begin{pmatrix} (1-t)\rho & 0 \\ 0 & t\sigma \end{pmatrix}$ and $K = \frac{1}{2} \begin{pmatrix} \rho & \rho \\ \rho & \lambda\sigma \end{pmatrix}$. Observe first that both $X$ and $K$ are positive semidefinite operators. In addition, recalling that $\alpha = \max\{\frac{3}{2(1-t)}, \frac{3\lambda}{2t}\}$, we get

$$ \alpha X - K = \begin{pmatrix} (\alpha(1-t) - \frac{1}{2})\rho & -\rho \\ -\rho & (\alpha t - \frac{\lambda}{2})\sigma \end{pmatrix} \ge \begin{pmatrix} \rho & -\rho \\ -\rho & \rho \end{pmatrix} \ge 0. \quad (71) $$

We are now ready to apply [62, Theorem 3.1]. Note that using the fact that $K \le \alpha X$, the operator monotonicity of the square root function, and the inequality $(a^2 + ax)^{1/2} \le a + \frac{x}{2}$ for $a, x \ge 0$, we get for any $s > 0$

$$ \begin{aligned} (X^2 + 2sK)^{1/2} - X &\le (X^2 + 2\alpha sX)^{1/2} - X \\ &= \begin{pmatrix} ((1-t)^2\rho^2 + 2\alpha s(1-t)\rho)^{1/2} - (1-t)\rho & 0 \\ 0 & (t^2\sigma^2 + 2\alpha s t\sigma)^{1/2} - t\sigma \end{pmatrix} \\ &\le \begin{pmatrix} \alpha sI & 0 \\ 0 & \alpha sI \end{pmatrix} \\ &\le \alpha sI. \end{aligned} \quad (72) $$

Thus we can apply [62, Theorem 3.1] and obtain the existence of a positive semidefinite operator $T$ on $H \oplus H$ with operator norm bounded by $\alpha$ that satisfies:

$$ XT + TX = 2K. \quad (73) $$

Writing $T = \begin{pmatrix} T_{11} & Z \\ Z^* & T_{22} \end{pmatrix}$, the equality (73) implies

$$ (1-t)\rho Z + tZ\sigma = \rho \quad (74) $$

and we have $\|Z\| \le \|T\| \le \alpha$.

## C The NPA hierarchy

In this section we briefly detail the NPA hierarchy which allows us to relax non-commutative polynomial optimization problems to a hierarchy of semidefinite programs. For a full exposition, we refer the reader to the original work [24].

We begin by recalling the basic definitions of noncommutative polynomial optimization problems. We consider noncommutative polynomials in the variables $x = (x_1, \dots, x_n)$ and their adjoints
$x^* = (x_1^*,...,x_n^*)$. A word corresponds to a monomial constructed from the variables above, i.e. products of the variables. We let $\mathcal{A}_k$ be the set of noncommutative polynomials that are linear combinations of words with $|w| \le k$. Similarly let $\mathcal{W}_k$ be the set of words with $|w| \le k$. We are interested in a noncommutative polynomial optimization problem of the form

$$ c_{\text{opt}} = \inf_{\substack{\mathcal{H}, \psi \in \mathcal{H}, \|\psi\|^2 = 1 \\ X_1, \dots, X_n \in \mathcal{B}(\mathcal{H})}} \langle \psi, p(X)\psi \rangle \\ \text{subject to} \quad \langle \psi, r_i(X)\psi \rangle \ge b_i \\ q_j(X) \ge 0 \tag{75} $$

Let $Q = \{q_j : j = 1, \dots, m\}$ be the set of positive polynomials in (75). The quadratic module $M_Q$ is the set of all polynomials of the form $\sum_{ij} f_i^* f_i + \sum_{ij} g_{ij}^* q_j g_{ij}$. If $M_Q$ contains the polynomial $C - \sum_{i=1}^n x_i x_i^* + x_i^* x_i$ for some $C > 0$ then we say the problem is Archimedean. Note that if each of the variables $x_i$ in the problem has an explicit bound on their operator norm, i.e. $x_i^* x_i \le C_i$ and $x_i x_i^* \le C_i$ then the problem becomes Archimedean if these constraints are added to the problem.

A moment relaxation of level $k$ is defined by a positive semidefinite linear functional $L: \mathcal{A}_{2k} \to \mathbb{C}$, i.e. for $f \in \mathcal{A}_{2k}$ with $f = \sum_{w \in \mathcal{W}_{2k}} f_w w$ for some $f_w \in \mathbb{C}$ then $L(f) = \sum_{w \in \mathcal{W}_{2k,2k}} f_w L(w)$. These positive semidefinite linear functionals are in one-to-one correspondence with so-called moment matrices of level $k$ $M_k$ which is a positive semidefinite matrix whose rows and columns are indexed by words in $\mathcal{W}_k$ and whose $(v, w)$ entry is given by

$$ M_k(v, w) = L(v^*w). \tag{76} $$

For each $q \in Q$ we also define the localizing moment matrix $M_{k-d_q}^q$ where $d_q = \lceil \deg(q)/2 \rceil$ as the matrix whose rows and columns are indexed by words in $\mathcal{W}_{k-d_q}$ and whose $(u, v)$ entry corresponds to

$$ M_{k-d_q}^q (v, w) = L(v^*qw). \tag{77} $$

We may then define the $k^{\text{th}}$ level relaxation of (75) as

$$ \begin{array}{l} c_k = \min \sum_{w} p_w L(w) \\ \text{s.t.} \quad M_k(1, 1) = 1 \\ \qquad \sum_{w} r_{i,w} L(w) \ge b_i \\ \qquad M_k \ge 0 \\ \qquad M_{k-d_q}^q \ge 0 \quad \text{for all } q \in Q \end{array} \tag{78} $$

If the problem is Archimedean then the authors of [24] showed that $\lim_{k \to \infty} c_k = c_{\text{opt}}$.

## D Additional plots

In this section we provide some additional plots that demonstrate the convergence of our technique. In particular we recover instances of other known tight analytical bounds [18, 19] in addition to [15] which was shown in Figure 1.

In Figure 5 we demonstrate that our technique can be used to recover the tight analytical bound [18] on $\inf H(A|X = 0, Q_E)$ for devices constrained to violate the Holz inequality [63] for three parties. Where the Holz inequality is formulated as follows: let Alice, Bob and Charlie each have binary input devices that give outputs in $\{1, -1\}$. Given a projective measurement $\{A_{1|x}, A_{-1,x}\}$ for Alice on input $x$ let $A_x = A_{1|x} - A_{-1|x}$ be the corresponding observable. Define observables $B_y$ and $C_z$ for Bob and Charlie in the same way. Finally, given a pair of observables $\{X_0, X_1\}$ define $X_{\pm} := \frac{1}{2}(X_0 \pm X_1)$. Then the Holz Bell-expression for three parties is given as

$$ \mathcal{B}_H = \langle A_1 B_+ C_+ \rangle - \langle A_0 (B_- + C_-) \rangle - \langle B_- C_- \rangle, \tag{79} $$
<figure><img src="image_5.png" /><figcaption>Figure 5: Recovering the tight analytical bound of [18]. Comparison of lower bounds on $H(A|X = 0, Q_E)$ for quantum devices that constrained to achieve some minimal violation of the Holz inequality for three parties [63]. Numerical bounds were computed at a relaxation level $2 + ABC$. A single SDP takes a few seconds to run at this level.</figcaption></figure>

it has a classical bound of 1 and a quantum mechanical bound of 3/2. In the figure we see that, in a similar fashion to Figure 1, we recover the tight analytical bound at a low NPA relaxation level and rapidly in the size of Gauss-Radau quadrature.

In Figure 6 we look at recovering instances of the family of tight analytical bounds on the asymmetric CHSH inequality [19]. The asymmetric CHSH inequality generalizes the standard CHSH inequality in the following way. Let the observables for Alice and Bob be defined in the same way as above, then the asymmetric CHSH expression with weight $\alpha \in \mathbb{R}$ is given as

$$ \mathcal{B}_{\alpha} = \alpha(\langle A_0 B_0 \rangle + \langle A_0 B_1 \rangle) + \langle A_1 B_0 \rangle - \langle A_1 B_1 \rangle. \qquad (80) $$

We see that for $\alpha = 1$ we recover the standard CHSH expression. The maximum classical value is given by $\max\{2, 2|\alpha|\}$ and the maximum quantum value is given by $2\sqrt{1+\alpha^2}$. In the figure we plot two instances of this family, $\alpha = 1.1$ and $\alpha = 0.9$. Similar to the Holz inequality we find that for both values of $\alpha$ we chose, we recover the analytical bound. However, note that in the case where $\alpha = 0.9$ it was slightly more challenging. If we use speedup (3) from Remark 2.6 then for the chosen relaxation level and Gauss-Radau quadratures, we do not recover part of the analytical bound. In particular, with the splitting of the objective function it appears to be unable to recover the tight rate curve. In particular, the construction of the analytical rate curve for $|\alpha| < 1$ in [19] is split into two parts, a convex function $g(s)$ is first defined which gives the rate curve for the larger violations $s$. Then the actual rate curve $r(s)$ is defined as

$$ r(s) := \begin{cases} g(s) & \text{if } |\alpha| \ge 1 \text{ or } s \ge s_* \\ g'(s_*)(|s| - 2) & \text{otherwise} \end{cases} \qquad (81) $$

where $s_* \in [2\sqrt{1+\alpha^2-\alpha^4}, 2\sqrt{1+\alpha^2}]$ is the unique point such that the tangent of $g$ gives a value of 0 at $s = 2$. Basically, the construction is to take the function $g$ to be the rate curve until the tangent of $g$ intersects the point $(2, 0)$ (moving from high scores to low scores) then the remaining rate curve is given by this tangent function. Numerical testing indicates that it is the flat part of the curve that the split objective (speedup (3) in Remark 2.6) cannot capture. Nevertheless, at the cost of more computational time, we can still run the numerics without speedup (3) and we're able to recover all of the rate curve in this way.
<figure><img src="image_6.png" /><figcaption>Figure 6: Recovery of the tight analytical bounds from [19]. Comparison of lower bounds on $H(A|X = 0, Q_E)$ for quantum devices that constrained to achieve some minimal violation of the asymmetric CHSH inequality for different values of the weighting parameter $\alpha$. For subfigure (a) the numerical bounds were computed at a relaxation level $2 + ABZ$ including all monomials present in the objective function and implemented using speedups (1) and (3) from Remark 2.6. For subfigure (b) we used speedups (1) and (2) from Remark 2.6.</figcaption></figure>
E Entropy accumulation and min-tradeoff functions

In order to compute the rates of finite round DI-RE or DI-QKD protocols, one can bound the total conditional smooth min-entropy $H_{\min}^{\epsilon}(\mathbf{AB}|\mathbf{XYE})$ accumulated by the devices during the protocol [64]. There are two tools, the entropy accumulation theorem [9, 10] and quantum probability estimation [11], that allow one to break down this total smooth min-entropy into smaller, easier to compute quantities. For instance, the entropy accumulation theorem roughly states that in an n-round protocol

$$H_{\min}^{\epsilon}(\mathbf{AB}|\mathbf{XYE}) \geq nf(q) - O(\sqrt{n}) \quad (82)$$

where $f(q)$ is a so-called min-tradeoff function defined as any function $f$ that maps probability distributions over a finite set $C$ (see the discussion around (14)) to $\mathbb{R}$ such that

$$f(q) \leq \inf H(AB|XYE) \quad (83)$$

where the infimum is a device-independent optimization over all strategies compatible with the average statistics $(C, q)$ observed in the protocol.

For a fixed distribution $q$ the rate calculations as performed in the main text will give a valid value for a potential min-tradeoff function. However, as we will now show, a single rate computation for some distribution $q'$ can also be used to construct a min-tradeoff function $f_{q'}$ that can then be evaluated for any distribution $q$. This construction is a consequence of weak duality for semidefinite programming problems [65].

Consider the following primal and dual pair of semidefinite programs

$$p(b) = \inf_{X} \mathrm{Tr}[CX] \\ \text{s.t.} \quad \mathrm{Tr}[F_i X] = b_i \quad \text{for all } i \\ X \ge 0 \quad (84)$$

and

$$d(b) = \sup_{\lambda, Y} \sum_{i} \lambda_i b_i \\ \text{s.t.} \quad C - \sum_{i} \lambda_i F_i - Y \ge 0 \quad (85) \\ Y \ge 0$$

where $b \in \mathbb{R}^m$ is some constraint vector and $C, F_i$ are real symmetric matrices. For the purposes of this exposition the primal problem should be thought of as the NPA moment matrix relaxation, where $X$ is the moment matrix.<fnref n="8" /> The constraints $\mathrm{Tr}[F_i X] = b_i$ can then be any general NPA constraints, as well as any statistical constraints that may have been imposed on the devices, e.g., for some statistical test $(C, q)$ we would impose constraints of the form $\sum_{(a,b,x,y):C(a,b,x,y)=c} p(a, b, x, y) = q(c)$. We want to show that for any fixed $\hat{b}$, a feasible point $(\hat{\lambda}, \hat{Y})$ to the dual problem parameterized by $\hat{b}$ defines a function $q(b) = \sum_i \hat{\lambda}_i b_i$ that is everywhere a lower bound on the optimal primal value $p(b)$. To see this note that for any feasible point $X$ of the primal problem parameterized by $b$ we have

$$\begin{aligned} 0 &\le \mathrm{Tr} \left[ X \left( C - \sum_i \hat{\lambda}_i F_i - \hat{Y} \right) \right] \\ &= \mathrm{Tr} [X C] - \sum_i \hat{\lambda}_i \mathrm{Tr} [X F_i] - \mathrm{Tr} [X \hat{Y}] \\ &= \mathrm{Tr} [X C] - \sum_i \hat{\lambda}_i b_i - \mathrm{Tr} [X \hat{Y}] \\ &\le \mathrm{Tr} [X C] - g(b) \end{aligned} \quad (86)$$

<footnote n="8">⁸Note that it is sufficient to consider real moment matrices when running our computations. This is due to the fact that for any feasible complex moment matrix $\Gamma$ the matrix $(\Gamma + \bar{\Gamma})/2$ is also a feasible moment matrix with the same objective value where $\bar{\Gamma}$ denotes the matrix corresponding to elementwise complex conjugation of $\Gamma$.</footnote>
where for both inequalities we used the fact that $\mathrm{Tr}[AB] \ge 0$ when $A$ and $B$ are both positive semidefinite. Rearranging and taking the infimum over all feasible $X$ we find that $g(b) \le p(b)$. Thus, from a single solution to the dual problem we can derive a function $g$ which is everywhere a lower bound on the optimal primal value. Moreover, by modifying the choice of constraint vector $\hat{b}$ with respect to which we solve the dual we can derive different lower bounding functions.

As mentioned above in the problems we consider the constraint vector will in general consist of constraints that are fixed and constraints that correspond to the statistical test $(C, q)$ which we vary as $q$ varies. As such we can split our constraint vector $b$ into $(b_{\mathrm{fixed}}, b_{\mathrm{vary}})$ where $b_{\mathrm{vary}}$ is precisely the distribution $q$ viewed as a vector. Then by defining $\alpha = \lambda_{\mathrm{fixed}} \cdot b_{\mathrm{fixed}}$ where $\lambda_{\mathrm{fixed}}$ is the part of the dual solution $\lambda$ corresponding to the $b_{\mathrm{fixed}}$ part of $b$, we can derive a function $g(q) = \alpha + \lambda_{\mathrm{vary}} \cdot q$ where we have written $q$ in place of $b_{\mathrm{vary}}$. Thus, for any primal problem corresponding to a lower bound on rate of a protocol satisfying some statistical test $(C, q)$, e.g., $\inf H(A|X=0, E)$, we immediately get an affine function $g(q)$ that satisfies

$$g(q) \le \inf H(A|X=0, E) \qquad (87)$$

which is in essence a min-tradeoff function. One could then use $g(q)$ together with the entropy accumulation theorem to get bounds on the rates of finite round protocols and prove security in exactly the same manner as in [37].

## F Improved bounds for noisy-preprocessing

For this section we assume that all Hilbert spaces are finite dimensional. When deriving the statement (5) we removed the $m^{\mathrm{th}}$ term in the summation $D_{F_1}$ by lower bounding it by a trivial value $\mathrm{Tr}[\rho] - \lambda \mathrm{Tr}[\rho]$. In the application to computing device-independent rates $\lambda = 1$ and the lower bound is 0. This bounding was done in order to guarantee that all the operators in the subsequent NPA hierarchy relaxations were explicitly bounded which can help improve numerical stability. However, in the application to DI-QKD where we consider noisy-preprocessing we can derive a tighter bound that depends on the bitflip probability $q \in [0, 1/2]$.

First note that if $\rho$ is a state, i.e $\mathrm{Tr}[\rho] = 1$, then

$$D_{F_1}(\rho \| \sigma) = 1 - \mathrm{Tr}[\rho^2 \sigma^{-1}]. \qquad (88)$$

Following the proof of Theorem 3.11 we need to lower bound this quantity or equivalently upper bound $\mathrm{Tr}[\rho^2 \sigma^{-1}]$. In the setting of DI-QKD considered in the main text we have that after measuring but before the noisy-preprocessing Alice and Eve's joint state is

$$\rho_{AE} = \sum_{a} |a\rangle\langle a| \otimes \rho_{E}(a) \qquad (89)$$

where $\rho_E(a) = \mathrm{Tr}_{Q_A Q_B} [\rho_{Q_A Q_B Q_E}(M_{a|x^*} \otimes I \otimes I)]$. After Alice applies noisy-preprocessing the state transforms to

$$\begin{aligned} \tilde{\rho}_{AE} &= \sum_{a} ((1-q)|a\rangle\langle a| + q|\bar{a}\rangle\langle\bar{a}|) \otimes \rho_E(a) \\ &= \sum_{a} |a\rangle\langle a| \otimes ((1-q)\rho_E(a) + q\rho_E(\bar{a})) \end{aligned} \qquad (90)$$

where we have defined $\bar{a} = a \oplus 1$ and we can also see that $\mathrm{Tr}_A[\tilde{\rho}_{AE}] = \rho_E$, i.e. Eve's marginal state is unchanged by Alice's local data processing. We are interested in upper bounds on $\mathrm{Tr}[\tilde{\rho}_{AE}^2(I \otimes \rho_E^{-1})]$ that are smaller than the previously obtained bound of 1. Taking the partial trace over the $A$ system we have

$$\mathrm{Tr}[\tilde{\rho}_{AE}^2(I \otimes \rho_E^{-1})] = \mathrm{Tr}\left[\sum_{a} ((1-q)\rho_E(a) + q\rho_E(\bar{a}))^2 \rho_E^{-1}\right]. \qquad (91)$$
Now by a direct calculation we have

$$
\begin{align}
\sum_a ((1-q)\rho_E(a) + q\rho_E(\bar{a}))^2 &= \sum_a (1-q)^2 \rho_E(a)^2 + q^2 \rho_E(\bar{a})^2 + q(1-q)(\rho_E(a)\rho_E(\bar{a}) + \rho_E(\bar{a})\rho_E(a)) \nonumber \\
&= (1-2q(1-q))(\rho_E(0)^2 + \rho_E(1)^2) + 2q(1-q)(\rho_E(0)\rho_E(1) + \rho_E(1)\rho_E(0)) \nonumber \\
&= (1-4q(1-q))(\rho_E(0)^2 + \rho_E(1)^2) + 2q(1-q)(\rho_E(0) + \rho_E(1))^2 \nonumber \\
&= (1-4q(1-q))(\rho_E(0)^2 + \rho_E(1)^2) + 2q(1-q)\rho_E^2. \tag{92}
\end{align}
$$

Applying this to the above expression we find

$$
\begin{align}
\mathrm{Tr} [\tilde{\rho}_{AE}^2(I \otimes \rho_E^{-1})] &= (1 - 4q(1 - q))\mathrm{Tr} [\rho_{AE}^2(I \otimes \rho_E^{-1})] + 2q(1 - q) \tag{93} \\
&\leq 1 - 2q(1 - q) \nonumber
\end{align}
$$

where we used the fact that Tr [ρ̃²ₐ₆(I ⊗ ρ̃⁻¹)] ≤ 1. We can see that this bound is tight at q = 1/2 as this implies the Petz entropy of order 2, H₂(A|E) = -log₂ Tr [ρ̃²ₐ₆ρ̃⁻¹], must be at least 1. This is exactly what we expect as Alice's system at q = 1/2 is just a uniformly random bit. Similarly if q = 0 we recover our original bound without preprocessing.

Putting everything together, we have

$$
D_{F_1}(\rho \| \sigma) \geq 1 - 2q(1 - q) \qquad (94)
$$

this implies that the constant term $c_m$ in Lemma 2.3 can be replaced by

$$
c_m = \frac{2q(1-q)}{m^2 \ln 2} + \sum_{i=1}^{m-1} \frac{w_i}{t_i \ln 2} \quad (95)
$$

when computing key-rates for DI-QKD protocols that include the noisy-preprocessing step.