Practical decoy-state method for twin-field quantum key distribution

Twin-Field (TF) quantum key distribution (QKD) represents a novel QKD approach whose principal merit is to beat the point-to-point private capacity of a lossy quantum channel, thanks to performing single-photon interference in an untrusted node. Indeed, recent security proofs of various TF-QKD type protocols have confirmed that the secret key rate of these schemes scales essentially as the square root of the transmittance of the channel. Here, we focus on the TF-QKD protocol introduced by Curty et al, whose secret key rate is nearly an order of magnitude higher than previous solutions. Its security relies on the estimation of the detection probabilities associated to various photon-number states through the decoy-state method. We derive analytical bounds on these quantities assuming that each party uses either two, three or four decoy intensity settings, and we investigate the protocol's performance in this scenario. Our simulations show that two decoy intensity settings are enough to beat the point-to-point private capacity of the channel, and that the use of four decoys is already basically optimal, in the sense that it almost reproduces the ideal scenario of infinite decoys. We also observe that the protocol seems to be quite robust against intensity fluctuations of the optical pulses prepared by the parties.

The last few decades have witnessed major advancements in the field of quantum communication [1,2], with quantum key distribution (QKD) [3][4][5][6][7][8][9][10][11][12][13] being its most developed application. Recent experiments over about 400 km of optical fibers [14,15] and over about 1000 km of satellite-to-ground links [16,17] demonstrated that QKD over long distances is possible. Despite such remarkable experimental achievements, the private capacity of point-to-point QKD is intrinsically limited by fundamental bounds [18,19]. These bounds state that in the high-loss regime the key rate scales basically linearly with the transmittance of the channel connecting the end-users Alice and Bob, i.e. it decreases exponentially with the total channel length. This imposes strict practical constraints on the possibility of achieving point-to-point QKD over arbitrary long distances. A way to overcome this limitation is to employ one or more intermediate nodes in the quantum channel connecting the parties. For instance, the use of quantum repeaters [20] yields a polynomial scaling of the communication efficiency with the distance [21]. Moreover, a quantum repeater scheme can be arbitrarily iterated along the quantum channel, thus increasing in principle the total communication distance between Alice and Bob as much as desired. Unfortunately, however, quantum repeaters are very challenging to build in practice with current technology: they either require quantum memories [20][21][22] or quantum error correction [23,24]. Of course, technology is improving, and quantum repeaters may become viable in the future. Other solutions, which attain a square-root improvement in the scaling of the key rate with respect to the transmittance of the channel, are obtained by placing a single untrusted relay between Alice and Bob. Such protocols include, for instance, Measurement-Device-Independent-QKD [6] (MDI-QKD) with quantum memories [25,26] and adaptive MDI-QKD featuring quantum non-demolition measurements [27]. The philosophy behind both types of protocols is that the central relay is able to adapt the pairings of photons received from Alice and Bob to the photon losses. In this way, for every signal sent by Alice and Bob to the central relay, just one of the two signals is required to arrive, leading to the mentioned square-root improvement in the key rate scaling. However, both protocols still require two-photon interference in the central node, as in the original MDI-QKD scheme [6]. More recently, [28] proposed the Twin-Field (TF) QKD protocol, still characterized by an untrusted central node, and conjectured a square-root improvement in the key rate scaling. This scaling has been later on confirmed in [29,30] for two variants of the original scheme. The advantage of TF-QKD lies in the fact that it is designed to generate key bits from single-photon interference in the central node, thus naturally retaining the scaling with the square-root of the transmittance without the need to adapt to photon losses via sophisticated devices. Since the original proposal, there has been an intense research activity to develop different versions of TF-QKD protocols equipped with their security proofs [29][30][31][32][33] as well as to investigate their experimental feasibility [34][35][36]. Among these protocols, the one that seems to deliver the higher secret ket rate [37] is that introduced in [33]. Its security relies on the ability to estimate the detection statistics (usually called yields) of various Fock states sent by Alice and Bob through the decoy-state method [38][39][40]. The key-rate simulations provided in [33] indeed exhibit an improved scaling with the loss, but the estimation of the yields is only carried out by means of numerical tools based on linear programming and considering only the case of three decoy intensity settings. In this paper, we derive analytical bounds on the yields which are required to evaluate the key rate formula of [33], assuming two, three and four decoy intensity settings. In so doing, we are able to show, for instance, that the use of two decoy intensity settings is already enough to beat the point-to-point private capacity bound reported in [19]. Also, we show that the use of four decoys is basically optimal in the sense that the resulting secret key rate is already very close to the ideal scenario which assumes infinite decoy intensity settings. Analytical bounds imply a fully-analytical expression for the protocol's secret key rate, which could be very convenient for performance optimization in scenarios where the number of parameters is high, like for instance in finite-key security analyses. In addition, we study how the performance of TF-QKD is affected under intensity fluctuations, which are inevitable in practice, and we demonstrate that the protocol in [33] seems to be actually quite robust against such fluctuations.
Like in [33], for simplicity, here we focus on the asymptotic-key rate scenario. However, we remark that by using the techniques reported in [41], it is cumbersome but straightforward to adapt our analytical methods also to the finite-key rate scenario, where, as mentioned above, it becomes particularly useful to have analytical bounds for the main quantities that enter the key rate formula. The article is structured as follows. In Sec. 1 we present the TF protocol from [33] and highlight the main yields that need to be bounded. In Sec. 2 we provide the analytical bounds on the yields for the case of two decoys (the cases of three and four decoys are treated in Appendix C and Appendix D, respectively). In Sec. 3 we provide simulations of the secret key rate versus the loss for a typical channel model (briefly described in Appendix A), and we also evaluate the effect of intensity fluctuations. We conclude the paper in Sec. 4.  As discussed above, we consider the TF-QKD protocol presented in [33] and sketched in figure 1. Alice and Bob establish a secret shared key by sending optical pulses to a central untrusted node, C. It is assumed that the node C shares a phase reference with Alice and Bob, which can be achieved by the transmission of strong optical pulses. The protocol is composed of the following five steps:

The TF-QKD protocol
(i) Alice (Bob) chooses the X-basis with probability p X and the Z-basis with probability p Z = 1 − p X .
Upon choosing the X-basis, Alice (Bob) prepares an optical pulse in a coherent state |α or |−α at random, corresponding to the key bit b A = 0 (b B = 0) or b A = 1 (b B = 1), respectively. Upon choosing the Z-basis, she (he) prepares an optical pulse in a phase-randomized coherent state: A n! |n n| (1.1) (ρ β B ) whose intensity β 2 A (β 2 B ) is drawn randomly from a set S A = {β 2 i } i (S B = {β 2 j } j ) of real nonnegative numbers.
(ii) Both parties send their optical pulses to the untrusted node C via optical channels in a synchronized manner.
(iii) The central node C applies a balanced beamsplitter to the incoming pulses and features two threshold detectors at its output ports. The detector placed at the output port associated to constructive (destructive) interference is denoted by D c (D d ).
(iv) The node C announces the measurement outcome k c (k d ) of detector D c (D d ), with k c = 0 and k c = 1 (k d = 0 and k d = 1) corresponding to a no-click and a click event, respectively.
(v) Alice and Bob form their raw keys with the bits b A and b B collected when both parties chose the X-basis and node C reported a click in only one detector (k c + k d = 1). Bob flips his bits b B for which the click occurred in D d .

