Practical decoy-state method for twin-field quantum key distribution
Federico Grasselli Marcos Curty

	 
	 
	 
	 
	 
	




^1Institut für Theoretische Physik III, Heinrich-Heine-Universität Düsseldorf, Universitätsstraße 1, D-40225 Düsseldorf, Germany 
 ^2Escuela de Ingeniería de Telecomunicación, Dept. of Signal Theory and Communications, University of Vigo, E-36310 Vigo, Spain
federico.grasselli@hhu.de 

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 <cit.>, with quantum key distribution (QKD) <cit.> being its most developed application. Recent experiments over about 400 km of optical fibers <cit.> and over about 1000 km of satellite-to-ground links <cit.> 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 <cit.>.
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 <cit.> yields a polynomial scaling of the communication efficiency with the distance <cit.>. 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 <cit.> or quantum error correction <cit.>. 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 <cit.> (MDI-QKD) with quantum memories <cit.> and adaptive MDI-QKD featuring quantum non-demolition measurements <cit.>. 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 <cit.>. More recently, <cit.> 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 <cit.> 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
<cit.> as well as to investigate their experimental feasibility <cit.>. Among these protocols, the one that seems to deliver the higher secret ket rate <cit.> is that introduced in <cit.>. 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 <cit.>. The key-rate simulations provided in <cit.> 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 <cit.>, 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 <cit.>. 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 <cit.> seems to be actually quite robust against such fluctuations.


Like in <cit.>, for simplicity, here we focus on the asymptotic-key rate scenario. However, we remark that by using the techniques reported in <cit.>, 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 <ref> we present the TF protocol from <cit.> and highlight the main yields that need to be bounded. In <ref> we provide the analytical bounds on the yields for the case of two decoys (the cases of three and four decoys are treated in <ref> and <ref>, respectively). In <ref> we provide simulations of the secret key rate versus the loss for a typical channel model (briefly described in <ref>), and we also evaluate the effect of intensity fluctuations. We conclude the paper in <ref>.



§ THE TF-QKD PROTOCOL
   

As discussed above, we consider the TF-QKD protocol presented in <cit.> and sketched in <ref>. 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:


  * 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=1/2π∫_0^2π dθ|β_A e^iθ⟩⟨β_A e^iθ|= ∑_n=0^∞(β_A^2)^n e^-β_A^2/n!|n⟩⟨n|

(ρ̂_β_B) whose intensity β_A^2 (β^2_B) is drawn randomly from a set 𝒮_A={β^2_i}_i (𝒮_B={β^2_j}_j) of real nonnegative numbers.

  * Both parties send their optical pulses to the untrusted node C via optical channels in a synchronized manner.

  * 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).

  * 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.

  * 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 <cit.> yields the following lower bound on the asymptotic key rate R:

    R ≥max{R_10,0 } + max{R_01,0 }  ,

where the terms R_k_c k_d, for (k_c,k_d) ∈{ (1,0), (0,1)}, are defined as:

    R_k_c k_d = p_X^2  p(k_c,k_d) [1-f  h(e_k_c k_d) -h(e^ph_k_c k_d) ]    ,

with h(x)=-x log_2 x  -(1-x) log_2 (1-x) 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:

    p(k_c,k_d)= ∑_b_A,b_B=0^1 p(b_A,b_B) p(k_c,k_d|b_A,b_B)   ,

where p(b_A,b_B) 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: p(b_A,b_B)=1/4   ∀  b_A,b_B. p(k_c,k_d|b_A,b_B) 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_k_c k_d and e^ph_k_c k_d in (<ref>) represent the bit-error rate in the X-basis and an upper bound on the phase-error rate, respectively. The former is defined as:

    e_10   = ∑_i,j=0 | i⊕ j=1^1 p(b_A=i,b_B=j)  p(k_c=1,k_d=0|b_A=i,b_B=j)/p(k_c=1,k_d=0)  ,
    
     e_01   = ∑_i=0^1 p(b_A=i,b_B=i)  p(k_c=0,k_d=1|b_A=i,b_B=i)/p(k_c=0,k_d=1)  ,

and the latter as:

    e^ph_k_c k_d = 1/p(k_c,k_d)[(∑_n,m=0^∞ c_2 n c_2 m√(Y^k_c,k_d_2n 2m))^2 + (∑_n,m=0^∞ c_2 n+1 c_2 m+1√(Y^k_c,k_d_2n+1 2m+1))^2 ]     ,

where the coefficients c_n are defined as c_n=e^-α^2/2α^n/√(n!) and the yields Y^k_c,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 (<ref>) are those Y^k_c,k_d_nm such that n+m is an even number.

The yields Y^k_c,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 <cit.>. 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^k_c,k_d_nm and thus to a higher key rate. However, a finite number of decoys only allows to derive non-trivial upper bounds[Every yield is a probability, thus it is trivially bounded by 1.] on a limited number of yields in (<ref>), 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 point-to-point private capacity bound (PLOB bound) <cit.> at high losses (see <ref>). 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) <cit.>, while the TF-QKD protocol considered here upper bounds the value of the phase-error rate (<ref>) by upper bounding several yields. In particular, we upper bound the yields Y^k_c,k_d_nm for (n,m)∈ℐ, where ℐ is a certain subset of {(n,m)|  n,m ∈ℕ_0} which depends on the number of decoys. Thanks to the derived upper bounds on the yields (which we shall denote by Y^U,  k_c,k_d_nm) we are able to estimate the phase error rate (<ref>) as follows:

    e^ph_k_c k_d≤1/p(k_c,k_d)   [(∑_(2n,2m)∈ℐ c_2 n c_2 m√(Y^U,  k_c,k_d_2n 2m) +∑_(2n,2m)∉ℐ c_2 n c_2 m)^2 . 
       . + (∑_(2n+1,2m+1)∈ℐ c_2 n+1 c_2 m+1√(Y^U,  k_c,k_d_2n+1 2m+1) +∑_(2n+1,2m+1)∉ℐ c_2 n+1 c_2 m+1)^2 ]     .





§ 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_β_A^2 (n) and P_β_B^2 (m) (see <ref>), 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 β_A^2 and β_B^2. 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^k_c,k_d_n m 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^k_c,k_d_n m by expressing the experimentally observed gains Q^β_A^2,β_B^2_k_c,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  β_A^2 and β_B^2, respectively– in terms of the yields:

    Q^β_A^2,β_B^2_k_c,k_d = ∑_n,m=0^∞ e^-β_A^2-β_B^2(β_A^2)^n(β_B^2)^m/n!m! Y^k_c,k_d_n m  .

As it is clear from (<ref>), to every distinct pair of decoy intensities (β_A^2,β_B^2) corresponds a new constraint on the set of infinite yields {Y^k_c,k_d_n m}_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 contribution[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.] in (<ref>) –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 (<ref>) 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 <ref> and <ref>, 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 β_A^2 and β_B^2 is the same for both parties <cit.> 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 (<ref>). Hence the yields are subjected to the following four equality constraints:

    Q̃^k,l≡ e^μ_k + μ_l Q^k,l =∑_n,m=0^∞Y_nm/n!m!μ_k^n μ_l^m    k,l ∈{0,1}  ,

and to the inequality constraints:

    0 ≤ Y_nm≤ 1   ∀  n,m    .

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:

    G_11   = Q̃^0,0 + Q̃^1,1 - (Q̃^0,1 +Q̃^1,0)  
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0 -μ^n_1) (μ^m_0 - μ^m_1)   .

The subscript in G_11 indicates the yield that is going to be bounded with this combination of gains. In (<ref>) the coefficients of the yields Y_0m and 
Y_n0, for any n and m, are identically zero. Thus (<ref>) can be rewritten as:

    G_11= Y_11 (μ_0 - μ_1)^2 + ∑_[n+m>2]n,m=1^∞Y_nm/n!m!(μ^n_0 -μ^n_1) (μ^m_0 - μ^m_1)
    	  .

We observe that the coefficients that multiply the yields Y_n m 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 (<ref>) holds. 
Since all the yields' coefficients carry the same sign in (<ref>) –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 (<ref>):

    Y^U_11 = G_11/(μ_0 - μ_1)^2  ,

where G_11 is defined in (<ref>).

We remark that by combining the gains as in (<ref>), 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 (<ref>) since it has the largest coefficient. All the yields' bounds presented in this work follow the same philosophy.



 §.§ Upper bound on Y_02

Consider the following combination of gains:

    G_02   = μ_1 Q̃^0,0 + μ_0 Q̃^1,1 - μ_1 Q̃^0,1 -μ_0 Q̃^1,0
       = ∑_n,m=0^∞Y_nm/n!m!(μ_1 μ^n_0 -μ_0 μ^n_1) (μ^m_0 - μ^m_1)   .

In (<ref>) the coefficients of the yields Y_n0 and Y_1m are identically zero. Thus (<ref>) can be rewritten as:

    G_02 =    -Y_01 (μ_0 - μ_1)^2 -Y_02/2(μ_0 + μ_1)(μ_0 - μ_1)^2 -∑_m=3^∞Y_0m/m!(μ_0 - μ_1)(μ_0^m - μ_1^m) 
    	
       +∑_[m=1]n=2^∞Y_nm/n!m!μ_0 μ_1 (μ^n-1_0 -μ^n-1_1) (μ^m_0 - μ^m_1)
    	  .

Like in the derivation of Y_11's bound given by (<ref>), a valid upper bound for Y_02 is obtained by considering the worst-case scenario for the remaining yields in (<ref>). 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 (<ref>), this means setting Y_01 and Y_0m to zero and Y_nm with n≥ 2 and m≥ 1, to 1 in (<ref>). In so doing, after rearranging the terms we obtain:

    Y^U_02 = 2/(μ_0 + μ_1)(μ_0 - μ_1)^2[-G_02 + (∑_m=1^∞μ_0^m/m!-μ_1^m/m!)
    	(∑_n=2^∞μ_1 μ_0^n/n!-μ_0 μ_1^n/n!) ]   ,

which leads to the following upper bound on Y_02:

    Y^U_02 = 2(e^μ_0-e^μ_1)(μ_0-μ_1 +μ_1 e^μ_0 -μ_0 e^μ_1) -2G_02/(μ_0 + μ_1)(μ_0 - μ_1)^2  .





 §.§ Upper bound on Y_20

Consider the following combination of gains:

    G_20   = μ_1 Q̃^0,0 + μ_0 Q̃^1,1 - μ_0 Q̃^0,1 -μ_1 Q̃^1,0
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0 - μ^n_1) (μ_1 μ^m_0 -μ_0 μ^m_1)    .

In (<ref>) the coefficients of the yields Y_n1 and Y_0m are identically zero. Thus (<ref>) can be rewritten as:

    G_20 =    -Y_10 (μ_0 - μ_1)^2 -Y_20/2(μ_0 + μ_1)(μ_0 - μ_1)^2 -∑_n=3^∞Y_n0/n!(μ_0 - μ_1)(μ_0^n - μ_1^n) 
    
       +∑_[m=2]n=1^∞Y_nm/n!m!μ_0 μ_1 (μ^n_0 -μ^n_1) (μ^m-1_0 - μ^m-1_1)
      .

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 (<ref>) 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:

    Y^U_20 = 2/(μ_0 + μ_1)(μ_0 - μ_1)^2[-G_20 + (∑_n=1^∞μ_0^n/n!-μ_1^n/n!)
    (∑_m=2^∞μ_1 μ_0^m/m!-μ_0 μ_1^m/m!) ]   ,

which leads to the following upper bound on Y_20:

    Y^U_20 = 2(e^μ_0-e^μ_1)(μ_0-μ_1 +μ_1 e^μ_0 -μ_0 e^μ_1) -2G_20/(μ_0 + μ_1)(μ_0 - μ_1)^2  .





 §.§ Upper bound on Y_00

Consider the following combination of gains:

    G_00   = μ^2_1 Q̃^0,0+μ^2_0 Q̃^1,1 -μ_0 μ_1 (Q̃^0,1 +Q̃^1,0)  
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0 μ_1 -μ_0 μ^n_1) (μ^m_0 μ_1 -μ_0 μ^m_1)   .

In (<ref>) the coefficients of the yields Y_1m and Y_n1, for any n and m, are identically zero. Thus (<ref>) can be rewritten as:

    G_00=    Y_00(μ_0-μ_1)^2 -μ_0 μ_1(μ_0 -μ_1)[∑_m=2^∞Y_0m/m!(μ^m-1_0-μ^m-1_1)+∑_n=2^∞Y_n0/n! 
    	(μ^n-1_0-μ^n-1_1) ] 
       + μ^2_0 μ^2_1 ∑_n,m=2^∞Y_nm/n! m! (μ^n-1_0-μ^n-1_1)(μ^m-1_0-μ^m-1_1)    .

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 (<ref>). 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 (<ref>). The upper bound on Y_00 thus satisfies:

    G_00 =    Y^U_00(μ_0-μ_1)^2 -μ_0 μ_1(μ_0 -μ_1)[(μ_0 -μ_1)/2(Y^U_02+Y^U_20) +(μ^2_0 -μ^2_1)/6(Y^U_03+Y^U_30) 
    	  .  
       . +(μ^3_0 -μ^3_1)/24(Y^U_04+Y^U_40) + 2∑_n=5^∞(μ^n-1_0-μ^n-1_1)/n!] 
    	+ μ^2_0 μ^2_1 (μ_0 -μ_1)^2/4 Y^L_22  .

