Reliable numerical key rates for quantum key distribution

In this work, we present a reliable, efficient, and tight numerical method for calculating key rates for finite-dimensional quantum key distribution (QKD) protocols. We illustrate our approach by finding higher key rates than those previously reported in the literature for several interesting scenarios (e.g., the Trojan-horse attack and the phase-coherent BB84 protocol). Our method will ultimately improve our ability to automate key rate calculations and, hence, to develop a user-friendly software package that could be used widely by QKD researchers.


Introduction
The possibility of large-scale quantum computers in the near future has spawned the field of quantum-safe cryptography [1]. This includes classical techniques based on computational hardness as well as information-theoretic approaches via physical assumptions. The latter invokes either assumptions about the physical channel [2], or less restrictive, the assumption only that the laws of quantum physics are correct, which is the basis for quantum key distribution (QKD). See Ref. [3] for a review of QKD and Ref. [4] for an update of recent progress. Both classical methods and QKD will likely play a role in quantum-safe cryptographic implementations.
The maturity of QKD technology is evidenced by a recent QKD satellite launch [5] as well as developments of fiber-based networks [6][7][8], suggesting that global networks are on the horizon. Still, there remains important open problems in QKD theory, such as (1) optimizing the practicality of QKD protocols to make them easily implementable, and (2) understanding the effects of device imperfections and side channels.
Solving these problems requires a robust theoretical method for evaluating a QKD protocol's performance, which involves a detailed security analysis. Performance is then quantified by the key rate -the number of bits of secret key obtained per exchange of quantum signal. Unfortunately, analytical methods for calculating the key rate are highly technical, are often limited in scope to particular protocols, and invoke inequalities that introduce looseness into the calculation.
We therefore focus our efforts on numerical methods, which are inherently more robust to both device imperfections and changes in protocol structure. Furthermore, numerics can be made user-friendly, such that the user needs only to define the specifications of the protocol and then the computer performs the key rate calculation. As an example, our group recently released a software package for this purpose. 1 The key rate calculation involves minimizing a convex function over all eavesdropping attacks that are consistent with the experimental data [9][10][11][12]. When employing numerics, one issue that arises is the efficiency of this optimization. This issue is particularly important for high-dimensional QKD protocols, or protocols with many signal states, since the relevant optimization involves many parameters. In such cases, the computational time can be very long -sometimes days -so it is crucial to implement a high efficiency algorithm.
Another issue with numerics is reliability. This is more subtle but also more important than the efficiency issue. Due to the inherent paranoia in cryptography, it is natural to ask whether numerically calculated key rates are trustworthy. After all, computers have finite numerical precision. Furthermore, optimization algorithms never truly reach the global optimum, as termination conditions always have some non-zero tolerance. Since QKD is now a serious, real-world technology, key rates must come with a security guarantee, and handwaving at these numerical issues will not suffice.
In this work, we present a numerical method that solves both the reliability and efficiency issues. Our method provides reliable lower bounds on the key rate with arbitrary tightness for finite-dimensional QKD protocols. Furthermore it is highly efficient and typically returns a key rate within seconds or less on one's personal computer. To illustrate our method, we apply it to three practically interesting scenarios. Namely we consider the Trojan-horse attack [13][14][15], the BB84 protocol with phase-coherent signal states [16,17], and the BB84 protocol with detector efficiency mismatch [18]. We improve upon literature key rates in all three cases.
Directly calculating the key rate involves a minimization problem (see Sec. 2). When solving this on a computer, the algorithm will terminate before reaching the global optimum and hence will return an upper bound on the true key rate. However, we are interested in reliable lower bounds, i.e., achievable key rates. In previous work [19], we noted this issue as motivation for transforming the optimization problem to the so-called dual problem [20]. This transforms the minimization problem into a maximization problem. Therefore, the dual problem will return a lower bound on the key rate, as desired. This method led to novel insights for particular protocols as discussed in [19]. However, in order to simplify the optimization in the dual problem, we invoked an inequality that in some cases introduces looseness into the key rate and, furthermore, makes the optimization problem non-convex. Ultimately the non-convexity reduces the efficiency of this approach, making it difficult to apply to protocols with large numbers of signal states.
We therefore present a method here that retains the efficiency of convex optimization, but also has the reliability of the dual problem. Our approach is to break up the calculation into two steps. The first step approximately minimizes the convex function, and hence finds an eavesdropping attack that is close to optimal. The second step takes this approximately optimal attack and converts it into a lower bound on the key rate. Breaking it up into these two steps adds flexibility to our method, in that any algorithm can be employed for the initial minimization of the convex function.
Our main technical result is to provide a recipe for performing the second step, i.e., for converting a nearoptimal attack into a tight lower bound on the key rate. At the technical level, we derive our main result first by linearizing the problem and then by transforming to the dual problem of the subsequent linearized problem. The idea is that, for a convex function, any linearization (about any point) will undercut (and hence lower bound) the curve. One obtains the tightest lower bound by this method if one linearizes about a point corresponding to the global minimum of the convex function.
We emphasize several differences between this work and Ref. [19]. Our method presented herein is tight, essentially giving the optimal key rate, whereas Ref. [19] invoked the Golden-Thompson inequality which for certain protocols introduces looseness into the calculated key rates. This difference can be seen below in our Figures 3,4, and 5, which show a clear gap between our key rates and those computed from the method in Ref. [19]. Another difference, as noted earlier, is the lack of convexity of the method in Ref. [19]. For a small number of optimization parameters, e.g., as encountered with entanglement-based QKD protocols, the non-convex method in Ref. [19] works quite well. On the other hand, for the method in Ref. [19], prepareand-measure QKD protocols require adding a number of optimization parameters that scale quadradically in the number of signal states, as we need to describe the pairwise quantum-mechanical overlaps between the signals (see Ref. [19] and the discussion below around Eq. (38)), which can lead to long computation times. In contrast, the convexity of the method presented herein leads to more efficient scaling of the computation time with the number of signal states for prepare-and-measure QKD, despite the larger amount of open parameters of this method.
In what follows we first give background on key rate calculations in the next section. Then we present our main result in Sec. 3. In Sec. 4 we describe how our approach applies to a general class of QKD protocols. We illustrate our method for three interesting example protocols in Sec. 5, and finally we conclude in Sec. 6. Technical details can be found in the Appendix.