Secret key rate formula
The security analysis performed in [33] yields the following lower bound on the asymptotic key rate R: where the terms R kck d , for (k c , k d ) ∈ {(1, 0), (0, 1)}, are defined as: being the binary entropy function, f the inefficiency function associated to error correction, and p(k c , k d ) the conditional probability that node C announces the outcome (k c , k d ) when both parties selected the X-basis. The probability p(k c , k d ) can be expressed as: is the joint probability of Alice and Bob preparing the coherent states |(−1) b A α and |(−1) b B α , respectively. According to the protocol description above, we have: instead denotes the conditional probability that node C announced (k c , k d ) given that Alice and Bob sent the coherent states |(−1) b A α and |(−1) b B α , respectively. Since we consider the asymptotic key-rate scenario, we assume that p(k c , k d |b A , b B ) coincides with the correspondent distribution observed by the parties. Finally, the terms e kck d and e ph kck d in (1.3) represent the bit-error rate in the X-basis and an upper bound on the phase-error rate, respectively. The former is defined as: 6) and the latter as: where the coefficients c n are defined as c n = e −α 2 2 α n / √ n! and the yields Y kc,k d nm are the conditional probabilities that node C announces the outcome (k c , k d ) given that Alice and Bob emitted an n-photon state and an m-photon state, respectively. Note that the only yields contributing to (1.7) are those Y kc,k d nm such that n + m is an even number. The yields Y kc,k d nm are quantities that are not directly observed by the parties, however they can be estimated either numerically or analytically with techniques based on the decoy-state method [38][39][40].
Here we consider the analytical approach. In particular, we assume that Alice and Bob have at their disposal either two, three or four decoy intensity settings when choosing the Z-basis. To each further decoy intensity correspond additional linear constraints on the yields, leading to tighter estimations of Y kc,k d nm and thus to a higher key rate. However, a finite number of decoys only allows to derive non-trivial upper bounds ‡ on a limited number of yields in (1.7), whereas the other yields are set to 1. Nevertheless, even bounding just four yields in a non-trivial way is enough for the secret key rate to beat the pointto-point private capacity bound (PLOB bound) [19] at high losses (see Sec. 3). Also, as we show below, with four decoy intensity settings one can already obtain a secret key rate very close to that achievable with infinite decoy intensity settings. We remark that standard decoy-state-based QKD protocols require to lower bound the value of a few yields (typically those associated to vacuum and single-photon pulses) [42], while the TF-QKD protocol considered here upper bounds the value of the phase-error rate (1.7) by upper bounding several yields. In particular, we upper bound the yields Y kc,k d nm for (n, m) ∈ I, where I is a certain subset of {(n, m)| n, m ∈ N 0 } which depends on the number of decoys. Thanks to the derived upper bounds on the yields (which we shall denote by Y U, kc,k d nm ) we are able to estimate the phase error rate (1.7) as follows:

Yields estimation
When both Alice and Bob choose the Z-basis in the first step of the TF-QKD protocol, they prepare phase-randomized coherent states with intensities β 2 A and β 2 B , respectively, and send them to C. From Eve's viewpoint, she cannot distinguish this scenario from the case in which the parties prepared number states |n and |m according to the Poissonian distributions P β 2 A (n) and P β 2 B (m) (see Eq. 1.1), where P µ (n) = e −µ µ n /n!. Therefore Eve's attack can only depend on the number states |n and |m but not on the signals' intensities β 2 A and β 2 B . As a consequence, the probability that Eve announces outcomes (k c , k d ) only depends on the number of photons (n, m) she received from Alice and Bob, i.e. the yields Y kc,k d nm are independent of the decoy intensities chosen by the parties. For this reason, one can derive a set of linear constraints on the yields Y kc,k d nm by expressing the experimentally observed gains Q kc,k d -which are defined as the conditional probabilities that node C announced the outcome (k c , k d ) given that Alice and Bob sent phase-randomized coherent states of intensities β 2 A and β 2 B , respectively-in terms of the yields: As it is clear from (2.1), to every distinct pair of decoy intensities (β 2 A , β 2 B ) corresponds a new constraint on the set of infinite yields {Y kc,k d nm } n,m , which leads to tighter upper bounds and thus to a higher secret key rate. On the other hand, having a large number of decoy intensities is experimentally demanding, hence the need to derive the tightest possible bounds on the yields with a limited number of decoys.
In this Section we present a simple analytical method to obtain tight bounds on the yields of largest ‡ Every yield is a probability, thus it is trivially bounded by 1. 4 contribution § in (1.7) -i.e. relative to the largest coefficients c n -when the parties use two intensity settings in the Z-basis. It is basically a Gaussian elimination-type technique but involving infinite-size coefficient matrices. In particular, the guiding principle that we use is to combine the constraints (2.1) so that in the resulting expression the yield to be bounded is the one with the largest coefficient, while the yields which had larger coefficients in the initial constraints have been removed in the combination. However, in some cases it turns out that is not possible to remove all the yields with larger coefficients than the one to be bounded, due to a lack of decoy intensity settings (i.e. constraints). In other cases, we manage to remove from the resulting expression even some yields which had a smaller coefficient than the one to be bounded. Such a procedure can be readily extended to the case of three and four decoy intensity settings. The results for these last two cases are presented in Appendix C and Appendix D, respectively. From now on, we assume that both optical channels linking the parties to the central node C have the same transmittance √ η. Therefore the set of optimal decoy intensities β 2 A and β 2 B is the same for both parties [43] and we define it as: {µ 0 , µ 1 }. In order to simplify the notation, we also omit the measurement outcome (k c , k d ) from the constraints given by (2.1). Hence the yields are subjected to the following four equality constraints: and to the inequality constraints: Below we derive upper bounds on the yields: Y 00 , Y 11 , Y 02 and Y 20 .