In this equation Y^U_ij are upper bounds and Y^L_ij are lower bounds. From (<ref>) we obtain the following upper bound on Y_00:

    Y^U_00 =     G_00/(μ_0-μ_1)^2 + μ_0 μ_1/μ_0-μ_1[(μ_0 -μ_1)/2(Y^U_02+Y^U_20) 
    	+(μ^2_0-μ^2_1)/6(Y^U_03+Y^U_30) + (μ^3_0 -μ^3_1)/24(Y^U_04+Y^U_40)] 
       + 2/μ_0 - μ_1[μ_1 (e^μ_0-1-μ_0^2/2-μ_0^3/6-μ_0^4/24) -
    	   μ_0 (e^μ_1-1-μ_1^2/2-μ_1^3/6-μ_1^4/24)]	-μ_0^2 μ_1^2/4 Y_22^L   .

where Y^U_02 and Y^U_20 are given in (<ref>) and (<ref>), respectively. The expressions for Y^U_03 and Y^U_04 in (<ref>) can be found by starting from the same expression (<ref>) that we used to derive Y^U_02, i.e.:

    G_02 =     -∑_m=1^∞Y_0m/m!(μ_0 - μ_1)(μ_0^m - μ_1^m)+∑_[m=1]n=2^∞Y_nm/n!m!μ_0 μ_1
    	(μ^n-1_0 -μ^n-1_1) (μ^m_0 - μ^m_1)    .

From this expression we can extract an upper bound on any generic Y_0m as follows:

    Y^U_0m = min{m!/(μ_0 - μ_1)(μ_0^m - μ_1^m)[-G_02 + (e^μ_0-e^μ_1)(μ_0 -μ_1 +μ_1 e^μ_0-μ_0 e^μ_1)]
    	   ,  1 }  ,

where we employ the constraint (<ref>). Similarly, the expressions for Y^U_30 and Y^U_40  are obtained starting from (<ref>) and deriving an upper bound on a generic Y_n0 as follows:

    Y^U_n0 = min{n!/(μ_0 - μ_1)(μ_0^n - μ_1^n)[-G_20 + (e^μ_0-e^μ_1)(μ_0 -μ_1 +μ_1 e^μ_0-μ_0 e^μ_1)]
      ,  1 }  .

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 (<ref>) we have that:

    G_11= ∑_n,m=1^∞Y_nm/n!m!(μ^n_0 -μ^n_1) (μ^m_0 - μ^m_1)   .

Then, by setting to 1 all the yields whose coefficient has equal sign to Y_22's we obtain:

    G_11= ∑_n,m=1^∞μ^n_0 -μ^n_1/n!μ^m_0 - μ^m_1/m! -(μ_0^2 -μ_1^2)^2/4 +(μ_0^2 -μ_1^2)^2/4 Y_22^L
    	  ,

which yields:

    Y_22^L =  max{4/(μ_0 -μ_1)^2 (μ_0 +μ_1)^2[G_11 - (e^μ_0-e^μ_1)^2] +1   ,  0}  .

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 (<ref>). 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 (<ref>) against the overall loss (-10log_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 Q^β_A^2,β_B^2_k_c,k_d– is given in <ref> <cit.>. 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 <cit.>, 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 optimization is carried out by using the built-in function “NMaximize” of the software Wolfram Mathematica 10.0.] 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 <cit.>, so we fix μ_1 to a reasonable small value from an experimental point of view, say μ_1=10^-5 <cit.>, 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 <ref> 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 <ref>), 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 <ref>.



 §.§ Two decoy intensity settings
  


In <ref> 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 <ref> 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 optimal values for α^2 and μ_0 are shown in <ref> and <ref>, respectively. The dashed lines are instead obtained by employing the exact expression of the yields[By “exact expression” we mean that if the experimental apparatus were accurately described by the channel model in <ref>, then the yields associated to that experimental setup would be precisely predicted by (<ref>).] which is given by (<ref>) 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 <ref> 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 (<ref>), as argued in <cit.>. 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 <ref> the PLOB bound on the secret key capacity <cit.>, which reads as follows in terms of the transmittance η :

    K(η)=-log_2 (1-η)    .

In <ref> 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 <ref>, and these yields are the most relevant terms in the phase-error rate formula given by (<ref>) <cit.>. 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 (<ref>).

However, <ref> 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 <ref> 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 <ref> 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 <ref>) and the decoy intensity μ_0 (see <ref>) are almost constant with the loss, for losses ≳ 20. 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 <ref>) and four decoy intensity settings (see <ref>).



 §.§ Three decoy intensity settings
 


In <ref> 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  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 <ref>, 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 <ref> and <ref>, respectively. The dashed lines are again obtained by employing the exact expression of the yields given by the channel model (<ref>) and coincide with those plotted in <ref>.

We observe in <ref> that the use of three decoys yields a significant improvement in the protocol's performance with respect to the two-decoys case (see <ref>). As a matter of fact, in <ref> 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 <ref>) is roughly double the value of the correspondent intensity when using two decoys (see <ref>). 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 the key rate formula given by (<ref>) 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 (<ref>) 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 <ref>.

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 <ref>) and three decoys (see <ref>). 



 §.§ Four decoy intensity settings
 


In <ref> 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 <ref>). 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 <ref> and <ref>, respectively. The dashed lines are the same as in <ref> and <ref>.

With four decoys (see <ref>) 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 <ref>), especially in the region of high losses.

Concerning the optimal signal intensity (see <ref>), we notice a slight increase with respect to the three-decoys case (see <ref>) 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 <ref>) looks quite different (with values above 1) from the optimal μ_0 plots for the cases of two and three decoys (see <ref> and <ref>) 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 <cit.>. 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. 



 §.§ 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 <cit.>. It cannot be used however to analyze set-ups where more than one laser source is used <cit.>. 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 α_opt^2 ≤α^2 ≤ 1.3 α_opt^2 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.

The results of this study are given in <ref>. Here we plot the original key rates –i.e. without fluctuations of the signal and decoy intensities– as dashed lines[The dashed lines of the key rates without fluctuations correspond to the solid lines in <ref>, <ref> and <ref>.] and the key rates affected by intensity fluctuations as solid lines. The plots are given for the same dark count rates and misalignments used in <ref>, in the case of two, three and four decoy intensity settings. The color of the solid lines becomes brighter for increasing fluctuation magnitude.