Background
The well-known asymptotic key rate formula [21] is given by the difference of two information-theoretic quantities associated, respectively, with privacy amplification (PA) and error correction (EC). These two terms appear in the following expression for the key rate per signal We explain this expression in more detail in Sec. 4. For now, we note that p pass refers to the probability for passing the post-selection (e.g., sifting) in the protocol, leak EC obs denotes the number of bits of information (per signal that passed post-selection) that Alice publicly reveals during error correction, and f (ρ) = p pass · f (ρ) is a function defined below in Eq. (5).
The first term in (1) is the PA term. Here, ρ is the density operator shared by Alice, Bob, and possibly other parties involved in the protocol. (Note that prepare-and-measure protocols can be recast as entanglement-based protocols and are described by same mathematics, see Sec. 4 for elaboration.) We assume that the state has an i.i.d. (independent, identically distributed) structure, and hence it makes sense to discuss the state ρ associated with a single round of quantum communication. This i.i.d. structure corresponds to Eve doing a so-called collective attack. However, the security of our derived asymptotic key rate also holds against the most general attacks (coherent attacks) if one imposes that the protocol involves a random permutation of the rounds (a symmetrization step) such that the Quantum de Finetti theorem [9,22] applies.
The density operator ρ is unknown, but the asymptotic experimental data gives linear constraints on it, of the form where the Γ i are Hermitian operators. Let S denote the set of states that satisfy these constraints where H + is the set of positive semidefinite operators. Also, we add the identity to the set {Γ i } to enforce that Tr(ρ) = 1, giving a total of n constraints.
The key rate calculation is an optimization problem, since we must consider the worst-case scenario (the most powerful eavesdropping attack) that is consistent with the experimental data. Hence Eq. (1) involves minimizing over all ρ ∈ S. Note that the error correction term is exactly determined by the observations, and hence we can pull it out of the optimization in (2). Similarly, p pass is known from the observations and is pulled into the optimization in (2).
As discussed in Sec. 4, f (ρ) can be written as where D(σ||τ ) := Tr(σ log σ)−Tr(σ log τ ) is the relative entropy, G is a completely positive (CP) map, and Z is a completely positive trace preserving (CPTP) map, more specifically a pinching quantum channel. (Sec. 4 discusses the meaning of G and Z, which are respectively related to the post-selection and the key map of the QKD protocol). Due to the joint convexity of the relative entropy and the fact that G and Z are linear maps, the function f (ρ) is convex in ρ. Furthermore the problem is a convex optimization problem since the set S is convex (see, e.g., Ref. [20]). While efficient numerical methods are known for such convex problems, the key rate calculation is unique compared to other convex problems, in that getting "close" to the optimal point is not good enough. One needs a reliable lower bound on the key rate, i.e., guaranteed security.

Main result 3.1 Reliable lower bound
We now show how to lower bound the minimization problem in (6). Our strategy is to break up the key rate calculation into two steps: • Step 1: Find an eavesdropping attack that is close to optimal, which gives an upper bound on the key rate. • Step 2: Convert this upper bound to a lower bound on the key rate.
With our approach, Step 1 does not need to be perfect -any eavesdropping attack may be used as an input for Step 2. However, if Step 1 returns the optimal attack, our lower bound calculated by Step 2 will be tight. Furthermore, our method for Step 2 is continuous around the optimal attack. Thus, finding a near-optimal attack produces a near-optimal lower bound.
Step 1 may be solved in various ways using convex optimization methods [20]. For concreteness, Sec. 3.5 presents one such method, which exploits the structure of our problem and is relatively fast.
On the other hand, our main result is a method for performing Step 2. We approach Step 2 via a sequence of theorems that successively improve the reliability and robustness of the lower bounds which they return. First, Theorem 1 presents the conceptual foundation for our lower bounding method. However, this theorem is stated under a restrictive assumption that may not hold in special cases. Therefore, we extend our result in Theorem 2. Finally, we improve our result once more in Theorem 3, which addresses numerical imprecision and is directly useful for numerical key rate calculations. Sec. 3.4 presents the argument that our method yields arbitrarily tight bounds on the key rate.
We now present our main result in its simplest conceptual form. To state this result, we first define the gradient of f at point ρ, whose representation in the standard basis {|j } is and σ jk := j|σ|k .
Theorem 1: Given any ρ ∈ S, if ∇f (ρ) exists, then Step 1 Step 2 Our lower bound L in ea ri za ti o n D u a l Figure 1: Illustration of our lower-bounding method.
Step 1 is any algorithm that takes an initial feasible point ρ0 and outputs another feasible point ρ, which may or may not be close to the optimal attack ρ * . Note that f (ρ) provides an upper bound on f (ρ * ) and, hence, on the key rate.
Step 2 converts this into a lower bound, by solving the dual problem of the linearization of f about point ρ. Since the linearization undercuts the curve f and since the dual problem is a maximization, our lower bound is reliable even if the numerical calculation does not reach the global optimum.
where α was defined in (6) and Here, the transpose T is taken in the same basis as that used to define the gradient in (7), and γ = {γ i } is the vector of expectation values from (3), which is of size n (the total number of constraints). Furthermore, equality in (8) holds if f (ρ) = α, i.e., if ρ corresponds to an optimal attack.
The proof is given in Appendix A. It involves linearizing the convex function f about point ρ and then transforming the subsequent linearized problem to its dual problem (see [20] for discussion of duality). Figure 1 illustrates the basic idea of Theorem 1. Theorem 1 takes any feasible eavesdropping attack ρ, which gives an upper bound f (ρ) on α, and converts it into a reliable lower bound on α. The fact that (9) involves a maximization is crucial for the reliability of the calculation. Since maximization involves approaching the solution from below, every number that the computer outputs is a lower bound on α, even if the computer does not reach the global maximum.

A more robust bound
The assumption that ∇f exists over the entirety of S does not necessarily hold. In particular, the gradient has the form which is derived using rules for derivatives of matrix traces [23]. (For example, the derivative of Tr(G(X)) is g(X) T where g(·) is the scalar derivative of G(·).) By inspection, we see that if G(ρ) is singular (i.e., not full rank), ∇f (ρ) may not exist, and hence Theorem 1 would not apply. Ideally we would like a bound that holds for all ρ ∈ S. Towards this end, it is helpful to introduce a slight perturbation that maps G(ρ) to its interior, i.e., points on the boundary become interior points. The key fact that makes this perturbation useful is that we can upper bound the difference between the function f evaluated on the perturbed and unperturbed inputs, and this ultimately leads to our result in Theorem 2 below. We introduce this perturbation as follows.
Let be such that 0 < < 1, and consider a perturbed channel where D is the depolarizing channel, where τ = 1/d is the maximally mixed state. Here (and for future reference) we define d and d , respectively, to be the dimension of ρ and G(ρ). Furthermore, we define a perturbed function which is just (5) with G replaced by G . The key point is the following lemma.