Upper bound on Y 11
Consider the following combination of gains: The subscript in G 11 indicates the yield that is going to be bounded with this combination of gains. In (2.4) the coefficients of the yields Y 0m and Y n0 , for any n and m, are identically zero. Thus (2.4) can be rewritten as: We observe that the coefficients that multiply the yields Y nm are always positive, being the product of two factors of equal sign. A valid upper bound for Y 11 is obtained considering the worst-case scenario for the other yields, taking into account that (2.3) holds. Since all the yields' coefficients carry the same sign in (2.5) -regardless of the relation between µ 0 and µ 1 -, the yield Y 11 is maximal when all the other yields are minimal. Thus the upper bound on Y 11 is extracted by setting all the other yields to zero in (2.5): where G 11 is defined in (2.4).
We remark that by combining the gains as in (2.4), we manage to obtain a closed expression for Y 11 in which the contribution of all the yields Y 0m and Y n0 is removed. Additionally, Y 11 is now the yield with the "highest weight" in (2.5) since it has the largest coefficient. All the yields' bounds presented in this work follow the same philosophy. § The same method can -in principle-be applied to any yield, however the limited number of decoy settings prevents from obtaining a non-trivial bound on every yield.

Upper bound on Y 02
Consider the following combination of gains: In (2.7) the coefficients of the yields Y n0 and Y 1m are identically zero. Thus (2.7) can be rewritten as: Like in the derivation of Y 11 's bound given by (2.6), a valid upper bound for Y 02 is obtained by considering the worst-case scenario for the remaining yields in (2.8). More specifically, Y 02 is maximal when the yields whose coefficient has the same sign as Y 02 's coefficient are minimal, and the yields whose coefficient has opposite sign to Y 02 's are maximal. Recalling constraint (2.3), this means setting Y 01 and Y 0m to zero and Y nm with n ≥ 2 and m ≥ 1, to 1 in (2.8). In so doing, after rearranging the terms we obtain: which leads to the following upper bound on Y 02 : (2.10)

Upper bound on Y 20
Consider the following combination of gains: In (2.11) the coefficients of the yields Y n1 and Y 0m are identically zero. Thus (2.11) can be rewritten as: A valid upper bound for Y 20 is obtained by setting to zero the yields whose coefficient has the same sign as Y 20 's coefficient, and by setting to 1 the yields whose coefficient has opposite sign to Y 20 's. In the case of (2.12) this means setting Y 10 and Y n0 to zero and Y nm with n ≥ 1 and m ≥ 2, to 1. In this way we obtain: which leads to the following upper bound on Y 20 : (2.14) 6

Upper bound on Y 00
Consider the following combination of gains: In (2.15) the coefficients of the yields Y 1m and Y n1 , for any n and m, are identically zero. Thus (2.15) can be rewritten as: As usual we extract an upper bound on Y 00 by setting to their lowest value the yields whose coefficient has the same sign as Y 00 's coefficient (which correspond to the Y nm with n, m ≥ 2), and by setting to their maximum value the yields whose coefficient has opposite sign to Y 00 's coefficient (which correspond to Y 0m and Y n0 ). We know that every yield is trivially bounded by (2.3). However, in order to derive a tighter bound on Y 00 , we employ non-trivial bounds for all the yields Y nm with n + m ≤ 4 in (2.16). The upper bound on Y 00 thus satisfies: In this equation Y U ij are upper bounds and Y L ij are lower bounds. From (2.17) we obtain the following upper bound on Y 00 : where Y U 02 and Y U 20 are given in (2.10) and (2.14), respectively. The expressions for Y U 03 and Y U 04 in (2.18) can be found by starting from the same expression (2.8) that we used to derive Y U 02 , i.e.: From this expression we can extract an upper bound on any generic Y 0m as follows: where we employ the constraint (2.3). Similarly, the expressions for Y U 30 and Y U 40 are obtained starting from (2.12) and deriving an upper bound on a generic Y n0 as follows: At last, the expression for Y L 22 can be derived from the same combination of yields which led to Y U 11 . In particular, from (2.5) we have that: Then, by setting to 1 all the yields whose coefficient has equal sign to Y 22 's we obtain: 7 which yields: (2.23) Note that the upper bounds derived on Y 04 and Y 40 in this Section could be used to improve the estimation of the phase error rate given by (1.8). However, the resulting improvement in the secret key rate would be extremely small in this case and we neglect it for simplicity.

Simulations
In this Section we provide plots of the secret key rate given by (1.2) against the overall loss (−10 log 10 η) measured in dB of the two optical channels linking Alice and Bob to node C. The channel model we use to simulate the quantities that would be observed experimentally -i.e. the gains p(k c , k d |b A , b B ) and kc,k d -is given in Appendix A [33]. It accounts for: the loss in the optical channels together with the non-unity detection efficiency of D c and D d (altogether described by the parameter η), the polarization and phase misalignments introduced by the channel and a dark count probability p d in each detector. For concreteness, in all the plots below we assume fixed polarization and phase misalignments of 2%, independently of the channel loss. Note that, as pointed out in [33], the TF-QKD protocol analyzed in this work is quite robust against phase mismatch. This is so because phase misalignment only affects the quantum bit error rate but not the phase error rate. For illustration purposes every plot is obtained for three different values of the dark count rate of the detectors, p d ∈ {10 −6 , 10 −7 , 10 −8 }. The plots are obtained by numerically optimizing the secret key rate -for every value of the loss-over the signal intensity (α 2 ) and over one decoy intensity, while for simplicity the other decoy intensities are fixed to near-to-optimal values for all values of the overall loss. More specifically, we preliminarily performed an optimization of the key rate over the whole set of intensity settings and noticed that most of the decoy intensities are roughly constant with the loss and tend to be as low as possible. For instance, if we consider the case with two decoy intensity settings (µ 0 and µ 1 , with µ 0 > µ 1 ), we observe that the optimal value for the weakest decoy µ 1 is basically the lowest possible for any value of the loss. In practice, however, it might be difficult to generate very weak signals due to the finite extinction ratio of a practical intensity modulator [44], so we fix µ 1 to a reasonable small value from an experimental point of view, say µ 1 = 10 −5 [34,36], while keeping the optimization over the remaining intensities. Similarly, if we consider the case with three decoy intensity settings (µ 0 , µ 1 and µ 2 , with µ 0 > µ 1 > µ 2 ), we find that the optimal values for the weakest decoys µ 1 and µ 2 are also the lowest possible for any value of the loss. Moreover, in this last case, we show in Appendix B that the system performance remains basically unchanged if one increases the value of the weakest intensity to say µ 2 = 10 −3 , which might be even easier to implement experimentally than 10 −5 . Thus, we fix µ 2 = 10 −3 and we differentiate it from µ 1 by, for example, one order of magnitude (i.e. we take µ 1 = 10 −2 ). The same argument holds in the case with four decoy intensity settings (see Appendix B), where we fix µ 2 = 10 −3 , µ 1 = 10 −2 , and µ 0 = 10 −1 . We remark, however, that our method is general in the sense that the analytical upper bounds on the yields can be evaluated with any desired combination of intensity settings, while we select these particular decoy intensity values only for illustration purposes. Also, let us emphasize that the optimal decoy intensity values in the finite-key regime might be different from the values mentioned above. The analysis of the finite-key regime is, however, beyond the scope of this paper. Importantly, it turns out that the resulting asymptotic secret key rates in these scenarios are almost indistinguishable from those obtained by optimizing the value of all the intensity settings.
The optimal values of the signal and decoy intensities which are optimized as a function of the loss are also plotted in this Section. In this regard, we also study how the key rate is affected when the intensities are subjected to fluctuations around their optimal values in Subsec. 3.4.