We observe that the performance of the protocol is considerably affected by intensity fluctuations in the case of two decoys, while the effect becomes almost negligible for three and four decoys. The reason for 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[By “maximum tolerated loss” we mean the loss threshold above which the protocol's key rate becomes roughly zero.]. 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 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 <cit.> 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 <cit.> 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 <cit.> 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 <cit.> 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 <cit.>. 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 <cit.>.



We thank Dagmar Bruß and Hermann Kampermann for helpful discussions, and an anonymous referee for very useful comments.
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 675662, and from the Spanish Ministry of Economy and Competitiveness (MINECO), the Fondo Europeo de Desarrollo Regional (FEDER) through grant TEC2017-88243-R.




§ REFERENCES

99
 quantum-communication-review N. Gisin and R. Thew. Nat. Photon. 1, 165–171 (2007).
 Kimble H. J. Kimble. Nature 453, 1023 (2008).
 BB84 C. H. Bennett and G. Brassard. Proc. IEEE Int. Conf. on Computers, Systems and Signal Processing, pp 175–9 (1984).
 E91 A. K. Ekert. Phys. Rev. Lett. 67, 661 (1991).
 Scarani-review V. Scarani, H. Pasquinucci, N. J. Cerf, M. Dušek, N. Lütkenhaus, and M. Peev. Rev. Mod. Phys. 81, 1301 (2009).
 Curty-MDIQKD H.-K. Lo, M. Curty, and B. Qi. Phys. Rev. Lett. 108, 130503 (2012).
 DIQKD1 U. Vazirani and T. Vidick. Phys. Rev. Lett. 113, 140501 (2014).
 Curty-review H.-K. Lo, M. Curty, and K. Tamaki. Nat. Photon. 8, 595–604 (2014).
 Diamanti-review E. Diamanti, H.-K. Lo, B. Qi and Z. Yuan. npj Quantum Inf. 2, 16025 (2016).
 Epping M. Epping, H. Kampermann, C. Macchiavello, and D. Bruß. New J. Phys. 19, 093012 (2017).
 DIQKD2 R. A. Friedman, F. Dupuis, O. Fawzi, R. Renner, and T. Vidick.  Nat. Commun. 9, 459 (2018).
 Ribeiro J. Ribeiro, G. Murta, and S. Wehner. Phys. Rev. A 97, 022307 (2018).
 Grasselli F. Grasselli, H. Kampermann and D. Bruß. New J. Phys. 20, 113014 (2018).
 404km H.-L. Yin et al.. Phys. Rev. Lett. 117, 190501 (2016).
 421km A. Boaron et al. Phys. Rev. Lett. 121, 190502 (2018).
 free-space-QKD1 S.-K. Liao et al. Nature 549, 43 (2017).
 free-space-QKD2 H. Takenaka et al. Nat. Photon. 11, 502 (2017). 
 Takeoka M. Takeoka, S. Guha, and M. M. Wilde. Nat. Comm. 5, 5235 (2014).
 PLOB S. Pirandola, R. Laurenza, C. Ottaviani, and L. Banchi. Nat. Comm. 8, 15043 (2017).
 review-quantum-repeaters N. Sangouard, C. Simon, N. de Riedmatten, and N. Gisin. Rev. Mod. Phys. 83, 33–80 (2011).
 DLCZ L.-M. Duan, M. D. Lukin, J. I. Cirac, and P. Zoller. Nature 414, 413–418 (2001).
 Grudka-quantumrep A. Grudka, et al.. Phys. Rev. A 90, 062311 (2014).
 No-quantum-memories W. J. Munro, A. M. Stephens, S. J. Devitt, K. A. Harrison, and K. Nemoto. Nat. Photon. 6, 777–781 (2012).
 all-photonic-repeaters K. Azuma, K. Tamaki, and H.-K. Lo. Nat. Comm. 6, 6787 (2015).
 Abruzzo-MDIQKD S. Abruzzo, H. Kampermann, and D. Bruß. Phys. Rev. A 89, 012301 (2014).
 mem-assisited-MDIQKD C. Panayi, M. Razavi, X. Ma, and N. Lütkenhaus. New J. Phys. 16, 043005 (2014).
 Azuma-intercityQKD K. Azuma, K. Tamaki, and W. J. Munro. Nat. Comm. 6, 10171 (2015).
 Lucamarini-TF M. Lucamarini, Z. L. Yuan, J. F. Dynes, and A. J. Shields. Nature 557, 400 (2018).
 Tamaki-security-TF K. Tamaki, H.-K. Lo, W. Wang, and M. Lucamarini. preprint arXiv:1805.05511.
 Ma-security-TF X. Ma, P. Zeng, and H. Zhou. Phys. Rev. X 8, 031043 (2018).
 Cui-security-TF C. Cui et al. preprint arXiv:1807.02334.
 Lutkenhaus-security-TF J. Lin and N. Lütkenhaus, Phys. Rev. A 98, 042332 (2018).
 Curty-security-TF M. Curty, K. Azuma, and H.-K. Lo. preprint arXiv:1807.07667.
 experiment-Toshiba M. Minder, M. Pittaluga, G. L. Roberts, M. Lucamarini, J. F. Dynes, Z. L. Yuan, and A. J. Shields. Nat. Photon. 13, 334-338 (2019).
 experiment-chinese Y. Liu et al. preprint arXiv:1902.06268. 
 experiment-Toronto X. Zhong, J. Hu, M. Curty, L. Qian and H.-K. Lo. preprint arXiv:1902.10209.
 comparison-TF M. Lucamarini. 8th International Conference on Quantum Cryptography (QCrypt’2018), Shanghai (China). Available at http://2018.qcrypt.net.
 decoy-state-method-Hwang W.-Y. Hwang. Phys. Rev. Lett. 91, 057901 (2003).
 decoy-state-method-Lo H.-K. Lo, X. Ma, and K. Chen. Phys. Rev. Lett. 94, 230504 (2005).
 decoy-state-method-Wang X.-B. Wang. Phys. Rev. Lett. 94, 230503 (2005).
 Curty-finitekey M. Curty, F. Xu, W. Cui, C. Lim, K. Tamaki, and H.-K. Lo. Nature Comm. 5, 3732 (2014).
 decoy-state-method-Ma X. Ma, B. Qi, Y. Zhao, and H.-K. Lo. Phys. Rev. A 72, 012326 (2005).
 asymmetric-MDI-QKD W. Wang, F. Xu, H.-K. Lo. preprint arXiv:1807.03466.
 ros D. Rosenberg  et al. New J. Phys. 11, 045009 (2009).
 




	

§ CHANNEL MODEL
  
	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_k_c,k_d)  is taken from <cit.>. 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: a^†_in→cosθ_A a^†_out- sinθ_A a^†_out⊥ (b^†_in→cosθ_B b^†_out- sinθ_B b^†_out⊥), 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:
	
    p(k_c,k_d|b_A,b_B)= (1-p_d)[p_d e^-2γ + q(k_c,k_d|b_A,b_B)]   ,

	where γ=√(η)α^2 (with α being the amplitude of the signal states) and
	
    q(k_c,k_d|b_A,b_B)= {[ e^-γ(1-cosϕcosθ)-e^-2γ                            k_c ⊕ b_A ⊕ b_B =1; e^-γ(1+cosϕcosθ)-e^-2γ                            k_c ⊕ b_A ⊕ b_B =0 ].

	with θ=θ_A -θ_B. Starting from (<ref>), one can readily compute the probability p(k_c,k_d) and the bit-error rate e_k_c,k_d by means of equations (<ref>) and (<ref>),(<ref>), respectively:
	
    p(k_c,k_d)    = 1/2(1-p_d)(e^-γcosϕcosθ+e^γcosϕcosθ)e^-γ -(1-p_d)^2 e^-2γ  , 
    
    		e_k_c,k_d   = e^-γcosϕcosθ-(1-p_d)e^-γ/e^-γcosϕcosθ+e^γcosϕcosθ -2(1-p_d)e^-γ  .

	The gains in the Z-basis instead read:
	
    Q^k,l_k_c,k_d= (1-p_d)[(p_d -1) e^-√(η)(μ_k+μ_l) + e^-√(η)(μ_k+μ_l)/2 I_0 (√(ημ_k μ_l)cosθ)]   ,

	where the function I(z)=1/2π i∮ e^(z/2)(t+1/t)t^-1dt is the modified Bessel function of first kind.

	In the simulations shown in <ref> 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:
	
    Y^k_c,k_d_nm= (1-p_d)[(p_d -1)(1-√(η))^n+m + y^k_c,k_d_nm]   ,

	where
	
    y^k_c,k_d_nm=∑_k=0^nn k∑_l=0^mm l√(η)^k+l(1-√(η))^n+m-k-l/2^k+l k!l!∑_r=0^kk r∑_p=0^ll  p∑_q=max(0,r+p-l)^min(k,r+p)k  q
       l  r+p-q (r+p)!(k+l-r-p)! cos^r+q(θ_A)cos^r+2p-q(θ_B) sin^2k-r-q (θ_A) sin^2l-r-2p+q (θ_B) 
    		  .

	To conclude, we remark that all the quantities entering the key rate formula (<ref>) –i.e. (<ref>),(<ref>) and the gains (<ref>) indirectly through the yields' bounds– are symmetric under the swap k_c ↔ k_d due to the symmetries of the channel model.

	In all the simulations shown in <ref> we fix both polarization and phase misalignments to 2%, which means that: θ_A=-θ_B=arcsin√(0.02) and δ=0.02.
	
	

§ STRONGER AND WEAKER DECOY INTENSITIES
 
	
	As explained in <ref>, 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 <ref>, 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 <ref> and <ref> 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.
	
	

§ YIELDS' BOUNDS WITH THREE DECOYS
 
	Here we derive analytical upper bounds on the yields appearing in (<ref>), following the same lines of <ref>. 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 √(η) <cit.>.

	The whole set of infinite yields is subjected to the following nine equality constraints:
	
    Q̃^k,l≡ e^μ_k + μ_l Q^k,l =∑_n,m=0^∞Y_nm/n!m!μ_k^n μ_l^m    k,l ∈{0,1,2}  ,

	and to the inequality constraints given by (<ref>).

	We derive bounds on the yields Y_00,Y_11,Y_02,Y_20,Y_22,Y_13,Y_31,Y_04 and Y_40.
	
	

 §.§ Upper bound on Y_22
 
	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):
	
    G_22^0,1   = μ^2_1 Q̃^0,0+μ^2_0 Q̃^1,1 -μ_0 μ_1 (Q̃^0,1 +Q̃^1,0) 
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0 μ_1 -μ_0 μ^n_1) (μ^m_0 μ_1 -μ_0 μ^m_1)  ;
    
    		G_22^0,2   = μ^2_2 Q̃^0,0+μ^2_0 Q̃^2,2 -μ_0 μ_2 (Q̃^0,2 +Q̃^2,0) 
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0 μ_2 -μ_0 μ^n_2) (μ^m_0 μ_2 -μ_0 μ^m_2)  ;
    
    		G_22^1,2   = μ^2_2 Q̃^1,1+μ^2_1 Q̃^2,2 -μ_1 μ_2 (Q̃^1,2 +Q̃^2,1) 
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_1 μ_2 -μ_1 μ^n_2) (μ^m_1 μ_2 -μ_1 μ^m_2)    ,

	where the superscripts in G_22^k,l indicate which intensities are involved, while the subscripts indicate the yield that is going to be bounded.

	We now combine G_22^0,1,G_22^0,2 and G_22^1,2 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:
	
    G_22^0,1+c_0  G_22^0,2+c_1  G_22^1,2   =  ∑_n,m=0^∞Y_nm/n!m![(μ^n_0 μ_1 -μ_0 μ^n_1) 
    		(μ^m_0 μ_1 -μ_0 μ^m_1) . 
       . +c_0 (μ^n_0 μ_2 -μ_0 μ^n_2) (μ^m_0 μ_2 -μ_0 μ^m_2) 
    		+c_1 (μ^n_1 μ_2 -μ_1 μ^n_2) (μ^m_1 μ_2 -μ_1 μ^m_2)]  .

	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_22^0,1,G_22^0,2 and G_22^1,2 one can always factor out the coefficient in front of G_22^0,1, 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_22^0,1, 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_22^0,2, and this one could be found in our case given by (<ref>).

	For Y_0m and Y_n0 to be removed in (<ref>) it suffices that:
	
    (μ_1-μ_0)(μ^m_0 μ_1 -μ_0 μ^m_1)+c_0 (μ_2-μ_0)(μ^m_0 μ_2 -μ_0 μ^m_2)+c_1(μ_2-μ_1)(μ^m_1 μ_2 -μ_1 μ^m_2) = 0   ∀  m ,

	which implies:
	
    μ_0^m[μ_1(μ_1-μ_0)+ c_0 μ_2(μ_2 -μ_0)] + μ_1^m[-μ_0(μ_1-μ_0)+ c_1 μ_2(μ_2 -μ_1)] 
       + μ_2^m[-c_0 μ_0(μ_2-μ_0) - c_1 μ_1 (μ_2 -μ_1)] = 0 	  ∀  m   .

	A sufficient condition for this is that every coefficient of μ_i^m is identically zero, which happens for:
	
    c_0    = -μ_1(μ_0-μ_1)/μ_2(μ_0-μ_2)  , 
    
    		c_1    =μ_0(μ_0-μ_1)/μ_2(μ_1-μ_2)  .

	Substituting (<ref>) and (<ref>) back into (<ref>) 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:
	
    μ_2 G_22^0,1 -μ_1(μ_0-μ_1)/(μ_0-μ_2) G_22^0,2 + μ_0(μ_0-μ_1)/(μ_1-μ_2) G_22^1,2 = 
    		  ∑_n,m=2^∞Y_nm/n!m![μ_2 (μ^n_0 μ_1 -μ_0 μ^n_1)(μ^m_0 μ_1 -μ_0 μ^m_1) 
       -μ_1(μ_0-μ_1)/(μ_0-μ_2)(μ^n_0 μ_2 -μ_0 μ^n_2) (μ^m_0 μ_2 -μ_0 μ^m_2)+
    		   μ_0(μ_0-μ_1)/(μ_1-μ_2)(μ^n_1 μ_2 -μ_1 μ^n_2) (μ^m_1 μ_2 -μ_1 μ^m_2)]   .

	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 (<ref>) may be recast as follows:
	
    Y_nm/n!m!μ_0 μ_1 μ_2[(μ^n-1_0 -μ^n-1_1)(μ^m_0 μ_1 -μ_0 μ^m_1) 
    		   -(μ_0-μ_1)/(μ_0-μ_2)(μ^n-1_0 -μ^n-1_2) (μ^m_0 μ_2 -μ_0 μ^m_2) 
       +(μ_0-μ_1)/(μ_1-μ_2)(μ^n-1_1 - μ^n-1_2) (μ^m_1 μ_2 -μ_1 μ^m_2) ]   ,

	or equivalently as:
	
    Y_nm/n!m!μ_0 μ_1 μ_2/(μ_0-μ_2)(μ_1-μ_2)  A_22(μ_0,μ_1,μ_2,m)· A_22(μ_0,μ_1,μ_2,n)   ,

	where
	
    A_22(μ_0,μ_1,μ_2,m)    ≡μ_1^m (μ_0 - μ_2) + μ_2^m (μ_1 - μ_0) + μ_0^m (μ_2 - μ_1)   .

	We can now rewrite factor A_22 as:
	
    A_22(μ_0,μ_1,μ_2,m) = μ_1 [μ_1^m-1(μ_0 - μ_2)-(μ_0^m - μ_2^m)]+ μ_0 μ_2 (μ_0^m-1-μ_2^m-1) 
       = μ_1 [μ_1^m-1(μ_0 - μ_2)-(μ_0 - μ_2)(∑_k=0^m-1μ_0^m-1-kμ_2^k)]+ 
    		μ_0 μ_2 (μ_0-μ_2)(∑_j=0^m-2μ_0^m-2-jμ_2^j) 
       = (μ_0 - μ_2) [μ_1^m -μ_1 ∑_k=0^m-1μ_0^m-1-kμ_2^k +μ_0 μ_2 ∑_j=0^m-2μ_0^m-2-jμ_2^j ] 
       = (μ_0 - μ_2) [μ_1^m + ∑_k=0^m-1μ_2^k (-μ_1 μ_0 ^m-1-k + μ_0μ_2μ_0^m-2-k) -μ_0μ_2μ_2^m-1/μ_0]
    		   
       = (μ_0 - μ_2) [-(μ_2^m - μ_1^m) + ∑_k=0^m-1μ_2^k μ_0^m-1-k (μ_2 - μ_1)] 
       = (μ_0 - μ_2)(μ_2 - μ_1)[∑_k=0^m-1μ_2^k μ_0^m-1-k - ∑_j=0^m-1μ_2^j μ_1^m-1-j]  
       = (μ_0 - μ_2)(μ_2 - μ_1)∑_k=0^m-1μ_2^k (μ_0^m-1-k-μ_1^m-1-k)    .

	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 (<ref>) into (<ref>), we get the final expression for each term of the sum in (<ref>):
	
    Y_nm/n!m!μ_0 μ_1 μ_2/(μ_0-μ_2)(μ_1-μ_2) (μ_0 - μ_2)^2(μ_2 - μ_1)^2 
       × [∑_k=0^m-1μ_2^k (μ_0^m-1-k-μ_1^m-1-k)] [∑_j=0^n-1μ_2^j (μ_0^n-1-j-μ_1^n-1-j)]   .

	That is, the sign of Y_nm's coefficient is independent of n and m and it is the same for all terms in (<ref>) (note that the product of the two sums in (<ref>) is always positive). Thus a valid upper bound for Y_22 is obtained by setting all the other yields to zero in (<ref>), except for Y_22. We obtain:
	
    μ_2 G_22^0,1 -μ_1(μ_0-μ_1)/(μ_0-μ_2) G_22^0,2 + μ_0(μ_0-μ_1)/(μ_1-μ_2) G_22^1,2 =
    		Y^U_22μ_0 μ_1 μ_2/4(μ_0-μ_2)(μ_1-μ_2)(μ_0 - μ_1)^2   ,

	which implies the following expression for the upper bound on Y_22:
	
    Y^U_22=4G_22^0,1/μ_0 μ_1 (μ_0 - μ_1) -G_22^0,2/μ_0 μ_2(μ_0-μ_2) +
    			 G_22^1,2/μ_1 μ_2(μ_1-μ_2)/(μ_0-μ_1)(μ_0-μ_2)(μ_1-μ_2)  .

	We remark that the bound given by (<ref>) 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 (<ref>) 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.    
	
	

 §.§ Upper bound on Y_11

	Consider the following combinations of gains in which all the terms Y_0m and Y_n0 are removed:
	
    G_11^0,1   = Q̃^0,0+ Q̃^1,1 - (Q̃^0,1 +Q̃^1,0) 
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0 - μ^n_1) (μ^m_0 -μ^m_1)  ;
    
    	G_11^0,2   = Q̃^0,0+ Q̃^2,2 - (Q̃^0,2 +Q̃^2,0) 
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0 - μ^n_2) (μ^m_0 - μ^m_2)  ; 
    
    	G_11^1,2   = Q̃^1,1+ Q̃^2,2 - (Q̃^1,2 +Q̃^2,1) 
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_1 - μ^n_2) (μ^m_1 - μ^m_2)    .

	We now combine G_11^0,1,G_11^0,2 and G_11^1,2 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:
	
    G_11^0,1   + c_0  G_11^0,2+c_1  G_11^1,2 = 
       ∑_n,m=0^∞Y_nm/n!m![(μ^n_0 - μ^n_1)(μ^m_0 - μ^m_1) +c_0 (μ^n_0 - μ^n_2)(μ^m_0 - μ^m_2)+ 
    	 c_1 (μ^n_1 - μ^n_2)(μ^m_1 - μ^m_2)]    .

	For Y_2m and Y_n2 to be removed it suffices:
	
    (μ^n_0 - μ^n_1)(μ^2_0 - μ^2_1) +c_0 (μ^n_0 - μ^n_2)(μ^2_0 - μ^2_2)+	c_1 (μ^n_1 - μ^n_2)(μ^2_1 - μ^2_2) =0   ∀  n ,

	which is fulfilled by:
	
    c_0    = -(μ_0^2-μ_1^2)/(μ^2_0-μ^2_2)  ,
    
    	c_1    =(μ^2_0-μ^2_1)/(μ^2_1-μ^2_2)  .

	Substituting these terms back into (<ref>) yields a combination of gains in which the terms Y_0m,Y_n0,Y_2m and Y_n2 are removed:
	
    G_11^0,1 -(μ_0^2-μ_1^2)/(μ^2_0-μ^2_2) G_11^0,2+(μ^2_0-μ^2_1)/(μ^2_1-μ^2_2) G_11^1,2 = 
       Y_11 (μ_0 - μ_1)[(μ_0 -μ_1) -(μ_0 +μ_1)/(μ_0+μ_2)(μ_0-μ_2) +(μ_0 +μ_1)/(μ_1+μ_2)(μ_1-μ_2) ] 
       +∑_m=3^∞Y_1m/m!(μ_0 -μ_1)[(μ^m_0 -μ^m_1) -(μ_0 +μ_1)/(μ_0+μ_2)(μ^m_0-μ^m_2) 
    		                                                           +(μ_0 +μ_1)/(μ_1+μ_2)(μ^m_1-μ^m_2) ] 
       +∑_n=3^∞Y_n1/n!(μ_0 -μ_1)[(μ^n_0 -μ^n_1) -(μ_0 +μ_1)/(μ_0+μ_2)(μ^n_0-μ^n_2) 
    	                                                              	+(μ_0+μ_1)/(μ_1+μ_2)(μ^n_1-μ^n_2)] 
       +∑_n,m=3^∞Y_nm/n!m![(μ^n_0 - μ^n_1)(μ^m_0 - μ^m_1) 
    	    -(μ_0^2-μ_1^2)/(μ^2_0-μ^2_2) (μ^n_0 - μ^n_2)(μ^m_0 - μ^m_2)+ 
    	    (μ^2_0-μ^2_1)/(μ^2_1-μ^2_2) (μ^n_1 - μ^n_2)(μ^m_1 - μ^m_2)]    .

	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 (<ref>) corresponding to the Y_nm, with n,m≥ 3, as follows:
	
    Y_nm/n!m!1/(μ^2_0-μ^2_2)(μ^2_1-μ^2_2) A_11(μ_0,μ_1,μ_2,m)·  A_11(μ_0,μ_1,μ_2,n)   ,

	where
	
    A_11(μ_0,μ_1,μ_2,m)    ≡μ_1^m (μ^2_0 - μ^2_2) + μ_2^m (μ^2_1 - μ^2_0) + μ_0^m (μ^2_2 - μ^2_1)   .

	The factor A_11 can be rewritten as:
	
    A_11(μ_0,μ_1,μ_2,m) = μ^2_1 [μ_1^m-2(μ^2_0 - μ^2_2)-(μ_0^m - μ_2^m)]+ μ^2_0 μ^2_2 (μ_0^m-2-μ_2^m-2) 
       = (μ_0 - μ_2) [μ_1^m(μ_0+μ_2)  -μ^2_1 ∑_k=0^m-1μ_0^m-1-kμ_2^k +μ^2_0 μ^2_2 ∑_j=0^m-3μ_0^m-3-jμ_2^j ] 
       = (μ_0 - μ_2) [μ_1^m(μ_0+μ_2) + ∑_k=0^m-1μ_2^k (-μ^2_1 μ_0 ^m-1-k + μ^2_0μ^2_2μ_0^m-3-k)
    	      -μ^2_0μ^2_2(μ_2^m-2/μ_0+μ_2^m-1/μ^2_0) ]  
       = (μ_0 - μ_2) [(μ_1^m-μ_2^m)(μ_0+μ_2) + ∑_k=0^m-1μ_2^k μ_0 ^m-1-k(μ_2^2 - μ_1^2)] 
       = (μ_0 - μ_2) [(μ_0+μ_2)(μ_1-μ_2)(∑_j=0^m-1μ_1^m-1-jμ_2^j) -(μ_1+μ_2)(μ_1-μ_2)
    	                      ∑_k=0^m-1μ_2^k μ_0 ^m-1-k] 
       = (μ_0 - μ_2)(μ_1-μ_2) ∑_k=0^m-1μ_2^k [(μ_0+μ_2)μ_1^m-1-k - (μ_1+μ_2)μ_0 ^m-1-k] 
       = (μ_0 - μ_2)(μ_1-μ_2){∑_k=0^m-3μ_2^k [μ_2(μ_1^m-1-k-μ_0 ^m-1-k) + μ_0μ_1(μ_1 ^m-2-k-μ_0 ^m-2-k) ] 
    		
        +μ_2^m-1(μ_0 -μ_1) + μ_2^m-1(μ_1-μ_0) }
       = (μ_0 - μ_2)(μ_1-μ_2)(μ_1-μ_0)∑_k=0^m-3μ_2^k [μ_2 ∑_j=0^m-2-kμ_1^m-2-k-jμ_0^j +μ_0μ_1 
    	      ∑_j=0^m-3-kμ_1^m-3-k-jμ_0^j]  
       = (μ_0 - μ_2)(μ_1-μ_2)(μ_1-μ_0)∑_k=0^m-3μ_2^k [(μ_2+μ_0)∑_j=0^m-3-kμ_1^m-2-k-jμ_0^j +μ_2μ_0^m-2-k] 
       ≡ (μ_0 - μ_2)(μ_1-μ_2)(μ_1-μ_0) F(m)      ,

	where the factor F(m)≥ 0  ,  ∀  m≥ 3. Substituting (<ref>) back into (<ref>), we recast each term of the sum in (<ref>) corresponding to the Y_nm, with n,m≥ 3, as:
	
    Y_nm/n!m!(μ_0 - μ_2)^2(μ_1-μ_2)^2(μ_1-μ_0)^2/(μ^2_0-μ^2_2)(μ^2_1-μ^2_2) F(n)F(m)    ,

	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 (<ref>) corresponding to the Y_1m, with m≥ 3, as:
	
    -Y_1m/m!(μ_0 - μ_1)^2(μ_1-μ_2)(μ_0-μ_2)/(μ_0+μ_2)(μ_1+μ_2) F(m)    ,

	thus deducing that this expression has opposite sign with respect to that given by (<ref>). Same holds for Y_n1, since it can be shown that its coefficient is exactly (<ref>) with the	substitution m→ n.

	Finally, by showing that the term corresponding to Y_11 in (<ref>) can be factorized as:
	
    Y_11(μ_0 - μ_1)^2(μ_1-μ_2)(μ_0-μ_2)/(μ_0+μ_2)(μ_1+μ_2)  ,

	one concludes that this expression has the same sign as that given by (<ref>).

	Putting together these considerations into (<ref>), 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 <ref> and <ref> we derive upper bounds on Y_13 and Y_31 (see <ref> and <ref>), we can employ them in (<ref>) instead of trivially bounding these yields with 1. In this way we obtain:
	
    G_11^0,1   -(μ_0^2-μ_1^2)/(μ^2_0-μ^2_2) G_11^0,2+(μ^2_0-μ^2_1)/(μ^2_1-μ^2_2) G_11^1,2 =
    		    Y^U_11(μ_0 - μ_1)^2(μ_1-μ_2)(μ_0-μ_2)/(μ_0+μ_2)(μ_1+μ_2)
       +(μ_0-μ_1)/6 (Y^U_13+Y^U_31)  [μ^3_0 -μ^3_1 -(μ_0 +μ_1)/(μ_0+μ_2)(μ^3_0-μ^3_2)+(μ_0+μ_1)/(μ_1+μ_2)(μ^3_1-μ^3_2)] 
       + 2(μ_0-μ_1) ∑_n=4^∞[(μ^n_0 -μ^n_1)/n! -(μ_0 +μ_1)/(μ_0+μ_2)(μ^n_0-μ^n_2)/n! 
    		    +(μ_0+μ_1)/(μ_1+μ_2)(μ^n_1-μ^n_2)/n!]   ,

	which leads to the following upper bound on Y_11:
	
    Y^U_11=     (μ_0+μ_2)(μ_1+μ_2)/(μ_0 - μ_1)^2(μ_1-μ_2)(μ_0-μ_2)[G_11^0,1 -(μ_0^2-μ_1^2)/(μ^2_0-μ^2_2) G_11^0,2+(μ^2_0-μ^2_1)/(μ^2_1-μ^2_2) G_11^1,2 - 2(μ_0-μ_1)E_11]  
       +(μ_1μ_2+μ_0μ_1+μ_0μ_2)/6(Y^U_13+Y^U_31)   ,
 
	where the term E_11 is defined as:
	
    E_11 =      e^μ_0-e^μ_1-(μ_0 -μ_1)(1+μ_0/2+μ_1/2+μ^2_0/6+μ^2_1/6+μ_0μ_1/6)  
       +μ_0+μ_1/μ_1+μ_2[e^μ_1-e^μ_2-(μ_1 -μ_2)(1+μ_1/2+μ_2/2+μ^2_1/6+μ^2_2/6+μ_1μ_2/6)] 
       -μ_0+μ_1/μ_0+μ_2[e^μ_0-e^μ_2-(μ_0 -μ_2)(1+μ_0/2+μ_2/2+μ^2_0/6+μ^2_2/6+μ_0μ_2/6)]   .

	
	

 §.§ Upper bound on Y_02 and Y_04
  
	Consider the following combinations of gains in which all the terms Y_1m and ,Y_n0 are removed:
	
    G_02^0,1   = μ_1Q̃^0,0+ μ_0Q̃^1,1 - μ_1Q̃^0,1 -μ_0Q̃^1,0
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0μ_1 - μ_0μ^n_1) (μ^m_0 -μ^m_1)  ;
    
    	G_02^0,2   = μ_2Q̃^0,0+ μ_0Q̃^2,2 - μ_2Q̃^0,2 -μ_0Q̃^2,0
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0μ_2 - μ_0μ^n_2) (μ^m_0 -μ^m_2)  ;
    
    	G_02^1,2   = μ_2Q̃^1,1+ μ_1Q̃^2,2 - μ_2Q̃^1,2 -μ_1Q̃^2,1
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_1μ_2 - μ_1μ^n_2) (μ^m_1 -μ^m_2)    .

	We now combine G_02^0,1,G_02^0,2 and G_02^1,2 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:
	
    G_02^0,1 + c_0  G_02^0,2+c_1  G_02^1,2 = 
       ∑_n,m=0^∞Y_nm/n!m![(μ^n_0μ_1 - μ_0μ^n_1)(μ^m_0 - μ^m_1) +c_0 (μ^n_0μ_2 - μ_0μ^n_2)(μ^m_0 - μ^m_2)+ 
    	c_1 (μ^n_1μ_2 - μ_1μ^n_2)(μ^m_1 - μ^m_2)]   .

	For Y_2m and Y_n1 to be removed the coefficients c_0 and c_1 must satisfy:
	
    {[            (μ^n_0μ_1 - μ_0μ^n_1)(μ_0 - μ_1) +c_0 (μ^n_0μ_2 - μ_0μ^n_2)(μ_0 - μ_2)+ c_1 (μ^n_1μ_2 - μ_1μ^n_2)(μ_1 - μ_2)                                                                                                                       =                                                                                                                0    ∀ n; (μ^2_0μ_1 - μ_0μ^2_1)(μ^m_0 - μ^m_1) +c_0 (μ^2_0μ_2 - μ_0μ^2_2)(μ^m_0 - μ^m_2)+c_1 (μ^2_1μ_2 - μ_1μ^2_2)(μ^m_1 - μ^m_2)                                                                                                                       =                                                                                                            0
    	     ∀  m ].

	or equivalently:
	
    2.2{[                           μ_0^n[μ_1(μ_0-μ_1)+c_0μ_2(μ_0-μ_2)]+μ_1^n[-μ_0(μ_0-μ_1)+c_1μ_2(μ_1-μ_2)];                                                 -μ_2^n[μ_0 c_0(μ_0-μ_2)+μ_1 c_1(μ_1-μ_2)] =0   ∀ n; μ_0^m[μ_1μ^2_0-μ_0μ^2_1+c_0(μ_2μ^2_0-μ_0μ^2_2)]+μ_1^m[-(μ_1μ^2_0-μ_0μ^2_1)+c_1(μ_2μ^2_1-μ_1μ^2_2)];                                  -μ_2^m[c_0(μ_2μ^2_0-μ_0μ^2_2)+c_1(μ_2μ^2_1-μ_1μ^2_2)] =0   ∀  m . ].

	A sufficient condition for this is that the coefficient of every μ_i^n and every μ_i^m 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:
	
    c_0    = -μ_1(μ_0-μ_1)/μ_2(μ_0-μ_2)  , 
    
    	c_1    = μ_0(μ_0-μ_1)/μ_2(μ_1-μ_2)  .

	Substituting these expressions back into (<ref>) 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:
	
    μ_2 G_02^0,1 -μ_1(μ_0-μ_1)/(μ_0-μ_2) G_02^0,2+μ_0(μ_0-μ_1)/(μ_1-μ_2) G_02^1,2 = 
       ∑_m=2^∞Y_0m/m!(μ_0 -μ_1)[-μ_2(μ^m_0 -μ^m_1) +μ_1(μ^m_0-μ^m_2) -μ_0 (μ^m_1-μ^m_2) ] 
       +∑_[m=2]n=3^∞Y_nm/n!m![μ_2 (μ^n_0μ_1 - μ_0μ^n_1)(μ^m_0 - μ^m_1) -μ_1(μ_0-μ_1)/(μ_0-μ_2)
    	 (μ^n_0μ_2 - μ_0μ^n_2)(μ^m_0 - μ^m_2) . 
       .+ μ_0(μ_0-μ_1)/(μ_1-μ_2)(μ^n_1μ_2 - μ_1μ^n_2)(μ^m_1 - μ^m_2) ]      .

	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 (<ref>) as follows:
	
    Y_nm/n!m!1/(μ_0-μ_2)(μ_2-μ_1) A_22(μ_0,μ_1,μ_2,m)· B_02(μ_0,μ_1,μ_2,n)   ,

	where
	
    B_02(μ_0,μ_1,μ_2,n)    ≡μ_1μ_2μ_0^n(μ_1-μ_2) + μ_0^2(μ_1μ_2^n-μ_2μ_1^n)+μ_0(μ_2^2μ_1^n-μ_1^2μ_2^n)

	and A_22 is the one found when bounding Y_22, thus we know from (<ref>) it can be recast as:
	
    A_22(μ_0,μ_1,μ_2,m)=(μ_0 - μ_2)(μ_2 - μ_1)∑_k=0^m-1μ_2^k (μ_0^m-1-k-μ_1^m-1-k)   .

	We can rewrite B_02 as:
	
    B_02(μ_0,μ_1,μ_2,n) =μ_0μ_1μ_2 (μ_1-μ_2)[μ_0^n-1 -μ_0 ∑_k=0^n-2μ_1^n-2-kμ_2^k +μ_1μ_2∑_j=0^n-3μ_1^n-3-jμ_2^j] 
       = μ_0μ_1μ_2 (μ_1-μ_2)[μ_0^n-1 + ∑_k=0^n-2μ_1^n-2-kμ_2^k(μ_2-μ_0) -μ_2^n-1] 
       = μ_0μ_1μ_2 (μ_1-μ_2)[μ_0^n-1-μ_2^n-1 - ∑_k=0^n-2μ_1^n-2-kμ_2^k(μ_0-μ_2)] 
       = μ_0μ_1μ_2 (μ_1-μ_2)(μ_0-μ_2)[∑_k=0^n-2μ_2^k μ_0^n-2-k - ∑_k=0^n-2μ_2^k μ_1^n-2-k] 
       = μ_0μ_1μ_2 (μ_1-μ_2)(μ_0-μ_2)∑_k=0^n-2μ_2^k (μ_0^n-2-k - μ_1^n-2-k)   .

	Employing (<ref>) and (<ref>) into (<ref>) we get:
	
    Y_nm/n!m!μ_0μ_1μ_2 (μ_1-μ_2)(μ_0-μ_2)[∑_k=0^n-2μ_2^k (μ_0^n-2-k - μ_1^n-2-k)]
    	  [∑_k=0^m-1μ_2^k (μ_0^m-1-k-μ_1^m-1-k)]    ,

	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 (<ref>) is always positive).

	Concerning the terms that appear in the sum in (<ref>) corresponding to the Y_0m, with m≥ 2, we have:
	
    Y_0m/m!(μ_1 -μ_0)[μ_2(μ^m_0 -μ^m_1) -μ_1(μ^m_0-μ^m_2) +μ_0 (μ^m_1-μ^m_2) ] 
       = Y_0m/m!(μ_1 -μ_0) A_22(μ_0,μ_1,μ_2,m)  
       = Y_0m/m!(μ_0 - μ_2)(μ_1 - μ_2)(μ_0 -μ_1)  ∑_k=0^m-1μ_2^k (μ_0^m-1-k-μ_1^m-1-k)   ,

	where we used (<ref>) in the first equality and (<ref>) in the second equality. Expression (<ref>) implies that its sign is always equal to the sign of the terms given by (<ref>), since it is determined by the same factor (μ_1-μ_2)(μ_0-μ_2) (note that the product of the last two factors in (<ref>) is always positive).

	A valid upper bound on Y_02 is thus obtained by setting all the other yields to zero in (<ref>). By doing so, we obtain:
	
    Y^U_02=2μ_2 G_02^0,1/μ_0-μ_1 -μ_1 G_02^0,2/μ_0-μ_2 +μ_0 G_02^1,2/μ_1-μ_2/(μ_0 - μ_2)(μ_1 - μ_2)(μ_0 -μ_1)  .

	One can do the same when bounding Y_04, i.e. setting all the other yields to zero except for Y_04, in (<ref>). We find that:
	
    Y^U_04=4!μ_2 G_02^0,1/μ_0-μ_1 -μ_1 G_02^0,2/μ_0-μ_2 +μ_0 G_02^1,2/μ_1-μ_2/μ_1(μ^4_0 - μ^4_2)-μ_0 (μ^4_1 - μ^4_2) -μ_2(μ^4_0 -μ^4_1)  .


	
	

 §.§ Upper bound on Y_20 and Y_40

	Consider the following combinations of gains in which all the terms Y_0m and Y_n1 are removed:
	
    G_20^0,1   = μ_1Q̃^0,0+ μ_0Q̃^1,1 - μ_0Q̃^0,1 -μ_1Q̃^1,0
       = ∑_n,m=0^∞Y_nm/n!m! (μ^n_0 -μ^n_1)(μ^m_0μ_1 - μ_0μ^m_1)   ;
    
    	G_20^0,2   = μ_2Q̃^0,0+ μ_0Q̃^2,2 - μ_0Q̃^0,2 -μ_2Q̃^2,0
       = ∑_n,m=0^∞Y_nm/n!m! (μ^n_0 -μ^n_2)(μ^m_0μ_2 - μ_0μ^m_2)  ;
    
    	G_20^1,2   = μ_2Q̃^1,1+ μ_1Q̃^2,2 - μ_1Q̃^1,2 -μ_2Q̃^2,1
       = ∑_n,m=0^∞Y_nm/n!m! (μ^n_1 -μ^n_2)(μ^m_1μ_2 - μ_1μ^m_2)    .

	We now combine G_20^0,1,G_20^0,2 and G_20^1,2 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:
	
    G_20^0,1 + c_0  G_20^0,2+c_1  G_20^1,2 = 
       ∑_n,m=0^∞Y_nm/n!m![(μ^n_0 -μ^n_1)(μ^m_0μ_1 - μ_0μ^m_1) +c_0 (μ^n_0 -μ^n_2)(μ^m_0μ_2 - μ_0μ^m_2)+ 
    	c_1 (μ^n_1 -μ^n_2)(μ^m_1μ_2 - μ_1μ^m_2)]     .

	For Y_1m and Y_n2 to be removed the coefficients c_0 and c_1 must satisfy:
	
    {[            (μ^m_0μ_1 - μ_0μ^m_1)(μ_0 - μ_1) +c_0 (μ^m_0μ_2 - μ_0μ^m_2)(μ_0 - μ_2)+ c_1 (μ^m_1μ_2 - μ_1μ^m_2)(μ_1 - μ_2)                                                                                                                       =                                                                                                                0    ∀ m; (μ^2_0μ_1 - μ_0μ^2_1)(μ^n_0 - μ^n_1) +c_0 (μ^2_0μ_2 - μ_0μ^2_2)(μ^n_0 - μ^n_2)+c_1 (μ^2_1μ_2 - μ_1μ^2_2)(μ^n_1 - μ^n_2)                                                                                                                       =                                                                                                          0
    		  ∀  n   . ].

	This system of linear equations coincides with the one given by (<ref>) that we found when bounding Y_02, thus the solution is given by (<ref>) for c_0 and by (<ref>) for c_1.
	Substituting these expressions back into (<ref>) 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:
	
    μ_2 G_20^0,1 -μ_1(μ_0-μ_1)/(μ_0-μ_2) G_20^0,2+μ_0(μ_0-μ_1)/(μ_1-μ_2) G_20^1,2 = 
       ∑_n=2^∞Y_n0/n!(μ_0 -μ_1)[-μ_2(μ^n_0 -μ^n_1) +μ_1(μ^n_0-μ^n_2) -μ_0 (μ^n_1-μ^n_2) ] 
       +∑_[m=3]n=2^∞Y_nm/n!m![μ_2 (μ^n_0 -μ^n_1)(μ^m_0μ_1 - μ_0μ^m_1) -μ_1(μ_0-μ_1)/(μ_0-μ_2)
    	(μ^n_0 -μ^n_2)(μ^m_0μ_2 - μ_0μ^m_2) . 
       .+ μ_0(μ_0-μ_1)/(μ_1-μ_2)(μ^n_1 -μ^n_2)(μ^m_1μ_2 - μ_1μ^m_2) ]      .

	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 <ref> to recast the terms that contain the Y_nm with n≥ 2 and m≥ 3. In particular, according to (<ref>), we obtain:
	
    Y_nm/n!m!μ_0μ_1μ_2 (μ_1-μ_2)(μ_0-μ_2)[∑_k=0^m-2μ_2^k (μ_0^m-2-k - μ_1^m-2-k)]
    	[∑_k=0^n-1μ_2^k (μ_0^n-1-k-μ_1^n-1-k)]    ,

	and according to (<ref>) the terms that contain the yields Y_n0 can be written as:
	
    Y_n0/n!(μ_0 - μ_2)(μ_1 - μ_2)(μ_0 -μ_1)  ∑_k=0^n-1μ_2^k (μ_0^n-1-k-μ_1^n-1-k)    .

	Like in the case of Y_02 (see <ref>), a valid upper bound on Y_20 is thus obtained setting all the other yields to zero in (<ref>). We obtain:
	
    Y^U_20=2μ_2 G_20^0,1/μ_0-μ_1 -μ_1 G_20^0,2/μ_0-μ_2 +μ_0 G_20^1,2/μ_1-μ_2/(μ_0 - μ_2)(μ_1 - μ_2)(μ_0 -μ_1)  .

	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:
	
    Y^U_04=4!μ_2 G_20^0,1/μ_0-μ_1 -μ_1 G_20^0,2/μ_0-μ_2 +μ_0 G_20^1,2/μ_1-μ_2/μ_1(μ^4_0 - μ^4_2)-μ_0 (μ^4_1 - μ^4_2) -μ_2(μ^4_0 -μ^4_1)  .


	
	

 §.§ Upper bound on Y_13
 
	We look for that combination of gains in which all the terms proportional to Y_n0,Y_n1,Y_0m and Y_2m are removed. In order to find it, we consider the most general combination of all gains:
	
    G_13=∑_i,j=0^2 c_i,jQ̃^i,j = ∑_n,m=0^∞Y_nm/n!m![∑_i,j=0^2 c_i,jμ_i^n μ_j^m ]   ,

	and impose proper conditions on the real coefficients c_i,j:
	
    Y_n0  ∑_i=0^2μ_i^n (∑_j=0^2 c_i,j) =0   ∀  n   ⇐   c_i,0+c_i,1+c_i,2=0 
    	          i=0,1,2  
       Y_n1  ∑_i=0^2μ_i^n (∑_j=0^2μ_j c_i,j) =0   ∀  n   ⇐  μ_0 c_i,0+ μ_1 c_i,1+μ_2 c_i,2=0      i=0,1,2  
       Y_0m  ∑_j=0^2μ_j^m (∑_i=0^2 c_i,j) =0   ∀  m   ⇐   c_0,j+c_1,j+c_2,j=0 
    	     j=0,1,2  
       Y_2m  ∑_j=0^2μ_j^m (∑_i=0^2μ_i^2 c_i,j) =0   ∀  m   ⇐  μ_0^2 c_0,j+μ_1^2 c_1,j+ μ_2^2 c_2,j=0      j=0,1,2   .

	The conditions given by equations (<ref>, <ref>, <ref>, <ref>) 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_0,0=1   ,
       c_0,1=-(μ_0-μ_2)/μ_1-μ_2  ,
       c_0,2=-1-c_0,1=μ_0-μ_1/μ_1 -μ_2  ,
       c_1,0=-(μ_0^2-μ_2^2)/μ_1^2-μ_2^2  ,
       c_1,1= c_1,0c_0,1 =(μ_0^2-μ_2^2)(μ_0-μ_2)/(μ_1^2-μ_2^2)(μ_1 -μ_2)  ,
       c_1,2= -c_1,0-c_1,1=c_1,0c_0,2=-(μ_0^2-μ_2^2)(μ_0-μ_1)/(μ_1^2-μ_2^2)(μ_1 -μ_2)  ,
       c_2,0=-1-c_1,0=μ_0^2-μ_1^2/μ_1^2-μ_2^2  ,
       c_2,1=-c_0,1-c_1,1=c_0,1c_2,0= (μ^2_1-μ^2_0)(μ_0-μ_2)/(μ_1^2-μ_2^2)(μ_1 -μ_2)  ,
       c_2,2=-c_2,0-c_2,1=(1+c_1,0)(1+c_0,1)=(μ_0^2-μ_1^2)(μ_0-μ_1)/(μ_1^2-μ_2^2)(μ_1-μ_2)  .

	By substituting (<ref>) back into (<ref>) we get an expression in which the terms Y_n0,Y_n1,Y_0m and Y_2m are removed:
	
    G_13=    ∑_m=2^∞Y_1m/m!(μ_0-μ_1)(μ_0-μ_2)/(μ_1-μ_2)(μ_1+μ_2)· A_22(μ_0,μ_1,μ_2,m) 
       + ∑_[m=2]n=3^∞Y_nm/n!m!A_22(μ_0,μ_1,μ_2,m)· A_11(μ_0,μ_1,μ_2,n)/(μ_1-μ_2)^2(μ_1+μ_2)  ,

	where A_22 is the factor given by (<ref>) also present in the bounds for Y_02 and Y_22, whereas A_11 is the factor given by (<ref>) which appears in the bound on Y_11. Note that this is somehow expected: when bounding Y_02 and Y_22 we removed the terms Y_n0 and Y_n1 as we just did for Y_13, and in bounding Y_11 we removed the terms Y_0m and Y_2m as we did here. Therefore, by exploiting the result given by (<ref>) we can recast each term of the sum corresponding to the Y_1m, with m≥ 2, in (<ref>) as:
	
    -Y_1m/m!(μ_0-μ_2)^2/(μ_1+μ_2)(μ_0-μ_1)∑_k=0^m-1μ_2^k (μ_0^m-1-k-μ_1^m-1-k)  ,

	and realize that it is always negative, regardless of the value of the intensities.

	By employing the results (<ref>, <ref>) we can recast each term of the sum corresponding to the Y_nm with n≥ 3 and m≥ 2, in (<ref>) as:
	
    Y_nm/n!m!(μ_0 - μ_2)^2(μ_1-μ_2)^2/(μ_1-μ_2)^2(μ_1+μ_2)(μ_0-μ_1)∑_k=0^m-1μ_2^k (μ_0^m-1-k-μ_1^m-1-k)F(n)  ,
 
	and realize that it is always positive[F(n) is defined in (<ref>).], 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 (<ref>). As a result we obtain:
	
    G_13=    -Y^U_13/3!(μ_0-μ_2)^2/(μ_1+μ_2)(μ_0-μ_1)[μ_0^2-μ_1^2+μ_2(μ_0-μ_1)] 
       + ∑_[m=2]n=3^∞[μ_1^m (μ_0 - μ_2) + μ_2^m (μ_1 - μ_0) + μ_0^m (μ_2 - μ_1)]·[μ_1^n (μ^2_0 - μ^2_2) + μ_2^n (μ^2_1 - μ^2_0) + μ_0^n (μ^2_2 - μ^2_1)]/n!m!(μ_1-μ_2)^2(μ_1+μ_2) ,

	which implies:
	
    Y^U_13/6(μ_0-μ_2)^2(μ_0-μ_1)^2(μ_0+μ_1+μ_2)/μ_1+μ_2 =-G_13
        + (e^μ_1-μ_1-1)(μ_0 - μ_2) + (e^μ_2-μ_2-1)(μ_1 - μ_0) + (e^μ_0-μ_0-1)(μ_2 - μ_1)/(μ_1-μ_2)^2(μ_1+μ_2)
       ×[(e^μ_1-μ_1^2/2-μ_1-1)(μ^2_0 - μ^2_2) + (e^μ_2-μ_2^2/2-μ_2-1)(μ^2_1 - μ^2_0) + 
    		(e^μ_0-μ_0^2/2-μ_0-1)(μ^2_2 - μ^2_1)]   .

	We thus obtain the following upper bound on Y_13:
	
    Y^U_13 = -6(μ_1+μ_2)G_13/(μ_0-μ_2)^2(μ_0-μ_1)^2(μ_0+μ_1+μ_2)+ 6/(μ_0-μ_2)^2(μ_1-μ_2)^2(μ_0-μ_1)^2(μ_0+μ_1+μ_2)
       ×[e^μ_2(μ_1-μ_0)+e^μ_1(μ_0-μ_2)+e^μ_0(μ_2-μ_1)]  
       ×[e^μ_2(μ^2_1-μ^2_0)+e^μ_1(μ^2_0-μ^2_2)+e^μ_0(μ^2_2-μ^2_1)-(μ_0-μ_1)(μ_1-μ_2)(μ_0-μ_2)]
    		  ,

	where G_13 is defined in (<ref>) and the coefficients of the combination of gains in (<ref>).
	
	
	

 §.§ 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:
	
    G_31=∑_i,j=0^2 c_i,jQ̃^i,j = ∑_n,m=0^∞Y_nm/n!m![∑_i,j=0^2 c_i,jμ_i^n μ_j^m ]   ,

	and impose proper conditions on the real coefficients c_i,j:
	
    Y_n0  ∑_i=0^2μ_i^n (∑_j=0^2 c_i,j) =0   ∀  n   ⇐   c_i,0+c_i,1+c_i,2=0 
    	     i=0,1,2  
       Y_n2  ∑_i=0^2μ_i^n (∑_j=0^2μ^2_j c_i,j) =0   ∀  n   ⇐  μ^2_0 c_i,0+ μ^2_1 c_i,1+μ^2_2 c_i,2=0      i=0,1,2  
       Y_0m  ∑_j=0^2μ_j^m (∑_i=0^2 c_i,j) =0   ∀  m   ⇐   c_0,j+c_1,j+c_2,j=0 
    	     j=0,1,2  
       Y_1m  ∑_j=0^2μ_j^m (∑_i=0^2μ_i c_i,j) =0   ∀  m   ⇐  μ_0 c_0,j+μ_1 c_1,j+ μ_2 c_2,j=0      j=0,1,2   .

	The conditions (<ref>, <ref>, <ref>, <ref>) 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_0,0=1   ,
       c_0,1=-(μ_0^2-μ_2^2)/μ_1^2-μ_2^2  ,
       c_0,2=-1-c_0,1=μ_0^2-μ_1^2/μ_1^2-μ_2^2  ,
       c_1,0=-(μ_0-μ_2)/μ_1-μ_2  ,
       c_1,1= c_1,0c_0,1 =(μ_0^2-μ_2^2)(μ_0-μ_2)/(μ_1^2-μ_2^2)(μ_1 -μ_2)  ,
       c_1,2= -c_1,0-c_1,1=c_1,0c_0,2=(μ^2_1-μ^2_0)(μ_0-μ_2)/(μ_1^2-μ_2^2)(μ_1 -μ_2)  ,
       c_2,0=-1-c_1,0=μ_0-μ_1/μ_1 -μ_2  ,
       c_2,1=-c_0,1-c_1,1=c_0,1c_2,0= -(μ_0^2-μ_2^2)(μ_0-μ_1)/(μ_1^2-μ_2^2)(μ_1 -μ_2)  ,
       c_2,2=-c_2,0-c_2,1=(1+c_1,0)(1+c_0,1)=(μ_0^2-μ_1^2)(μ_0-μ_1)/(μ_1^2-μ_2^2)(μ_1-μ_2)  .

	By substituting (<ref>) back into (<ref>) we get an expression in which the terms Y_n0,Y_n2,Y_0m and Y_1m are removed:
	
    G_31=    ∑_n=2^∞Y_n1/n!(μ_0-μ_1)(μ_0-μ_2)/(μ_1-μ_2)(μ_1+μ_2)· A_22(μ_0,μ_1,μ_2,n) + 
       ∑_[m=3]n=2^∞Y_nm/n!m!A_22(μ_0,μ_1,μ_2,n)· A_11(μ_0,μ_1,μ_2,m)/(μ_1-μ_2)^2(μ_1+μ_2)  ,

	where A_22 and A_11 are again the factors from Y_22 and Y_11 bounds given by equations (<ref>,<ref>), similarly to what happens when bounding Y_13 (see <ref>). Therefore the analysis of the coefficients' sign is the same as in <ref>. Hence a valid upper bound on Y_31 is obtained by setting Y_n1→ 0 (except for Y_31) and Y_nm→ 1 in (<ref>) for all n≥ 2 and m≥ 3 in (<ref>).
	Analogous steps to those in <ref> lead to the following upper bound:
	
    Y^U_31 = -6(μ_1+μ_2)G_31/(μ_0-μ_2)^2(μ_0-μ_1)^2(μ_0+μ_1+μ_2)+ 6/(μ_0-μ_2)^2(μ_1-μ_2)^2(μ_0-μ_1)^2(μ_0+μ_1+μ_2)
       ×[e^μ_2(μ_1-μ_0)+e^μ_1(μ_0-μ_2)+e^μ_0(μ_2-μ_1)]  
       ×[e^μ_2(μ^2_1-μ^2_0)+e^μ_1(μ^2_0-μ^2_2)+e^μ_0(μ^2_2-μ^2_1)-(μ_0-μ_1)(μ_1-μ_2)(μ_0-μ_2)]
    	  ,

	where G_31 is defined in (<ref>) and the coefficients of the combination of gains in (<ref>).
	
	
	

 §.§ Upper bound on Y_00

	Consider the following combinations of gains in which all the terms Y_1m and Y_n1 are removed:
	
    G_00^0,1   = μ^2_1 Q̃^0,0+μ^2_0 Q̃^1,1 -μ_0 μ_1 (Q̃^0,1 +Q̃^1,0) 
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0 μ_1 -μ_0 μ^n_1) (μ^m_0 μ_1 -μ_0 μ^m_1)  ;
    
    		G_00^0,2   = μ^2_2 Q̃^0,0+μ^2_0 Q̃^2,2 -μ_0 μ_2 (Q̃^0,2 +Q̃^2,0) 
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_0 μ_2 -μ_0 μ^n_2) (μ^m_0 μ_2 -μ_0 μ^m_2)  ;
    
    		G_00^1,2   = μ^2_2 Q̃^1,1+μ^2_1 Q̃^2,2 -μ_1 μ_2 (Q̃^1,2 +Q̃^2,1) 
       = ∑_n,m=0^∞Y_nm/n!m!(μ^n_1 μ_2 -μ_1 μ^n_2) (μ^m_1 μ_2 -μ_1 μ^m_2)   .

	We now combine G_00^0,1,G_00^0,2 and G_00^1,2 with arbitrary real coefficients c_0 and c_1 and impose that the terms Y_2m and Y_n2 are also removed in the resulting expression:
	
    G_00^0,1+c_0  G_00^0,2+c_1  G_00^1,2 =  ∑_n,m=0^∞Y_nm/n!m![(μ^n_0 μ_1 -μ_0 μ^n_1) 
    		(μ^m_0 μ_1 -μ_0 μ^m_1) . 
       . +c_0 (μ^n_0 μ_2 -μ_0 μ^n_2) (μ^m_0 μ_2 -μ_0 μ^m_2) 
    		+c_1 (μ^n_1 μ_2 -μ_1 μ^n_2) (μ^m_1 μ_2 -μ_1 μ^m_2)]    .

	For Y_2m and Y_n2 to be removed it suffices that for every m it holds:
	
    (μ_0^2μ_1-μ_0μ^2_1)(μ^m_0 μ_1 -μ_0 μ^m_1)+c_0 (μ_0^2μ_2-μ_0μ^2_2)(μ^m_0 μ_2 -μ_0 μ^m_2)+c_1(μ_1^2μ_2-μ_1μ^2_2)(μ^m_1 μ_2 -μ_1 μ^m_2) = 0  ,

	which is fulfilled by:
	
    c_0    = -μ^2_1(μ_0-μ_1)/μ^2_2(μ_0-μ_2)  , 
    
    		c_1    =μ^2_0(μ_0-μ_1)/μ^2_2(μ_1-μ_2)  .

	Substituting (<ref>) and (<ref>) back into (<ref>) and multiplying both sides by μ^2_2, we get an expression where all the terms Y_0m,Y_2m,Y_n0 and Y_n2 are removed and where the term Y_00 gives the largest contribution. More precisely, we find that:
	
    μ^2_2 G_00^0,1 -μ^2_1(μ_0-μ_1)/(μ_0-μ_2) G_00^0,2 + μ^2_0(μ_0-μ_1)/(μ_1-μ_2) G_00^1,2 = 
       Y_00[μ_2^2 (μ_0 -μ_1)^2 -μ_1^2 (μ_0 -μ_1)(μ_0 -μ_2) +μ_0^2 (μ_0 -μ_1)(μ_1 -μ_2)]  
       +∑_m=3^∞Y_0m/m![μ_2^2 (μ_1 -μ_0)(μ^m_0 μ_1 -μ_0 μ^m_1) + μ_1^2(μ_0 -μ_1)(μ^m_0 μ_2 -μ_0 μ^m_2) -
    		        μ_0^2 (μ_0-μ_1)(μ^m_1 μ_2 -μ_1 μ^m_2)]  
       +∑_n=3^∞Y_n0/n![μ_2^2 (μ_1 -μ_0)(μ^n_0 μ_1 -μ_0 μ^n_1) + μ_1^2(μ_0 -μ_1)(μ^n_0 μ_2 -μ_0 μ^n_2) -
    		  	    μ_0^2 (μ_0-μ_1)(μ^n_1 μ_2 -μ_1 μ^n_2)]  
       +∑_n,m=3^∞Y_nm/n!m!μ^2_0μ^2_1μ^2_2 [(μ^n-1_0 -μ^n-1_1)(μ^m-1_0 -μ^m-1_1) 
       - (μ_0-μ_1)/(μ_0-μ_2)(μ^n-1_0 -μ^n-1_2)(μ^m-1_0 -μ^m-1_2) +(μ_0-μ_1)/(μ_1-μ_2)(μ^n-1_1 -μ^n-1_2)(μ^m-1_1 -μ^m-1_2)]   .

	In order to extract an upper bound on Y_00 we need to study the sign of the yields' coefficients. We start by recasting the term corresponding to Y_00 as:
	
    Y_00 (μ_0 -μ_1)[μ_2^2(μ_0 -μ_1)-μ_1^2 (μ_0 -μ_2) +μ_0^2 (μ_1 -μ_2)]  
       = Y_00 (μ_0 -μ_1)^2 (μ_1-μ_2)(μ_0-μ_2)   .

	We observe that the sign of this expression is determined by the factors (μ_1-μ_2)(μ_0-μ_2).

	We then proceed by recasting each term of the sum corresponding to the Y_nm, with n,m≥ 3 in (<ref>) as:
	
    Y_nm/n!m!A_00(μ_0,μ_1,μ_2,m)· A_00(μ_0,μ_1,μ_2,n)/(μ_0-μ_2)(μ_1-μ_2)  ,

	where
	
    A_00(μ_0,μ_1,μ_2,m)    ≡μ_1^m (μ^2_2μ_0 - μ_2μ^2_0) + μ_2^m (μ^2_0μ_1 - μ_0μ^2_1) + μ_0^m (μ^2_1μ_2 - μ_1μ_2^2)   .

	This factor can be rewritten as:
	
    A_00(μ_0,μ_1,μ_2,m)    = μ_0μ_1 μ_2[μ_1^m-1 (μ_2 - μ_0) + μ_2^m-1 (μ_0 - μ_1) + μ_0^m-1 (μ_1 - μ_2)] 
       = -μ_0 μ_1 μ_2   A_22(μ_0,μ_1,μ_2,m-1)   ,

	where A_22 is defined as (<ref>) in <ref>. Thus we can use the result (<ref>) obtained in <ref> to directly recast A_00 as:
	
    A_00(μ_0,μ_1,μ_2,m)=  μ_0μ_1 μ_2 (μ_0 - μ_2)(μ_1 - μ_2)∑_k=0^m-2μ_2^k (μ_0^m-2-k-μ_1^m-2-k)    .

	By substituting (<ref>) back into (<ref>), we get the final expression for each term of the sum corresponding to the Y_nm, with n,m≥ 3 in (<ref>):
	
    Y_nm/n!m!μ^2_0 μ^2_1 μ^2_2(μ_0-μ_2)(μ_1-μ_2)[∑_k=0^m-2μ_2^k (μ_0^m-2-k-μ_1^m-2-k)][∑_k=0^n-2μ_2^k (μ_0^n-2-k-μ_1^n-2-k)]
    		     .

	which has manifestly the same sign as the expression given by (<ref>), for any value of the intensities (the product of the last two factors is always positive).

	Finally, we recast the Y_0m's terms (Y_n0's terms are identical under the replacement m→ n) as:
	
    Y_0m/m!μ_0μ_1 μ_2(μ_0-μ_1) [μ_2 (μ^m-1_1  - μ^m-1_0) + μ_1(μ^m-1_0 - μ^m-1_2) - μ_0(μ^m-1_1 - μ^m-1_2)]  
       =Y_0m/m!(μ_0-μ_1) A_00(μ_0,μ_1,μ_2,m) 
       = Y_0m/m!μ_0μ_1 μ_2 (μ_0-μ_1) (μ_0 - μ_2)(μ_1 - μ_2)∑_k=0^m-2μ_2^k (μ_0^m-2-k-μ_1^m-2-k)   ,

	where we employed (<ref>) in the first equality and (<ref>) in the second one. We observe that the sign of the Y_0m's terms is again determined by the factors (μ_0 - μ_2)(μ_1 - μ_2).

	We conclude that the coefficients of Y_0m, Y_n0 and Y_nm, with n,m≥ 3, carry the same sign as Y_00's, which implies that a valid upper bound on Y_00 is obtained by setting all the other yields to zero in (<ref>). In so doing, we find that:
	
    Y^U_00=μ_2^2 G_00^0,1/μ_0 - μ_1 -μ_1^2 G_00^0,2/μ_0-μ_2 + μ_0^2 G_00^1,2/μ_1-μ_2/(μ_0-μ_1)(μ_0-μ_2)(μ_1-μ_2)  .

	
	
	

