Entanglement Enhanced Estimation of a Parameter Embedded in Multiple Phases

Quantum-enhanced sensing promises to improve the performance of sensing tasks using non-classical probes and measurements that require far fewer scene-modulated photons than the best classical schemes, thereby granting previously-inaccessible information about a wide range of physical systems. We propose a generalized distributed sensing framework that uses an entangled quantum probe to estimate a scene-parameter encoded within an array of phases, with a functional dependence on that parameter determined by the physics of the actual system. The receiver uses a laser light source enhanced by quantum-entangled multi-partite squeezed-vacuum light to probe the phases and thereby estimate the desired scene-parameter. The entanglement suppresses the collective quantum vacuum noise across the phase array. We report simple analytical expressions for the Cram\'er Rao bound that depend only on the optical probes and the physical model of the measured system, and we show that our structured receiver asymptotically saturates the quantum Cram\'er-Rao bound in the lossless case. Our approach enables Heisenberg limited precision in estimating a scene-parameter with respect to total probe energy, as well as with respect to the number of modulated phases. Furthermore, we study the impact of uniform loss in our system and examine the behavior of both the quantum and the classical Cram\'er-Rao bounds. We apply our framework to examples as diverse as radio-frequency phased-array directional radar, beam-displacement tracking for atomic-force microscopy, and fiber-based temperature gradiometry.


I. INTRODUCTION
Quantum phenomena are now known to be powerful and viable tools to enhance estimation precision in diverse fields, e.g., astronomy [1], general relativity [2][3][4][5], models for quantum-to-classical transition [6], microscopy [7], and optical imaging [8][9][10]. Quantumenhanced estimation in sensing applications, which arguably comprises the nearest-term realizable quantum technologies of practical importance, promises an improvement the sensitivity of estimating an unknown parameter of the physical system being probed. In an idealized sensing context, this improvement takes the form of an improved scaling of estimation variance with probe power, known as Heisenberg scaling. Moreover, this Heisenberg scaling for sensitivity can be obtained using Gaussian quantum states of light (that can be generated using lasers, linear optics, and squeezed light, e.g., produced using parametric amplifiers) and Gaussian measurements (i.e., homodyne and heterodyne detection). This is an especially interesting point, since Gaussian resources alone do not suffice to perform a variety of other quantum information tasks [11][12][13][14][15][16][17][18][19].
In the context of distributed, or networked, quantum sensing, entangled Gaussian states have been shown to yield additional, significant advantages over separable quantum-enabled probes, for which the states of individual optical modes can be independently characterized [20][21][22][23][24][25][26]. These entangled probes have the advantage of being generated from separable Gaussian states and Gaussian unitary operations [20], both of which can be readily analyzed [27,28] and realized experimentally [29,30]. * michaelgrace@email.arizona.edu One widely applicable scenario for distributed quantum sensing is as follows: a quantitative parameter of interest x modulates a series of M optical phase delays in M optical modes with non-symmetric but known functional relationships to the parameter, and the parameter must be estimated using user-controlled probes and measurements. This problem statement can be used to model various practical photonic sensing tasks, some specific examples of which have been recently studied. These include the estimation of the angle of incidence, or other attributes of a radio-frequency (RF) wave upon an array of sensor pixels where each pixel is a phase modulator that is read out optically [31]; estimating a small transverse displacement of an optical beam in an atomic-force microscope (AFM) [32]; estimation of a small angular velocity of rotation of a Sagnac-based fiber-optical gyroscope (FOG) [25]; and estimating material defects with a fiber-based temperature gradiometer. General results for variants of this scenario have been found in recent work focusing on measurements that perfectly achieve the ultimate bounds on estimation precision [33] and handle imperfections in prior knowledge of the parameter [34,35].
Here, we propose a framework for quantum sensing that combines the fundamental Heisenberg scaling advantage offered by entangled Gaussian probe states with practicality in analysis and implementation for real-world applications. The two estimation theory tools we use are the quantum Fisher information (QFI) H x and the classical Fisher information (CFI) I x . Both give lower bounds to the mean squared error (MSE) (x − x) 2 for unbiased estimatorsx (i.e., x = x), via the classical and quantum Cramér-Rao bounds: (x − x) 2 ≥ I −1 x ≥ H −1 x . We consider the lossless case to derive the fundamental optimal quantum performance and we also study the effects of uniform pure loss. Applications of our model reach be- Each twomode MZI is probed with a coherent state and one mode of an M -mode-entangled squeezed vacuum state. F , the Fourier gate, is an M -mode linear-optical interferometer, and the detector at the output of the circuit is a homodyne detector. L denotes the pure loss channel which is modeled as a beam splitter with transmissivity τ whose lower input mode is set to vacuum. A local oscillator mode for homodyne detection is implied but not shown for simplicity.
yond photonic sensors, and could serve as the foundation for more general distributed quantum sensing tasks with applications to quantum process tomography of quantum computers and quantum network tomography.

