Covert sensing using floodlight illumination
Christos N. Gagatsos Boulat A. Bash Animesh Datta Zheshen Zhang Saikat Guha




College of Optical Sciences, University of Arizona, 1630 E. University Blvd., Tucson, Arizona 85719, USA

Department of Electrical and Computer Engineering, University of Arizona, 1230 E Speedway Blvd., Tucson, Arizona 85719, USA

Department of Physics, University of Warwick, Coventry, CV4 7AL, UK

Department of Materials Science and Engineering, University of Arizona, 1235 E. James E. Rogers Way, Tucson, Arizona 85721, USA
College of Optical Sciences, University of Arizona, 1630 E. University Blvd., Tucson, Arizona 85719, USA

College of Optical Sciences, University of Arizona, 1630 E. University Blvd., Tucson, Arizona 85719, USA
Department of Electrical and Computer Engineering, University of Arizona, 1230 E Speedway Blvd., Tucson, Arizona 85719, USA




We propose a scheme for covert active sensing using floodlight illumination from a THz-bandwidth amplified spontaneous emission (ASE) source and heterodyne detection. We evaluate the quantum-estimation-theoretic performance limit of covert sensing, wherein a transmitter's attempt to sense a target phase is kept undetectable to a quantum-equipped passive adversary, by hiding the signal photons under the thermal noise floor. Despite the quantum state of each mode of the ASE source being mixed (thermal), and hence inferior compared to the pure coherent state of a laser mode, the thousand-times higher optical bandwidth of the ASE source results in achieving a substantially superior performance compared to a narrowband laser source by allowing the probe light to be spread over many more orthogonal temporal modes within a given integration time. Even though our analysis is restricted to single-mode phase sensing, this system could be applicable extendible for various practical optical sensing applications.

Covert sensing using floodlight illumination
    Saikat Guha
    Received; accepted
============================================



  
Introduction. Quantum technologies and in particular quantum and quantum-enhanced sensing <cit.> have high-impact near-term applications that include microscopy <cit.>, vibrometry <cit.>, ranging <cit.>, astronomy <cit.>, and medical imaging <cit.>. Moreover, they have contributed to fundamental discoveries such as gravitational waves detection <cit.>. In particular, the problem of high-precision estimation of a small unknown phase modulation has attracted a lot of attention in the quantum sensing literature, because of its applicability to many of the aforesaid real-life optical imaging and sensing problems. The quantum-physical modeling of a typical sensor results in the quantum state ρ(θ) of a single mode of the target-modulated light encoding a parameter of interest θ. We will be interested in scenarios where n ≈ WT independent temporal modes are modulated identically by a static target resulting in a target-return state of the form, ρ(θ)^⊗ n, where W (Hz) is the optical bandwidth of the optical source, and T (s) is the interrogation time. One desires to find the optimal measurement on this target-return light that generates an estimate θ̂_n with the minimum possible deviation from the true value, often quantified as the variance, or the mean squared error (MSE), ⟨ (θ-θ̂_n)^2. The quantum Cramér-Rao bound (QCRB) <cit.> gives a lower bound on the MSE attainable over all measurements generating an unbiased estimate of θ, which is given by: ⟨ (θ-θ̂_n)^2⟩≥1/(nF), where F is the quantum Fisher information (QFI) of θ, calculated on ρ(θ). The bound gets asymptotically tight as n →∞. In this Letter, we discuss the notion of quantum-secure covert sensing, inspired by recent work on covert communications <cit.>. The task is to minimize the MSE of sensing an unknown carrier-phase modulation θ by using active illumination, while ensuring that an all-powerful non-line-of-sight passive adversary Willie who is able to observe the entire electromagnetic background around the sensor's field of operation equipped with all measurements allowed by quantum physics, and who has knowledge of the exact probing interval and the true value of θ, cannot reliably distinguish between the two equally-likely hypotheses: H_0 (the target is not being probed), and H_1 (the target is being probed). A sensing scheme is ϵ-covert if we can ensure Willie's detection error probability to satisfy, 1/2 - ϵ≤ℙ_e^(det)≤ 1/2, ϵ>0, i.e., his ability to distinguish the two hypotheses is ϵ-close to that of flipping a coin. Our goal is to find—with the optimal choice of an optical transmitter and receiver—the minimum MSE that can be attained as a function of n (which is proportional to the integration time T for a given optical bandwidth W), while ensuring ϵ-covert operation. We note that passive sensing—i.e., sensing a self-luminous or naturally-illuminated scene—is by definition the most covert form of sensing. The quintessential example of passive sensing is human vision. Passive sensing can however be impractical when the illumination levels are so weak that the signal-to-noise ratio at the receiver is insufficient to obtain the desired accuracy within the integration time, or if the target is hidden from direct sight. In this paper, we focus on covert sensing using active illumination <cit.> and consider the problem of phase estimation in the presence of loss and noise.


Problem setup. Figure <ref> is a schematic of the problem setup. Alice (the sensor) generates an n-mode optical probe, each mode of which traverses the forward propagation path to the target, modeled by a single-mode bosonic channel C_1(η_1, n̅_B_1) of transmissivity η_1 and additive thermal noise of mean photon number n̅_B_1 per mode. Each mode of the output light illuminates the same (unknown) target phase θ. The target-modulated light propagates through the return channel C_2(η_2, n̅_B_2) back to Alice. All the lost photons on both the forward and return path channels are given to Willie, who is allowed to make an arbitrary measurement on this 2n-mode optical state in order to do his hypothesis test to determine whether Alice is sensing or not. The problem is to quantify the minimum MSE attainable as a function of n, with optimal choices of Alice's probe state, receiver, and the transmitted mean photon number per mode n̅_S so as to ensure ϵ-covertness no matter what measurement Willie uses. We assume the forward and return path transmissivities η_1 and η_2 to be known to Alice, and different in general, to account for bi-static sensor configurations with unequal path lengths. We also allow for different thermal noise temperatures for the forward and return path channels, to account for scenarios such as one of the two paths being actively jammed. 

Main results. For our achievability results, we will consider two sensor probes: (1) a laser-light (coherent-state) transmitter and heterodyne detection, and (2) a broadband amplified spontaneous emission (ASE) noise source (see Fig. <ref>), whose output is split on a highly unbalanced beam splitter which produces a phase-insensitive classically-correlated two-mode output. The weaker of the two outputs is used to probe the target, while the brighter beam is kept locally at the transmitter to be used as local oscillator to heterodyne detect the target-return light. Our two main results are:

Alice can achieve an MSE, ⟨ (θ - θ̂_n)^2 ⟩ = 𝒪(1/√(n)). Attempting to have the MSE diminish any faster as a function of n must result in detection by Willie with high probability. 


The MSE when Alice uses a coherent state probe state scales as ⟨ (θ - θ̂_n)^2 ⟩ = c_coh/(ϵ√(n)), ϵ>0, while when she uses the ASE source ⟨ (θ - θ̂_n)^2 ⟩ = c_ASE/(ϵ√(n)), ϵ>0. The per-mode performance of the coherent state probe is better than the ASE source probe, i.e., c_coh≤ c_ASE. However, the ASE source can achieve a much lower MSE for a given integration time T, since n ≈ WT and optical bandwidth W of the ASE source is larger by 2-3 orders of magnitude.