Lemma 1:
The gradient ∇f (ρ), which is given by always exists for all ρ 0.
As elaborated in Appendix B, Lemma 1 follows from the fact that G (ρ) > 0 and similarly Z(G (ρ)) > 0. Hence the logarithm of these two matrices in (15) is well-defined.
The above lemma motivates the following theorem, which is a more robust version of Theorem 1.
, where e is the base of the natural logarithm, then where .
The proof is given in Appendix C.3. The basic idea of Theorem 2 is that it is essentially just Theorem 1, but applied to a slightly perturbed function f . Alternatively, one can think of Theorem 2 as being Theorem 1 applied to a perturbed state, G (ρ) instead of G(ρ), where the perturbation maps states on the boundary to the interior.
Note that Theorem 2 generalizes Theorem 1, since we have that assuming that ∇f (ρ) exists. However, Theorem 2 is more robust than Theorem 1, because it holds for any ρ ∈ S, even if ∇f (ρ) does not exist.

Finite precision computation
Finding a lower bound with our method requires a computer program for any nontrivial QKD protocol. The finite precision inherent to computational methods introduces errors that threaten the reliability of the calculated lower bounds. With this threat in mind, we identify all possible sources of errors in evaluating Theorem 1. First, the variables will not be precise: ρ will not exactly satisfy the constraints, and {Γ i } and {γ i } will not be exact. Second, function evaluations will not be precise: every function evaluation introduces errors.
The aforementioned errors fall into broader categories. Variable imprecision is implementation independent since we can characterize the imprecision in a universal manner. Function-evaluation imprecision is implementation dependent since in general, its characterization varies widely with the particular algorithm (particularly for evaluating nontrivial functions such as the matrix logarithm appearing in our objective function). Since the latter kind of errors depend on the implementation, a universal treatment of them is not possible. In principle, it is possible to bound the effect of such errors for a particular implementation [24], although it is beyond the scope of this article. On the other hand, we give a full treatment of implementationindependent errors in what follows.
From a computational perspective, it is virtually impossible to find an element strictly in S. Furthermore, for many problems, the computer representation {Γ i } and {γ i } do not equal {Γ i } and {γ i } by the nature of their numerical construction. This lack of strictly constraint-satisfying density matrices motivates the need for a relaxed theorem.
We show (See Appendix D.1) that both strict set membership and imprecise constraints can be described by the inequalities where > 0 is an appropriately chosen number. Determining depends heavily on how {Γ i } and {γ i } are constructed, so we leave out a precise analysis. In general, given a suitable high-precision computing environment, may be made arbitrarily small. (For example, for our calculations in Sec. 5, we choose < 10 −12 .) The bound on the unknown constraint violations motivates the introduction of a relaxed set So long as is larger than the constraint violations, the following relation holds (See Appendix D.1): With this new set, we now present a relaxed version of Theorem 2, with the proof given in Appendix D.3.
where ζ was defined in (19) and Here, the setS * (σ) is defined bỹ Notice that as goes to zero the last term in (25) vanishes and hence which implies that Theorem 3 generalizes Theorem 2. The idea is that Theorem 3 provides a method that is robust to constraint violation due to numerical imprecision. Hence, Theorem 3 is directly useful for numerical key rate calculations. Although it is more complicated than Theorems 1 and 2, Theorem 3 is what we employ in practice for our key rate calculations. For example, the calculations presented in Sec. 5 use this theorem.

Convergence and tightness
Theorem 3 is directly useful in key rate calculations, and as such we discuss some of its convergence properties.
In practice, Step 1 does not return a density matrix ρ * that minimizes f over S. Instead it finds a matrix ρ ∈ S that approximately minimizes f overS . Hence, we answer the natural question of how close the lower bound produced from ρ will be to α = f (ρ * ).
Note that by introducing the perturbed function f , it follows from Lemma 1 that the lower bound produced by Theorem 3 will be continous with respect to its argument. Thus, if ρ * minimizes f overS and ρ is close to ρ * under some norm, it follows that the lower bounds produced by applying Theorem 3 to ρ and ρ * will be close. Concretely, we have showing that our lower bounding method converges (i.e. the optimal lower bound of ρ * is approachable).
In Appendix E we show that the minimizer ρ * of f over S and the minimizer ρ * of f overS satisfy Therefore, the lower bound produced by Theorem 3 is tight when applied to ρ * . Combining (28) and (29), it follows that the lower bound produced by applying Theorem 3 to ρ will be arbitrarily close to f (ρ * ) provided that , are small and ρ is close to ρ * . So provided a suitable computer implementation, our method will produce lower bounds on the true key rate that are arbitrarily tight.