II. PHYSICAL MODEL AND ESTIMATION TASK
In this work, we consider a scenario shown in Fig. 1: a single scalar scene-parameter x modulates phases θ m (x), 1 ≤ m ≤ M , in an array of M Mach-Zehnder Interferometers (MZIs), each with a pair of input and output modes. The state preparation circuit and the phase modulation can be described by the passive unitary transformations where F signifies a linear-optic Fourier gate, D θ(x) is an M × M diagonal matrix with entries θ m (x), and I is an identity matrix. In addition, the 2M optical modes of the circuit each undergo identical pure loss channels with optical transmissivity 0 ≤ τ ≤ 1, imposing modesymmetric loss on the probe state before modulation by the parameter-dependent phases (e.g., induced during propagation to spatially distributed sensors). Given a total mean photon number budget N , one is tasked to design a 2M -mode (quantum) optical probe and an associated receiver to minimize the variance of the estimate of x. The probe we consider consists of a single-mode squeezed vacuum (SV) |0; r with squeezing parameter r > 0 and a (laser light) coherent state |α with complex amplitude α ∈ C. Each of these two inputs is equally mixed with M − 1 vacuum modes on balanced Fourier gates. The output of the first Fourier gate is an M -mode entangled continuous-variable (CV) state [22], and the output of the second is M identical (product) coherent states, |α/ √ M . Our scheme applies a linearoptical receiver circuit to the M modulated modes followed by single-mode homodyne detection. We will show that when the parameter x is known to be close to some value x 0 , an appropriate recombining operation involves a set of phase shifts U x0 , the inverse unitary U † I , and a single-mode phase shift U H .