Covert sensing with a coherent state probe. Let us summarize the main results from <cit.>. The model is depicted in Fig. <ref> with η_eff=η and n̅_B_eff = n̅_B. Over n modes, Alice transmits a mean photon number of ⟨ N_S⟩=n n̅_S, where n̅_S is the mean transmitted photon number per mode. We make the assumption that the total variance, ⟨Δ N_S^2 ⟩=𝒪(n). Let us begin by stating the converse: Alice's sensing attempt is either detected by Willie with arbitrarily low detection probability, or the MSE of Alice's estimation is tightly lower bounded by 1/√(n), i.e.,  ⟨ (θ - θ̂_n)^2 ⟩ = Ω(1/√(n)) <cit.>. This notation implies that if Alice tries to get her MSE to surpass the 1/√(n) scaling (e.g., by increasing her transmit power), Willie will detect her sensing attempt with high probability. Now, we state the achievability: Suppose that Willie can perform any measurement permitted by physics on the 2n optical modes he collects. Alice can lower bound Willie's detection error probability ℙ_e^(det)≥ 1/2-ϵ, for any ϵ>0 while achieving the MSE ⟨ (θ - θ̂_n)^2 ⟩ = 𝒪(1/√(n)) using an n-mode coherent-state probe. For a coherent probe state ⊗_i=1^n|α_i⟩, where each complex amplitude α_i is drawn from the distribution p(α)=e^-|α|^2/n̅_S/πn̅_S, Willie's detection error probability is shown to satisfy <cit.>:

    ℙ_e^(det)≥1/2-(1-η)n̅_S√(n)/4√(ηn̅_B(1+η)).

Alice can choose

    n̅_S = 4ϵ√(ηn̅_B(1+ηn̅_B))/√(n)(1-η)

to ensure ϵ-covertness <cit.>. Using the aforesaid coherent state probe and ideal heterodyne detection measurement, Alice can achieve the MSE <cit.>:

    ⟨ (θ - θ̂_n)^2 ⟩≈c_het/ϵ√(n),

where

    c_het = (1-η)(1+n̅_B(1-η))/8 η√(ηn̅_B(1+ηn̅_B)).

Using a coherent state probe, but no restriction on the receiver, the MSE achieved is at best two times lower than that is achieved with a heterodyne receiver <cit.>. In other words, with Alice's optimal receiver,

    ⟨ (θ - θ̂_n)^2 ⟩≥c_coh/ϵ√(n),

where

    c_coh = (1-η)(1+2n̅_B(1-η))/16 η√(ηn̅_B(1+ηn̅_B)).


Covert sensing with ASE floodlight illumination. Let us first consider the achievability. Since all states and channels involved are Gaussian, we will work with the symplectic formalism. Further, all the first moments of the states involved are zero. The ASE source, after being split on an un-balanced beam-splitter, emits a two-mode zero-mean Gaussian state with covariance matrix (CM) given by <cit.>:

    V_ASE = 
    [  A_ASE 0_2× 2; 0_2× 2  A_ASE ],

where

    A_ASE = 
    [     n̅_S+1/2 √(n̅_Sn̅_LO); √(n̅_Sn̅_LO)    n̅_LO+1/2 ].

Here, n̅_S is the mean transmit photon number per mode, and n̅_LO≫n̅_S is the mean photon number per mode of Alice's local oscillator (LO) for her heterodyne detection receiver. If Alice is probing the target, Willie gets a state with the following CM (see App. Sec.<ref>),

    V_W =
    [   B_W   C_W; C_W^T   B_W ],

where

    B_W =
    [      w_11 -w_12cosθ; -w_12cosθ      w_22 ], and


    C_W =
    [         0  w_12sinθ; -w_12sinθ         0 ],

where w_11 = (1-η_2)η_1 n̅_S + (1-η_1)(1-η_2)n̅_B_1+η_2 n̅_B_2 + 1/2, 
w_12 = √((1-η_2)η_1(1-η_1)) (n̅_B_1-n̅_S), w_22 =  η_1 n̅_B_1 + (1-η_1) n̅_S + 1/2.
When Alice is not sensing, Willie's CM is given by (see App. Sec. <ref>) V_W^(0)=V_W(n̅_S=0), where w_ij^(0)=w_ij(n̅_S=0). 
We assume n̅_S ≪n̅_B_1 and n̅_S≪n̅_B_2, as this is intuitively needed for remaining covert. Willie must discriminate between the state ρ̂_1 with CM V_W (Alice is sensing) and the state ρ̂_0 with CM V_W^(0) (Alice is not sensing). After transmission of n modes over the probing interval, Willie's average probability of discrimination error is

    ℙ_e^(w)   ≥1/2[1 - 1/2ρ̂_0^⊗ n - ρ̂_1^⊗ n_1].

See also App. Sec. <ref> for a brief discussion on detectability.
The trace distance ρ̂_0-ρ̂_1_1 between states 
ρ̂_0 and ρ̂_1 is upper-bounded by the quantum relative entropy (QRE) D(ρ̂_0ρ̂_1)=tr(ρ̂_0lnρ̂_0 )-tr(ρ̂_0lnρ̂_1 ) using quantum Pinsker's inequality, i.e.,
ρ̂_0-ρ̂_1_1≤√(2D(ρ̂_0ρ̂_1))
which implies that

    ℙ_e^(w)   ≥1/2-√(1/8D(ρ̂_0^⊗ nρ̂_1^⊗ n)).

Thus, ensuring that D(ρ̂_0^⊗ nρ̂_1^⊗ n)≤ 8ϵ^2 ensures that ℙ_e^(w)≥1/2-ϵ over n channel uses. QRE is additive for tensor product states, therefore
D(ρ̂_0^⊗ nρ̂_1^⊗ n)=nD(ρ̂_0ρ̂_1).
By expanding the QRE in Taylor series with respect to n̅_S at n̅_S=0 and by using Taylor's theorem with the remainder we can upper bound D(ρ̂_0ρ̂_1) by its second order term c_2=f(η_1,η_2,n̅_B_1,n̅_B_2)>0. For the aforesaid analysis on the QRE, see App. Sec. <ref>. Consequently Willie's error probability will be lower bounded as,

    ℙ_e^(w)≥1/2-√(c_2)/4√(n) n̅_S.

Thus, if Alice sets

    n̅_S =4 ϵ/√(c_2)1/√(n),

she can ensure ϵ-covertness, i.e., ℙ_e^(w)≥1/2-ϵ. As shown in Fig. <ref>, Alice sees one effective thermal channel (see App. Sec. <ref>), with η_eff=η_1 η_2 and n̅_B_eff = [n̅_B_1(1-η_1)η_2+n̅_B_2(1-η_2)]/(1-η_1 η_2), which simplifies the analysis on Alice's measurement. We consider heterodyne detection and in App. Sec. <ref> we prove that the heterodyne output comprises of a pair of quadrature measurements with Gaussian distributions 𝒩(μ_i,σ_i^2), i=1,2, of means μ_i and variances σ_i^2. 𝒩(cosθ,[1+n̅_B_eff(1-η_eff)]/(2η_effn̅_S)+cos^2θ) and 𝒩(sinθ,[1+n̅_B_eff(1-η_eff)]/(2η_effn̅_S)+sin^2θ).
The variances σ_1,2^2 are taken in the limit n̅_LO→∞, i.e., a very bright local oscillator. By considering n modes, taking into account that independent Gaussian variables are additive, substituting n̅_S from Eq. (<ref>) and assuming that n ≫ 1, we get σ_1^2 ≈σ_2^2=σ^2_het (App. Sec. <ref>) where

    σ^2_het = c̃_het/ϵ√(n),