Finding a near-optimal attack
Thus far we focused on Step 2 of the key rate calculation procedure. However, Step 1 needs to be addressed since obtaining a reasonable lower bound from Step 2 requires a ρ ∈ S that is sufficiently close to the optimal solution ρ * . Here we present an algorithm that has proved effective in practice.
Note that in what follows we do not distinguish between exact and inexact representations of {Γ i } or {γ i }, nor do we distinguish between exact and approximate set membership. The validity of this approach is justified by observing that Step 1 is decoupled from Step 2. Specifically, we can apply Theorem 3 to any positive semidefinite matrix that is approximately consistent with the constraints, and still obtain a reliable lower bound.
By applying the Gram-Schmidt process to the original set of observables {Γ i }, we obtain a set {Γ i } of Hermitian operators that are orthogonal under the Hilbert-Schmidt norm. The expectation value of each of these operators is denotedγ for states ρ ∈ S. We extend this set to an orthonormal basis {Γ i } ∪ {Ω j } for the complete Hermitian operator space with suitable operators {Ω j }. With this basis, we may rewrite (4) as where m is the number of free parameters. This perspective on S divides the operator space into "fixed" and "free" subspaces. From a practical optimization pointof-view, the benefit of this representation is that it reduces the original equality constrained problem (which is cumbersome to work with) to a constrained minimization subject to a single semidefinite constraint. We now adapt the Frank-Wolfe method [25] to problem (6). As discussed below, the minimization in Algorithm 1 is a semidefinite program and hence is easier to perform than the original minimization problem in (6).
There are two concrete reasons why adopting the subspace perspective, as in (31), is useful for solving Algorithm 1. First, finding a ρ 0 ∈ S becomes simple. We only need to find ω so that Second, ∆ρ has the form Calculating ω = {ω 1 , ..., ω m } requires solving the standard linear semidefinite program (SDP) Problems of this form have been extensively studied (e.g., see [20]) and efficient SDP solvers are widely available.
4 General framework for protocols Having stated our main result in abstract form, we now connect our result to concrete QKD protocols. Our goal is to present a framework that applies to the known discrete-variable (DV) QKD protocols in the literature. This includes both entanglement-based (EB) and prepare-and-measure (PM) protocols. Examples of protocols that fall under our framework include the BB84 [26], B92 [27], SARG [3,28], six-state [29], and decoy-state protocols [30]. Rather than show how our approach applies to each of these examples, we instead construct a generic protocol that encompasses these examples. We call this the "prototypical protocol", shown in Fig. 2. (In Fig. 2, the distinction between the EB and PM scenarios is depicted by a box around the source with a dashed outline, indicating that the source may or may not be located inside Alice's lab.) We will now show how to apply our approach to this prototypical protocol. We remark that our framework can be easily extended to cover protocols with a central node between Alice and Bob, such as the MDI (measurement-device independent) [31] and the simplified trusted node [32] protocols. We direct the reader to Ref. [19] for a discussion of this extension.
Let us now describe the basic steps involved in our prototypical protocol: 1. In the EB scenario, Alice and Bob each receive a quantum signal (A and B, respectively) from a source and they measure the signal according to the respective POVMs P A = {P A j } and P B = {P B j }, producing the raw data. In the PM scenario, the same mathematical description applies (via the socalled source-replacement scheme [33,34]) since one can think of Alice's prepared states as resulting from Alice performing a measurement P A on a register system A. Figure 2: A prototypical QKD protocol that we use to illustrate our framework. The source sends systems A and B to Alice and Bob, respectively. The dashed line around the source indicates that it may or may not be inside of Alice's laboratory, respectively treating the PM and EB scenarios. Alice and Bob respectively measure POVMs P A and P B . Then they publicly announce some information related to their measurement outcomes. These announcements are used to perform postselection, where some announcement outcomes are discarded. Alice then implements a key map that maps her information (raw data + announcements) to a key variable. Finally, Alice leaks some information (about the result of the key map) to Bob for error correction purposes, and then Bob forms his key.
2. Alice and Bob make a public announcement, announcing some aspect of their measurement outcomes.
3. Alice and Bob perform post-selection based on these announcements.

Alice implements a key map.
The key map is a function that maps Alice's raw data and the announcements to a key symbol, chosen from {0, 1, ..., N − 1} where N is the number of key symbols.
5. Alice performs one-way error correction, leaking some information to Bob, and Bob forms his key.
6. Alice performs privacy amplification, typically by applying a random universal hash function and then communicating the choice of hash function to Bob.

Mathematical model for prototypical protocol
The quantum state shared by Alice and Bob (prior to their measurements) can be written as ρ AB . In the worst-case scenario, Eve possesses a purification of ρ AB denoted as system E, and we assume this worst-case scenario to be as pessimistic as possible. In the following discussion of the protocol, we will note that Eve obtains access to additional systems due to public announcements made by Alice and Bob. (In total, by the end of the protocol, Eve will have access to E A B where A and B are respectively the registers that store Alice's and Bob's public announcements.) Consider the experimental constraints on the state ρ AB . These constraints have the form: For prepare-and-measure (PM) protocols, we add additional constraints, as follows. We employ the sourcereplacement scheme [33,34], which treats system A as a register that stores the information about which state Alice prepared. This corresponds to Alice preparing the bipartite state where {|φ i } are the signal states and {p i } are their associated probabilities. Eve's attack maps system A to system B, producing the state ρ AB . On the other hand, system A is inaccessible to Eve, and hence is fixed, independent of Eve's attack. To fix ρ A we add in constraints of the form where {Θ j } is a finite set of tomographically complete observables on A. Hence, (36) and (39) together represent the constraints that Alice and Bob have on their state.
Step 2 of the above protocol involves Alice and Bob each making an announcement based on their measurement results. In this case, it is helpful to group together the POVM elements into sets associated with particular announcements. Let us write Alice's POVM as Here the first index denotes the announcement, and the second index denotes a particular element associated with that announcement. Bob's announcement is given by a quantum channel with Kraus operators We remark that K B b expands the Hilbert space by introducing the registers B and B, which is why the operators on these subsystems appear as kets in (40). Similarly, Alice's announcement is given by a quantum channel with Kraus operators Here, A and B are registers that store Alice's and Bob's announcements, respectively. Also, A and B are registers that store Alice's and Bob's measurement outcomes, for a given announcement. So, the state after making these announcements becomes where A is a CPTP map. We remark that the form of (42) is such that A and B are classical registers, meaning that the purifying system has a copy of the registers. This is the way one models a public announcement.
Step 3 is post-selection. Here, Alice and Bob select some announcements to keep and some to discard. Let A be the set of all announcements that are kept. Then define the projector with identity acting on the other systems. The postselection is modeled by projecting with this projector to obtain the state where p pass = Tr(ρ A AAB BB Π).
In step 4, Alice chooses a key map. A key map is a function g whose arguments include the outcome of Alice's measurements (a, α a ) and Bob's announcement b. The function outputs a value in {0, 1, ..., N −1} where N is the number of key symbols. Hence, we write the key map as the function g(a, α a , b). We define an isometry V that stores the key information in a register system R, as follows We first act with this isometry on the state ρ (3) in (45), which stores the key information in the standard basis {|j R } on R. Then we decohere R in this basis, which turns R into a classical register denoted Z R , giving the final state where Z is a pinching quantum channel, whose action is given by Z(σ) = j (|j j| R ⊗ 1)σ(|j j| R ⊗ 1).

Key rate
Finally, Alice performs error correction, which gives leak EC obs number of bits about the key map results to Eve, followed by privacy amplification. This gives the following formula for the key rate [21]: Here, H(A|B) ρ = H(ρ AB ) − H(ρ B ) denotes the conditional von Neumann entropy with H(σ) = −Tr(σ log σ).
Note that the first term in (49) refers to the state ρ Z R A AAB BB defined in (48). Also, as noted above, E is a purifying system of ρ AB .
More precisely, the expression in (49) must be minimized over all density operators ρ AB that satisfy the constraints in (36) and (39). Hence we write: where S has the general form in (4), with the constraints given by (36) and (39).
Having now expressed the key rate in explicit form, we can now relate (50) back to the discussion in Sec. 2. This involves simplifying the notation. Namely, for the constraints in (36) and (39), we rewrite them as Tr(Γ i ρ) = γ i , as in (3). Note that we drop the subsystem labels on the state ρ = ρ AB . Finally, we write the optimization problem in (50) as where = p pass · D ρ (4) RA AAB BB ||ρ (5) Note that (53) removes the dependence on Eve's system E and is derived using Theorem 1 from [35]. Equation (54) is derived from the previous line using the property D(cσ||cτ ) = cD(σ||τ ) for any constant c > 0. Furthermore we define G such that its action on an operator σ is given by Note that (54) has the same form as Eq. (5). In summary, to apply our numerical approach to a given protocol, one formulates G via (55) and the constraints via (36) and (39). With these objects defined, one then applies our numerical method outlined in Sec. 3.

Examples
In this section we consider three practically important scenarios that show the power of our approach. In particular, we consider (1) the BB84 protocol with detector efficiency mismatch, (2) the Trojan-horse attack on the BB84 protocol, and (3) the BB84 protocol with phasecoherent signal states. Each of these three scenarios involves a BB84-style protocol but with some "imperfection" accounted for (detector inefficiency, existence of a side channel, and lack of phase randomization). For each scenario, our approach yields significantly higher key rates than those previously obtained in the literature. We remark that the robustness of our numerical approach could allow us to investigate all three imperfections in a single protocol, although for simplicity we consider them separately.
For illustration purposes, for the following examples we assume that error correction is performed at the Shannon limit. This means that the error correction term in (50) can be written as a conditional entropy, with the state ρ (5) defined in (48). Here, Z B can be viewed as the classical register that one would obtain from measuring in the standard basis on B.
Let us briefly remark about the gaps between our upper and lower bounds on the key rate, where the upper and lower bounds are respectively obtained from our "Step 1" and "Step 2" from Sec. 3. In all of our figures below, we only plot our lower bounds, since they are guaranteed to be reliable. But it is worth noting that, for the examples considered below, the gaps between our upper and lower bounds are extremely small. Namely we observed the gap to be ∼ 10 −6 times smaller than the key rate value. Hence for these examples our upper and lower bounds essentially coincide. In these cases, one can think of our Step 2 as simply verifying that Step 1 worked well, i.e., that Step 1 found the optimal attack. We remark that the examples considered below are low dimensional, and it is possible that the gap between the upper and lower bounds could grow with the number of signal states, which could cause Step 1 to be less reliable.
For each example, further details about the constraints used in our calculations can be found in Appendix F.

Efficiency mismatch
Consider a polarization-encoded BB84 protocol, where Bob actively choses his detection basis setting. In this case, Bob's measurement involves two detectors, D 1 and Our method Ref. [  D 2 , that are associated with the two polarization states for a given basis. In practice, it is likely that D 1 and D 2 have different efficiencies, commonly referred to as efficiency mismatch. This mismatch can be further enhanced (by Eve) by manipulating the spatial mode of the incoming light [36]. When efficiency mismatch is large enough, successful hacking strategies on QKD systems have been demonstrated [37]. Furthermore, even with a small amount of efficiency mismatch, the security analysis of QKD becomes difficult to perform.
This motivates the application of our numerical approach the case of detector efficiency mismatch. For simplicity, we consider single-photon signal states, and we assume that no multiple photons arrive at Bob's measurement apparatus. (Our numerical approach can handle multi-photon signals and multi-photon detection events; however, we leave a detailed discussion of the general case for future work.) The single-photon case was previously treated analytically by Fung et al. [18], and hence it provides an opportunity to compare our numerics to the literature.
To further simplify the analysis, we assume one of Bob's detectors is perfect while the other detector has an efficiency η. This allows us to plot the key rate as a function of η, as shown in Fig. 3 Figure 4: Key rate vs error rate for the single-photon BB84 protocol under a Trojan-horse attack. The key rate is plotted for different values of µout. Our numerical method improves on our previous approach in Ref. [19], which in turn gives higher key rates than the analytical method of Ref. [15].
channel is the actual attack; we simply use this channel to generate artificial data for Alice and Bob.) The plot shows that our new numerical method outperforms our previous numerical method based on the dual problem [19], which in turn outperforms the analytical method from Ref. [18].