Two decoy intensity settings
In figure 2 we plot the secret key rate against the overall loss for the case where Alice and Bob use two decoy intensity settings each. The solid lines are obtained by bounding from above the yields Y 00 , Y 02 , Y 20 and Y 11 by means of the expressions derived in Sec. 2 and by optimizing the rate over the signal intensity α 2 and the decoy intensity µ 0 , while the other decoy intensity is fixed to µ 1 = 10 −5 as explained above.
The optimization is carried out by using the built-in function "NMaximize" of the software Wolfram Mathematica 10.0. Loss Alice-Bob (dB)

PLOB bound
Secret key rate Secret key rate in logarithmic scale as a function of the overall loss in the channels Alice-C and Bob-C for three different dark count rates (10 −6 green, 10 −7 red, 10 −8 blue). The solid lines correspond to the case where the yields Y00, Y02, Y20 and Y11 are estimated by means of two decoy intensity settings through the bounds presented in Sec. 2 and the key rate is optimized over the signal intensity α 2 (see figure 3a) and the decoy intensity µ0 (see figure 3b). The other decoy intensity, µ1, is fixed to µ1 = 10 −5 . The dashed lines assume that all the yields are known from the channel model and the secret key rate is optimized over α 2 . That is, these lines show the maximum value of the secret key rate which could be achieved with an infinite number of decoy intensity settings and the security analysis reported in [33]. The solid magenta line illustrates the PLOB bound [19]. The plot shows that in the presence of a dark count rate of at most about p d = 10 −7 the protocol can beat the PLOB bound even with just two decoy intensity settings.  The optimal values for α 2 and µ 0 are shown in figure 3a and figure 3b, respectively. The dashed lines are instead obtained by employing the exact expression of the yields ¶ which is given by (A.6) for the channel model considered. This represents the ideal scenario in which the parties have an infinite number of decoys through which they can estimate all the yields precisely. Note that in order to obtain the dashed lines in figure 2 we use the exact expression of the yields Y nm only for n, m ≤ 12 while we set the other yields to 1. This is enough to basically reproduce the behavior of the secret key rate when all the infinite number of yields are computed via the channel model's formula given by (A.6), as argued in [33]. The dashed lines are only optimized over the signal intensity, since the yields are directly given by the channel model. Finally, we also insert in figure 2 the PLOB bound on the secret key capacity [19], which reads as follows in terms of the transmittance η : In figure 2 we observe that even by means of just two decoy intensity settings the key rate can beat the PLOB bound, provided that the dark count rate is p d 10 −7 . This happens because with two decoys the parties can already non-trivially estimate the yields Y nm with n + m ≤ 2 as we showed in Sec. 2, and these yields are the most relevant terms in the phase-error rate formula given by (1.7) [33]. Note that we did not estimate the yields Y 01 and Y 10 since only the yields Y nm with n + m an even number contribute to the phase-error rate (1.7).
However, figure 2 also shows that there is a sensible gap between the rates where the yields are estimated with two decoys (solid lines) and the best possible rates one could achieve (dashed lines) if all the yields were known. This clearly indicates that, although two decoys allow to estimate the yields of largest contribution in the phase-error rate, such estimations are not sufficiently tight and the ability to estimate a larger number of yields would increase the performance of the protocol. By considering figure 3 and the fixed value of the decoy intensity µ 1 , one notices that the optimal intensities are rather small and thus, in a real experimental implementation, intensity fluctuations might be an issue. In Subsec. 3.4 we address this problem by studying how the key rate is affected under intensity fluctuations and show that for fluctuations up to about 40% the change in the key rate performance is minimal. Also, we notice that the optimal values of the signal intensity α 2 (see figure 3a) and the decoy intensity µ 0 (see figure 3b) are almost constant with the loss, for losses 20 dB. This means that in a scenario where the loss in the quantum channels varies dynamically with time within a reasonable interval, one could still fix the signal intensity and both decoy intensities to constant values which happen to be close to the optimal ones. This argument also holds in the case of three (see Subsec. 3.2) and four decoy intensity settings (see Subsec. 3.3).