III. QUANTUM FISHER INFORMATION
The QFI H x quantifies optimal precision in estimating x by any receiver. As H x is directly a function of ρ x , it is dependent on the choices made for the input state |ψ 0 and the preparation circuit U I . In the quantum metrology literature [36], it has been shown for various settings that in the absence of loss and noise, quantum optical probes can attain the so-called Heisenberg scaling H x = O(N 2 ), where N is the total photon-unit energy in the probe field. In contrast, sensing schemes using classical probes can only achieve H x = O(N ) at best. If the CFI I x for a specific quantum probe and a receiver scales as O(N 2 ), that system design achieves Heisenberg scaling with probe energy, while a CFI I x = H x indicates a quantum optimal receiver for a given probe.
The QFI leading to a quantum Crámer-Rao bound on the estimation of the parameter x can be expressed as follows (∂ x ≡ ∂/∂x): where H kl are the indexed matrix elements of the quantum Fisher information matrix for multiple parameter estimation of the M phases θ m (x). Since QFI is agnostic of the chosen receiver, we compute the H kl from the intermediate, generally mixed, state ρ x that is found immediately after modulation by the parameter-dependent phases.
Since the Gaussian input state |ψ 0 is transformed to ρ x solely by Gaussian transformations, the state ρ x will also be Gaussian. The QFI matrix elements H kl can then be computed from the displacement vector d x and covariance matrix (CM) V x of the state ρ x . The displacement vector and CM of the probe state shown in Fig. 1 are [37] where α = (q 0 + ip 0 )/ √ 2 and In general, the displacement vector and CM of a Gaussian state transform evolve through a Gaussian unitary U as d out = S d in and V out = SV in S T , where S is the symplectic matrix satisfying with Ω = antidiag(1, −1). Furthermore, a generic displacement vector and CM evolve through a symmetric pure loss channel via where X τ = √ τ I, Y τ = [(1 − τ )/2]I, I is the 4M × 4M identity matrix, and 0 ≤ τ ≤ 1 is the transmittance the channel acting identically on each mode.
After some calculations (Appendix A), we find that where the mean thermal photon number and the reduced squeezing parameter s = 1 4 ln 1 + (e 2r − 1)τ 1 + (e −2r − 1)τ , are given as functions of the input squeezing r and the symmetric transmittance τ , and where Therefore, the QFI for the parameter x under uniform pure loss is found by utilizing Eq. (12) and inserting Eq. (9) into Eq. (3) and takes the form, where ∂θ( By choosing r > 0 ⇒ s > 0, the choice of α that maximizes the QFI of Eq. (13) is α = α * . This can be seen by noting that the only dependence of Eq. (13) on an α which is not in a modulo, is in Eq. (12). There for s > 0, the prefactor of α 2 + α * is positive and therefore the choice α = α * is optimal.
In the absence of loss (i.e., τ = 1), we find that as a function of the mean photon numbers N s = sinh 2 r and N v = |α| 2 of the input SV state and coherent state, the QFI of Eq. (13) reduces to where N = N v +N s . The lossless QFI is clearly optimized by choosing α ∈ R, a property inherited from Eq. (13) which is valid for all 0 ≤ τ ≤ 1 . By inspection, the first two terms of Eq. (14) scale quadratically with input energy as N v 1 and N s 1, providing the Heisenberg scaling advantage available with quantum probes; the first term depends on the energy contribution from both input sources, while the second depends on the energy from the SV source.

IV. PERFORMANCE EVALUATION OF PROPOSED RECEIVER
Let us assume that x is known to be close to a value x 0 , i.e., |x − x 0 | 1, due to either a priori information or a preliminary estimate of x. Under this condition, applying the conjugate phases −θ m (x 0 ) to the state ρ x followed by a sequence of 50-50 beamsplitters (i.e., the second beamsplitters of the MZIs) and inverse Fourier gates F † efficiently recombines the information-bearing light to one desired output modeb 1 (Fig. 1). The phase φ H controls which field quadrature is measured via homodyne detection. The phases −θ m (x 0 ) could be unitarily evolved to a different set of phases −θ m (x 0 ) that are applied after the second set of beamsplitters; this configuration maintains the physical sensor's natural mathematical description as an array of M MZIs [25,31,32] while not changing the CFI. If prior knowledge is not available of the functional form of the phases θ m (x) [38] or of the parameter x itself [34,35], the linear recombination of light must be determined via an adaptive strategy.
We begin by assuming lossless optical sensing (i.e., τ = 1). The 2M -mode output state ρ H is a pure state and is determined by the input state |ψ 0 and the full system unitary U (x) = U H U † I U † x0 U x U I , where U H = diag(e iφH , 1, . . . , 1). Only two matrix elements in the system unitary U (x) will be necessary to calculate the CFI in this case; these elements can be evaluated as Since |ψ 0 and U (x) are Gaussian, the output of a realquadrature homodyne measurement on modeb 1 is characterized by the mean d H,1 (x) and variance V H,1,1 (x) of the first mode of the output state ρ H [27,28]. These can be read off from the first moment vector d where the displacement vector and CM of the input state are given by Eqs. (4) and (5) and the symplectic matrix S(x) can be found from U (x) using Eq. (6). Recognizing the optimality of α ∈ R for the QFI, we set α = q 0 (i.e., p 0 = 0). Therefore, the real quadrature of modê a M +1 has the only non-zero first moment among the input quadratures, and we have where To evolve the CM of the probe state [Eq. (5)], we use the fact S(x)S T (x) = I for the symplectic matrix of any passive unitary to find that the variance of the measured quadrature is where from Eq. (6) we have S 1,1 (x) = Re{U 1,1 (x)} and S 1,2M +1 (x) = Im{U 1,1 (x)}. We now consider symmetric pure loss acting on the optical modes of Fig. 1. We show in Appendix B that a symmetric pure loss channel across all modes of a system can be commuted with any passive unitary transformation. As a result, we can conceptually commute the pure loss channels to the very end of the optical circuit just before the homodyne measurement, such that only the effect of loss on the modeb 1 will affect the CFI. Taking into account a transmittance of τ ≤ 1 and using Eqs. (7) and (8), the displacement vector and covariance matrix of the measured mode become and V H,1,1 (x; τ ) = 1 2 1 + τ S 1,1 (x) 2 (e 2r − 1) For estimating a parameter x from a Gaussian random variable with mean d H,1 (x; τ ) and variance V H,1,1 (x; τ ), the CFI is known to be given by I x = I x,d + I x,V , where the two terms are given by The phase φ H is a free parameter for the receiver design, and three modes of operation can be identified. First, setting φ H = π/2 maximizes the sensitivity of d H,1 (x; τ ) to changes in the parameter x, such that This result is consistent with previously reported works on distributed sensing that involve a coherent state and SV probe and uniform loss but only consider equal phases on the M modes [22,25]. This modality is preferred in the realistic scenario where the photon-unit energy of the laser source dominates that of the squeezed light and also has the advantage of a fixed homodyning angle that does not depend on the characterization of the system. In particular, in the absence of loss (τ = 1) we have which exhibits Heisenberg scaling with probe energy. Second, the sensitivity of the covariance matrix V H,1,1 (x; τ ) can be optimized, yielding the optimal phase φ H = (1/2) arccos τ sinh(2r)/[1 + 2τ sinh(r) 2 ] . In this case, the CFI becomes Eq. (27) represents the best Fisher information that can be achieved with solely a squeezed vacuum input; in the absence of loss, this Heisenberg-limited term takes the form which resembles the optimal performance achievable with Gaussian probes [33]. Finally, the total CFI can be optimized by maximizing Eqs. (21) and (22) with respect to s ≡ sin 2 (φ H ), which yields the optimal value where g = 1 2 |α| 2 [csch(2r) +τ tanh(r) +τ ]. This optimization gives which reduces tõ Clearly, 0 ≤ s ≤ 1 for any operationally valid homodyne phase φ H , and since it can be shown that I x0 is monotonically increasing with respect to 0 ≤ s ≤ s opt , the optimal receiver configuration uses s = s opt when s opt < 1 and s = 1 (i.e., φ H = π/2) when s opt ≥ 1. The resulting CFI is In the lossless case, we havẽ and I (1) x0 is given by Eq. (25).

V. OPTICAL SENSING APPLICATIONS
Our framework applies to any situation in which the phases of M orthogonal optical modes are modulated by a common physical signal, one unknown parameter of which is to be estimated.

A. RF Signal Estimation with a Photonic Receiver
One example is the estimation of the angle of incidence x ≡ φ of an incident RF field using an M -pixel sensor array in an RF-photonic receiver antenna (Fig. 2), an application for which CV entanglement was recently shown to improve upon classical estimation precision [31].
Each phase element in Fig. 2 is optically read out by an integrated-photonic MZI circuit. The optical-frequency continuous-wave (cw) field in the waveguide mode in one arm of the m th MZI is E m (t) = Ee i[ωt+θm(φ)] , with where A is the RF-photonic amplitude-phase modulation efficiency, Ω is the center-frequency of the RF field, mb is the relative position of the m th sensor, and c = 3 × 10 8 m/s is the speed of light. We assume |φ − φ 0 | 1, e.g., when the RF wave is known to arrive at close to normal incidence (φ 0 = 0). Eq. (35) can be used to compute the prefactors ∂θ(φ) and ∂θ(φ) 2 in the QFI and CFI calculations. Fig. 3 plots the CFI I φ0 for our fully optimized sensor design compared with the QFI H φ0 evaluated at φ = φ 0 , given total probe energy N ∝ M and varying degrees of symmetric loss. In each case, the energy allocation between the SV and coherent probe states is optimized to maximize the QFI or CFI (see insets). The QFI and CFI for classical sensors, where r = 0, are also shown for reference. In the absence of loss (Fig. 3a), the Heisenberg scaling H φ0 ∝ I This scaling advantage is only accessed by the use of the quantum SV probe; indeed, in the absence of loss, it is optimal to allocate all of the input energy to the SV state, in which case I φ0 = 2 ∂θ(φ 0 ) 2 N s (N s + 1), (36) indicating that our receiver asymptotically achieves the quantum limit for distributed angle-of-incidence estima- (yellow diamonds), we optimize over input energy allocations for squeezing and displacement such that each point corresponds to a different numerically determined squeezing parameter. In the remaining two cases all energy is designated to the coherent state, such that r = 0. For all plots we set A = 0.1, Ω = 30 kHz, b = 10 m, and φ0 = 0 in Eq. (35). Insets: optimal fraction of total probe energy allocated to squeezed vacuum to maximize H φ 0 (blue circles) and I φ 0 (yellow diamonds).
tion in the lossless case. Importantly, our entanglementenhanced scheme enables an increasingly large advantage over classical designs for large RF sensor arrays.
While the use of the complementary coherent state probe adds no benefit for our distributed sensor design in the absence of loss, it contributes to improved distributed sensing performance in the presence of loss. As loss is introduced (Figs. 3b-d), we observe several changes to the performance of our receiver. First, the Heisenberg scaling effect is immediately lost with the introduction of loss, which is a well-known phenomenon for quantum sensing. On the other hand, the QFI of the classical sensor (r = 0) draws nearer to the QFI of the entangled sensor as the severity of the loss increases, reflecting the fact that energy is increasingly shifted away from the SV state to the coherent state in the optimization of the QFI H φ0 . The coherent state is even more quickly prioritized as loss increases in the optimization of the CFI φ0 , which falls gradually further from the QFI H φ0 as loss is increased, though it never becomes sub-optimal by more than a factor of ∼ 5. In the high loss regime, it is easy to see that s opt > 1, so I φ0 , and with most of the energy allocated to the coherent state the performance of the entangled sensor converges to that of the classical sensor. These results demonstrate the value of including a coherent state in the probe design for lossy entangled sensing, as does Fig. 4, which compares the fully optimized CFI against the QFI as the energy allocation ratio is traded off between the SV and coherent state. Allocating energy to the coherent state allows our design to achieve optimal estimation precision, especially with higher loss and larger distributed sensing networks, and also allows for a practically feasible increase of probe energy. In the conclusion section we briefly describe how to design a structured receiver that exactly achieves the QFI in all parameter regimes.

B. Optical Beam Displacement Tracking
Our framework allows us to reproduce results from sensing tasks that can be unraveled into M MZIs. One example is the use of spatial-mode entanglement to estimate a small lateral displacement δ of an optical beam [32], for example to implement a quantum-enhanced AFM. The input modesâ m are a set of 2M spatially overlapping, mutually-orthogonal (e.g., Hermite-Gauss) modes with negligible loss in the near-field regime. The output modesb m are vacuum-propagation normal modes extracted by the receiver. Beam displacement causes modal-energy crosstalk quantified by a matrix Γ. This 2M -mode crosstalk can be unitarily converted into a set of M MZIs with phases 2λ m δ, where the λ m depend on the eigenvalues of Γ. We thus can find the CFI evaluated at δ 0 = 0, using Eq. (25), where the prefactor becomes Notably, since it was proven that M m=1 λ m ∝ M 3/2 [32], our analysis recovers the linear dependence of the CFI prefactor on M , which arises because the spatial-mode crosstalk becomes progressively more sensitive to the beam displacement δ as the mode order m increases [32].

C. Optical Temperature Gradiometry
Our framework is also equipped to describe temporally entangled optical probes and dynamic physical systems, for which the time-bandwidth product M = W T gives the number of orthogonal temporal modes for an optical source bandwidth W and integration time T . For example, consider the estimation of the thermal conductivity k of a uniform, dielectric rod with density ρ and specific heat c p , where the assumption |k − k 0 | 1 could stem from knowledge that k diverges slightly from that of a known material due to physical impurities. We embed one branch of an optical fiber-based MZI at a po-sition y = y 0 along the rod. If the rod is heated to an initial temperature distribution u(y, 0) and allowed to relax to steady-state, the sensor could probe the kdependent optical phases induced by the temperature u(y 0 , t) at times t = m/W using an M -temporal-modeentangled CV state and M (product) coherent states. The functional form θ m (k) of the phases in the M effective MZIs will depend both on the solution to the heat equation ∂ t u(y, t) = (k/ρc p )∂ 2 y u(y, t) and on the temperature-dependent Sellmeier equation for the optical refractive index of the fiber material. As long as the first derivatives ∂ k θ m (k) can be computed analytically or numerically, Eq. (33) gives the fully optimized CFI for the temporally entangled sensor, which can be compared with the QFI of Eq. (13) in the presence of fiber loss. In addition, with constant-power sources, the photonunit probe energy N will naturally scale linearly with M . The Heisenberg scaling I k = O(N 2 ) therefore extends to I k = O(M 2 ) under a lossless approximation. For low to moderate levels of loss, we expect a significant advantage from entanglement for long integration times (Fig. 3b-d).

VI. DISCUSSION
Many photonic sensing tasks can be reduced to estimating a scalar parameter x that modulates phases θ m (x) in M MZIs. We discussed examples which include passive sensors whose receivers use quantum-enhanced computing, e.g., a quantum-enhanced RF photonic receiver, as well as active sensors that probe a scene with a non-classical illumination, e.g., beam tracking for AFM. We proved a Heisenberg limited scaling of the Fisher Information I x = O(N 2 ) in estimating the parameter in terms of the total photon-unit energy N employed by the sensor. Furthermore, we argued that under certain circumstances, we see Heisenberg scaling in M as well, i.e., I x = O(M 2 ). The latter, which is true for two of our examples, becomes significant when M is large.
Additional constant factor improvements in I x are possible over our simple receiver design, e.g., by optimizing U I for known functions θ m (x). This will involve arbitrary tuning of O(M 2 ) phases in a programmable linear-optic circuit [39,40]. Under optical loss and noise, the Heisenberg scalings with respect to N and M will disappear, and I x = O(M N ) will prevail. However, there will be a constant factor improvement in I x over a classical sensor in the long integration time limit, which can be significant for moderate losses and if M is large [22,25,26].
We note that if one desires to find the optimal receiver, a tractable method could be based on the symmetric logarithmic derivative (SLD) whose eigenvectors provide an optimal measurement. Indeed, the SLD is a Gaussian polynomial [41,42] (i.e. at most quadratic inâ i andâ † i or equivalently inq i andp i , where index i counts the modes). Therefore, one could follow standard methods to diagonalize the Gibbs matrix corresponding to the SLD operators, i.e., similarly to [43] where single FIG. 5. A schematic of the optimal structured receiver design that attains the QFI for estimating the single scalar parameter x for the problem studied in this paper, as sketched in Fig. 1. UL and UR are multimode passive linear-optical transformations, and S1, . . . , S2M are single-mode squeezers. The homodyne detection receiver strategy depicted in Fig. 1 is strictly suboptimal in the presence of losses. mode systems have been examined. If the SLD operators for each parameter commute [44] so that there exists a measurement which attains the QFI, a similar method to [26] could work in a multi-parameter and multi-mode setting potentially revealing a rich measurement structure. For estimating a single scalar parameter embedded in a multimode Gaussian state, which is relevant to the problem studied in this paper, the optimal receiver design takes the form shown in Fig. 5. The optimal receiver takes the form of a multi-mode Gaussian unitary transformation of the modulated lossy multimode Gaussian state ρ H followed by mode-resolved photon number resolving (PNR) on all 2M modes. The general Gaussian unitary U G preceding the PNR detector array can be further decomposed into a linear optical unitary U L followed by single-mode squeezers S 1 , . . . , S 2M , and another linear optical unitary U R , as shown [45]. We will explore the optimal receiver design problem for general Gaussian multi-parameter estimation problems in future work.
To calculate the QFI for a parameter x which is embedded in a set of phases {θ 1 (x), . . . , θ M (x)}, we need first to calculate the QFI matrix for the multiple phase estimation problem whose matrix elements H kl are, where E n and |e n are respectively the eigenvalues and eigenvectors of the 2M -dimensional state ρ x just after modulation by the parameter-dependent phases. Then, we apply the Jacobian transformation and the QFI for parameter x is given by Eq. (3). Here we focus on calculating H kl . The final state (and the state in any stage of our setup) is Gaussian. Therefore, we can exploit the phase-space toolbox to find the eigenvalues and eigenvectors of the final state. We also note that since we assume uniform loss, we can commute the single-mode single channels with unitary transformations (see Appendix B). We choose for our derivation to conceptually commute the pure loss channels to act on the input state, i.e., before the interferometer. Also, the eigenvalues of the state will not change after the action of pure loss as subsequently the state is transformed under unitary operations.
The 4M × 4M CM V τ of the state just after the pure loss channels is found by transforming the input CM V 0 [Eq. (5)] according to Eq. (8). By calculating the (usual) eigenvalues of |iΩV τ |, we find that the symplectic diagonalization of V τ is where D τ = diag (ν 1 , 1/2, . . . , 1/2), and where ν 1 = τ (1 − τ ) sinh 2 r + 1/4 (we will examine the S τ after we find the eigenvalues). We write the eigenvalues compactly as, whereN 1 (10) is the mean thermal photon number in the first mode arising from the action of the pure loss channel on squeezed vacuum. As discussed before, Eq. (A3) are the eigenvalues of ρ x as well. Now, let us find the eigenvectors of the final state. It is easy to verify the following equations, S τ = diag e s , 1 . . . , 1, e −s , 1, . . . , 1 (A4) where S τ is a 4M × 4M matrix. Equations (A5), (A6), and (A7) guarantee that S τ is a symplectic matrix and, as per Eqs. (A2) and (A4), the symplectic matrix S τ diagonalizes (in the symplectic sense) the CM V τ . From the structure of S τ we deduct that it corresponds to a single mode squeezer on the first mode and identity on the rest of the modes. In the absence of loss (i.e., τ = 1), Eq. (11) gives s = r. We must not forget the displacement of the state prior to loss in the M + 1 mode. After loss is applied [Eq. (7)], the displacement of the state is given by d τ = (0, . . . , 0, √ τ q 0 , 0, . . . , 0, √ τ p 0 , 0, . . . 0), describing a phase-space displacement by √ τ α on the M + 1 mode. Therefore, the diagonal form of final state is where the eigenvalues are given in Eq. (A3), the eigenvectors |e n are where U x is the unitary which imprints the phases onto the state, U I is the passive unitary interferometer, D( √ τ α) (M +1) is a displacement operator that acts only on mode M + 1, U (1) sq (s) is a single-mode squeezer acting on the first mode, and | n = |n 1 , . . . , n 2M represents a product of 2M Fock states.
Coming back to the QFI of Eq. (A1) we observe that our eigenvalues [(A3)] do not depend on the parameter, as expected since we our parameters are imprinted by a unitary operation. Therefore, Eq. (A1) simplifies to e n |n k |e m e m |n l |e n , (A10) where we took under consideration the fact that the QFI elements for our case are always real. Taking into account the form of Eq. (A3), Eq. (A10) can be simplified further as the Kronecker deltas will set all Fock number to zero, except of n 1 and n M +1 . We get, where, A (k) n1,m1 = ψ(n 1 )|U (1) † sq (s)U † I n k U I U (1) sq (s)|ψ(m 1 ) (A13) B (k,l) n1 = ψ(n 1 )|U (1) † sq (s)U † I n k n l U I U (1) sq (s)|ψ(n 1 ) (A14) and where the state |ψ(n) = |n, 0, . . . , 0, √ τ α, 0, . . . , 0 is a product state between a Fock state |n in mode 1, a coherent state | √ τ α in mode M + 1, and vacua in the other 2M − 2 modes. For the specific U I considered in Eq. (1), working in the Heisenberg picture, and converging all summations involved in Eq. (A11), it is straightforward to calculate the amplitudes of Eqs. (A13) and (A14) to arrive at Eq. (9). Given a specific set of functions θ m (x), one can straightforwardly apply Eq. (3) to obtain the QFI for the parameter x.

Appendix B: Commutation of a Symmetric Pure Loss Channel with Gaussian Operations
Here we show that a symmetric set of M -mode pure loss channels with transmission coefficient τ on each mode commutes with any passive (i.e., energypreserving) Gaussian unitary transformation when the state being acted upon is Gaussian. Consider an arbitrary M -mode Gaussian state with displacement vector d and covariance matrix V . An arbitrary passive unitary matrix U B , which can in general be decomposed into twomode beamsplitters and phase shifts [46], is described by its symplectic matrix S B from Eq. (6) with S T B S B = I. In addition, if the state is subjected to symmetric pure loss, the resulting displacement vector d τ and CM V τ are calculated according to Eq. (7) and Eq. (8). Clearly, S B commutes with both X τ and Y τ .
Acting the unitary U B on the state at the output of the loss channel, we find a displacement vector where d B = S B d is the displacement vector of the initial state transformed by U B . The corresponding covariance matrix is where V B = S B V S T B is the transformed covariance matrix of the initial state subjected to U B . Taken together, Eqs. (B1) and (B2) prove the commutation of symmetric pure loss and passive Gaussian unitaries for Gaussian states.