Trojan-horse attack
Consider the phase-encoded BB84 protocol. Here, Alice's light source produces a pulse that passes through an interferometer, one arm of which applies a variable phase θ chosen from the set {0, π/2, π, 3π/2} to encode the information. Bob decodes this phase information with an interferometer in his lab.
There is a simple hacking attack on this protocol that exploits a side channel in Alice's encoder (i.e., a channel by which Eve can obtain additional information, beyond the direct channel from Alice to Bob). The attack involves Eve sending a bright pulse of light into Alice's lab [13]. Some fraction of this pulse reaches Alice's phase encoder and is encoded with the same information that Alice is attempting to send to Bob. A portion of this light is reflected back to Eve, who can then decode some of the phase information, potentially compromising the protocol's security. This is called the Trojan-horse attack (THA) [14], since it involves a "malicious gift" from Eve.
For simplicity, let us restrict our attention here to the case where Alice's light source outputs only a single photon per signal. Following the approach of Lucamarini  et al. [15], we describe Eve's input pulse and backreflected pulse as coherent states denoted by | √ µ in and |e iθ √ µ out , respectively. Here, θ stores Alice's phase encoding setting, and the input and output intensities typically satisfy µ out µ in . The signal states emerging from Alice's source are: where Here, |n L (|n M ) is the n-photon state of the long (short) arm of the interferometer. Figure 4 shows the results of our numerics for the THA. Here we plot key rate versus error rate (assuming the z and x error rates are identical, for simplicity) for various values of µ out . The plot shows that our new numerical method gives higher key rates than both the analytical method from Ref. [15] as well as the numerical method from Ref. [19].

BB84 protocol with phase-coherent signal states
Here we consider a protocol proposed by Huttner et al. [16] and analyzed by Lo and Preskill [17]. This is a phase-encoded BB84 protocol, but using coherent states instead of single-photon states. This is quite practical, since one can use attenuated laser pulses from modelocked lasers to generate these signal states. In addition the protocol is practical because the experimenter does not need to do phase randomization. So it is worth investigating the key rate of this protocol.
The signal states prepared by Alice are [17] |φ z+ = | + α S ⊗ |α S (64) where α is the amplitude of the coherent state, S is the signal mode and S is the reference mode. When Bob receives the signal, he performs a polarization measurement (in one of two complementary bases), discarding no-click events, and assigning a random bit value to double-click events. Lo and Preskill [17] gave an analytical lower bound on the key rate for this protocol, as a function of transmission probability η and amplitude α. Their theoretical curves are shown as dotted lines in Fig. 5, for several values of η. In the same plot, we show the result of our numerical optimization as solid lines, with the key rates obtained from the method in Ref. [19] shown as dasheddotted lines. Interestingly our numerics give higher key rates than the previous literature over the entire parameter range. This is an important result due to the practicality of this protocol.