and c̃_het=[1+n̅_B_eff(1-η_eff)]√(c_2)/(8η_eff). Therefore,

    ⟨ (θ-θ̂_het,n)^2⟩≈c̃_het/ϵ√(n)=𝒪(1/√(n)).

For a justification of the approximation in Eq. (<ref>), refer to <cit.> (Appendix C therein) and App. Sec. <ref>. In the same limit n̅_LO→∞, the QFI for the parameter θ for Alice is given by (see App. Sec. <ref>):

    F_A = 4 n̅_Sη_eff/1+2n̅_B_eff(1-η_eff).

Therefore, by imposing the covertness condition Eq. (<ref>) we get the QCRB for the ASE source,

    ⟨ (θ-θ̂_n)^2⟩≥c_ASE/ϵ√(n)

where c_ASE = [1+2n̅_B_eff(1-η_eff)]√(c_2)/(16η_eff).
Therefore, with the ASE source, from Eqs. (<ref>) and (<ref>), heterodyne detection yields an MSE that is at most twice compared to what is attainable by the optimal quantum receiver.

Let us discuss briefly the converse, i.e, the statement that regardless of the n-mode probe state employed by Alice, the sensing attempt is either detected with arbitrarily low detection probability, or the MSE of Alice's estimator is tightly lower bounded as  ⟨ (θ - θ̂_n)^2 ⟩ = Ω(1/√(n)). Since for the converse a specific probe state is not fixed and having proved that Alice sees only one effective thermal channel with (η_eff, n̅_B_eff), the converse proved in <cit.> using the methods and results from <cit.> is still valid.

Laser light vs. ASE source probe. Let us compare the quantum limits corresponding to the ASE and the coherent state sources. In <cit.>, it was found that for a coherent state probe, ⟨ (θ-θ̂_n)^2⟩≥ F_coh^-1 =  c_coh/(ϵ√(n)), where c_coh=(1-η_eff)[1+2n̅_B_eff(1-η_eff)]/[16η_eff√(η_effn̅_B_eff(1+η_effn̅_B_eff))] and F_coh is the QFI for coherent state probe and for the effective single thermal channel Alice sees. When the ratio μ=F_A^-1/F_coh^-1 is μ<1, the ASE source can outperform the coherent state probe. For n=WT channel uses, where the optical bandwidth W depends on the source and probing time T is fixed, the ratio is μ=μ_c/√(μ_w), where μ_c=c_ASE/c_coh and μ_w=W_ASE/W_coh. 
If n̅_B_1=n̅_B_2, as expected in most situations since thermal equilibrium across the forward and return paths is not an onerous assumption, we get μ_c=1 which is its minimum value (App. Sec. <ref>), i.e., the per-mode performance for the two sources is identical. Therefore, we have μ=√(W_coh/W_ASE)<1 since the bandwidth of a typical laser source is in the GHz regime while that of an ASE source is in the THz regime. We thus conclude that the ASE source would outperform a coherent state probe because of its larger bandwidth, under identical thermal environments.

Performance evaluation. Let us consider sensing over a mono-static free-space channel of line-of-sight range L, with an ASE source of center wavelength λ. We employ a single spatial mode (the fundamental Gaussian mode) probe. Assuming vacuum propagation (ignoring turbulence and atmospheric extinction), and Gaussian-attenuation apertures (for simplifying the expressions without losing the essence of the problem), we get the transmissivities of both the forward and return paths as η_1=η_2=η = (A_t A_T)/ (λ L)^2, where A_t = π r_t^2 is the area of the radius-r_t exit pupil of Alice's transmitter aperture, which we take to be equal to the area of the entrance pupil of the aperture of her receiver telescope. A_T = π r_T^2 is the effective (radius-r_T circular) area of the target cross-section interrogated. The cw ASE source of optical bandwidth W illuminates the target for T seconds, which results in n ≈ WT temporal modes to probe the target.

In Fig. <ref>, we plot c_ASE of Eq. (<ref>) as a function of the center frequency f = hc/λ of the source, where h is Planck constant and c is the speed of light. We assume Planck law limited thermal noise, with n̅_B_1 = n̅_B_2 = n̅_B = [exp(hc/(λ k_B T_0))-1]^-1, where k_B is the Boltzmann constant and T_0 = 300K is the ambient temperature. Longer wavelengths (smaller frequencies) have more noise (larger n̅_B) to hide the probe, which tends to make c_ ASE smaller. On the other hand, longer wavelengths (smaller frequencies) have more diffraction limited loss (smaller η), giving more photons to Willie, which tends to make c_ ASE higher. So, there is a trade-off in the choice of the center frequency, and we see that the minimum of c_ ASE happens in the mid to long wave infrared (ML-IR) region. We numerically evaluate the optimal wavelength λ which minimizes c_ASE. For ϵ=10^-3, probing time T=1s, and for a source with bandwidth W=3THz we find the right hand side of Eq. (<ref>), which we denote as B, the MSE lower bound. For L=1km: λ=9.40 μm, B=0.00322, for L=3km: λ=6.35μm, B=0.09927, and for L=5km: λ=5.38μm, B=0.81438.       

Sources in the ML-IR region with W=3THz and wavelength 3μm≤λ≤ 4.4 μm and λ=4,4.6,8.7,10.2μm such as <cit.><cit.>, for L=1km and T=1s give B=0.00327 and B=0.08423 for wavelengths 8.7μm and 3μm respectively, which correspond to the extrema of the performance for the aforementioned parameters, wavelength range and values.



Conclusions. Our work introduces to the greater communities of optical sensing, strategic security, and fundamental quantum optics, the concept of quantum-secure covert (low probability of detect) sensing. Our analytical analysis has shown that using floodlight illumination from an ASE source, covert phase sensing is possible and experimentally feasible. In fact, the ASE source enhances covert sensing performance by exploiting its large optical bandwidth W, i.e., the number of available orthogonal modes n=WT can significantly increase for a given time T. An open question is to consider sensing multiple parameters, e.g., covert sensing of phase and loss (ranging) or multiple phases that could be the pixels of an image. We leave these questions for future work.
Our setup is experimentally feasible. The ASE source is divided into two arms by a highly reflective fiber coupler, so that  n̅_S≪ 1 and n̅_LO≫ 1. The sensing arm is frequency shifted prior to being sent to a very noisy environment and then is phase shifted by a weakly reflecting object. The LO is stored locally in a spool of optical fibers for later heterodyne measurement. Note that for sources with smaller wavelengths than those in <cit.>, one could inject artificial noise to get a reasonable MSE lower bound for a proof of concept experiment. Our present work indicates covert sensing is experimentally feasible and applicable to practical optical sensing tasks. A thorough analysis of modeling, experimental demonstration, and application development, are left for follow up works.