Three decoy intensity settings
In figure 4 we plot the secret key rate against the overall loss for the case where Alice and Bob use three decoy intensity settings each. The solid lines are obtained by bounding from above the relevant yields Y nm such that n + m ≤ 4 (i.e. we upper bound the yields Y 00 , Y 02 , Y 20 , Y 11 , Y 13 , Y 31 , Y 04 , Y 40 and Y 22 ). The exact expressions for the different upper bounds on the yields can be found in Appendix C, and we omit them here for simplicity. The solid lines are optimized over the signal intensity α 2 and the decoy intensity µ 0 , while the weakest decoy intensities are fixed for simplicity to µ 1 = 10 −2 and µ 2 = 10 −3 . As explained above, the resulting secret key rate in this scenario is almost indistinguishable from that obtained by optimizing over all the intensity settings. The optimal values for α 2 and µ 0 are shown in figure 5a and figure 5b, respectively. The dashed lines are again obtained by employing the exact expression of the yields given by the channel model (A.6) and coincide with those plotted in figure 2.
We observe in figure 4 that the use of three decoys yields a significant improvement in the protocol's performance with respect to the two-decoys case (see figure 2). As a matter of fact, in figure 4 the solid lines are almost overlapping the dashed lines for most values of the channel loss. This is due to the fact that with three decoys the parties constrain the yields with nine independent equations (instead of four equations as in the two-decoys case), which enable a tighter estimation of Y 00 , Y 02 , Y 20 and Y 11 and the non-trivial estimation of five additional yields. Moreover, in the case of three decoys the optimal signal intensity α 2 (see figure 5a) is roughly double the value of the correspondent intensity when using two decoys (see figure 3a). The reason for this is connected to the role of α 2 in the protocol's key rate. In fact, the prefactor p(k c , k d ) with k c + k d = 1 of PLOB bound Loss Alice-Bob (dB) Secret key rate Fig. 4: Secret key rate in logarithmic scale as a function of the overall loss in the channels Alice-C and Bob-C for three different dark count rates (10 −6 green, 10 −7 red, 10 −8 blue). The solid lines correspond to the case where the yields Y00, Y02, Y20, Y11, Y13, Y31, Y04, Y40 and Y22 are estimated by means of three decoy intensity settings through the bounds presented in Appendix C and the key rate is optimized over the signal intensity α 2 (see figure 5a) and the decoy intensity µ0 (see figure 5b). The other decoy intensities are fixed to µ1 = 10 −2 and µ2 = 10 −3 . The dashed lines assume that all the yields are known from the channel model and the secret key rate is optimized over α 2 . That is, these lines show the maximum value of the secret key rate which could be achieved with an infinite number of decoy intensity settings and the security analysis reported in [33]. The solid magenta line illustrates the PLOB bound [19]. The plot shows that already with three decoy intensity settings the key rate (solid lines) is sensibly close to the ideal one in which all the yields are known (dashed lines), meaning that the contribution of the other yields trivially bounded by 1 in the phase error rate is minimal.
the key rate formula given by (1.3) increases for increasing α 2 : the higher the mean number of photons sent by the parties (within certain limits) the higher the probability of having a click in one of the two detectors. On the other hand, increasing α 2 excessively also affects the phase-error rate. Note that by setting some yields to 1 in the phase error rate formula given by (1.7) we give rise to addends like c 2n c 2m and c 2n+1 c 2m+1 which increase for increasing α 2 , leading to an overall increase of the phase-error rate and thus decrease of the key rate. The optimal value of α 2 is thus given by the trade-off between the effect of the prefactor p(k c , k d ) and that of the terms c 2n c 2m and c 2n+1 c 2m+1 . Now, by noting that the contribution of the therms c 2n c 2m and c 2n+1 c 2m+1 decreases for increasing n, m, we understand that their negative effect on the key rate is diminished in the case of three decoys since we non-trivially estimate more yields, i.e. a lower number of yields is set to 1. This allows α 2 to acquire higher values with respect to the two-decoys case, as we observed in figure 5a. Finally we point out that such an argument does not apply to the discussion about the optimal value of the decoy intensity µ 0 in the case of two and three decoys. As a matter of fact, the key rate does not depend on the decoy intensities in the same way as on the signal intensity: the decoy intensities only appear in the yields' bounds inserted in the phase-error rate. Additionally, the analytical bounds on the yields when using two or three decoys cannot be compared in a straightforward way. Nonetheless we observe a similar behavior of the optimal µ 0 for two (see figure 3b) and three decoys (see figure 5b).

Four decoy intensity settings
In figure 6 we plot the secret key rate against the overall loss for the case where Alice and Bob use four decoy intensity settings each. Like in the three-decoys case, the solid lines are obtained by bounding from above the yields Y 00 , Y 02 , Y 20 , Y 11 , Y 13 , Y 31 , Y 04 , Y 40 and Y 22 by means of four decoys. In particular, for the  yields Y 00 , Y 02 , Y 20 , Y 11 and Y 22 we use the exact same analytical bounds derived with three decoys since they are tight enough, and the use of a fourth decoy intensity would just make them more cumbersome without providing a significant improvement of the resulting secret key rate. For the remaining four yields we instead derived tighter bounds with the help of the fourth intensity µ 3 (see Appendix D). The solid lines are obtained by optimizing the rate over the signal intensity α 2 and the fourth decoy intensity µ 3 . It turns out that the optimal values for the other decoy intensities are basically the lowest possible for any value of the loss, so, as explained above, for simplicity we fix the smallest one to an experimentally reasonable small value (say µ 2 = 10 −3 ), and then we differentiate it from the other two decoys, µ 1 and µ 0 , by one order of magnitude, i.e., we take µ 1 = 10 −2 and µ 0 = 10 −1 . Importantly, this decision has a neglectable effect on the resulting secret key rate, when compared to that obtained by optimizing over all intensity settings. The optimal values for α 2 and µ 3 are shown in figure 7a and figure 7b, respectively. The dashed lines are the same as in figure 2 and figure 4. With four decoys (see figure 6) the key rates basically reproduces the ideal ones (dashed lines) in which all the yields are known, with the gap being at maximum of 1 dB at the very end of the plot lines (i.e. in the very high loss regime). This demonstrates that there is no need to bound further yields than the nine yields we bounded in the cases of three and four decoys. Of course, the tighter estimation of the yields Y 13 , Y 31 , Y 04 and Y 40 achieved with four decoys results in an improvement of the key rate with respect to the case of three decoys (see figure 4), especially in the region of high losses. Concerning the optimal signal intensity (see figure 7a), we notice a slight increase with respect to the three-decoys case (see figure 5a) due to the tighter estimation of some yields in the phase-error rate formula, which allows their correspondent coefficients to acquire a slightly higher value under an increase of α 2 . Finally, the reason why the optimal µ 3 plot (see figure 7b) looks quite different (with values above 1) from the optimal µ 0 plots for the cases of two and three decoys (see figure 3b and figure 5b) is the following.
In the TF-QKD protocol considered, the most important yields (i.e., those with a bigger impact on the resulting phase error rate) are those associated to pairs of pulses with zero or with a very low number of photons. It is therefore very important to be able to estimate these yields as tightly as possible. For this, we have that the optimal intensities µ 0 and µ 1 (µ 0 , µ 1 and µ 2 ) for the case with two (three) decoys are well below 1, just like in standard decoy-state QKD protocols [39,40]. However, as explained above, here we use the intensity µ 3 to improve the upper bounds for the yields Y 13 , Y 31 , Y 04 and Y 40 . That is, the intensity µ 3 is only used to estimate yields associated to pairs of pulses with a total number of photons equal to four. Thus, it is natural that the optimal value of µ 3 is not too low and greater than 1. Loss Alice-Bob (dB)

PLOB bound
Secret key rate Fig. 6: Secret key rate in logarithmic scale as a function of the overall loss in the channels Alice-C and Bob-C for three different dark count rates (10 −6 green, 10 −7 red, 10 −8 blue). The solid lines correspond to the case where the yields Y00, Y02, Y20, Y11 and Y22 are estimated by means of three decoys through the bounds presented in Appendix C (i.e., for simplicity here we disregard the information provided by the additional fourth decoy intensity setting) and Y13, Y31, Y04 and Y40 are estimated with four decoys via the bounds in Appendix D.
The key rate is optimized over the signal intensity α 2 (see figure 7a) and the decoy intensity µ3 (see figure 7b), while the other decoy intensities are fixed to µ0 = 10 −1 , µ1 = 10 −2 and µ2 = 10 −3 . The dashed lines are optimized over α 2 and assume that all the yields are known from the channel model. They correspond to the maximum value of the secret key rate which could be achieved with an infinite number of decoy intensity settings. The solid magenta line illustrates the PLOB bound [19]. The plot indicates that the tighter estimation of the yields Y13, Y31, Y04 and Y40 with respect to the case of three decoy intensity settings is enough to basically reproduce the ideal scenario in which all the yields are known (dashed lines).