Conclusions
In conclusion, we presented a new numerical approach for calculating key rates for QKD. For concreteness, we name our approach the "reliable primal method" or the "reliable primal problem". As discussed in Sec. 3, Step 1 of our method is simply the primal optimization problem.
Step 2 of our method converts the output of Step 1 (a nearly optimal eavesdropping attack) into a reliable lower bound on the key rate. We presented an efficient method for Step 1 in Sec. 3.5. Our main contribution is a method for Step 2, which is presented in Theorems 1, 2, and 3.
Reliability is the most important issue with numerical key rate calculations, since key rates must come with a security guarantee. In this work, we highlighted the various issues associated with numerical key rate calculations, such as constraint violation and inexact variable storage by computers. Furthermore, we showed how to address these issues. Our most robust result, Theorem 3, allows one to lower bound the key rate despite numerical imprecision.
We discussed that our method is arbitrarily tight in Sec 3.4. This allowed us to make significant improvements over previous literature key rates for three interesting examples in Sec. 5. Furthermore, the tightness of our approach implies that the solid curves that we plotted in Figs. 3, 4, and 5 are essentially unbeatable, i.e., they cannot be improved upon. Eliminating looseness from key rate calculations is a major advance for the field of QKD research.
Future applications of our work include investigating device imperfections, side channels, multi-photon detection events, decoy-state protocols with partial phase randomization, measurement-device independent protocols, differential-phase shift protocols, and coherent one-way protocols. Perhaps more importantly, our approach can allow researchers to explore and evaluate novel protocol ideas that have yet to be discovered.
As noted in the Introduction, our group released a user-friendly software for key rate calculations based on the dual problem from Ref. [19]. Interestingly, our reliable primal method presented here improves on the approach of Ref. [19] in terms of both speed and tightness. Therefore, we plan to improve our publicly-available software in the future by incorporating the reliable primal method. We believe this software has the potential to be used throughout the QKD community, both in industry and academia.

A Proof of Theorem 1 A.1 Standard form for semidefinite programs
To prove Theorem 1, we will make use of the so-called standard form for semidefinite programs (e.g., see page 265 of [20]), given as follows: subject to: subject to: Here, the inner product is of the Hilbert-Schmidt form, A, B := Tr(A † B), and the vector notation a· b is shorthand for m j=1 a j b j . If µ and ν are, respectively, the solutions of the primal and dual problems, then in general one has which is known as weak duality [20]. Furthermore, under certain conditions, (69) turns into an equality µ = ν, known as strong duality. A sufficient condition for strong duality to hold is Slater's condition (e.g., page 265 of [20]), which can be stated as follows for the standard-form SDP in (68). Slater's condition: Let X = {X 0 : B j , X = b j , ∀j} be the set of feasible X for the primal problem in (68). Suppose that X = ∅. Furthermore, suppose that there exists a y ∈ R n such that A > y · B. Then µ = ν . (70)

A.2 Inequality in (8)
In what follows, we first prove the inequality in (8), and in the next subsection we establish the equality condition. For some matrix σ, let σ := vec(σ) denote the vectorization obtained by stacking the columns of σ. For two matrices σ and τ , note that the inner product between their vectorizations can be written as where T is the transpose taken in the basis used to define the vectorization. Now consider the gradient matrix ∇f (ρ) defined in (7) and its vectorization ∇f (ρ). Note that the quantity quantifies how much the function f changes whenever one moves from point ρ to point σ.
[Note that Eq. (74) rewrote g(ρ, σ) in matrix notation, where T is the transpose in the same basis that is used to represent the gradient matrix.] More precisely, g(ρ, σ) quantifies how much the linearization L ρ of f changes, where L ρ is the linearization about point ρ. Specifically, the linearization is given by The assumption of Thm. 1 is that f is differentiable about the point ρ of interest. Since f is a convex, differentiable function over a convex set S, the linearization L ρ always lies below the curve f (e.g. page 69 of [20]). Hence, for any two points ρ, σ ∈ S we have Now let ρ * ∈ S minimize f over S. Then where (78) exploits the fact that ρ * ∈ S. Hence finding a lower bound on α reduces to the minimization problem This is a linear semidefinite program (SDP) and we may apply duality theory to obtain the dual SDP. In particular, our problem is essentially the standard SDP form given in Sec. A.1, which gives the following dual problem where Weak duality from (69) implies that Inserting (84) into (79) gives the desired lower bound in (8).

A.3 Equality in (8)
With the inequality in Theorem 1 proven, we now show that equality in (8) holds if ρ satisfies f (ρ) = f (ρ * ) = α. First, consider the inequality in (84). Slater's condition stated in Sec. A.1 provides a sufficient criteria for strong duality to hold and, hence, for (84) to be an equality. To satisfy this condition it is adequate to show that S = ∅ and there exists y ∈ R n such that i y i Γ T i < ∇f (ρ). Since the set of constraints {Γ i } correspond to a valid density matrix, it immediately follows that S = ∅. Since density matrices are constrained to have trace one, without loss of generality we may take Γ 1 = 1 and γ 1 = 1. Thus, if λ min is the smallest eigenvalue of ∇f (ρ), it follows that (λ min − 1)Γ T 1 < ∇f (ρ). So y = (λ min − 1, 0, ..., 0) T satisfies i y i Γ T i < ∇f (ρ). With Slater's condition satisfied, strong duality holds and min σ∈S Tr(σ T ∇f (ρ)) = max We now show that if f (ρ) = f (ρ * ) then equality in (8) holds. Suppose that f (ρ) = f (ρ * ), then (79) yields Next, we state a lemma that provides a bound in the opposite direction of (86).

Lemma 2:
For a point ρ that minimizes f over S, Proof. The Karush-Kuhn-Tucker (KKT) conditions (e.g. page 243 of [20]) provide necessary conditions for the optimality of ρ. For our problem they say that if ρ is optimal, then there exists a pair ( λ, Z) ∈ R n × H where H denotes the set of d × d Hermitian matrices such that Let σ ∈ S then where we have used the definition of S and (89) in the last equality. Since σ and Z T are both positive semidefinite it follows that the trace of their product is nonnegative. Hence, and since σ is arbitrary, the desired result follows.