Acknowledgments. CG was supported on a General Dynamics research and development contract to the University of Arizona. SG, BB and ZZ acknowledge the Office of Naval Research program Communications and Networking with Quantum Operationally-Secure Technology for Maritime Deployment (CONQUEST), awarded under Raytheon BBN Technologies prime contract number N00014-16-C-2069, and subcontract to University of Arizona. AD was partly supported by the UK EPSRC (EP/K04057X/2), and the UK National Quantum Technologies Programme (EP/M013243/1, EP/M01326X/1). The authors thank Darlene Hart (General Dynamics), William Clark (General Dynamics), and Bob Shine (DRS Daylight Solutions) for useful discussions on light sources.







§ APPENDIX


Aequation


We give information on the derivations used in the main paper. Specifically, we discuss detectability, we derive the covariance matrices for the global output, for Alice, and for Willie. We calculate the quantum relative entropy (QRE) for Willie's part between the states which correspond to Alice is sensing and Alice is not sensing scenarios and we discuss the Taylor expansion of the QRE. We prove that Alice sees only one effective thermal channel. We give explicit calculations for the heterodyne detection. Finally, we derive the quantum Fisher information for when Alice is sensing the unknown phase. 



 §.§ Detectability

The adversary performs a binary hypothesis test on his sample to determine
whether the target is being interrogated or not.
Performance of the hypothesis test is typically measured by its
detection error probability denoted as ℙ_ e^ (det) here since we speak generally, but it is ℙ_ e^ (w) in the main paper when we talk about the ASE source. We have, 
ℙ_ e^ (det)=ℙ_ FA+ℙ_ MD/2,
where equal prior probabilities on sensor's interrogation state are
assumed, ℙ_ FA is the probability of false alarm and
ℙ_ MD is the probability of missed detection.
The sensor desires to remain covert by
ensuring that ℙ_ e^ (det)≥1/2-ϵ for
an arbitrary small ϵ>0 regardless of adversary's measurement
choice (since ℙ_ e^ (det)=1/2 for a random 
guess).
By decreasing the power used in a probe, the sensor can decrease
the effectiveness of the adversary's hypothesis test at the expense of
the increased mean squared error (MSE) of the estimate.




 §.§ Alice's and Willie's covariance matrices

We give the full expressions of: the final covariance matrix (CM) of the global state, i.e, the final state that contains both Alice's and Willie's modes (4 modes in total, therefore the dimensions of the CM is 8× 8), and the CM of each one Alice and Willie. We remind the reader that the first moments are zero. We also provide a guide on how these CM were derived.

We use the symplectic formalism <cit.>, where under a symplectic matrix S, a CM V_0 is evolved to the CM V=SV_0S^T. All CM are in the qqpp representation, meaning that for a 2N × 2N CM the upper left block contains information only on position, the lower right block concerns only momentum, and the off-diagonal blocks contain information on position-momentum correlations. The phase space vector is r=(r_1,… r_2N)^T=(q_1,… q_N, p_1 … p_N )^T. In our set up, we assume that the first two modes belong to Willie, while the last two modes belong to Alice. We have worked with ħ=1, therefore the CM for vacuum is I/2, where I is the identity matrix.     

Since initially Alice's and Willie's systems are uncorrelated, in the qqpp representation the initial, global CM is

    V_AW_0=
    	[   n̅_B_2+1/2            0            0            0            0            0            0            0;            0   n̅_B_1+1/2            0            0            0            0            0            0;            0            0     n̅_S+1/2 √(n̅_Sn̅_LO)            0            0            0            0;            0            0 √(n̅_Sn̅_LO)  n̅_LO + 1/2            0            0            0            0;            0            0            0            0   n̅_B_2+1/2            0            0            0;            0            0            0            0            0   n̅_B_1+1/2            0            0;            0            0            0            0            0            0     n̅_S+1/2 √(n̅_Sn̅_LO);            0            0            0            0            0            0 √(n̅_Sn̅_LO)  n̅_LO + 1/2 ]




Our setup consists of three symplectic transformations, which are depicted in Fig. <ref> and in order of action are the following. A beam splitter B_23(η_1) with transmissivity η_1 acting on modes 2 and 3, a phase shift R_3(θ) acting on mode 3, and a beam splitter B_13(η_2) with transmissivity η_2 acting on modes 1 and 3. By noticing on which modes each symplectic matrix acts, to find a the final CM becomes a straightforward job. Let us demonstrate an example. The beam splitter B_23(η_1) will act on the sub-matrices V^(pos_1)_AW_0 and V^(mom_1)_AW_0 of V_AW_0 with elements V_AW_0^ij, i,j=2,3 and V_AW_0^ij, i,j=6,7 respectively, i.e., to the momentum and position coordinates with indices  2,3. The symplectic matrix B_23(η_1) is,

    B_23(η_1) =
    	[ B_23^(pos)(η_1)          0_2× 2;          0_2× 2 B_23^(mom)(η_1) ],

where

    B_23^(pos)(η_1) = B_23^(mom)(η_1) =
    	[    √(η_1)  √(1-η_1); -√(1-η_1)    √(η_1) ].

Therefore we apply the transformation V_1=B_23(η_1) ( V^(pos_1)_AW_0⊕ V^(mom_1)_AW_0) B_23^T(η_1) where all involved matrices are block diagonal and therefore the resulting matrix will also be block-diagonal. Then in V_AW_0, we substitute the V^(pos_1)_AW_0 and V^(mom_1)_AW_0 blocks with the resulting blocks of V_1. The rest of V_AW_0 elements remain unchanged.

Subsequently, we apply the phase shift which corresponds to the symplectic matrix,

    R_3(θ)=
    	[     1     0     0     0     0     0     0     0;     0     1     0     0     0     0     0     0;     0     0  cosθ     0     0     0 -sinθ     0;     0     0     0     1     0     0     0     0;     0     0     0     0     1     0     0     0;     0     0     0     0     0     1     0     0;     0     0  sinθ     0     0     0  cosθ     0;     0     0     0     0     0     0     0     1 ]

or one can again simplify calculations by noticing which modes the phase shift affects. Finally, the beam splitter B_13(η_2) is applied on the resulting CM. The final, global CM is

    V=[   B   C; C^T   B ]

where,

    B=
    	[  B_W B_AW; B_AW  B_A ]

and

    C=
    	[     C_W    C_AW; -C_AW^T     C_A ]

with

    B_W=[ (1-η_2)η_1 n̅_S + (1-η_1)(1-η_2)n̅_B_1+η_2 n̅_B_2 + 1/2                 √((1-η_2)η_1(1-η_1)) (n̅_B_1-n̅_S) cosθ;                 √((1-η_2)η_1(1-η_1)) (n̅_B_1-n̅_S) cosθ                         η_1 n̅_B_1 + (1-η_1) n̅_S + 1/2 ],


    B_AW=[ √(η_2(1-η_2))[η_1 n̅_S+(1-η_1)n̅_B_1-n̅_B_2]                  √((1-η_2)η_1 n̅_Sn̅_LO)cosθ;          √(η_2 η_1 (1-η_1))(n̅_B_1-n̅_S)cosθ                          √((1-η_1)n̅_S N_LO) ],


    B_A=[ (1-η_1)η_2 n̅_B_1 + (1-η_2)n̅_B_2+η_1η_2 n̅_S +1/2                            √(η_1η_2 n̅_Sn̅_LO)cosθ;                            √(η_1η_2 n̅_Sn̅_LO)cosθ                                          n̅_LO+1/2 ],


    C_W=[                                        0 -√((1-η_2)η_1(1-η_1)) (n̅_B_1-n̅_S) sinθ;  √((1-η_2)η_1(1-η_1)) (n̅_B_1-n̅_S) sinθ                                        0 ],


    C_AW=[                                   0        -√((1-η_2)η_1 n̅_Sn̅_LO)sinθ; √(η_2η_1(1-η_1)) (n̅_B_1-n̅_S) sinθ                                   0 ],


    C_A=[                        0 -√(η_2η_1 n̅_Sn̅_LO)sinθ;  √(η_2η_1 n̅_Sn̅_LO)sinθ                        0 ].