Intensity fluctuations
Here we investigate the robustness of the TF-QKD protocol against intensity fluctuations that may occur in the preparation of the pulses sent by Alice and Bob. This is motivated by the fact that the optimal signal and decoy intensities that the parties should adopt in order to maximize the key rate for a given loss are quite small, thus the effect of intensity fluctuations might be an issue in practice. On the other hand, we also note that the optimal value of a given decoy or signal intensity is either constant or varies very moderately with the loss.
Here we consider the simple scenario in which the intensity fluctuations are symmetric, i.e., we assume that the intensity of Alice's signal matches perfectly with the intensity of Bob's signal. Or, to put it in other words, we consider that Alice's and Bob's signals suffer from the same intensity fluctuations and thus their intensities are equal. This means that such analysis is only valid to evaluate auto-compensating TF-QKD set-ups like, for instance, the one introduced in [36]. It cannot be used however to analyze set-ups where more than one laser source is used [34,35]. Although we do not expect a dramatic change of our results when asymmetric intensity fluctuations are considered in the latter case, specially if they are not too large. Also, we assume that the signal and all the decoy intensities suffer from a fluctuation of magnitude 30%, 40% or 50% around their optimal value. This means for example that, for a fluctuation say of 30%, the signal intensity α 2 and all the decoy intensities µ k fluctuate in the intervals: 0.7 α 2 opt ≤ α 2 ≤ 1.3 α 2 opt and 0.7 µ opt k ≤ µ k ≤ 1.3 µ opt k , respectively, where α 2 opt and µ opt k represent the optimal values. We then account for the worst-case scenario by numerically minimizing the key rate over all the intensities constrained in their respective fluctuation interval. Only in this way we can still guarantee that the resulting key rate is associated to a secure protocol.   3) with the secret key rates affected by increasing intensity fluctuations (solid lines): 30%,40% and 50% (brighter colors; right to left). We assume that the fluctuations affect each decoy intensity and the signal intensity as well. The plots show that the TF-QKD protocol is quite robust against intensity fluctuations, and that its robustness increases with the number of decoys.
this lies in the fact that the tightness of the yields' bounds has a stronger dependence on the value of the decoy intensities when the number of decoys -and thus constraints on the yields-is low. In other words, if the parties have at their disposal a larger number of decoys, they can properly combine the numerous constraints on the yields and obtain inherently tight bounds, i.e. bounds that are tight regardless of the actual values of the intensities involved. If, instead, the parties have few decoys, say two, then the bounds they derive on the yields can be tight or loose depending on the values assigned to the decoy intensities, since the constraints on the yields are fewer.
In conclusion, in the case of two decoys the parties can tolerate intensity fluctuations up to 40%, which correspond to a decrease in the protocol's key rate especially in the high-loss region, quantified by a reduction of about 5 to 6 dB of the maximum tolerated loss * . Remarkably, with three decoys the decrease of the maximum tolerated loss would be under 5 dB for fluctuations up to 50%. Finally, for four decoys the protocol's performance remains almost the same for fluctuations up to about 50% around * By "maximum tolerated loss" we mean the loss threshold above which the protocol's key rate becomes roughly zero.
the optimal values (except when the dark count probability is the smallest considered: p d = 10 −8 ).
We deduce that the TF-QKD protocol introduced in [33] seems to be quite robust against intensity fluctuations.