B Proof of Lemma 1
As noted in (15), the gradient of interest has the form To show that this expression is well-defined, we simply need to show that logarithms can be evaluated, which would be true if G (ρ) and Z(G (ρ)) are full-rank matrices. Consider the first of these matrices Since G(ρ) 0, it follows that In other words, G (ρ) is full rank. Now consider the second of these matrices where the second line notes that Z is a pinching quantum channel and hence Z(1) = 1. Since Z(G(ρ)) 0, we So both G (ρ) and Z(G (ρ)) are full rank, implying that (99) is well-defined.

C Proof of Theorem 2 C.1 Continuity in ρ
Here we state that the objective function f (ρ) is continuous, which will be useful for our proof of Theorem 2. First we state two helpful lemmas for the trace distance (or trace norm). The first lemma was proved in Ref. [41].

Lemma 3:
Let ρ 0 and σ 0 be positive semidefinite matrices. Then where The next lemma is a corollary of Lemma 3.

Lemma 4:
Let ρ 0 and σ 0 be positive semidefinite matrices. Let E be a completely positive trace-nonincreasing (CPTNI) map. Then where Proof. Let ρ − σ = Q − S, where Q 0 and S 0, and Q and S have orthogonal support. Hence Here, (109) used the fact that E is trace-nonincreasing. Also, (110) follows directly from (105).
Next we state a lemma, known as Fannes' inequality, that entropy is continuous. The proof can be found, e.g., in Ref. [42].

C.2 Continuity in
Here we show another sort of continuity, namely, the continuity of f (ρ) in . First we need the following lemma. Then we have Proof. Note that (121) follows from (120) due to Lemma 4. So we just need to prove (120). Towards this end, we define with ∆ := G(ρ) − 1/d . Now note that 0 G(ρ) 1, which implies This means that the d eigenvalues {λ j } of ∆ each satisfy |λ j | (1 − 1/d ). Hence Combining (127) with (124) proves (120). Now we prove our desired continuity statement.
where ζ := 2 (d − 1) log d (d −1) . Proof. We begin by expanding the expression of interest as follows Next we exploit Lemma 7 and Lemma 5. In particular, applying Lemma 5 with κ = (d − 1) and n = d gives .
Finally, we obtain

C.3 Proof of Theorem 2
We first define a perturbed optimization problem. Let ∈ R such that 0 < < 1/(de). Then consider where f was defined in (14) as Next we relate α( ) to α with the following lemma.
where ζ : Proof. Note that (136) is equivalent to Now let α = f (ρ * ) and α( ) = f (ρ * ), where ρ * ∈ S and ρ * ∈ S. We prove (137) in two steps. First we have where the last line follows from Lemma 8. Second we have where again the last line follows from Lemma 8. Taken together, (140) and (143) give the desired result.
Lemma 9 connects α to α( ). So now we can proceed to lower bound α( ) to ultimately get a lower bound on α.
The key observation is that α( ) corresponds precisely to α provided that one makes the substitution G → G . Since this change simply involves substituting a different quantum channel for G, the lower bound that we previously found on α in Theorem 1 directly applies, with the appropriate substitution made. Recall that in Theorem 1 we found that, for any ρ ∈ S such that ∇f (ρ) exists, where This result can be directly applied find a similar lower bound on α( ). One simply needs to substitute G for G, which is equivalent to substituting f for f . Hence we obtain where It is crucial to note that, while (144) applies only if ∇f (ρ) exists, (146) applies to all ρ ∈ S since ∇f (ρ) always exists. Furthermore, we remark that the argument in Sec. A.3 directly applies to α( ), so that we have where ρ * is a state that satisfies α( ) = f (ρ * ). Finally, we combine (146) with Lemma 9, specifically (143). This gives the following bound for any ρ ∈ S: But as illustrated in Fig. 6,S does not coincide with the true set of interest S. To address this, we will expand the setS until it encompasses S, as the set S does in Fig. 6. We begin by relating the approximate and exact representations bỹ where ||δΓ i || HS < 1 and |δγ i | < 2 , for all i. Here, the Hilbert-Schmidt norm is defined by ||A|| HS := Tr(A † A). Determining the constants 1 and 2 is rather technical and depends heavily on how the approximate observables and expectation values are computed. However, if the mantissa of the underlying representation is increased in length (or arbitrary precision arithmetic is used) and appropriate numerical algorithms are applied, the approximate representations will become arbitrarily accurate. Hence, 1 and 2 may be made as small as needed.
We define the quantity which measures our overall uncertainty in the variable representation. We now prove a lemma that motivates our treatment of the imprecision. Namely, it allows us to upper bound the numerical constraint violation.
Proof. We can apply the triangle inequality, the Cauchy-Schwarz inequality and the fact that ρ is a density matrix to obtain Motivated by Lemma 10, we consider the set where it is clear that lim →0S =S. Furthermore, if we choose = rep , we obtain the set which follows directly from Lemma 10. Hence,S rep encompasses both S andS. So we could useS rep as the basis for our optimization. However, we have not yet accounted imprecision in the numerical solver. As discussed in the next subsection, we may need to expandS rep further to allow for solver imprecision.

D.2 Imprecise solvers
Up to this point in our analysis, we have neglected the fact that no numerical solver is exact. In reality, the matrix ρ returned by the solver may not be positive semidefinite or satisfy the approximate constraints defined by {Γ i } and {γ i }. Here we present a method for addressing this situation. First, we add a bit of the identity toρ to obtain a positive semidefinite matrix. Let λ min denote the smallest eigenvalue ofρ, thenρ is positive semidefinite. Second, we expand the set over which we optimize such that it encompassesρ , as illustrated in Fig. 6. Let sol be a positive real number such that ρ ρ * Figure 6: Illustration of our method for handling numerical imprecision. Ideally we want to represent the set S in the computer. However, since the computer has finite precision, it storesS that approximates S, whereS is defined by (151). Furthermore, when we attempt to optimize overS, the answerρ returned by the computer may not even be inS. We growS into a new setS that contains both S andρ , whereS is defined in (161). (Note that lim →0S =S.) With this expanded setS , we may then find a reliable lower bound by applying the technique discussed in the main text (Sec. 3.3).

S SS
for all i. The quantity sol describes how closeρ is to satisfying the approximate constraints provided that it is chosen to be as small as possible. The closerρ is to satisfying the constraints, the smaller sol will be.
We now have two quantities: rep describes the representation precision and sol describes the solver precision. We define the quantity The reason for defining is that it specifies how large of a set we need to contain both the imprecise solution returned by the solver and the exact solution. If we can construct such a set, then we may relax both Theorem 1 and Theorem 2 to obtain the computationally feasible Theorem 3.
Recall the relaxed set of approximate density matrices defined in (161), rewritten here for convenience: From Lemma 10, it follows that S ⊆S , and (164) implies thatρ ∈S . This is illustrated in Fig. 6. In summary, our method for handling numerical imprecision involves minimizing the objective function f (ρ) over S . Because S ⊆S , our method will yield a reliable lower bound on α = min ρ∈S f (ρ).