By tracing out, i.e, by keeping the corresponding sub-matrices, we can find the covariance matrices V_A for Alice and V_W for Willie,

    V_A=
    	[   B_A   C_A; C_A^T   B_A ]=
    	[      a_11 -a_12cosθ         0  a_12sinθ; -a_12cosθ      a_22 -a_12sinθ         0;         0 -a_12sinθ      a_11 -a_12cosθ;  a_12sinθ         0 -a_12cosθ      a_22 ],


    V_W=
    	[   B_W   C_W; C_W^T   B_W ]=
    	[      w_11 -w_12cosθ         0  w_12sinθ; -w_12cosθ      w_22 -w_12sinθ         0;         0 -w_12sinθ      w_11 -w_12cosθ;  w_12sinθ         0 -w_12cosθ      w_22 ],

where

    a_11   =    (1-η_1)η_2 n̅_B_1 + (1-η_2)n̅_B_2+η_1η_2 n̅_S +1/2,
     a_12   =    -√(η_1η_2 n̅_Sn̅_LO)cosθ,
     a_22   =   n̅_LO+1/2,
    
    	w_11   =    (1-η_2)η_1 n̅_S + (1-η_1)(1-η_2)n̅_B_1+η_2 n̅_B_2 +1/2,
    
    	w_12   =   √((1-η_2)η_1(1-η_1)) (n̅_B_1-n̅_S),
    
    	w_22   =   η_1 n̅_B_1 + (1-η_1) n̅_S + 1/2.
  




 §.§ The quantum relative entropy

To calculate the quantum relative entropy (QRE) D(ρ̂_0ρ̂_1)=tr(ρ̂_0lnρ̂_0 )-tr(ρ̂_0lnρ̂_1 ) between two Gaussian states ρ̂_0 and ρ̂_1, we need to find the symplectic eigenvalues and symplectic eigenvectors of one of the the two covariance matrices and the symplectic eigenvalues of the other <cit.> (section Methods therein). Let ρ̂_0, ρ̂_1 to have a covariance matrix V_0, V_1 respectively. The covariance matrix V_1 is for the case where Alice is sensing, i.e, (n̅_S≠ 0), while the covariance matrix V_0 is for the case where Alice is not sensing, i.e., (n̅_S = 0), therefore V_0=V_W(n̅_S=0) and V_1=V_W(n̅_S≠ 0), where V_W is Willie's covariance matrix given in Eq. (<ref>). That means Willie's CM matrix elements for n̅_S=0, i.e, the matrix elements of V_0, are

    w_11^(0)   =    (1-η_1)(1-η_2)n̅_B_1+η_2 n̅_B_2 +1/2,
     w_12^(0)   =   √((1-η_2)η_1(1-η_1))n̅_B_1,
     w_22^(0)   =   η_1 n̅_B_1 + 1/2

and V_0 will have the same structure as V_W in Eq. (<ref>). The QRE is,

    D(ρ̂_0ρ̂_1) = 	-Σ(V_0,V_0)+Σ(V_0,V_1).

In Eq. (<ref>), the first term is the von Neumann entropy (with a minus in front) of a Gaussian state with CM V_0 and zero first moments, which for our case is,

    Σ(V_0,V_0) = 1/2∑_k=1^2[ (1+2u^(0)_k ) ln(u^(0)_k+1/2) + (1-2u^(0)_k ) ln(u^(0)_k-1/2) ],

where u^(0)_k, k=1,2 are the symplectic eigenvalues of V_0 given by

    u^(0)_1    =   1/2(w_11^(0)+w_22^(0) +√(4 w_12^(0)2+(w_11^(0)-w_22^(0))^2))
    	u^(0)_2    =   1/2(w_11^(0)+w_22^(0) -√(4 w_12^(0)2+(w_11^(0)-w_22^(0))^2)).

The second term in Eq. (<ref>), for the case at hand is given by

    Σ(V_0,V_1) = 1/2∑_k=1^2[ (1+2d_k ) ln(u_k+1/2) + (1-2d_k ) ln(u_k-1/2) ],

where u_k are the symplectic eigenvalues of V_1 and d_k are the diagonal elements (which are doubly degenerate in our case, just as the symplectic spectrum of a CM) of V_0'=M^T V_0 M, where M is the symplectic eigenvectors' matrix of V_1.
The symplectic eigenvalues of V_1 are found to be

    u_1    =   1/2(w_11+w_22 +√(4 w_12^2+(w_11-w_22)^2))
    	u_2    =   1/2(w_11+w_22 -√(4 w_12^2+(w_11-w_22)^2))

and the symplectic eigenvectors' matrix of V_1 is a beam splitter followed by a 2-mode phase rotation, 

    M= R(θ/2+π,-θ/2)B(τ)=
    	[   -√(τ)cosθ/2  √(1-τ)sinθ/2    √(τ)sinθ/2 -√(1-τ)sinθ/2;  √(1-τ)cosθ/2    √(τ)cosθ/2  √(1-τ)sinθ/2    √(τ)sinθ/2;   -√(τ)sinθ/2  √(1-τ)sinθ/2   -√(τ)cosθ/2  √(1-τ)cosθ/2; -√(1-τ)sinθ/2   -√(τ)sinθ/2  √(1-τ)cosθ/2    √(τ)cosθ/2 ],

where

    τ = w_11-w_22/2√(4 w_12^2+(w_11-w_22)^2)+1/2.

The matrix M can be easily checked to satisfy MV_1M^T=diag(u_0,u_1,u_0,u_1) and MΩ M^T=Ω, 
where,

    Ω=
    	[  0_2× 2  I_2× 2; -I_2× 2  0_2× 2 ],

meaning that M is a valid symplectic matrix that diagonalizes (in the symplectic sense) the matrix V_1.  
The diagonal elements of V_0' are,

    d_1   =   w_11^(0)+w_22^(0)/2+1/2√(4 w_12^2+(w_11-w_22)^2)(4w_12w_12^(0)+(w_11-w_22)(w_11^(0)-w_22^(0)))
    	d_2   =   w_11^(0)+w_22^(0)/2-1/2√(4 w_12^2+(w_11-w_22)^2)(4w_12w_12^(0)+(w_11-w_22)(w_11^(0)-w_22^(0))).
 
Given the expressions (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), the relative entropy is now known as a function of the V_0 and V_1 matrix elements through Eqs. (<ref>), (<ref>), and (<ref>) (see also Mathematica file QRE.nb).




 §.§ Taylor expansion of the relative entropy