§ YIELDS' BOUNDS WITH FOUR DECOYS
 
	Here we derive analytical upper bounds on the yields appearing in (<ref>), following the same lines of <ref>. In this case we assume that Alice and Bob can prepare their phase-randomized coherent pulses with four different intensity settings: {μ_0,μ_1,μ_2,μ_3}, 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 √(η) <cit.>.

	The whole set of infinite yields is subjected to the following sixteen equality constraints:
	
    Q̃^k,l≡ e^μ_k + μ_l Q^k,l =∑_n,m=0^∞Y_nm/n!m!μ_k^n μ_l^m    k,l ∈{0,1,2,3}  ,

	and to the same inequality constraints given by (<ref>).

	In this appendix we only obtain bounds on the yields Y_13,Y_31,Y_04 and Y_40 since the bounds derived on the yields Y_00, Y_11, Y_02, Y_20 and Y_22 in <ref> are already good enough, i.e bounding them with one additional decoy intensity would not result in a significant improvement of the performance of the protocol.
	
	

 §.§ Upper bound on Y_04
 
	Consider the following combinations of gains in which all the terms Y_1m and Y_n0 are removed:
	
    G_04^i,j   = μ_jQ̃^i,i+ μ_iQ̃^j,j - μ_jQ̃^i,j -μ_iQ̃^j,i= ∑_n,m=0^∞Y_nm/n!m!(μ_j μ_i^n -μ_i μ_j^n)(μ_i^m -μ_j^m)    ,

	where i,j ∈{0,1,2,3}. Since G_04^i,i=0 and G_04^i,j=G_04^j,i, we only have six distinct combinations that read (for j>i): G_04^0,1,G_04^0,2,G_04^0,3,G_04^1,2,G_04^1,3,G_04^2,3.

	We now take the linear combination of the G_04^i,j such that even the yields Y_2m,Y_3m,Y_n1 and Y_n2 are removed:
	
    ∑_j>i c_i,j G_04^i,j = ∑_n,m=0^∞Y_nm/n!m!∑_j>i c_i,j (μ_j μ_i^n -μ_i μ_j^n)(μ_i^m -μ_j^m)    ,

	where we implicitly assume that both indexes i,j run over the set {0,1,2,3}.
	For Y_2m,Y_3m,Y_n1 and Y_n2 to be removed, the real coefficients c_i,j must satisfy:
	
    {[ ∑_j>i c_i,j (μ_j μ_i^2 -μ_i μ_j^2)(μ_i^m -μ_j^m)                                                =                                         0    ∀ m; ∑_j>i c_i,j (μ_j μ_i^3 -μ_i μ_j^3)(μ_i^m -μ_j^m)                                                =                                         0    ∀ m;     ∑_j>i c_i,j (μ_j μ_i^n -μ_i μ_j^n)(μ_i -μ_j)                                                =                                         0    ∀ n; ∑_j>i c_i,j (μ_j μ_i^n -μ_i μ_j^n)(μ_i^2 -μ_j^2)                                                =                                         0    ∀ n;                                                  ].

	In order to solve system (<ref>), we look for those coefficients c_i,j such that the multiplicative factors of μ_i^m and μ_i^n (for i=0,1,2,3) are all set to zero. This corresponds to imposing sixteen conditions on the six coefficients c_i,j. These conditions are not all independent, and a solution can be found even when we require (for simplicity) that c_0,1=1:
	
    c_0,1=1   ,
       c_0,2= -(μ_0-μ_1)μ_1(μ_1-μ_3)/(μ_0-μ_2)μ_2 (μ_2-μ_3)  , 
       c_0,3= (μ_0-μ_1)μ_1(μ_1-μ_2)/(μ_0-μ_3)μ_3 (μ_2-μ_3)  ,
       c_1,2= (μ_0-μ_1)μ_0(μ_0-μ_3)/(μ_1-μ_2)μ_2 (μ_2-μ_3)  ,
       c_1,3= -(μ_0-μ_1)μ_0(μ_0-μ_2)/(μ_1-μ_3)μ_3 (μ_2-μ_3)  ,
       c_2,3= μ_0μ_1(μ_0-μ_1)^2/μ_2μ_3 (μ_2-μ_3)^2  .

	By substituting the solution for the coefficients given by (<ref>) back into (<ref>), one gets:
	
    ∑_j>i c_i,j G_04^i,j = ∑_m=3^∞Y_0m/m! A_04(μ_0,μ_1,μ_2,μ_3,m) + ∑_[m=3]n=4^∞Y_nm/n!m! 
    	B_04(μ_0,μ_1,μ_2,μ_3,n,m)  ,

	where:
	
    A_04(μ_0,μ_1,μ_2,μ_3,m) = -(μ_0 -μ_1)/μ_2 μ_3 (μ_2 - μ_3)[μ_0^m (μ_1 -μ_2)(μ_1 -μ_3)(μ_2 -μ_3)-μ_1^m (μ_0 -μ_2)(μ_0 -μ_3)(μ_2 -μ_3)  . 
       .+μ_2^m (μ_0 -μ_1)(μ_0 -μ_3)(μ_1 -μ_3)- μ_3^m (μ_0 -μ_1)(μ_0 -μ_2)(μ_1 -μ_2) ]   
       = -(μ_0 -μ_1)^2 (μ_0 -μ_2)(μ_1 -μ_2) (μ_0 -μ_3)(μ_1 -μ_3)/μ_2 μ_3(∑_i_1 ≤ i_2 ≤…≤ i_m-3μ_i_1μ_i_2·…·μ_i_m-3)   ,

	and
	
    B_04(μ_0,μ_1,μ_2,μ_3,n,m) = -μ_0 μ_1/(μ_0 -μ_2)(μ_1 -μ_2)(μ_1 -μ_3)(μ_0 -μ_3)(μ_2 -μ_3)^2[ μ_0^m (μ_1 -μ_2)(μ_1 -μ_3)(μ_2 -μ_3). 
       .-μ_1^m (μ_0 -μ_2)(μ_0 -μ_3)(μ_2 -μ_3) +μ_2^m (μ_0 -μ_1)(μ_0 -μ_3)(μ_1 -μ_3)- μ_3^m (μ_0 -μ_1)(μ_0 -μ_2)(μ_1 -μ_2) ] 
       ×[- μ_0^n-1 (μ_1 -μ_2)(μ_1 -μ_3)(μ_2 -μ_3)+μ_1^n-1 (μ_0 -μ_2)(μ_0 -μ_3)(μ_2 -μ_3) . 
       .-μ_2^n-1 (μ_0 -μ_1)(μ_0 -μ_3)(μ_1 -μ_3)+μ_3^n-1 (μ_0 -μ_1)(μ_0 -μ_2)(μ_1 -μ_2)] 
       =  -μ_0μ_1 μ_2 μ_3  A_04(μ_0,μ_1,μ_2,μ_3,m)·(∑_i_1 ≤ i_2 ≤…≤ i_n-4μ_i_1μ_i_2·…·μ_i_n-4)   .

	In (<ref>,<ref>) we again assume that the indexes in the sums run over the set {0,1,2,3} and we define ∑_i_1 ≤ i_2 ≤…≤ i_m-3μ_i_1μ_i_2·…·μ_i_m-3|_m=3 =1. From (<ref>) we deduce that the sign of Y_0m's coefficient is independent of m, while from (<ref>) we deduce that Y_nm's coefficient has always opposite sign to that of Y_0m. Therefore a valid upper bound on Y_04 is obtained by setting to zero all the other yields Y_0m and to 1 the yields Y_nm with n≥ 4 and m≥ 3 in (<ref>). We thus obtain:
	
    ∑_j>i c_i,j G_04^i,j = Y^U_04/4! A_04(μ_0,μ_1,μ_2,μ_3,4) + ∑_[m=3]n=4^∞B_04(μ_0,μ_1,μ_2,μ_3,n,m)/n!m!  ,

	which implies the following upper bound on Y_04:
	
    Y^U_04 = 4!/A_04(μ_0,μ_1,μ_2,μ_3,4)[∑_j>i c_i,j G_04^i,j - ∑_[m=3]n=4^∞B_04(μ_0,μ_1,μ_2,μ_3,n,m)/n!m!]    ,

	where c_i,j are given in (<ref>), G_04^i,j is defined in (<ref>), the coefficient A_04 reads:
	
    A_04(μ_0,μ_1,μ_2,μ_3,4) = -(μ_0 -μ_1)^2 (μ_0 -μ_2)(μ_1 -μ_2) (μ_0 -μ_3)(μ_1 -μ_3)(μ_0+μ_1+μ_2+μ_3)/μ_2 μ_3  ,

	and the sum over the coefficient B_04 reads:
	
    ∑_[m=3]n=4^∞B_04(μ_0,μ_1,μ_2,μ_3,n,m)/n!m! = μ_0 μ_1/(μ_0 -μ_2)(μ_1 -μ_2)(μ_1 -μ_3)(μ_0 -μ_3)(μ_2 -μ_3)^2  
       ×[ (e^μ_0-1-μ_0 -μ_0^2/2) (μ_1 -μ_2)(μ_1 -μ_3)(μ_2 -μ_3)-(e^μ_1-1-μ_1 -μ_1^2/2) (μ_0 -μ_2)(μ_0 -μ_3)(μ_2 -μ_3) . 
       . +(e^μ_2-1-μ_2 -μ_2^2/2) (μ_0 -μ_1)(μ_0 -μ_3)(μ_1 -μ_3)- (e^μ_3-1-μ_3 -μ_3^2/2) (μ_0 -μ_1)(μ_0 -μ_2)(μ_1 -μ_2) ]^2    .


	
	

 §.§ Upper bound on Y_40
 
	Consider the following combinations of gains in which all the terms Y_0m and Y_n1 are removed:
	
    G_40^i,j   = μ_jQ̃^i,i+ μ_iQ̃^j,j - μ_iQ̃^i,j -μ_jQ̃^j,i= ∑_n,m=0^∞Y_nm/n!m!(μ_i^n -μ_j^n)(μ_j μ_i^m -μ_i μ_j^m)   ,

	where i,j ∈{0,1,2,3}. Since G_40^i,i=0 and G_40^i,j=G_40^j,i, we only have six distinct combinations that read (for j>i): G_40^0,1,G_40^0,2,G_40^0,3,G_40^1,2,G_40^1,3,G_40^2,3.

	We now take the linear combination of the G_40^i,j such that even the yields Y_1m,Y_2m,Y_n2 and Y_n3 are removed:
	
    ∑_j>i c_i,j G_40^i,j = ∑_n,m=0^∞Y_nm/n!m!∑_j>i c_i,j (μ_i^n -μ_j^n)(μ_j μ_i^m -μ_i μ_j^m)     ,

	where we implicitly assume that both indexes i,j run over the set {0,1,2,3}.
	For Y_1m,Y_2m,Y_n2 and Y_n3 to be removed, the real coefficients c_i,j must satisfy:
	
    {[ ∑_j>i c_i,j (μ_i^n -μ_j^n)(μ_j μ_i^2 -μ_i μ_j^2)                                                =                                         0    ∀ n; ∑_j>i c_i,j (μ_i^n -μ_j^n)(μ_j μ_i^3 -μ_i μ_j^3)                                                =                                         0    ∀ n;     ∑_j>i c_i,j (μ_i -μ_j)(μ_j μ_i^m -μ_i μ_j^m)                                                =                                         0    ∀ m; ∑_j>i c_i,j (μ_i^2 -μ_j^2)(μ_j μ_i^m -μ_i μ_j^m)                                                =                                       0    ∀ m .;                                                  ].

	We now notice that the system (<ref>) is exactly the same system solved in <ref> while bounding Y_04, thus the solution for the coefficients c_i,j is given in (<ref>).
	By substituting the solution (<ref>) back into (<ref>), one gets:
	
    ∑_j>i c_i,j G_40^i,j = ∑_n=3^∞Y_n0/n! A_04(μ_0,μ_1,μ_2,μ_3,n) + ∑_[m=4]n=3^∞Y_nm/n!m! B_04(μ_0,μ_1,μ_2,μ_3,m,n)      ,

	where A_04 and B_04 are the coefficients defined in (<ref>,<ref>) while bounding Y_04. Hence we can adopt the observations made on the sign of A_04 and B_04 from <ref> and conclude that a valid upper bound on Y_40 is obtained by setting to zero all the other yields Y_n0 and to 1 the yields Y_nm with n≥ 3 and m≥ 4 in (<ref>). The upper bound on Y_40 then reads:
	
    Y^U_40 = 4!/A_04(μ_0,μ_1,μ_2,μ_3,4)[∑_j>i c_i,j G_40^i,j - ∑_[m=4]n=3^∞B_04(μ_0,μ_1,μ_2,μ_3,m,n)/n!m!]    ,

	where c_i,j, G_40^i,j, A_04(μ_0,μ_1,μ_2,μ_3,4) and the sum over B_04 are given in (<ref>), (<ref>), (<ref>) and (<ref>), respectively.
	
	
	

 §.§ Upper bound on Y_13
 
	We consider the most general combination of all sixteen gains:
	
    ∑_i,j=0^3 c_i,jQ̃^i,j = ∑_n,m=0^∞Y_nm/n!m![∑_i,j=0^3 c_i,jμ_i^n μ_j^m ]  ,

	and require that the terms Y_n0,Y_n1,Y_n2,Y_0m,Y_2m and Y_3m are removed, by imposing proper conditions on the real coefficients c_i,j:
	
    Y_n0  ∑_i,j=0^3 c_i,jμ_i^n =0   ∀  n   ⇐  ∑_j=0^3 c_i,j=0       i=0,1,2,3  
       Y_n1  ∑_i,j=0^3 c_i,jμ_i^n μ_j =0   ∀  n   ⇐  ∑_j=0^3 c_i,jμ_j =0      i=0,1,2,3  
       Y_n2  ∑_i,j=0^3 c_i,jμ_i^n μ_j^2 =0   ∀  n   ⇐  ∑_j=0^3 c_i,jμ_j^2 =0      i=0,1,2,3  
       Y_0m  ∑_i,j=0^3 c_i,jμ_j^m=0   ∀  m   ⇐  ∑_i=0^3 c_i,j=0      j=0,1,2,3  
       Y_2m  ∑_i,j=0^3 c_i,jμ_i^2 μ_j^m=0   ∀  m   ⇐  ∑_i=0^3 c_i,jμ_i^2 =0       j=0,1,2,3  
       Y_3m  ∑_i,j=0^3 c_i,jμ_i^3 μ_j^m=0   ∀  m   ⇐  ∑_i=0^3 c_i,jμ_i^3 =0       j=0,1,2,3    .

	The twenty-four conditions given by (<ref> - <ref>) form an over-determined system of equations for the sixteen 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_0,0=1   ,
       c_0,1=(μ_0-μ_2) (μ_0-μ_3)/(μ_2-μ_1) (μ_1-μ_3)  ,
       c_0,2=(μ_0-μ_1) (μ_0-μ_3)/(μ_1-μ_2) (μ_2-μ_3)  , 
       c_0,3=(μ_0-μ_1) (μ_0-μ_2)/(μ_1-μ_3) (μ_3-μ_2)  ,
       c_1,0=-(μ_0-μ_2) (μ_0-μ_3) [μ_0 (μ_2+μ_3)+μ_2 μ_3]/(μ_1-μ_2)(μ_1-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_1,1=(μ_0-μ_2)^2 (μ_0-μ_3)^2 [μ_0 (μ_2+μ_3)+μ_2 μ_3]/(μ_1-μ_2)^2 (μ_1-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_1,2=-(μ_0-μ_1) (μ_0-μ_2) (μ_0-μ_3)^2 [μ_0 (μ_2+μ_3)+μ_2 μ_3]/(μ_1-μ_2)^2 (μ_1-μ_3) (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_1,3=(μ_0-μ_1) (μ_0-μ_2)^2 (μ_0-μ_3) [μ_0 (μ_2+μ_3)+μ_2 μ_3]/(μ_1-μ_2) (μ_1-μ_3)^2 (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_2,0= (μ_0-μ_1) (μ_0-μ_3) [μ_0 (μ_1+μ_3)+μ_1 μ_3]/(μ_1-μ_2) (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_2,1= -(μ_0-μ_1) (μ_0-μ_2) (μ_0-μ_3)^2 [μ_0 (μ_1+μ_3)+μ_1 μ_3]/(μ_1-μ_2)^2 (μ_1-μ_3) (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_2,2= (μ_0-μ_1)^2 (μ_0-μ_3)^2 [μ_0 (μ_1+μ_3)+μ_1 μ_3]/(μ_1-μ_2)^2 (μ_2-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_2,3=-(μ_0-μ_1)^2 (μ_0-μ_2) (μ_0-μ_3) [μ_0 (μ_1+μ_3)+μ_1 μ_3]/(μ_1-μ_2) (μ_1-μ_3) (μ_2-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_3,0=(μ_0-μ_1) (μ_0-μ_2) [μ_0 (μ_1+μ_2)+μ_1 μ_2]/(μ_1-μ_3) (μ_3-μ_2) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_3,1=(μ_0-μ_1) (μ_0-μ_2)^2 (μ_0-μ_3) [μ_0 (μ_1+μ_2)+μ_1 μ_2]/(μ_1-μ_2) (μ_1-μ_3)^2 (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_3,2=-(μ_0-μ_1)^2 (μ_0-μ_2) (μ_0-μ_3) [μ_0 (μ_1+μ_2)+μ_1 μ_2]/(μ_1-μ_2) (μ_1-μ_3) (μ_2-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_3,3=(μ_0-μ_1)^2 (μ_0-μ_2)^2 [μ_0 (μ_1+μ_2)+μ_1 μ_2]/(μ_1-μ_3)^2  (μ_2-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  .

	By substituting these expressions back into (<ref>) and by making some simplifications, one gets:
	
    ∑_i,j=0^3 c_i,jQ̃^i,j = ∑_m=3^∞Y_1m/m! A_13(μ_0,μ_1,μ_2,μ_3,m) + ∑_[m=3]n=4^∞Y_nm/n!m! A_13(μ_0,μ_1,μ_2,μ_3,m)· C_n     ,

	where:
	
    A_13(μ_0,μ_1,μ_2,μ_3,m)    = (μ_0 -μ_1)^2 (μ_0 -μ_2)^2(μ_0 -μ_3)^2 /μ_2 μ_3+ μ_1 μ_2+μ_1 μ_3(∑_i_1 ≤ i_2 ≤…≤ i_m-3μ_i_1μ_i_2·…·μ_i_m-3)   ,

	and C_n (n ≥ 5) is defined recursively as:
	
    {[                                                                                                            C_n                                                                                                              = [∑_j=1^n-4 (μ^j_0+μ^j_1+μ^j_2+μ^j_3) C_n-j -μ_0 μ_1 μ_2 μ_3 (∑_i_1 ≤ i_2 ≤…≤ i_n-5μ_i_1μ_i_2·…·μ_i_n-5)]/(n-4);                                                                                                            C_4                                                                                                              =                                                       μ_0 μ_1 μ_2 +μ_0 μ_1 μ_3 +μ_0 μ_2 μ_3 + μ_1 μ_2 μ_3    . ].

	In (<ref>,<ref>) we assume that the indexes i_j in the sums run over the set {0,1,2,3} and we define ∑_i_1 ≤ i_2 ≤…≤ i_m-3μ_i_1μ_i_2·…·μ_i_m-3|_m=3 =1. From (<ref>) we deduce that the sign of Y_1m's coefficient is always positive, while from (<ref>) we deduce that Y_nm's coefficient has always equal sign to that of Y_1m, since C_n is always a positive quantity. Therefore a valid upper bound on Y_13 is obtained by setting to zero all the other yields  in (<ref>). The upper bound on Y_13 then reads:
	
    Y_13^U = 6/A_13(μ_0,μ_1,μ_2,μ_3,3)(	∑_i,j=0^3 c_i,jQ̃^i,j)  ,

	where c_i,j are defined in (<ref>) and A_13(μ_0,μ_1,μ_2,μ_3,3) reads:
	
    A_13(μ_0,μ_1,μ_2,μ_3,3)=  (μ_0 -μ_1)^2 (μ_0 -μ_2)^2(μ_0 -μ_3)^2 /μ_2 μ_3+ μ_1 μ_2+μ_1 μ_3  .

	
	
	

 §.§ Upper bound on Y_31
 
	We consider the most general combination of all sixteen gains:
	
    ∑_i,j=0^3 c_i,jQ̃^i,j = ∑_n,m=0^∞Y_nm/n!m![∑_i,j=0^3 c_i,jμ_i^n μ_j^m ]   ,

	and require that the terms Y_n0,Y_n2,Y_n3,Y_0m,Y_1m and Y_2m are removed, by imposing proper conditions on the real coefficients c_i,j:
	
    Y_n0  ∑_i,j=0^3 c_i,jμ_i^n =0   ∀  n   ⇐  ∑_j=0^3 c_i,j=0       i=0,1,2,3  
       Y_n2  ∑_i,j=0^3 c_i,jμ_i^n μ_j^2 =0   ∀  n   ⇐  ∑_j=0^3 c_i,jμ_j^2 =0      i=0,1,2,3  
       Y_n3  ∑_i,j=0^3 c_i,jμ_i^n μ_j^3 =0   ∀  n   ⇐  ∑_j=0^3 c_i,jμ_j^3 =0      i=0,1,2,3  
       Y_0m  ∑_i,j=0^3 c_i,jμ_j^m=0   ∀  m   ⇐  ∑_i=0^3 c_i,j=0      j=0,1,2,3  
       Y_1m  ∑_i,j=0^3 c_i,jμ_i μ_j^m=0   ∀  m   ⇐  ∑_i=0^3 c_i,jμ_i =0       j=0,1,2,3  
       Y_2m  ∑_i,j=0^3 c_i,jμ_i^2 μ_j^m=0   ∀  m   ⇐  ∑_i=0^3 c_i,jμ_i^2 =0       j=0,1,2,3    .

	The twenty-four conditions (<ref> - <ref>) form an over-determined system of equations for the sixteen 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_0,0=1   ,
       c_0,1=-(μ_0-μ_2) (μ_0-μ_3) [μ_0 (μ_2+μ_3)+μ_2 μ_3]/(μ_1-μ_2)(μ_1-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_0,2=(μ_0-μ_1) (μ_0-μ_3) [μ_0 (μ_1+μ_3)+μ_1 μ_3]/(μ_1-μ_2) (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_0,3=(μ_0-μ_1) (μ_0-μ_2) [μ_0 (μ_1+μ_2)+μ_1 μ_2]/(μ_1-μ_3) (μ_3-μ_2) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_1,0=(μ_0-μ_2) (μ_0-μ_3)/(μ_2-μ_1) (μ_1-μ_3)  ,
       c_1,1=(μ_0-μ_2)^2 (μ_0-μ_3)^2 [μ_0 (μ_2+μ_3)+μ_2 μ_3]/(μ_1-μ_2)^2 (μ_1-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_1,2=-(μ_0-μ_1) (μ_0-μ_2) (μ_0-μ_3)^2 [μ_0 (μ_1+μ_3)+μ_1 μ_3]/(μ_1-μ_2)^2 (μ_1-μ_3) (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_1,3=(μ_0-μ_1) (μ_0-μ_2)^2 (μ_0-μ_3) [μ_0 (μ_1+μ_2)+μ_1 μ_2]/(μ_1-μ_2) (μ_1-μ_3)^2 (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_2,0=(μ_0-μ_1) (μ_0-μ_3)/(μ_1-μ_2) (μ_2-μ_3)  ,
       c_2,1=-(μ_0-μ_1) (μ_0-μ_2) (μ_0-μ_3)^2 [μ_0 (μ_2+μ_3)+μ_2 μ_3]/(μ_1-μ_2)^2 (μ_1-μ_3) (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_2,2= (μ_0-μ_1)^2 (μ_0-μ_3)^2 [μ_0 (μ_1+μ_3)+μ_1 μ_3]/(μ_1-μ_2)^2 (μ_2-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_2,3=-(μ_0-μ_1)^2 (μ_0-μ_2) (μ_0-μ_3) [μ_0 (μ_1+μ_2)+μ_1 μ_2]/(μ_1-μ_2) (μ_1-μ_3) (μ_2-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_3,0=(μ_0-μ_1) (μ_0-μ_2)/(μ_1-μ_3) (μ_3-μ_2)  ,
       c_3,1=(μ_0-μ_1) (μ_0-μ_2)^2 (μ_0-μ_3) [μ_0 (μ_2+μ_3)+μ_2 μ_3]/(μ_1-μ_2) (μ_1-μ_3)^2 (μ_2-μ_3) [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_3,2=-(μ_0-μ_1)^2 (μ_0-μ_2) (μ_0-μ_3) [μ_0 (μ_1+μ_3)+μ_1 μ_3]/(μ_1-μ_2) (μ_1-μ_3) (μ_2-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  ,
       c_3,3=(μ_0-μ_1)^2 (μ_0-μ_2)^2 [μ_0 (μ_1+μ_2)+μ_1 μ_2]/(μ_1-μ_3)^2  (μ_2-μ_3)^2 [μ_1 (μ_2+μ_3)+μ_2 μ_3]  .

	By substituting these expressions back into (<ref>) and by making some simplifications, one gets:
	
    ∑_i,j=0^3 c_i,jQ̃^i,j = ∑_n=3^∞Y_n1/n! A_13(μ_0,μ_1,μ_2,μ_3,n) + ∑_[m=4]n=3^∞Y_nm/n!m! A_13(μ_0,μ_1,μ_2,μ_3,n) · C_m     ,

	where A_13 and C_m also appear in <ref> when bounding Y_13 and are defined as (<ref>) and (<ref>), respectively. Thus, following the same lines of <ref>, we conclude that all yields in (<ref>) are multiplied by a positive factor. A valid upper bound on Y_31 is then obtained by setting to zero all the other yields in (<ref>). We obtain:
	
    Y_31^U = 6/A_13(μ_0,μ_1,μ_2,μ_3,3)(	∑_i,j=0^3 c_i,jQ̃^i,j)   ,

	where c_i,j and A_13(μ_0,μ_1,μ_2,μ_3,3) are defined in (<ref>) and (<ref>), respectively.
	