Conclusions
In this paper we have investigated in detail the performance of the Twin-Field quantum key distribution (TF-QKD) protocol presented in [33] in the realistic scenario of a finite number of decoy intensity settings at the parties' disposal. Indeed, the protocol requires that Alice and Bob use the decoy-state method [38][39][40] to estimate the phase-error rate by upper bounding certain yields. Unlike most QKD protocols which employ such method, in this case the protocol's key rate depends -in principle-on infinitely many yields and it is essential to upper bound (rather than lower bound) their values. Clearly, the more yields the parties tightly upper bound, the better the protocol's performance is. We have introduced an analytical method to perform such estimation when Alice and Bob use two, three or four decoy intensity settings each. The yields' analytical bounds provided in this work imply a fully-analytical expression for the protocol's secret key rate, which is very convenient for performance optimization (e.g. in the finite-key scenario). Also, we remark that the secret key rates obtained with our analytical bounds basically overlap those achievable with numerical tools like linear programming for most values of the overall loss, which confirms that the analytical approach is actually quite tight.
In so doing, we have shown that the TF-QKD protocol can beat the PLOB bound [19] even with just two decoys for reasonable values of the setup parameters, which include: the loss, the dark count rate, the polarization misalignment and the phase mismatch. Furthermore the plots assuming four decoys demonstrate that one can approximately achieve the best possible performance by tightly estimating only nine yields. The optimization of the key rate over the signal and decoy intensities indicates that their optimal values are all either constant or weakly-dependent on the loss of the channel. This means that the protocol is particularly suitable for contexts where the channel loss varies in time, for instance in the scalable MDI-QKD networks conceived in [43]. Finally we have investigated the scenario where the intensities of the optical states prepared by Alice and Bob are affected by fluctuations and observed that the protocol seems to be very robust against such phenomena. A natural continuation of this work would take into account the finite-key effects due to the finite number of pulses sent by the parties to the central relay. This could be done by combining the results presented in this paper with the finite-keys estimation techniques used in [41].
The channel model that we employ to simulate the gains that would be observed experimentally in the X-basis (i.e. the probabilities p(k c , k d |b A , b B )) and Z-basis (i.e. the probabilities Q k,l kc,k d ) is taken from [33]. In all the expressions of this Section we assume k c + k d = 1. In particular, a beam splitter of transmittance √ η accounts for the loss in the quantum channel linking Alice (Bob) to node C and for the non-unity detection efficiency of detectors D c and D d . The polarization misalignment introduced by the channel Alice-C (Bob-C) is modeled with a unitary operation mapping the polarization input modes a † in (b † in ) to the orthogonal polarization output modes a † out and a † out⊥ (b † out and b † out⊥ ) according to: , for an angle θ A (θ B ). Moreover, the phase mismatch between Alice and Bob's signals arriving at node C is modeled by shifting the phase of Bob's signals by an angle φ = δπ, for a certain parameter δ. Finally the model considers that both detectors are affected by a dark count probability p d , which is independent of the signals received and has the same value for both detectors. With this setup, the gains in the X-basis can be written as: where γ = √ ηα 2 (with α being the amplitude of the signal states) and The gains in the Z-basis instead read: where the function I(z) = 1 2πi e (z/2)(t+1/t) t −1 dt is the modified Bessel function of first kind. In the simulations shown in Sec. 3 we compare the key rate computed with our analytical bounds on the yields with the key rate evaluated with the exact expressions of the yields, i.e. the expressions obtained directly from the channel model. According to the above channel model, the yields read: To conclude, we remark that all the quantities entering the key rate formula (1. Secret key rate (a) Optimal key rate as a function of the overall loss when the parties use three decoy intensity settings, for three different values of the dark count rate (p d ).
The solid lines are obtained by fixing the weaker decoy intensities to µ2 = 10 −3 and µ1 = 10 −2 , while the dot-dashed lines are obtained by fixing the same intensities to µ2 = 10 −5 and µ1 = 10 −4 . The dashed lines assume that all the yields are known from the channel model and the magenta line is the PLOB bound [19]. Note that the green dot-dashed lines and green solid lines (p d = 10 −6 ) are almost perfectly overlapping.  As explained in Sec. 3, the optimal key rates are basically not affected if their optimization is only performed over the signal intensity (α) and over one decoy intensity, while having the remaining weaker decoy intensities fixed to near-to-optimal values for all losses. In figure B1, we compare the optimal key rate that the parties can achieve when fixing their weaker decoy intensities to substantially different values, in the case of three (left) and four (right) decoy intensity settings. In particular, the solid lines are the same plotted in figure 4 and figure 6 for the three-and four-decoys case, respectively, i.e. they are obtained by fixing the weaker decoy intensities to µ 2 = 10 −3 and µ 1 = 10 −2 (three decoy intensity settings) and to µ 2 = 10 −3 , µ 1 = 10 −2 and µ 0 = 10 −1 (four decoy intensity settings). The dot-dashed lines, instead, are obtained by fixing the weaker intensities to values which are two orders of magnitude lower, that is µ 2 = 10 −5 and µ 1 = 10 −4 in the case of three decoy intensity settings and µ 2 = 10 −5 , µ 1 = 10 −4 and µ 0 = 10 −3 in the case of four decoy intensity settings. Clearly, the optimal key rates are basically not affected by employing relatively stronger pulses (those with µ 2 = 10 −3 as the weakest intensity) for the weaker decoy intensity settings. Such stronger pulses could be more easily implemented experimentally and, for this, have been chosen in our simulations.

Appendix C. Yields' bounds with three decoys
Here we derive analytical upper bounds on the yields appearing in (1.7), following the same lines of Sec. 2. In this case we assume that Alice and Bob can prepare their phase-randomized coherent pulses with three different intensity settings: {µ 0 , µ 1 , µ 2 }, which are the same for both parties. This choice is optimal since we assumed that the two optical channels linking the parties to the central node C have equal transmittance √ η [43].
The whole set of infinite yields is subjected to the following nine equality constraints: and to the inequality constraints given by (2.3 Consider the following combinations of gains in which all the terms Y 1m and Y n1 are removed (i.e. their coefficients are equal to zero): where the superscripts in G k,l 22 indicate which intensities are involved, while the subscripts indicate the yield that is going to be bounded. We now combine G 0, 1 22 , G 0,2 22 and G 1,2 22 with arbitrary real coefficients c 0 and c 1 and impose that the resulting expression has the yields Y 0m and Y n0 removed as well: 3) Note that the linear combination above is already the most general for our needs. As a matter of fact, for every linear combination of G 0, 1 22 , G 0,2 22 and G 1,2 22 one can always factor out the coefficient in front of G 0, 1 22 , as far as it is not zero. However, if the particular combination of gains which removes the terms Y 0m and Y n0 has a null coefficient in front of G 0, 1 22 , for symmetry reasons there would also exist another combination -that also removes the yields Y 0m and Y n0 -with a null coefficient in front of say G 0, 2 22 , and this one could be found in our case given by (C.3). For Y 0m and Y n0 to be removed in (C.3) it suffices that: A sufficient condition for this is that every coefficient of µ m i is identically zero, which happens for: Substituting (C.6) and (C.7) back into (C.3) and multiplying both sides by µ 2 , we get an expression where all the terms Y 0m , Y 1m , Y n0 and Y n1 are removed and where the term Y 22 gives the largest contribution: In order to extract a bound for Y 22 we need to recast the yields' coefficients in such a way that their sign becomes manifest. Each term of the sum in (C.8) may be recast as follows: or equivalently as: We can now rewrite factor A 22 as: ) .

(C.12)
Of course we can employ this expression also for A 22 (µ 0 , µ 1 , µ 2 , n), under the substitution m → n. We will apply this consideration from now on to similar scenarios. By substituting (C.12) into (C.10), we get the final expression for each term of the sum in (C.8): That is, the sign of Y nm 's coefficient is independent of n and m and it is the same for all terms in (C.8) (note that the product of the two sums in (C.13) is always positive). Thus a valid upper bound for Y 22 is obtained by setting all the other yields to zero in (C.8), except for Y 22 . We obtain: which implies the following expression for the upper bound on Y 22 : We remark that the bound given by (C.15) is not valid when any of the intensities µ 0 , µ 1 or µ 2 is equal to zero. As a matter of fact, in any of these cases the starting expression given by (C.8) becomes trivial. However, in most practical situations, due to the finite extinction ratio of amplitude modulators, none of the decoy intensities is actually equal to zero.

Appendix C.2. Upper bound on Y 11
Consider the following combinations of gains in which all the terms Y 0m and Y n0 are removed: We now combine G 0,1 11 , G 0,2 11 and G 1,2 11 with arbitrary real coefficients c 0 and c 1 and impose that the resulting expression has the yields Y 2m and Y n2 also removed: For Y 2m and Y n2 to be removed it suffices: , (C.19) . (C.20)

21
Substituting these terms back into (C.17) yields a combination of gains in which the terms Y 0m , Y n0 , Y 2m and Y n2 are removed: In order to get a valid upper bound for Y 11 we need to determine the signs of the coefficients of the remaining yields. We start by recasting each term of the sum in (C.21) corresponding to the Y nm , with n, m ≥ 3, as follows: The factor A 11 can be rewritten as: where the factor F (m) ≥ 0 , ∀ m ≥ 3. Substituting (C.24) back into (C.22), we recast each term of the sum in (C.21) corresponding to the Y nm , with n, m ≥ 3, as: 25) so that its sign is manifestly dependent on the factor (µ 0 − µ 2 )(µ 1 − µ 2 ). In a similar fashion, one can rewrite each term of the sum in (C.21) corresponding to the Y 1m , with m ≥ 3, as: 26) thus deducing that this expression has opposite sign with respect to that given by (C.25). Same holds for Y n1 , since it can be shown that its coefficient is exactly (C.26) with the substitution m → n.
Finally, by showing that the term corresponding to Y 11 in (C.21) can be factorized as: one concludes that this expression has the same sign as that given by (C.25).
Putting together these considerations into (C.21), a valid upper bound on Y 11 is obtained when the yields Y nm , with n, m ≥ 3, are set to zero and the yields Y 1m and Y n1 are set to their maximum allowed value.
Since in Appendix C.5 and Appendix C.6 we derive upper bounds on Y 13 and Y 31 (see C.65 and C.73), we can employ them in (C.21) instead of trivially bounding these yields with 1. In this way we obtain: which leads to the following upper bound on Y 11 : where the term E 11 is defined as: Consider the following combinations of gains in which all the terms Y 1m and , Y n0 are removed: We now combine G 0,1 02 , G 0,2 02 and G 1,2 02 with arbitrary real coefficients c 0 and c 1 and impose that the resulting expression has the yields Y 2m and Y n1 also removed: For Y 2m and Y n1 to be removed the coefficients c 0 and c 1 must satisfy: 33) or equivalently: A sufficient condition for this is that the coefficient of every µ n i and every µ m i is identically zero. This imposes six conditions on c 0 and c 1 , however thanks to the inherent symmetries of the system a solution exists, and reads: Substituting these expressions back into (C.32) and multiplying both sides by µ 2 , yields a combination of gains in which the terms Y n0 , Y n1 , Y 1m and Y 2m are removed. In particular, we obtain: In order to get a valid upper bound for Y 02 and Y 04 we need to study the sign of the coefficients of the remaining yields. We start by recasting each term of the sum corresponding to the Y nm , with n ≥ 3 and m ≥ 2, in (C.37) as follows: and A 22 is the one found when bounding Y 22 , thus we know from (C.12) it can be recast as: ) .
(C.40) 24 We can rewrite B 02 as: Employing (C.40) and (C.41) into (C.38) we get: which means that the sign of this expression is fully determined by the factor (µ 1 − µ 2 )(µ 0 − µ 2 ) (note that the product of the two sums in (C.42) is always positive).
Concerning the terms that appear in the sum in (C.37) corresponding to the Y 0m , with m ≥ 2, we have: where we used (C.11) in the first equality and (C.40) in the second equality. Expression (C.43) implies that its sign is always equal to the sign of the terms given by (C.42), since it is determined by the same factor (µ 1 − µ 2 )(µ 0 − µ 2 ) (note that the product of the last two factors in (C.43) is always positive).
A valid upper bound on Y 02 is thus obtained by setting all the other yields to zero in (C.37). By doing so, we obtain: One can do the same when bounding Y 04 , i.e. setting all the other yields to zero except for Y 04 , in (C.37). We find that: .
Consider the following combinations of gains in which all the terms Y 0m and Y n1 are removed: We now combine G 0,1 20 , G 0,2 20 and G 1,2 20 with arbitrary real coefficients c 0 and c 1 and impose that the resulting expression has the yields Y 1m and Y n2 also removed: For Y 1m and Y n2 to be removed the coefficients c 0 and c 1 must satisfy: (C.48) This system of linear equations coincides with the one given by (C.33) that we found when bounding Y 02 , thus the solution is given by (C.35) for c 0 and by (C.36) for c 1 . Substituting these expressions back into (C.47) and multiplying both sides by µ 2 , yields a combination of gains in which the terms Y n1 , Y n2 , Y 0m and Y 1m are removed: Since the coefficients of Y n0 and Y nm coincide with those found when bounding Y 02 if one exchanges m ←→ n, we can directly use the results obtained in Appendix C.3 to recast the terms that contain the Y nm with n ≥ 2 and m ≥ 3. In particular, according to (C.42), we obtain: ) , (C.50) and according to (C.43) the terms that contain the yields Y n0 can be written as: Like in the case of Y 02 (see Appendix C.3), a valid upper bound on Y 20 is thus obtained setting all the other yields to zero in (C.49). We obtain: One can do the same to bound Y 40 , i.e. to set all the other yields to zero, except for Y 40 . In this case we obtain: .
(C.53) removed the terms Y 0m and Y 2m as we did here. Therefore, by exploiting the result given by (C.12) we can recast each term of the sum corresponding to the Y 1m , with m ≥ 2, in (C.60) as: and realize that it is always negative, regardless of the value of the intensities. By employing the results (C.12, C.24) we can recast each term of the sum corresponding to the Y nm with n ≥ 3 and m ≥ 2, in (C.60) as: and realize that it is always positive , regardless of the intensities.
A valid upper bound on Y 13 is then obtained by setting Y 1m → 0 (except for Y 13 ) and Y nm → 1 for all n ≥ 3 and m ≥ 2 in (C.60). As a result we obtain: which implies: We thus obtain the following upper bound on Y 13 : where G 13 is defined in (C.54) and the coefficients of the combination of gains in (C.59).
Appendix C.6. Upper bound on Y 31 We look for that combination of gains in which all the terms proportional to Y n0 , Y n2 , Y 0m and Y 1m are removed. In order to find it, we proceed like in the previous case. That is, we consider the most general combination of all gains: and impose proper conditions on the real coefficients c i,j : Y n0 removed: µ i c i,j = 0 ∀ m ⇐ µ 0 c 0,j + µ 1 c 1,j + µ 2 c 2,j = 0 for j = 0, 1, 2 . (C.70) The conditions (C.67, C.68, C.69, C.70) form an overdetermined system of equations for the nine variables c i,j . However, thanks to the symmetries of the problem, a unique solution for c i,j exists and reads (we rescale every coefficient by requiring c 0,0 = 1): . (C.71) By substituting (C.71) back into (C.66) we get an expression in which the terms Y n0 , Y n2 , Y 0m and Y 1m are removed: Y nm n! m! A 22 (µ 0 , µ 1 , µ 2 , n) · A 11 (µ 0 , µ 1 , µ 2 , m) (µ 1 − µ 2 ) 2 (µ 1 + µ 2 ) , (C.72) where A 22 and A 11 are again the factors from Y 22 and Y 11 bounds given by equations (C.11,C.23), similarly to what happens when bounding Y 13 (see Appendix C.5). Therefore the analysis of the coefficients' sign is the same as in Appendix C.5. Hence a valid upper bound on Y 31 is obtained by setting Y n1 → 0 (except for Y 31 ) and Y nm → 1 in (C.72) for all n ≥ 2 and m ≥ 3 in (C.72). Analogous steps to those in Appendix C.5 lead to the following upper bound: where G 31 is defined in (C.66) and the coefficients of the combination of gains in (C.71).