We expand the QRE (<ref>) in Taylor series with respect to n̅_S at n̅_S=0

    D(ρ̂_0ρ̂_1) ≡ D(n̅_S) = D(0) + ∂ D(n̅_S)/∂n̅_S|_n̅_S=0 +1/2!∂^2 D(n̅_S)/∂n̅_S^2|_n̅_S=0 + 1/3!∂^3 D(n̅_S)/∂n̅_S^3|_n̅_S=0+…

The 0^th term in (<ref>) is zero because D(0)=D(ρ̂_0ρ̂_0) = 0. The 1^st order term is zero, as the QRE has a minimum at D(ρ̂_0ρ̂_1) = 0, and the 2^nd order term is non-negative as the QRE D(ρ̂_0ρ̂_1) ≥ 0, i.e., it has a global minimum at zero. Note that Σ(V_0,V_0) of the QRE (<ref>) has not a dependence on n̅_S, therefore one is concerned only with the term Σ(V_0,V_1) when taking derivatives on the QRE with respect to n̅_S.

Having expressed the QRE as a function of the CM elements and consequently as a function of (η_1,η_2,n̅_B_1,n̅_B_2), it is straightforward to derive the 2^nd and 3^rd order derivatives. Below, we give a guide on how to derive that the 2^nd derivative (at n̅_S=0) is non-negative, that μ_c ≥ 1, and that the 3^rd derivative (at n̅_S=0) is negative. 
The main tool is to upper and lower bound the ln(.) <cit.> with

    2x/x+2≤ln(1+x) ≤x/2x+2/x+1, x≥ 0.
  
The expression of the 2^nd derivative has the form

    c_2=∂^2 D(n̅_S)/∂n̅_S^2|_n̅_S=0= T_1+T_2 ln(1+x),
 
where T_1, T_2>0, and x>0 are functions of (η_1,η_2,n̅_B_1,n̅_B_2). Since T_2>0, we use the lower bound of (<ref>), c_2 ≥ T_1+T_22x/x+2, from which expression it can shown that c_2≥ 0.