D.3 Proof of Theorem 3
Our proof strategy for Theorem 3 will be similar to that for Theorem 2, in that we will consider a slightly perturbed optimization problem to ensure a well-defined gradient, and we will invoke continuity to argue that the perturbed problem is close to the problem of interest. The only difference is that now the problem of interest will allow for small violations of the constraints.
Of course, we are ultimately interested in lower bounding the quantity defined in (6), But due to numerical imprecision (discussed above), we must consider a slightly different problem, whereS was defined in (22). Fortunately, we have that S ⊆S (see previous subsection), which implies that Hence, any lower bound on α is also a lower bound on α. So let us now find a lower bound on α .
As noted above, we will in fact consider a further perturbation on the optimization problem, to ensure a welldefined gradient. Again, this is analogous to what we did in Theorem 2. This perturbed optimization problem has the form where f was defined in (14). Suppose that ρ * ∈S achieves this optimization, i.e., α ( ) = f (ρ * ). Next, we consider any ρ ∈S , and we invoke an argument identical to that used in Eqs. (77)-(79) to write Next we rewrite the minimization in (171) as a maximization via duality.
Lemma 11: Let ρ ∈S , then where v( γ, ) := ( γ + ) ⊕ (− γ + ) is a 2n-dimensional vector, and Here, 0 is the n × n matrix of zeros, and i = |i i| is the projector onto the i-th element of the standard basis in n-dimensions.
Proof. First we rewrite the minimization problem as min σ∈H Tr(σ T ∇f (ρ)) (177) Next, we utilize slack variables (e.g., see page 131 of [20]); the idea is to replace each inequality constraint with an equality constraint, and a nonnegativity constraint. Letting a and b denote slack variables, the problem becomes Next we recast the problem so that there is again one positive semidefinite variable. First we define the following n × n diagonal matrices: where 0 is the n-dimensional zero vector, and |i i| is the projector onto the i-th element of the standard ndimensional basis. Next we define the block-diagonal matrices The optimization problem then becomes This is a semidefinite program of a form identical to the one that appears in the proof of Theorem 1. Duality theory, as discussed briefly in Sec. A.1, yields the dual problem where v( γ, ) := ( γ + ) ⊕ (− γ + ) is a 2n-dimensional vector.
Next we rewrite the optimization problem in a simpler way, as follows.
whereS * (σ) := ( y, z) ∈ (R n , R n ) | − z y z, Proof. First, note that we can rewrite the setŜ * (ρ) as followŝ This formula is derived by breaking up the constraints in (173) associated with different blocks. The first block gives the constraint The second block gives the constraint n i x i |i i| 0, and the third block gives the constraint n i x i+n |i i| 0. These latter two constraints imply that x i 0 for all i ∈ [1, 2n], or in other words, that x 0.
and in (200) we used the fact that lim →0 ρ * = ρ * . Now by (148) This gives Finally, Lemma 9 implies that lim , →0 Thus, given suitably small , it follows that the bound produced by Theorem 3 applied to a near-optimal state is arbitrarily tight.

F.1 Efficiency Mismatch
For the BB84 protocol with detector efficiency mismatch, we model it as an entanglement-based protocol. We write Alice's POVM as where {|0 , |1 } is the z-basis on a qubit, and |± = 1/ √ 2(|0 ± |1 ). Here p z denotes the probability for Alice to measure in the z-basis. For our numerics, we chose p z ≈ 1, corresponding to using the z-basis most of the time.
We model Bob's system as a qutrit, where the one-photon subspace is modeled as a qubit subspace, and the third dimension is the vacuum. This third dimension is incorporated because of detector inefficiency, which may cause a no-click event. Bob's POVM elements associated with detecting a photon are given by where the direct sum is used here to embed qubit operators inside a qutrit Hilbert space. Note that P B 2 and P B 4 have a factor of η due to detector inefficiency. (We assume one detector has perfect efficiency, while the other has efficiency η.) The last of Bob's POVM elements corresponds to a no-click event, To generate the constraints in (36), we simulate the data using a depolarizing channel with depolarizing probability p, We consider the bipartite state generated from applying this channel to half of a maximally entangled state |Φ , and we emphasize that this state is only used to simulate experimental data. To obtain the constraints in (36), we compute p jk = Tr((P A j ⊗ P B k )ρ sim AB ) .
Since Alice's density operator is fixed, we add the additional constraints specified by (39). Now we consider the announcements made by Alice and Bob. Alice announces her choice of basis, so (41) becomes (We remark that, in this case, introducing the additional register system A is redundant since the key information can be read off directly from system A, but we do it here for completeness.) Similarly, Bob announces his choice of basis, so (40) becomes For the post-selection, Alice and Bob discard events where they measure in different bases. So (44) becomes Finally, consider the isometry V associated with the key map defined in (46). We can define the key map such that Alice stores 0 (1) in her key when she obtains outcome P A 1 or P A 3 (P A 2 or P A 4 ). This gives with identity acting on all other subsystems. The above expressions allow one to define G in (55), and hence define the optimization problem.

F.3 BB84 protocol with phase-coherent signal states
We model the BB84 protocol with phase-coherent signal states as a prepare-and-measure protocol with sifting, similar to how we modeled the Trojan-horse attack above. We apply the source-replacement scheme as described in Sec. 4, with the state |ψ AA from (37) given by where {|φ z± , |φ x± } are specified in (64)-(67). Here, p z denotes the probability of Alice preparing a state in the z-basis, and it is biased to be close to one. Alice's POVM acts on the register system A with the standard basis elements, namely By applying a squashing model [44], we model Bob's system as a qutrit, where the one-photon subspace is modeled as a qubit subspace, and the third dimension is the vacuum. This third dimension is incorporated because of channel loss, which may cause a no-click event. Bob's POVM elements are then Next we consider data simulation for the purpose of formulating the constraints in (36). We model the channel between Alice and Bob as a lossy channel E loss (ρ) with transmission probability η. Note that the action of this channel on a coherent state is |α → | √ ηα . We can apply the channel to the state |ψ AA to obtain the bipartite state ρ sim AB = (I ⊗ E loss )(|ψ ψ| AA ) .
The remainder of the model is identical to that of Sec. F.2 from (223) onwards. This defines the optimization problem.