To prove that μ_c=c_ASE/c_coh≥ 1 ⇒ c_ASE-c_coh≥ 0, where c_ASE = [1+2n̅_B_eff(1-η_eff)]√(c_2)/(16η_eff), we use a similar technique. We eliminate the ln(.) in c_2 using again the lower bound of (<ref>, a technique which gives a sufficient, but not necessary, condition. However, we get indeed μ_c ≥ 1.

The 3^rd derivative of the QRE has the form,

    c_3=∂^3 D(n̅_S)/∂n̅_S^3|_n̅_S=0= P_1+P_2 +P_3ln(1+x),
 
where P_1, P_2, P_3, and x>0 are functions of (η_1,η_2,n̅_B_1,n̅_B_2). In this case, P_3 can be positive or negative. Therefore, one has to take cases when eliminating the ln(.), which again gives a sufficient, but not necessary, condition. When P_3>0, we use the upper bound of (<ref>), c_3 ≤ P_1+P_2+P_3x/2x+2/x+1, while when P_3<0, we use the lower bound of (<ref>), c_3 ≤ P_1+P_2+P_32x/x+2. For both cases it can shown that c_3≤ 0.

The case where n̅_B_1=n̅_B_2=n̅_B yields simple results,

    D(ρ̂_0ρ̂_1)_(n̅_B)   =   η_effn̅_Bln[η_effn̅_B/1+η_effn̅_B1+(1-η_eff)n̅_S+η_effn̅_B/(1-η_eff)n̅_S+η_effn̅_B]+ln[1+(1-η_eff)n̅_S+η_effn̅_B/(1-η_eff)n̅_S+η_effn̅_B],
    
    	c_2 (n̅_B)   =   (1-η_eff)^2/η_effn̅_B(1+η_effn̅_B)>0,
    λ_c    =   1,
    
    	c_3 (n̅_B)   =    -2(1-η_eff)^3(1+2η_effn̅_B)/η_eff^2n̅_B^2(1+η_effn̅_B)^2<0.
  
Details can be found in Mathematica file QRE.nb.



 §.§ Effective channel

A covariance matrix V_1 under a thermal channel C(n̅_B,η) <cit.> is transformed as

    V_2 = X V_1 X^T +Y,
   
where,

    X   =   √(η) I
    	Y   =   (1-η)(n̅_B-1/2) I.

Alice's sensing arm goes through a thermal channel C_1(n̅_B_1,η_1), a phase shift R(θ) which is a symplectic orthogonal matrix, and another thermal channel C_2(n̅_B_2,η_2). Therefore, the sensing arm will see the transformation,

    X̃   =    X_2 R(θ) X_1
    Ỹ   =    X_2 R(θ) Y_1 R^T(θ) X_2^T + Y_2,

Where X_1, Y_1 and X_2, Y_2 correspond to C_1 and C_2 respectively. Using Eqs. (<ref>), (<ref>), and since the matrices X_1, X_2, Y_1, Y_2 are proportional to the identity matrix I they commute with R(θ), Eqs. (<ref>) and (<ref>) give,

    X̃   =   √(η_eff) I
    Ỹ   =    (1-η_eff) (n̅_B_eff+1/2)I,

where,

    η_eff   =   η_1η_2
    n̅_B_eff   =   (1-η_1)η_2 n̅_B_1+(1-η_2)n̅_B_2/1-η_1η_2.

Therefore, Alice's sensing arm sees only one effective thermal channel with transmittance η_eff, mean thermal photon number n̅_B_eff, and a phase shift R(θ). The phase shift R(θ) can be put before or after the thermal channel as it commutes with X̃ and Ỹ and R(θ) is orthogonal. The matrix elements a_11 and a_12 in Alice's CM can now be written,

    a_11   =    (1-η_eff) n̅_B_eff + η_effn̅_S +1/2,
     a_12   =    -√(η_effn̅_Sn̅_LO)cosθ.
 





 §.§ Heterodyne detection

We will consider dual homodyne detection on Alice's system as in Fig. <ref>, which is an implementation of heterodyne detection. We want to calculate the mean values μ_13, μ_42 and the variances σ^2_13, σ^2_42 of the Gaussian distributions 𝒩(μ_13,σ^2_13), 𝒩(μ_42,σ^2_42) respectively. The aforementioned Gaussian distributions correspond to the observables n̂_1-n̂_3 and n̂_4-n̂_2 in the following manner, 

    μ_13   =   ⟨n̂_1-n̂_3 ⟩
    μ_42   =   ⟨n̂_4-n̂_2 ⟩
    σ^2_13   =   ⟨ (n̂_1-n̂_3)^2 ⟩ = ⟨n̂_1^2 ⟩ - ⟨n̂_1 ⟩^2 + ⟨n̂_3^2 ⟩ - ⟨n̂_3 ⟩^2 - 2 (⟨n̂_1n̂_3⟩-⟨n̂_1⟩⟨n̂_3⟩) 
    σ^2_42   =   ⟨ (n̂_4-n̂_4)^2 ⟩ = ⟨n̂_4^2 ⟩ - ⟨n̂_4 ⟩^2 + ⟨n̂_2^2 ⟩ - ⟨n̂_2 ⟩^2 - 2 (⟨n̂_4n̂_2⟩-⟨n̂_4⟩⟨n̂_2⟩).

Therefore, we need to calculate the expectation values in Eqs. (<ref>), (<ref>), (<ref>), and (<ref>). To this end, we calculate the evolution of the CM for the setup of Fig. <ref>. The CM V_in just before the heterodyne detection is given by the Eq. (<ref>) and by taking into account that the CM for vacuum is I/2,

    V_in=
    	[       1/2         0         0         0         0         0         0         0;         0      a_11 -a_12cosθ         0         0         0  a_12sinθ         0;         0 -a_12cosθ      a_22         0         0 -a_12sinθ         0         0;         0         0         0       1/2         0         0         0         0;         0         0         0         0       1/2         0         0         0;         0         0 -a_12sinθ         0         0      a_11 -a_12cosθ         0;         0  a_12sinθ         0         0         0 -a_12cosθ      a_22         0;         0         0         0         0         0         0         0       1/2 ].

After the application of all symplectic matrices as described in Fig. <ref>, the CM just before the photon number resolution detectors is,

    V_het   =   [   A   B; B^T   A;     ],
    
    	A   =   [ A_11 A_12 A_13 A_14; A_12 A_22 A_23 A_24; A_13 A_23 A_33 A_34; A_14 A_24 A_34 A_44 ],


    B   =   [     0  B_12  B_13  B_14; -B_12     0  B_23  B_24; -B_13 -B_23     0  B_34; -B_14 -B_24 -B_34     0 ],


    A_11   =   1/4 (1 + a_11 + a_22 - 2 a_12cosθ)
    
    	A_12   =   1/8 [1 - 2 a_11 + 2 a_12(cosθ - sinθ)]
    
    	A_13   =   a_11-a_22/4
    
    	A_14   =   1/8 [1 - 2 a_11 + 2 a_12(cosθ + sinθ)]
    
    	A_22   =   1/4 (1 + a_11 + a_22 + 2 a_12sinθ)
    
    	A_23   =   1/8 [1 - 2 a_11 - 2 a_12 (cosθ + sinθ)]
    
    	A_24   =   a_11-a_22/4
    
    	A_33   =   1/4 (1 + a_11 + a_22 + 2 a_12cosθ)
    
    	A_34   =   1/8 [1 - 2 a_11 - 2 a_12( cosθ - sinθ)] 
    
    	A_44   =   1/4 (1 + a_11 + a_22 - 2 a_12sinθ),


    B_12   =    -1/8 [1 - 2 a_22 + 2 a_12 (cosθ - sinθ)]
    
    	B_13   =    -1/2 a_12sinθ
    
    	B_14   =   1/8 [1 - 2 a_22 + 2 a_12 (cosθ + sinθ)]
    
    	B_23   =    -1/8 [1 - 2 a_22 - 2 a_12 (cosθ + sinθ)]
    
    	B_24   =    -1/2 a_12cosθ
    
    	B_34   =    -1/8 [1 - 2 a_22 - 2 a_12 (cosθ - sinθ)].


Using the the CM of Eq. (<ref>), we can construct the Wigner function,

    W(r) = 1/(2π)^4√( V_het)∫ d^2rexp(-1/2r^T V_het^-1r)

where r=(q_1,…,q_4,p_1,…,p_4) is the phase space vector (q_i,p_i refer to position and momentum respectively of the i-th mode, also note that we work with ħ=1), and d^2r = dq_1… dq_4 dp_1… dp_4. Utilizing the Wigner function we can calculate mean values in the symmetric ordering,

    ⟨ : n̂_i n̂_j : ⟩ = ∫ d^2r W(r) (q_i^2+p_i^2)(q_j^2+p_j^2)/4.

Utilizing Wick-Isserlis theorem and the commutation relation [â_i,â_j^†]=δ_ij, for example, 

    ⟨ : n̂_i n̂_j : ⟩ = ⟨n̂_i n̂_j  ⟩+1/2⟨n̂_i⟩+1/2⟨n̂_j⟩-1/4,

we find

    ⟨n̂_1⟩   =   1/4 (a_11 + a_22 - 2 a_12cosθ-1)
    ⟨n̂_2⟩   =   1/4 (a_11 + a_22 + 2 a_12sinθ-1)
    ⟨n̂_3⟩   =   1/4 (a_11 + a_22 + 2 a_12cosθ-1)
    ⟨n̂_4⟩   =   1/4 (a_11 + a_22 - 2 a_12sinθ-1)
    ⟨n̂_1^2 ⟩   =   1/8 (a_11 + a_22 - 2 a_12cosθ-1) ( a_11 + a_22 - 2 a_12cosθ+1)
    ⟨n̂_2^2 ⟩   =   1/8 (a_11 + a_22 + 2 a_12sinθ-1) ( a_11 + a_22 + 2 a_12sinθ+1)
    ⟨n̂_3^2 ⟩   =   1/8 (a_11 + a_22 + 2 a_12cosθ-1) ( a_11 + a_22 + 2 a_12cosθ+1) 
    ⟨n̂_4^2 ⟩   =   1/8 (a_11 + a_22 - 2 a_12sinθ-1) ( a_11 + a_22 - 2 a_12sinθ+1)
    ⟨n̂_1 n̂_3 ⟩   =   1/16 [1 + 2 (a_11-1) a_11 + 2 (a_22-1) a_22 - 4 a_12^2 cos 2θ] 
    ⟨n̂_4 n̂_2 ⟩   =   1/16 [1 + 2 (a_11-1) a_11 + 2 (a_22-1) a_22 + 4 a_12^2 cos 2θ].


From Eqs. (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) to (<ref>), we find

    μ_13   =   √(η_effn̅_Sn̅_LO)cosθ
    μ_42   =   √(η_effn̅_Sn̅_LO)sinθ
    σ_13^2    =   1/2 [n̅_LO + n̅_eff (1 - η_eff) + n̅_eff (1 - η_eff) n̅_LO + n̅_Sη_eff] + n̅_LOn̅_Sη_effcos^2θ
    σ_42^2    =   1/2 [n̅_LO + n̅_eff (1 - η_eff) + n̅_eff (1 - η_eff) n̅_LO + n̅_Sη_eff] + n̅_LOn̅_Sη_effsin^2θ.

We normalize the mean values with k=√(η_effn̅_Sn̅_LO) and consequently the variances with k^2, and by taking the limit of a very strong local oscillator we get,

    μ_1   =   μ_13/k = cosθ
    μ_2   =   μ_42/k = sinθ
    σ_1^2    =   lim_n̅_LO→∞σ_13/k^2 =σ^2+cos^2θ
    σ_2^2    =   lim_n̅_LO→∞σ_42/k^2 =σ^2+sin^2θ,

where

    σ^2 = 1+(1-η_eff)n̅_B_eff/2 η_effn̅_S.

Alice is trying to determine the mean values of the distributions 𝒩(cosθ, σ^2+cos^2θ) and 𝒩(sinθ, σ^2+sin^2θ) by using the channel n times, i.e, by performing n measurements. Then the distributions of the mean values will be  𝒩(cosθ, σ^2/n+cos^2θ/n) and 𝒩(sinθ, σ^2/n+sin^2θ/n), were we have used the fact that independent Gaussian variables are additive. We set,

    σ^2_het=σ^2/n=1+(1-η_eff)n̅_B_eff/2 n η_effn̅_S

and by imposing the covertness condition n̅_S =(4 ϵ)/[√(c_2)√(n)] we get,

    σ^2_het = c̃_het/ϵ√(n),

where c̃_het=[1+n̅_B_eff(1-η_eff)]√(c_2)/(8η_eff).
Therefore the two distributions now are  𝒩(cosθ,  c̃_het/(ϵ√(n))+cos^2θ/n) and 𝒩(sinθ,  c̃_het/(ϵ√(n))+sin^2θ/n).
In the limit of n ≫ 1 (a limit which is justified by the large bandwidth of the ASE source), since sin^2θ≤ 1 and cos^2θ≤ 1 the term that scales as 1/n will be insignificant, therefore the approach in Sec.<ref> and in <cit.> (Appendix C therein) applies here as well and we conclude that the mean square error of the estimator θ̂_het,n=arctan((sinθ+σ^2_het)/(cosθ+σ^2_het)) scales as ⟨ (θ-θ̂_het,n)^2⟩ =𝒪(1/√(n)).



 §.§ Proof that ⟨ (θ - θ̂_ het)^2 ⟩≤σ^2_ het+𝒪(1/n) 

Consider the following estimator for θ,

    θ̂_ het   =   tan^-1(sin(θ)+1/n∑_i=1^n Z^ (Q,N)_i/cos(θ)+1/n∑_i=1^n Z^ (I,N)_i)
       =   tan^-1(sin(θ)+Z^ (Q)/cos(θ)+Z^ (I)),

where the random variables representing in-phase (I) and quadrature (Q) noise Z^ (I), Z^ (Q) respectively, obey the Gaussian distributions Z^ (I)∼𝒩(0,σ^2_ het) and Z^ (Q)∼𝒩(0,σ^2_ het), Z^ (I,N)∼𝒩(0,σ^2+cos^2θ) and Z^ (Q,N)∼𝒩(0,σ^2+sin^2θ). The distributions are indeed Gaussian for n≫ 1 and by applying the central limit theorem. The variance σ^2_ het is given in Eq. (<ref>). The MSE is:

    ⟨ (θ - θ̂_ het)^2 ⟩   =⟨(θ-tan^-1(sin(θ)+Z^ (Q)/cos(θ)+Z^ (I)))^2 ⟩
       =⟨(θ-tan^-1(sin(θ)+Rcos(φ)/cos(θ)+Rsin(φ)))^2 ⟩

where in (<ref>) we use circular symmetry of the two-dimensional added white Gaussian noise (AWGN) to change from the rectangular to polar coordinate system.
Thus, the radius is distributed as a Rayleigh random variable
R∼Rayleigh(σ^2_ het)
while the angle is distributed uniformly φ∼𝒰([0,2π]).
Now, the Taylor series expansion of
tan^-1(sin(θ)+rcos(φ)/cos(θ)+rsin(φ)) around r=0 is:

    tan^-1(sin(θ)+rcos(φ)/cos(θ)+rsin(φ))   =θ+rcos(θ+φ)-r^2/2sin(2(θ+φ))-r^3/3cos(3(θ+φ))+r^4/4sin(4(θ+φ))
       =+r^5/5cos(5(θ+φ))-r^6/6sin(6(θ+φ))-r^7/7cos(7(θ+φ))+r^8/8sin(8(θ+φ))
       =+…
       ≤θ+rcos(θ+φ)+∑_i=2^∞r^i/i
       =θ+rcos(θ+φ)-(log(1-r)+r) provided 0≤ r <1,

where the upper bound in (<ref>) is because 
sin(x),cos(x)∈[-1,1].
While this demonstrates the convergence of the Taylor series for 
r<1, the n^th root test shows that the Taylor series in 
(<ref>) does 
not converge for r>1 (the series converges for r=1 by the alternating 
series test, however, this is a zero-probability event).
However, since tan^-1(x)∈[-π/2,π/2] and 
θ∈(-π/2,π/2), for any r and 
φ, 

    |tan^-1(sin(θ)+rcos(φ)/cos(θ)+rsin(φ))-θ|≤π.

Therefore, using the Taylor series expansion of log(1-x) around x=0 in
(<ref>), and (<ref>), it is 
straightforward to show there exist constants a∈(0,1) and b>0 such that:

    |tan^-1(sin(θ)+rcos(φ)/cos(θ)+rsin(φ))-θ|≤{[ rcos(θ+φ)+br^2        if r≤ a;              π      otherwise ].

We can use (<ref>) to upper bound the MSE,

    ⟨ (θ - θ̂_ het)^2 ⟩   ≤1/2π∫_0^2π(∫_0^a(rcos(θ+φ)+br^2)^2re^-x^2/2σ^2_ het/σ^2_ hetdr+∫_a^∞π^2re^-x^2/2σ^2_ het/σ^2_ hetdr)dφ
       =σ^2_ het+8b^2σ^4_ het-1/2e^-a^2/2σ^2_ het(2a^4b^2+a^2(1+8c^2σ^2_ het)+2σ^2_ het(1+8c^2σ^2_ het)-2π^2)
       =σ^2_ het+𝒪(σ^4_ het)
       =σ^2_ het+𝒪(1/n).





 §.§ Quantum Fisher information

The final CM describing Alice's state is

    V_A =
    	[      a_11 -a_12cosθ         0  a_12sinθ; -a_12cosθ      a_22 -a_12sinθ         0;         0 -a_12sinθ      a_11 -a_12cosθ;  a_12sinθ         0 -a_12cosθ      a_22 ]

The QFI for the parameter θ can be calculated using the formalism developed in <cit.>. One first must derive the quantum fidelity

    ℱ(ω)=1/√(√(Γ)+√(Λ)-√((√(Γ)+√(Λ))^2-Δ))

between the two-mode states with CM V_A ≡ V_A(θ) and V_A(θ+ω) and zero first moments. Equation (<ref>) is based on the symplectic invariants,

    Δ   =   (V_A(θ)+V_A(θ+ω)) ≥ 1
    Γ   =   16 (Ω V_A(θ) Ω V_A(θ+ω)-I_4× 4/4) ≥Δ
    Λ   =   16 (V_A(θ)+i/2Ω) (V_A(θ+ω)+i/2Ω) ≥ 0

where,

    Ω=
    	[  0_2× 2  I_2× 2; -I_2× 2  0_2× 2 ].

The QFI for the parameter θ F_A' is given by the second derivative,

    F_A'=-4 d^2 ℱ(ω)/dω^2|_ω=0= 8 a_12^2/4(a_11a_22-a_12^2)-1=8 a_12^2/4√( V_A)-1.

Substituting Eqs. (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) into Eq. (<ref>) we get

    F_A' = 4 n̅_LOn̅_Sη_eff/n̅_LO+(1-η_eff)n̅_B_eff(1+2n̅_LO)+η_effn̅_S

and by taking the limit n̅_LO→∞ we get,

    F_A = lim_n̅_LO→∞ F_A' = 4 n̅_Sη_eff/1+2n̅_B_eff(1-η_eff).







