# Distributed Quantum Sensing with Squeezed-Vacuum Light
in a Configurable Network of Mach-Zehnder Interferometers

Marco Malitesta  QSTAR, INO-CNR and LENS, Largo Enrico Fermi 2, 50125 Firenze,
Italy  Università degli Studi di Napoli ”Federico II”, Via Cinthia 21, 80126
Napoli, Italy  Augusto Smerzi  QSTAR, INO-CNR and LENS, Largo Enrico Fermi 2,
50125 Firenze, Italy  Luca Pezzè  QSTAR, INO-CNR and LENS, Largo Enrico Fermi
2, 50125 Firenze, Italy

###### Abstract

We study a sensor network of distributed Mach-Zehnder interferometers (MZIs)
for the parallel (simultaneous) estimation of an arbitrary number $d\geq 1$ of
phase shifts. The scheme uses a squeezed-vacuum state that is split between
$d$ modes by a quantum circuit (QC). Each output mode of the QC is the input
of one of $d$ MZIs, the other input of each MZI being a coherent state. In
particular, i) we identify the optimal configuration of the sensor network
that allows the estimation of any linear combination of the $d$ phases with
maximal sensitivity. The protocol overcomes the shot-noise limit and reaches
Heisenberg scalings with respect to the total average number of particles in
the overall probe state, the multiphase estimation only requiring local
photocounting. Furthermore, the parallel multiphase estimation overcomes
optimal separable strategies for the estimation of any linear combination of
the phases: the sensitivity gain being a factor $d$, at most. Viceversa, ii)
given a specific QC, we identify the optimal linear combination of the phases
that maximizes the sensitivity and show that results are robust against random
choices of the QC. Our scheme paves the ways to a variety of applications in
distributed quantum sensing.

##  I Introduction

Optical interferometry exploiting squeezed light LoudonJMO1987 ; ScullyBOOK ;
BreitbachNATURE1997 ; AndersenPS2016  has been – since the pioneer 40-years-
old proposal by Caves CavesPRD1981  – a cornerstone of theoretical
ParisPLA1995 ; BarnettEPJD2003 ; PezzePRL2008 ; LangPRL2013 ; Ruo-
BercheraPRA2015 ; SparaciariPRA2016 ; BondurantPRD1984  and experimental
WuPRL1986 ; XiaoPRL1987 ; GrangierPRL1987 ; PolzikPRL1992 ; GodaNATPHYS2008
photonic quantum sensing SchnabelPR2017 ; PirandolaNATPHOT2018 ;
LawrieACSP2019 ; PolinoARXIV . A Mach-Zehnder interferometer (MZI) with a
high-power coherent state in one input port and a low-intensity squeezed-
vacuum light in the other input can reach a phase estimation uncertainty
$\Delta^{2}\theta=e^{-2r}/\bar{n}_{T}$ CavesPRD1981 , where $r\geq 0$ is the
squeeze parameter, $\bar{n}_{T}$ is the total average number of photons and
$\theta$ is the relative phase shift between the two arms of the
interferometer. This scheme can overcome the shot-noise (SN) limit
$\Delta^{2}\theta_{\rm SN}=1/\bar{n}_{T}$ by an amount depending on the
squeezing strength $r$. Currently, squeeze factors of more than 10 dB have
been observed in several experiments VahlbruchPRL2016 ; SchonbeckOPTLETT2018 ;
SchnabelPR2017 ; LawrieACSP2019 . Furthermore, when the coherent and the
squeezed-vacuum input states have approximately the same intensity, the MZI
can achieve PezzePRL2008  the Heisenberg limit (HL) $\Delta^{2}\theta_{\rm
HL}=1/\bar{n}_{T}^{2}$. This prediction has been associated HofmanPRA2007 ;
PezzePRL2008  to the onset of NOON states after the first beam splitter of the
MZI, as verified experimentally AfekSCIENCE2010 . Enhancing the sensitivity by
replacing the normally-empty input port with squeezed-vacuum light is relevant
when there are constraints limiting the total light intensity inside the
interferometer. Gravitational wave detection is an important application
SchnabelNATCOMM2010 ; RafalPRA2013 ; ChuaCQG2014 ; AbadieNATPHYS2011 ;
AasiNATPHOT2013 ; TsePRL2019 ; AcernesePRL2019 , where the quantum-enhancement
offered by squeezing allows to boost substantially the expected rate of
detectable events. Quantum imaging MoreauNPR2019 ; Ruo-BercheraMETROLOGIA2019
, microscopy CasacioNATURE2021  and the probing of biological samples
TaylorNATPHOT2013  are other relevant applications that require high
resolution but low probe power TaylorPHYSREP2016 . Squeezed-vacuum states can
be also generated via spin-changing collisions in a Bose-Einstein condensate
GrossNATURE2011 ; HamleyNATPHYS2012 ; PeiseNATCOMM2015  and used to enhance
the sensitivity of atomic MZIs KrusePRL2016 ; PezzeRMP2018 .

Figure 1:  A distributed quantum sensor for the estimation of $d$ relative
phases $\theta_{1},...,\theta_{d}$ can follow an entangled (a) or a separable
(b) strategy. Panel (a) shows the Mach-Zehnder sensor network composed by $d$
MZIs. The input of the $j$th MZI is given by a coherent state
$|\alpha_{j}\rangle$ in the mode $a_{j}$, while the other input $b_{j}$ is one
of the outputs of a linear $d$-mode quantum circuit. The input of the quantum
circuit is a single squeezed-vacuum state $|\xi\rangle$, which is mixed with
$d-1$ vacuum states $|0\rangle$. The output state of the $d$-mode quantum
circuit is mode entangled: we thus identify the scheme of panel (a) as an
entangled multiphase estimation strategy. In panel (b) the different MZIs are
independent. The $j$th MZI has a coherent state $|\alpha_{j}^{\prime}\rangle$
in the input mode $a_{j}$ and a squeezed-vacuum state $|\xi^{\prime}\rangle$
in mode $b_{j}$. This scheme realizes a separable multiphase estimation
strategy. The table in panel (c) summarizes the main results of the manuscript
for what concerns the gain $\mathcal{G}$ of the entangled strategy over the
separable one, under different constraints $\mathcal{C}$ (see text for
details). Here, $\bar{n}_{T}$ is the total average number of particles,
$\bar{n}_{s}$ is the total number of particles in the squeezed state(s), and
prime symbols refer to the separable strategy.

To date, Mach-Zehnder interferometry using squeezed states has focused on the
estimation of a single phase shift SchnabelPR2017 ; LawrieACSP2019 ;
RafalPROGOPT2015 ; PezzeRMP2018 . Yet, several applications require the
estimation of multiple phases encoded simultaneously in a network of
spatially-separated sensors AlbarelliPLA2020 . Distributed quantum sensing
HumphreysPRL2013 ; LiuJPA2016 ; CiampiniSCIREP2016 ; ProctorPRL2018 ;
GePRL2018 ; ElderidgePRA2018 ; GessnerPRL2018 ; ZhuangPRA2018 ; NicholsPRA2018
; OhPRR2020 ; ZhuangNJP2020 ; GessnerNATCOMM2020 ; TriggianiARXIV ;
GebhartARXIV  is thus attracting increasing interest. Recently, different
schemes have been experimentally realized using squeezed light GuoNATPHYS2020
; XiaPRL2020 , photonics Greenberger–Horne–Zeilinger (GHZ) LiuNATPHOT2021 ,
Bell ZhaoPRX2021 , multimode NOON HongNATCOMM2021  and single-photon Fock
PolinoOPTICA2019 ; ValeriNPJ2020  states. In Ref. GuoNATPHYS2020  a displaced
squeezed state is split among four spatial modes, which then undergo a phase
shift and are finally measured by homodyne detection. Reference XiaPRL2020
has instead reported sensing of phase-space displacement using a squeezed
vacuum state split among three modes ZhuangPRA2018 . The sensor network of
Refs. LiuNATPHOT2021 ; ZhaoPRX2021  is based on the polarization rotation of
photonic qubits. In particular, Ref. ZhaoPRX2021  has realized distributed
sensing over large spatial distances. Finally, Refs. PolinoOPTICA2019 ;
ValeriNPJ2020  have considered a multimode interferometer on a photonic chip,
where single photons are split among many modes by a generalized beam
splitter, phase shifted and recombined by a second multimode beam splitter
CiampiniSCIREP2016 .

Surprisingly, so far, little effort has been devoted to study multiphase
estimation in a network of distributed MZIs. This system is relevant in
quantum optics, as well as in atom interferometry, with possible applications
in arrays of quantum clocks KomarNATPHYS2014 ; PolzikPRA2016  and spatial
magnetometry BaumgratzPRL2016 ; AltenburgPRA2017 ; ApellanizPRA2018 ;
HouPRL2020 . A key aspect of the MZI is that phase sensitivity bounds are well
defined and quantified in term of the total average number of particles in the
input state PezzePRL2007 , without ambiguities related to the resource cost
necessary to establish a phase-reference JarzynaPRA2012 ; PezzePRA2015 ;
GoldbergPRA2020  (e.g. for homodyne detection). Total resource counting is
crucial to quantify improved performances of quantum devices over classical
strategies and to claim sub-SN sensitivities.

In this manuscript, we study a distributed-sensing scheme that generalizes the
single MZI with coherent$\otimes$squeezed-vacuum input light CavesPRD1981 ;
ParisPLA1995 ; BarnettEPJD2003 ; PezzePRL2008 ; LangPRL2013  to a network of
$d\geq 1$ spatially-distributed MZIs, see Fig. 1(a). In our setup, a single
squeezed-vacuum state $|\xi\rangle$ of squeeze parameter $r$ is first split by
a quantum circuit (QC) consisting of a $d$-mode beam splitter GuoNATPHYS2020 ;
XiaPRL2020 . The QC is identified by a unitary transformation $\bm{U}$ that
can be realized, in practice, by a sequence of two-mode linear operations
ReckPRL1994 ; NokkalaNJP2018 . The sensor network of Fig. 1(a) can thus be
scaled to an arbitrary number, $d$, of interferometers. The linear splitting
of the squeezed-vacuum state generates entanglement among the $d$ output modes
$b_{1},...,b_{d}$ of the QC. We thus identify the scheme of Fig. 1(a) as an
entangled multiphase estimation strategy. Each mode $b_{j}$ is used as the
input mode of a MZI, the other input $a_{j}$ being a coherent state. It should
be notices that the $d$ coherent states $|\alpha_{j}\rangle$ are phase locked
with the squeezed-vacuum state $|\xi\rangle$. In the $j$th MZI ($j=1,...,d$),
the two input modes mix at a balanced beam splitter, encode a relative phases
$\theta_{j}$ and are detected by photocounting after a final balanced beam
splitter. Here, we estimate arbitrary linear combinations
$\bm{v}\cdot\bm{\theta}=\sum_{j=1}^{d}v_{j}\theta_{j}$ of the $d$ phase shifts
$\bm{\theta}=\\{\theta_{1},...,\theta_{d}\\}$, where
$\bm{v}=\\{v_{1},...,v_{d}\\}$ is a real vector. The estimation method is
based on a multimode moment-matrix approach GessnerNATCOMM2020  and the
corresponding sensitivity is compared to the multiparameter quantum Cramér-Rao
bound HelstromBOOK ; HolevoBOOK ; ParrisIJQI2009 .

We identify different regimes – depending on the relative intensity of the
squeezed-vacuum and the coherent states – and predict sub-shot noise
sensitivities, up to the Heisenberg limit, with respect to the total number
$\bar{n}_{T}$ of photons used. In particular, i) given a specific linear
combination $\bm{v}\cdot\bm{\theta}$ that one wants to estimate, we identify
the optimal configuration of the sensor network of Fig. 1(a) that minimizes
the uncertainty $\Delta^{2}(\bm{v}\cdot\bm{\theta})$. The analytical
optimization of $\bm{U}$, $|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}$ and $r$ is
supported by numerical calculations. In particular, we identify conditions for
which the strategy of Fig. 1(a) is never surpassed (for any
$\bm{v}\cdot\bm{\theta}$) by the separable strategy of Fig. 1(b), which uses
independent MZIs for the estimation of each $\theta_{j}$. It should be also
noticed that the separable strategy uses $d$ squeezed-vacuum states (one for
each MZI), while the entangled strategy uses a single squeezed state. The
table in Fig. 1(c) summarizes our results concerning the maximum gain
$\mathcal{G}$ of the entangled strategy over the separable one, under
different constraints $\mathcal{C}$. There, $\mathcal{G}$ is understood as
optimized over all $\bm{v}$. When $\bar{n}_{s}=\sinh^{2}r\gg d$, we find a
maximum gain given by a factor $d$. In the opposite limit, for $\bar{n}_{s}\ll
d$ – in particular for $d\to\infty$ – we find a finite gain $e^{2r}$. The
later result is surprising since, in this case, the squeezed vacuum is mixes,
at the QC, with a diverging number of vacuum states. Viceversa, ii) for any
given QC transformation $\bm{U}$, we identify optimal and orthogonal linear
combinations $\bm{v}\cdot\bm{\theta}$ that can be estimated with the highest
possible sensitivity. This optimization problem leads us to introduce the
useful concepts of Fisher and squeezing spectra. In particular, we show that
results are robust against random choices of $\bm{U}$. Our findings pave the
way toward distributed multi-phase estimation in a network of MZIs, using
quantum states and detection capabilities that are common to many
laboratories.

The paper is structured as it follows. Secion II introduces basic and general
notions of distributed sensing. In particular, we recall the multimode moment-
matrix approach used in this work and the multiparameter Cramér-Rao bound.
Section III illustrates in details the Mach-Zehnder sensor network of Fig.
1(a). Section IV presents the optimization for arbitrary $\bm{v}$. First, we
provide, in Sec. IV.1, upper and lower bounds to the sensitivity and present
different sensitivity scalings with respect to the total average number of
particles. We then compare, in Sec. IV.2, optimal entangled and separable
strategies under different constraints, providing a detailed discussion of the
results presented in the table of Fig. 1(c). Section V studies the optimal
linear combination of phases that can be estimated for a given configuration
of the Mach-Zehnder sensor network. We first introduce the notion of Fisher
and squeezing spectra, in Sec. V.1, and then apply this formalism to calculate
the sensitivity achievable for random choices of the QC, in Sec. V.2. We
finally, in Sec. VI, compare our finding with the literature and conclude.

##  II Distributed quantum sensing

###  II.1 Definition and strategies

In a distributed quantum sensing problem, $d$ unknown parameters
$\theta_{1},...,\theta_{d}$, are encoded in independent (e.g. spatially-
separated) modes or interferometers. Parameter encoding is described by
commuting transformations. In the ideal noiseless scenario, this is given by
the unitary evolution
$e^{-i\bm{\theta}\cdot\hat{\bm{H}}}=\otimes_{j=1}^{d}e^{-i\hat{H}_{j}\theta_{j}}$,
where $\hat{\bm{H}}=\\{\hat{H}_{1},...,\hat{H}_{d}\\}$ is a set of commuting
Hermitian operators, $[\hat{H}_{i},\hat{H}_{j}]=0$ for $i,j=1,...,d$, and
$\bm{\theta}\cdot\hat{\bm{H}}=\sum_{j=1}^{d}\theta_{j}\hat{H}_{j}$.

The sensing scheme can follow an entangled (also indicated as parallel, or
global, in the literature) or a separable (sequential, or local) strategy
ProctorPRL2018 ; GePRL2018 ; GessnerPRL2018 ; KnottPRA2016 ; GuoNATPHYS2020 ;
XiaPRL2020 . In an entangled strategy, the overall probe state $\hat{\rho}$ of
the sensor network is prepared in a mode-entangled state. In contrast, a
separable strategy uses the product state
$\hat{\rho}=\bigotimes_{j=1}^{d}\hat{\rho}_{j}$, where $\hat{\rho}_{j}$ is the
probe state of the $j$th sensor: the different sensors are thus uncorrelated
and the parameters $\theta_{1},...,\theta_{d}$ are estimated independently.
Notice that classical correlations among the different $\hat{\rho}_{j}$ are
not useful to increase the multiparameter sensitivity, in general
GessnerPRL2018 . Finally, an interesting possibility is to consider local
measurements at each sensor, without requiring a mode entangled measures
(although also distributed sensing scheme based on a final recombination of
parameter-sensing modes have been considered CiampiniSCIREP2016 ; GePRL2018 ;
TriggianiARXIV ; PolinoOPTICA2019 ; ValeriNPJ2020 ; OhPRR2020 ): local
measurements are advantageous when the sensing modes are spatially delocalized
ZhaoPRX2021 .

###  II.2 Figure of merit

One of the goals of multiparameter estimation, in general, is to infer linear
combinations $\bm{v}\cdot\bm{\theta}=\sum_{j=1}^{d}v_{j}\theta_{j}$ of $d$
parameters $\bm{\theta}$ encoded in the quantum device GePRL2018 ;
ProctorPRL2018 ; GessnerPRL2018 ; XiaPRL2020 ; GuoNATPHYS2020 ; LiuNATPHOT2021
; ZhaoPRX2021 ; RubioJPA2020 ; GrossJPA2021 ; QianPRA2019 . In the following
we take $v_{j}$ real (either positive or negative) and $v_{j}\neq 0$ for all
$j=1,...,d$ to guarantee an irreducible $d$-parameter problem. We also
consider the normalization $|\bm{v}|^{2}=\sum_{j=1}^{d}v_{j}^{2}=1/d$, without
loss of generality. An example of linear combination of parameters is the
average $\bm{v}\cdot\bm{\theta}=(\theta_{1}+\theta_{2}+...+\theta_{d})/d$,
corresponding to $v_{j}=1/d$.

The method of moments is a feasible approach to multi-parameter estimation
GessnerNATCOMM2020 . Here, it is based on a set of $d$ Hermitian and commuting
measurement operators $\hat{X}_{j}$ whose mean
$\big{\langle}\hat{X}_{j}\big{\rangle}$ is a monotonic function of
$\theta_{j}$ only. The estimation method consists of repeating the measurement
of the local observable $\hat{X}_{j}$ several times. Taking the average value
$\bar{X}_{j}$ and inverting the equation
$\big{\langle}\hat{X}_{j}\big{\rangle}=\bar{X}_{j}$ provides an estimate of
$\theta_{j}$. For the entangled multiparameter scenario, the method achieves
an uncertainty GessnerNATCOMM2020

$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}=\bm{v}^{T}\bm{\mathcal{M}}^{-1}\bm{v},$ (1)

where $\bm{\mathcal{M}}=\bm{G}^{T}\bm{\Gamma}^{-1}\bm{G}$,
$\bm{G}_{ij}=\partial\langle\hat{X}_{i}\rangle/\partial\theta_{j}$, and
$\bm{\Gamma}_{ij}=\langle\hat{X}_{i}\hat{X}_{j}\rangle-\langle\hat{X}_{i}\rangle\langle\hat{X}_{j}\rangle$
are $d\times d$ matrices ($\bm{G}$ being diagonal in this case), and the
expectation values are calculated with respect to the joint output state of
the $d$ sensors,
$e^{-i\hat{\bm{H}}\cdot\bm{\theta}}\hat{\rho}e^{i\hat{\bm{H}}\cdot\bm{\theta}}$.
The covariance matrix $\bm{\Gamma}$ expresses correlations between measurement
observables. These correlations are directly linked to the entanglement in the
probe state $\hat{\rho}$ and can be engineered to enhanced the sensitivity in
the estimation of certain combinations $\bm{v}\cdot\bm{\theta}$. In Ref.
GessnerNATCOMM2020  the moment matrix $\bm{\mathcal{M}}$ has been also used to
characterize and detect metrological multimode squeezing. In separable
strategies, $\bm{\Gamma}$ is diagonal and Eq. (1) becomes

$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}=\sum_{j=1}^{d}\frac{v_{j}^{2}\Delta^{2}\hat{X}_{j}}{(d\big{\langle}\hat{X}_{j}\big{\rangle}/d\theta_{j})^{2}},$
(2)

where
$\Delta^{2}\hat{X}_{j}=\big{\langle}\hat{X}_{j}^{2}\big{\rangle}-\big{\langle}\hat{X}_{j}\big{\rangle}^{2}$,
with the expectation values calculated on the output state of the $j$th
sensor, $e^{-i\hat{H}_{j}\theta_{j}}\hat{\rho}_{j}e^{i\hat{H}_{j}\theta_{j}}$.

The ultimate sensitivity limit in the estimation of $\bm{v}\cdot\bm{\theta}$
is provided by the quantum Cramer-Rao bound HelstromBOOK ; HolevoBOOK ;
ParrisIJQI2009 . In the entangled setting, we have
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\geq\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm eQCR}$, where

$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm eQCR}=\bm{v}^{T}\bm{\mathcal{F}_{\rm
Q}}^{-1}\bm{v},$ (3)

and $\bm{\mathcal{F}_{\rm Q}}$ is the $d\times d$ quantum Fisher information
matrix (QFIM) notaQFIM . If the overall probe state is pure,
$\hat{\rho}=|\psi\rangle\langle\psi|$, then the QFIM is $(\bm{\mathcal{F}_{\rm
Q}})_{ij}=4(\langle\psi|\hat{H}_{i}\hat{H}_{j}|\psi\rangle-\langle\psi|\hat{H}_{i}|\psi\rangle\langle\psi|\hat{H}_{j}|\psi\rangle)$
and Eq. (3) can be saturated by optimal measurements and estimators
MatsumotoJPA2002 ; PezzePRL2017 . In the sequential setting, the quantum
Cramer-Rao bound is $\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}\geq\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm sQCR}$, where

$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
sQCR}=\sum_{j=1}^{d}\frac{v_{j}^{2}}{\mathcal{F}_{j}},$ (4)

and $\mathcal{F}_{j}$ is the (scalar) quantum Fisher information HelstromBOOK
; HolevoBOOK ; BraunsteinPRL1994 . For pure states, Eq. (4) is obtained from
Eq. (3) when taking the product
$|\psi\rangle=\bigotimes_{j=1}^{d}|\psi_{j}\rangle$ such that
$\bm{\mathcal{F}_{\rm Q}}$ becomes diagonal with entries
$\mathcal{F}_{j}=(\bm{\mathcal{F}_{\rm
Q}})_{jj}=4(\langle\psi_{j}|\hat{H}_{j}^{2}|\psi_{j}\rangle-\langle\psi_{j}|\hat{H}_{j}|\psi_{j}\rangle^{2})$.

We recall that the different terms on the right-hand side of Eqs. (1)-(4) are
understood as divided by the number of repeated independent measurements $m$
used for the estimation. In particular, Eqs. (1)-(4) can be saturated, in
general, for $m\gg 1$. To simplify the notation, we neglect the factor $m$
here and in the following (see Ref. GebhartARXIV  for a multi-parameter
Bayesian estimation analysis including $m$ as a resource).

##  III Mach-Zehnder sensor network

The quantum distributed sensing scheme considered in this manuscript consists
of a network of $d$ MZIs, see Fig. 1. The entangled strategy is shown in Fig.
1(a), while the separable one in Fig. 1(b). In both cases, the $j$th
interferometer is described by the unitary phase encoding transformation
$e^{-i\theta_{j}\hat{H}_{j}}$, where
$\hat{H}_{j}=(\hat{a}_{j}^{\dagger}\hat{b}_{j}-\hat{b}_{j}^{\dagger}\hat{a}_{j})/2i$,
$\hat{a}_{j}$ and $\hat{b}_{j}$ ($\hat{a}_{j}^{\dagger}$ and
$\hat{b}_{j}^{\dagger}$) are bosonic annihilation (creation) operators for the
two interferometer modes, respectively, and $\theta_{j}$ is a relative phase
shift between the two interferometer arms. The $d$ MZIs are independent
(namely,
$[\hat{a}_{i},\hat{a}_{j}^{\dagger}]=[\hat{a}_{i},\hat{b}_{j}^{\dagger}]=[\hat{a}_{i},\hat{b}_{j}]=0$
for $i\neq j$), which guarantees $[\hat{H}_{i},\hat{H}_{j}]=0$. Furthermore,
we take local measurement observables
$\hat{X}_{j}=(\hat{a}_{j}^{\dagger}\hat{a}_{j}-\hat{b}_{j}^{\dagger}\hat{b}_{j})_{\rm
out}/2$, counting the relative number of photons at the output port of the
$j$th MZI, for $j=1,...,d$.

In the entangled strategy of Fig. 1(a), the input mode $a_{j}$ of the $j$th
MZI is fed with a coherent state $|\alpha_{j}\rangle$, where
$\alpha_{j}=|\alpha_{j}|e^{i\phi_{j}}$. The other input mode $b_{j}$ is fed
with the state obtained by the multi-mode splitting of a single squeezed-
vacuum state $|\xi\rangle$. Here, $\xi=re^{i\varphi}$, $r$ is the squeeze
parameter, $\bar{n}_{s}=\sinh^{2}r$ is the mean number of photons, and
$\varphi$ is the phase of $|\xi\rangle$. The multi-mode splitting corresponds
to a QC described by a unitary $d\times d$ matrix $\bm{U}$. Denoting as
$(\hat{b}_{\textrm{in}})_{j}$ the annihilation operators associated to the
$j$th input mode, see Fig. 1(a), we have
$\hat{b}_{j}=\sum_{k}\bm{U}^{\dagger}_{jk}(\hat{b}_{\textrm{in}})_{k}$. The
squeezed vacuum is inserted in one input port of the network, that we indicate
as port $D$ (ranging from 1 to $d$), while all the other input modes are in
the vacuum state. The total average number of photons in the full Mach-Zehnder
sensor network of Fig. 1(a) is given by
$\bar{n}_{T}=\sum_{j=1}^{d}|\alpha_{j}|^{2}+\bar{n}_{s}$.

The inverse moment matrix $\bm{\mathcal{M}}^{-1}$ can be calculated
analytically, see Appendix A. Below, we report the explicit expression by
assuming the optimal phase-matching conditions $\Im(e^{i\chi_{j}}u_{j})=0$ for
$j=1,...,d$. Here $\Re(x)$ and $\Im(x)$ indicate the real and imaginary part
of $x$, respectively, and $u_{j}=\bm{U}_{Dj}$ ($\sum_{j=1}^{d}|u_{j}|^{2}=1$
being $\bm{U}$ unitary). The condition $\Im(e^{i\chi_{j}}u_{j})=0$ can be
fulfilled by adjusting the relative phase $\chi_{j}=\phi_{j}-\varphi/2$
between the coherent state in mode $a_{j}$ and the squeezed-vacuum. In other
words, the optimal sensing condition is obtained by matching the phases of
each $|\alpha\rangle_{j}$ relative to that of $|\xi\rangle$. At the optimal
working point $\theta_{j}=\pi/2$, we have

$(\bm{\mathcal{M}}^{-1})_{ij}=\frac{|\alpha_{i}|(e^{-2r}-1)|\alpha_{j}|\tilde{u}_{i}\tilde{u}_{j}}{(|\alpha_{i}|^{2}-\tilde{u}_{i}^{2}\bar{n}_{s})(|\alpha_{j}|^{2}-\tilde{u}_{j}^{2}\bar{n}_{s})}+\frac{|\alpha_{j}|^{2}+\tilde{u}_{j}^{2}\bar{n}_{s}}{(|\alpha_{j}|^{2}-\tilde{u}_{j}^{2}\bar{n}_{s})^{2}}\delta_{ij},$
(5)

where $\delta_{ij}$ is the Dirac delta function and
$\tilde{u}_{j}=\Re(e^{i\chi_{j}}u_{j})=\pm|u_{j}|$. We also calculate
analytically the QFIM (see Appendix A):

$(\bm{\mathcal{F}_{\rm
Q}})_{ij}=|\alpha_{i}|(e^{2r}-1)|\alpha_{j}|\tilde{u}_{i}\tilde{u}_{j}+(|\alpha_{j}|^{2}+\tilde{u}_{i}^{2}\bar{n}_{s})\delta_{ij},$
(6)

which is independent from $\bm{\theta}$. We finally notice that Eqs. (5) and
(6) do not depend on the QC transformation $\bm{U}$ as a whole, but only on
the vector $\bm{\tilde{u}}=\\{\tilde{u}_{1},\dots,\tilde{u}_{d}\\}$. Given a
specific QC transformation $\bm{U}$ and the actual input port $D$ in which
$|\xi\rangle$ is inserted, we can identify $2^{d}$ different (non-orthogonal)
vectors $\\{\tilde{u}_{i}\\}_{i=1,\dots,d}$, depending on the choice of sign
for each $\tilde{u}_{i}=\pm|\bm{U}_{Di}|$. We recall that $\bm{U}_{Di}$ is the
element of the matrix $\bm{U}$ at row $D$ and column $i$.

In the separable strategy of Fig. 1(b), the $j$th MZI is fed with a coherent
state $|\alpha_{j}^{\prime}\rangle$ in mode $a_{j}$ and a squeezed vacuum
state $|\xi^{\prime}\rangle$ in mode $b_{j}$, where
$\alpha_{j}^{\prime}=|\alpha_{j}^{\prime}|e^{i\phi_{j}^{\prime}}$ and
$\xi_{j}^{\prime}=r_{j}^{\prime}e^{i\varphi_{j}^{\prime}}$. Under the optimal
condition $\chi_{j}^{\prime}=\phi_{j}^{\prime}-\varphi_{j}^{\prime}=0$
(requiring phase locking between $|\alpha_{j}^{\prime}\rangle$ and
$|\xi^{\prime}\rangle$), we find CavesPRD1981 ; ParisPLA1995 ; PezzePRL2008

$\frac{\Delta^{2}\hat{X}_{j}}{(d\big{\langle}\hat{X}_{j}\big{\rangle}/d\theta_{j})^{2}}=\frac{|\alpha_{j}^{\prime}|^{2}e^{-2r_{j}^{\prime}}+(\bar{n}_{s}^{\prime})_{j}}{[|\alpha_{j}^{\prime}|^{2}-(\bar{n}_{s}^{\prime})_{j}]^{2}},$
(7)

and the quantum Fisher information PezzePRL2008 ; LangPRL2013

$\mathcal{F}_{j}=|\alpha_{j}^{\prime}|^{2}e^{2r_{j}^{\prime}}+(\bar{n}_{s}^{\prime})_{j},$
(8)

where $(\bar{n}_{s}^{\prime})_{j}=\sinh^{2}r_{j}^{\prime}$ is the mean number
of photon in the state $|\xi^{\prime}\rangle$. Differently from the entangled
strategy, the separable strategy uses $d$ squeezed-vacuum states. The total
average number of particles in the separable sensor is thus
$\bar{n}_{T}^{\prime}=\sum_{j=1}^{d}|\alpha_{j}^{\prime}|^{2}+\bar{n}_{s}^{\prime}$,
where $\bar{n}_{s}^{\prime}=\sinh^{2}r_{j}^{\prime}$ is the total average
number of photons in the $d$ squeezed states. When $d=1$, Eqs. (5) and (6)
agree with Eqs. (7) and (8), respectively, and recover a single MZI with
coherent$\otimes$squeezed-vacuum input state CavesPRD1981 ; ParisPLA1995 ;
BarnettEPJD2003 ; PezzePRL2008 ; LangPRL2013 ; Ruo-BercheraPRA2015 ;
SparaciariPRA2016 .

##  IV Optimal sensing configuration

In this section, we study the following problem: given a linear combination of
parameters $\bm{v}\cdot\bm{\theta}$, we want to find the optimal configuration
of the Mach-Zehnder sensor network of Fig. 1(a) that minimizes the phase
uncertainty $\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$, for instance when
using the method of moment as estimation strategy. In other words, for a given
$\bm{v}\in\mathbb{R}^{d}$, we search for

$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},r}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}.$ (9)

We recall that the minimization over the QC transformation $\bm{U}$
corresponds to a minimization over the vector $\bm{\tilde{u}}$ defined above.
Here, we approach Eq. (9) by a direct calculation of the inverse moment matrix
which, in some limits, assumes a convenient form suitable for analytical
optimization. Clearly, Eq. (9) can be also generalized to the quantum Cramer-
Rao bound, namely

$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},r}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}.$ (10)

In Sec. IV.1, we derive upper and lower bounds to Eqs. (9) and (10) that hold
for every $\bm{v}$. We then compare, in Sec. IV.2, the optimized entangled and
separable strategies, under different constraints. In particular, numerical
studies for $d=2$ and $d=3$ show that the optimal parallel strategy overcomes
a corresponding optimal sequential strategy for every $\bm{v}$.

###  IV.1 Bounds and scalings

For fixed $\bar{n}_{T}$ and $\bar{n}_{s}$, it is possible to find upper and
lower bounds to Eq. (9), for every $\bm{v}$:

$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\geq\frac{e^{-2r}}{d\bar{n}_{T}}+\frac{\bar{n}_{s}}{d\bar{n}_{T}^{2}},$
(11)

that holds for $\bar{n}_{T}\gg\bar{n}_{s}$, and

$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\leq\frac{e^{-2r}}{\bar{n}_{T}}+\frac{\bar{n}_{s}\mathcal{W}}{\bar{n}_{T}^{2}},$
(12)

that holds for $\bar{n}_{T}\gg(d+1)\bar{n}_{s}$, where
$\mathcal{W}=d^{3}\sum_{j=1}^{d}v_{j}^{4}$. The above inequalities are derived
in Appendix B. As shown below, the upper bound is tight for
$\bm{v}=\bm{v}_{\rm ave}=(\pm 1,\pm 1,...,\pm 1)/d$ such that $\bm{v}_{\rm
ave}\cdot\bm{\theta}=\sum_{j=1}^{d}\pm\theta_{j}/d$ (that, for brevity, we
indicate as generalized average). The lower bound is tight in the trivial case
when $v_{j}=1/\sqrt{d}$ and $v_{i\neq j}=0$ such that
$\bm{v}\cdot\bm{\theta}=\theta_{j}/\sqrt{d}$ (the factor $\sqrt{d}$ is due to
consistency with the normalization $|\bm{v}|^{2}=1/d$). In this case, the
optimal scheme consists of a single MZI with a squeezed vacuum in one port and
a coherent state of $\bar{n}_{T}-\bar{n}_{s}$ particles in the other port. The
bounds (11) and (12) are characterized by different regimes and scalings.

  * •

Regime $\bar{n}_{T}\gg\bar{n}_{s}e^{2r}\mathcal{W}$. Noticing that
$\mathcal{W}\geq 1$, this regime also implies
$\bar{n}_{T}\gg\bar{n}_{s}e^{2r}$. In this case, the first term in both Eqs.
(11) and (12) dominates over the second one, giving notaEq13

$\frac{e^{-2r}}{d\bar{n}_{T}}\leq\min_{\bm{U},|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\leq\frac{e^{-2r}}{\bar{n}_{T}}.$ (13)

Both terms correspond to a sub-SN uncertainty with prefactor related to the
squeeze parameter $r$.

  * •

Optimal squeezing. We minimize Eqs. (11) and (12) with respect to
$\bar{n}_{s}$, for a fixed $\bar{n}_{T}$. Considering $\bar{n}_{s}\gg 1$ (such
that $e^{2r}\approx 4\bar{n}_{s}$) and taking the derivative with respect to
$\bar{n}_{s}$, one finds

$\frac{1}{d\bar{n}_{T}^{3/2}}\leq\min_{U,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},r}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\leq\frac{\sqrt{\mathcal{W}}}{\bar{n}_{T}^{3/2}}.$ (14)

The optimal value of $\bar{n}_{s}$ minimizing the upper (lower) bound
$\bar{n}_{s}\approx\sqrt{\bar{n}_{T}/4\mathcal{W}}$
($\bar{n}_{s}\approx\sqrt{\bar{n}_{T}/4}$). These values are consistent with
the validity conditions of Eq. (11), namely $\bar{n}_{T}\gg\bar{n}_{s}$.
Furthermore, the upper bound in Eq. (14) holds under the additional condition
$\bar{n}_{s}\gg(d+1)/(4\mathcal{W})$. Overall, both bounds in Eq. (14)
corresponds to a scaling of phase uncertainty faster than the SN.

  * •

Transient Heisenberg scaling for
$(d+1)\bar{n}_{s}\ll\bar{n}_{T}\ll\bar{n}_{s}e^{2r}$. In this regime, which
requires $e^{2r}\gg d+1$, both Eq. (11) and Eq. (12) show a transient
Heisenberg scaling notaEq15 :

$\frac{\bar{n}_{s}}{d\bar{n}_{T}^{2}}\leq\min_{U,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\leq\frac{\bar{n}_{s}\mathcal{W}}{\bar{n}_{T}^{2}}.$ (15)

Equation (15) is understood as $1/\bar{n}_{T}^{2}$ scaling that holds as a
function of $\bar{n}_{T}$ in a restricted regime, and for a fixed value of
$\bar{n}_{s}$.

In analogy to Eqs. (11) and (12), we can find an upper and a lower bound to
Eq. (10):

$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\geq\frac{1}{d[\bar{n}_{T}e^{2r}-\bar{n}_{s}(e^{2r}-1)]},$ (16)

and

$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\leq\frac{e^{-2r}}{\bar{n}_{T}-\bar{n}_{s}}.$ (17)

The bounds holds for every $\bm{v}$ and do not require additional conditions
on $\bar{n}_{T}$ and $\bar{n}_{s}$. The demonstration of the inequalities (17)
and (16) is detailed in Appendix C. Also in this case, the upper bound is
tight for the estimation of $\bm{v}_{\rm ave}\cdot\bm{\theta}$, while the
lower bound is tight for the estimation of a single $\theta_{j}$. We
distinguish different regimes:

  * •

For $\bar{n}_{T}\gg\bar{n}_{s}$, Eqs. (16) and (17) simplify and we have

$\frac{e^{-2r}}{d\bar{n}_{T}}\leq\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\leq\frac{e^{-2r}}{\bar{n}_{T}}.$ (18)

Equation (18) holds, in particular, also in the regime
$\bar{n}_{T}\gg\bar{n}_{s}e^{2r}\mathcal{W}$, where the upper and lower bounds
to Eq. (18) coincide with that of Eq. (13).

  * •

Heisenberg limit for $\bar{n}_{T}\approx 2\bar{n}_{s}$. Taking $\bar{n}_{s}\gg
1$ (so that $\bar{n}_{s}\approx e^{2r}/4$), we can immediately optimize Eqs.
(16) and (17) with respect to $\bar{n}_{s}$, and for fixed $\bar{n}_{T}$:

$\frac{1}{d\bar{n}_{T}^{2}}\leq\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},r}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\leq\frac{1}{\bar{n}_{T}^{2}}.$ (19)

It should be noticed that the value of the squeeze parameter that minimizes
both bounds (16) and (17) is $\bar{n}_{s}=\bar{n}_{T}/2$. This is different
from the value that minimizes the bounds in Eq. (14). The different
optimizations correspond to different scalings, $O(\bar{n}_{T}^{-3/2})$ and
$O(\bar{n}_{T}^{-2})$, respectively.

Figure 2: (a) Optimized sensitivity for the estimation of the generalized
average $\bm{v}_{\rm
ave}\cdot\bm{\theta}=(\pm\theta_{1}\pm\theta_{2}...\pm\theta_{d})/d$, as a
function of the total average number of particles $\bar{n}_{T}$, for fixed
$\bar{n}_{s}$. The thin solid green line is
$\min\Delta^{2}(\bm{\theta}\cdot\bm{v}_{\rm ave})_{\rm emom}$, while the thin
solid red line $\min\Delta^{2}(\bm{\theta}\cdot\bm{v}_{\rm ave})_{\rm eQCR}$.
For simplicity here $\min$ indicates the minimum over $\bm{U}$, and
$|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}$: the minimization if performed
numerically for $d=2$. These numerics are compared with Eq. (20) (thick green
line) and Eq. (21) (thick red line). The thin red line is barely visible due
to the perfect superposition with the thick red line. The dotted vertical line
highlights the point $\bar{n}_{T}=\bar{n}_{s}e^{2r}$. Black lines correspond
to different analytical limit: the dot-dashed line is $e^{-2r}/\bar{n}_{T}$,
the dashed line is $1/\bar{n}_{T}^{3/2}$ and the dotted line is
$\bar{n}_{s}/\bar{n}_{T}^{2}$. We also indicate the SN ($1/\bar{n}_{T}$) and
the HL ($1/\bar{n}_{T}^{2}$). Here $\bar{n}_{s}=100$. Panels (b) and (c) show
$\min_{r}\Delta(\bm{\theta}\cdot\bm{v}_{\rm ave})_{\rm emom}$ and
$\min_{r}\Delta(\bm{\theta}\cdot\bm{v}_{\rm ave})_{\rm eQCR}$ as a function of
$\bar{n}_{T}$. We use $\min_{r}$ as shorthand notation to indicate the
minimization over $\bm{U}$, $|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}$ and over
the squeeze parameter $r$. Dots are numerical results, while lines are
expected analytical behaviours: $\bar{n}_{T}^{-3/2}$ [in panel (b)] and
$\bar{n}_{T}^{-2}$ [in panel (c)].

To summarize, through a set of bounds, we have identified different scalings
and behaviours that characterize Eqs. (9) and (10) in different regimes of
parameters. Below, we demonstrate the saturation of the upper bounds (12) and
(17).

For symmetry reasons, the best estimation of the generalized average phase
$\bm{v}_{\rm
ave}\cdot\bm{\theta}=(\pm\theta_{1}\pm\theta_{2}\pm...\pm\theta_{d})/d$
(corresponding to $|v_{{\rm ave},j}|=1/d$) is obtained when all the coherent
states have the same intensity, namely $|\alpha_{j}|^{2}=\bar{n}_{c},\forall
j$. Moreover, in the regime $\bar{n}_{T}\gg(d+1)\bar{n}_{s}$ it is possible to
prove that the QC satisfying the condition $\bm{\tilde{u}}=\sqrt{d}\bm{v}_{\rm
ave}$ is optimal (see Appendix D). We have

$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}_{\rm
ave}\cdot\bm{\theta})_{\rm
emom}=\frac{e^{-2r}}{\bar{n}_{T}}+\frac{\bar{n}_{s}}{\bar{n}_{T}^{2}},$ (20)

that coincides with the upper bound (12), when noticing that
$\mathcal{W}=d^{3}\sum_{j=1}^{2}v_{{\rm ave},j}^{4}=1$. Equation (20) holds
under the same condition of Eq. (12), namely $\bar{n}_{T}\gg(d+1)\bar{n}_{s}$.
Furthermore, for the uniform QC considered here, the optimization of the
quantum Cramer-Rao bound leads to

$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}_{\rm
ave}\cdot\bm{\theta})_{\rm
eQCR}=\frac{1}{\bar{n}_{T}e^{2r}-\bar{n}_{s}\left(e^{2r}-1\right)},$ (21)

that saturates Eq. (17) for $r\gg 1$. Equation (21) is proved in Appendix D
for $\bar{n}_{T}\gg(d+1)\bar{n}_{s}$, showing that the condition
$\bm{\tilde{u}}=\sqrt{d}\bm{v}_{\rm ave}$ is optimal in this regime. In Fig.
2(a) we plot numerical results for
$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}_{\rm
ave}\cdot\bm{\theta})_{\rm emom}$ (thin green line) and
$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}_{\rm
ave}\cdot\bm{\theta})_{\rm eQCR}$ (thin red line) as a function of
$\bar{n}_{T}$. These are compared with Eq. (20) and (21), shown as thick green
and red lines, respectively. The perfect superposition between the red lines
shows that Eq. (21) holds in all regimes. Equation (20) instead holds for
sufficiently large values of $\bar{n}_{T}$, as expected. In Fig. 2 we also
plot the analytical behaviours Eqs. (13)-(15) in the corresponding different
regimes. In particular, Eqs. (20) and (21) coincide for
$\bar{n}_{T}\gg\bar{n}_{s}e^{2r}$, indicating that the method of moments is an
optimal estimation strategy in that regime.

In Fig. 2(b) and (c) we plot, respectively,
$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},r}\Delta^{2}(\bm{v}_{\rm
ave}\cdot\bm{\theta})_{\rm emom}$ and
$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},r}\Delta^{2}(\bm{v}_{\rm
ave}\cdot\bm{\theta})_{\rm eQCR}$, as a function of $\bar{n}_{T}$. Dots are
numerical results. The solid line in panel (b) [panel (c)] is obtained by
minimizing Eq. (20) [Eq. (21)] with respect to $r$, predicting a sensitivity
$1/\bar{n}_{T}^{3/2}$ [$1/\bar{n}_{T}^{2}$].

Figure 3: Results of a numerical optimization of the sensor network of Fig.
1(a) for $d=2$ (upper row panels) and $d=3$ (lower row), under different
constraints (as detailed in the text). Panels (a) and (b) corresponds to the
optimization under the constraint $\mathcal{C}_{1}$, see Sec. IV.2.1. Panel
(a) plots
$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},\,r}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}$ (solid black line) and
$\min_{(|\alpha_{1}|^{2},r_{1}),...,(|\alpha_{d}|^{2},r_{d})}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}$ (dashed blue line). The dot-dashed orange lines are the upper and lower
bounds of Eq. (14). The gain Eq. (23) is shown in the inset of panel (a) for
$d=2$ and in panel (b) for $d=3$. Panels (c) and (d) corresponds to the
optimization under the constraint $\mathcal{C}_{2}$, see Sec. IV.2.2. Panel
(c) plots $\min_{\bm{U}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$ (solid
black line) and
$\min_{r_{1}^{\prime},...,r_{d}^{\prime}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}$ (dashed blue line). The corresponding gain, Eq. (25), is shown as solid
line in the inset, together with the analytical prediction Eq. (26) (dashed
line). Panel (d) shows $\mathcal{G}_{2}(\bm{v})$ for $d=3$. In panels (c) and
(d) we have set
$|\alpha_{j}|^{2}=|\alpha_{j}^{\prime}|^{2}=(\bar{n}_{T}-\bar{n}_{s})/d$ for
all $j$. Panels (e) and (f) corresponds to the optimization under the
constraint $\mathcal{C}_{3}$, see Sec. IV.2.3. In particular, we consider
$\bar{n}_{s}/\bar{n}_{T}=10^{-4}$ for the entangled strategy and
$(\bar{n}_{s}^{\prime})_{j}/(\bar{n}_{T}^{\prime})_{j}=10^{-4}$ ($j=1,...,d$)
for the separable strategy. Panel (e) plots
$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}$ (solid black line) and
$\min_{|\alpha_{1}^{\prime}|^{2},...,|\alpha_{d}^{\prime}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}$ (dashed blue line). The gain Eq. (33) is shown in the inset of panel
(c) for $d=2$ and in panel (d) for $d=3$. The dot-dashed orange lines are the
lower and upper bounds Eqs. (11) and (12), respectively. In all panels
$\bar{n}_{T}=\bar{n}_{T}^{\prime}=10^{6}$. In panels (c)-(f)
$\bar{n}_{s}=\bar{n}_{s}^{\prime}=10^{2}$.

###  IV.2 Comparison between optimal entangled and separable strategies

In the following, we compare optimal entangled and separable strategies for an
arbitrary linear combination of parameters $\bm{v}\cdot\bm{\theta}$. Let us
define the gain factor

$\mathcal{G}_{\mathcal{C}}=\frac{\min_{\\{|\alpha_{1}^{\prime}|^{2},r_{1}^{\prime},...,|\alpha_{d}^{\prime}|^{2},r_{d}^{\prime}\\}\in\mathcal{C}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}}{\min_{\\{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},r\\}\in\mathcal{C}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}},$ (22)

where $\mathcal{C}$ indicates a common constraint on resources for both cases.
Different constraints are discussed below. We recall that the uncertainties
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm smom}$ and
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$ in Eq. (22) are already
optimized with respect to the relative phase between the squeezed-vacuum
state(s) and the $d$ coherent states.

The minimizations in Eq. (22) are performed numerically for $d=2$ and $d=3$ by
using a variational approach, see Fig. 3. Analytical predictions for any $d$
can be derived in interesting cases. For $d=2$, the QC consists of a
generalized beam-splitter that implements the mode transformation
$\bm{\hat{b}}={U}^{\dagger}\bm{\hat{b}_{\rm in}}$, where
$U=\begin{pmatrix}\tau&\varrho\\\ -\varrho^{*}&\tau^{*}\end{pmatrix}$. Without
loss of generality, we consider $\tau$ and $\varrho$ being real numbers with
$\tau^{2}+\varrho^{2}=1$. The sensor network is given by two MZIs and
$\bm{v}^{T}=(\cos\phi_{v},\sin\phi_{v})/\sqrt{2}$ is expressed as a function
of $\phi_{v}\in[0,2\pi]$. For $d=3$, the sensing scheme consists of three MZIs
and the goal is to estimate a linear combinations of three relative phases,
$\theta_{1}$, $\theta_{2}$ and $\theta_{3}$. The QC is implemented as a
general $3\times 3$ orthogonal transformation. Any such transformation can be
realized, in general, by a sequence of three two-mode beam splitters
ReckPRL1994 . For $d=3$, the vector $\bm{v}$ is parametrized as
$\bm{v}^{T}=(\sin\theta_{v}\cos\phi_{v},\sin\theta_{v}\sin\phi_{v},\cos\phi_{v})/\sqrt{3}$,
as a function of $\phi_{v}\in[0,2\pi]$ and $\theta_{v}\in[0,\pi]$.

Upon imposing different constraints, we observe
$\mathcal{G}_{\mathcal{C}}(\bm{v})\geq 1$ for every $\bm{v}$, see Fig. 3. In
particular, we obtain $\mathcal{G}_{\mathcal{C}}(\bm{v})=1$ when the problem
reduces to the estimation of a single phase [e.g. for $\phi_{v}=0,\pi/2$ and
$d=2$, corresponding to $\bm{v}=(1/\sqrt{d},0)$ and $(0,1/\sqrt{d})$,
respectively]. The maximum gain $\mathcal{G}_{\mathcal{C}}(\bm{v})$ is
obtained for the estimation of the generalized average phase $\bm{v}_{\rm
ave}\cdot\bm{\theta}=\sum_{j=1}^{d}\pm\theta_{d}/d$, that is
$\phi_{v}=\pm\pi/4$ for $d=2$ and $\phi_{v}=\pm\pi/4$,
$\theta_{v}=\arccos(1/\sqrt{3})$, for $d=3$ (notice that the figure shows the
case $v_{j}\geq 0$, while $\mathcal{G}_{\mathcal{C}}(\bm{v})$ is symmetric
under $v_{j}\to-v_{j}$).

####  IV.2.1 $\mathcal{C}_{1}$: same total average number of particles

We impose here the only constraint of having the same total average number of
particles for both strategies: $\bar{n}_{T}=\bar{n}_{T}^{\prime}$. The
entangled strategy is thus optimized over the QC transformation $\bm{U}$, the
average number of particles of each coherent state, $|\alpha_{j}|^{2}$, and
the squeeze parameter $r$. The sequential strategy is instead optimized over
the input states of each independent MZI, namely over each $|\alpha_{j}|^{2}$
and squeeze parameter $r_{j}$, for $j=1,...,d$. The gain factor is thus

$\mathcal{G}_{1}(\bm{v})=\frac{\min_{|\alpha_{1}^{\prime}|^{2},r_{1}^{\prime},...,|\alpha_{d}^{\prime}|^{2},r_{d}^{\prime}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}}{\min_{\bm{U},|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},r}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}}.$ (23)

Numerical results for $d=2$ and $d=3$ are shown in Fig. 3(a) and (b),
respectively.

In the case $\bm{v}=\bm{v}_{\rm ave}$, the optimized entangled strategy is
$\min\Delta^{2}(\bm{v_{\rm ave}}\cdot\bm{\theta})_{\rm
emom}=\bar{n}_{T}^{-3/2}$. The optimized separable strategy is instead
$\min\Delta^{2}(\bm{v_{\rm ave}}\cdot\bm{\theta})_{\rm
smom}=1/d\times(\bar{n}_{T}/d)^{-3/2}=\sqrt{d}\times(\bar{n}_{T})^{-3/2}$,
where $\bar{n}_{T}/d$ is the total average number of photons injected in each
MZI. The highest gain in this case equals

$\mathcal{G}_{1}(\bm{v_{\rm ave}})=\sqrt{d}.$ (24)

A maximum gain $\sqrt{2}$ and $\sqrt{3}$ can be seen in the inset of panels
(a) and (b) of Fig. 3, respectively.

####  IV.2.2 $\mathcal{C}_{2}$: same total average number of particles and
same coherent state intensities

Here, we consider the same total average number of particles for both
strategies, $\bar{n}_{T}=\bar{n}_{T}^{\prime}$, and the same coherent state
intensities, $|\alpha_{j}|^{2}=|\alpha_{j}^{\prime}|^{2}$ for all $j$. These
constraints fix the total average number of photons in the squeezed state(s),
namely
$\bar{n}_{s}=\bar{n}_{T}-\sum_{j=1}^{d}|\alpha_{j}|^{2}=\bar{n}_{s}^{\prime}$.
The entangled strategy uses a single squeezed-vacuum state and we optimize the
multimode splitting in the QC. The separable strategy is instead optimized
with respect to the $d$ squeeze parameters $r_{j}^{\prime}$ ($j=1,\dots,d$)
with $\sum_{j=1}^{d}\sinh^{2}r_{j}^{\prime}=\bar{n}_{s}^{\prime}=\bar{n}_{s}$.
The gain factor writes

$\mathcal{G}_{2}(\bm{v})=\frac{\min_{r_{1}^{\prime},...,r_{d}^{\prime}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}}{\min_{\bm{U}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}}.$ (25)

In the following, for simplicity, we can take $|\alpha_{j}|^{2}=\bar{n}_{c}$
for all $j$. In Fig. 5(c) and (d), we show the results of a numerical
optimizations for $d=2$ and $d=3$, respectively. They agree well with the
analytical prediction

$\mathcal{G}_{2}(\bm{v})=d\,\Bigg{(}\sum_{i=1}^{d}|v_{i}|\Bigg{)}^{2},$ (26)

that can be derived under the conditions $r,r_{j}^{\prime}\gg 1$ and in the
regime $\bar{n}_{c}\gg\bar{n}_{s}e^{2r}$, see Appendix E. We observe
$\mathcal{G}_{2}(\bm{v})\geq 1$ for every $\bm{v}$, with maximum gain equal
for $\bm{v}=\bm{v}_{\rm ave}$.

Let us now focus the discussion on the maximum gain point $\bm{v}=\bm{v}_{\rm
ave}$ and $\bar{n}_{c}\gg\bar{n}_{s}$. In this case, the optimal entangled
strategy is obtained for $\tilde{u}_{j}=1/\sqrt{d}$ and we have

$\min_{\bm{U}}\Delta^{2}(\bm{v}_{\rm ave}\cdot\bm{\theta})_{\rm
emom}=\frac{e^{-2r}\bar{n}_{c}+\bar{n}_{s}/d}{d\bar{n}_{c}^{2}},$ (27)

see Eq. (20) and demonstration in Appendix D. The optimal separable strategy
is obtained using $d$ squeezed-vacuum states having the same squeeze parameter
$r_{j}^{\prime}=r^{\prime}={\rm arcsinh}\sqrt{\bar{n}_{s}/d}$ [since
$\bar{n}_{s}$ is fixed, the intensity of the squeezed-vacuum state in each MZI
is $(\bar{n}_{s}^{\prime})_{j}=\bar{n}_{s}/d$]. As shown in Appendix F, for
$\bar{n}_{c}\gg\bar{n}_{s}$, we obtain

$\min_{r^{\prime}_{1},...,r^{\prime}_{d}}\Delta^{2}(\bm{v}_{\rm
ave}\cdot\bm{\theta})_{\rm
smom}=\frac{\bar{n}_{c}e^{-2r^{\prime}}+\bar{n}_{s}/d}{d\bar{n}_{c}^{2}}.$
(28)

Finally, the gain factor

$\mathcal{G}_{2}(\bm{v}_{\rm
ave})=\frac{\bar{n}_{c}e^{-2r^{\prime}}+\bar{n}_{s}/d}{\bar{n}_{c}e^{-2r}+\bar{n}_{s}/d}$
(29)

is obtained by taking the ratio between Eqs. (27) and (28). It is interesting
to consider different limits of Eq. (29), see also Fig. 4.

Let us consider $d$ fixed, $\bar{n}_{s}/d\gg 1$, such that $r,r^{\prime}\gg
1$, and $d\bar{n}_{c}\gg\bar{n}_{s}e^{2r}$ (which also implies
$d\bar{n}_{c}\gg\bar{n}_{s}e^{2r^{\prime}}$). Notice that these conditions are
fulfilled in Fig. 3(c) and (d). In this case, Eq. (29) simplifies to

$\mathcal{G}_{2}(\bm{v}_{\rm ave})=\frac{e^{-2r^{\prime}}}{e^{-2r}}=d.$ (30)

The gain factor $d$ is obtained taking into account that, for $r,r^{\prime}\gg
1$, we have $e^{-2r^{\prime}}/4\approx\bar{n}_{s}^{\prime}/d$,
$e^{-2r}/4\approx\bar{n}_{s}$. A maximum gain close to $2$ and $3$ can be seen
in panels (c) and (d) of Fig. 3, respectively.

Figure 4: Gain factor $\mathcal{G}_{2}(\bm{v}_{\rm ave})$, Eq. (29), as a
function of $d$, for $\bar{n}_{c}=10^{8}$ and $\bar{n}_{s}=10^{2}$ (dots). The
solid line is $\mathcal{G}_{2}(\bm{v}_{\rm ave})=d$, which is expected for
$d\ll\bar{n}_{s}$; the dashed line is $\mathcal{G}_{2}(\bm{v}_{\rm
ave})=e^{2r}$, for $d\gg\bar{n}_{s}$.

Let us consider the opposite limit, $\bar{n}_{s}/d\ll 1$. The discussion
includes the limit $d\to\infty$ that is peculiar of the multiparameter
problem. First, we notice that, for any fixed value of $\bar{n}_{s}$, in the
limit $d\to\infty$, we have $r^{\prime}={\rm arcsinh}\sqrt{\bar{n}_{s}/d}\to
0$. The separable sensing scheme thus reduces to $d$ MZIs fed with coherent
state in one port and vacuum in the other port. As shown in Eq. (28), we
recover the shot noise limit,
$\min_{r^{\prime}_{1},...,r^{\prime}_{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}\to 1/(d\bar{n}_{c})=1/\bar{n}_{T}$. It is easy to show that the same
sensitivity is also achieved for the optimized quantum Cramér-Rao bound
$\min_{r^{\prime}_{1},...,r^{\prime}_{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
sQCR}$. The situation is completely different for the entangled scheme. In
this case, $r$ remains finite in the limit $d\to\infty$ and, according to Eq.
(27), we have $\min_{\bm{U}}\Delta^{2}(\bm{v}_{\rm ave}\cdot\bm{\theta})_{\rm
emom}=e^{-2r}/(d\bar{n}_{c})=e^{-2r}/\bar{n}_{T}$. The entangled strategy
still achieves a sub-SN sensitivity, with a gain factor

$\mathcal{G}_{2}(\bm{v}_{\rm ave})\to e^{2r},\qquad{\rm for}\,\,d\to\infty.$
(31)

Surprisingly, a finite gain is obtained when a single squeezed-vacuum state is
mixed, by the QC, with a diverging number ($d-1$) of vacuum states
$|0\rangle$. The physical reason for the finite gain Eq. (31) is due to
quantum correlations in the covariance matrix $\bm{\Gamma}$. Under the above
conditions $\bm{\Gamma}$ reads (see Appendix A for the general expression)

$4\bm{\Gamma}=\frac{\bar{n}_{c}(e^{2r}-1)}{d}\begin{pmatrix}1&\ldots&1\\\
\vdots&\ddots&\vdots\\\
1&\ldots&1\end{pmatrix}+\Big{(}\bar{n}_{c}+\frac{\bar{n}_{s}}{d}\Big{)}\begin{pmatrix}1&\ldots&0\\\
\vdots&\ddots&\vdots\\\ 0&\ldots&1\end{pmatrix}.$ (32)

In the limit $d\to\infty$, the prefactor of the all-ones matrix tends to zero
while the prefactor of the identity matrix remains finite. However, the all-
ones matrix is characterized by perfect correlations between all its $d^{2}$
elements, while the identity matrix does not feature any correlation. When
calculating
$\tfrac{(1,1,...,1)^{T}}{\sqrt{d}}4\bm{\Gamma}\tfrac{(1,1,...,1)}{\sqrt{d}}=\bar{n}_{c}e^{2r}+\bar{n}_{s}/d\approx\bar{n}_{c}e^{2r}$
[notice that the vector $(1,1,...,1)/\sqrt{d}$ is normalized to one] the
contributions coming from the two terms in Eq. (32) have the same magnitude
with respect to $d$. The finite mean value of the correlation matrix is
responsible for the finite phase sensitivity and gain.

####  IV.2.3 $\mathcal{C}_{3}$: same total average number of particles and
same total squeezed-vacuum intensities

We consider here $\bar{n}_{T}=\bar{n}_{T}^{\prime}$ and the same total average
number of particles in the squeezed state(s), namely
$\bar{n}_{s}=\bar{n}_{s}^{\prime}=\sum_{j=1}^{d}(\bar{n}_{s}^{\prime})_{j}$.
For the entangled strategy, we fix the ratio $\bar{n}_{s}/\bar{n}_{T}$ [equal
to $10^{-4}$ in Fig. 3(e) and (f)] and thus optimize
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$ over the QC tranformation
$\bm{U}$ and the $d$ coherent state intensities $|\alpha_{1}|^{2}$, …,
$|\alpha_{d}|^{2}$. For the separable strategy, we fix the ratio
$(\bar{n}_{s}^{\prime})_{j}/(\bar{n}_{T}^{\prime})_{j}$ [also equal to
$10^{-4}$ for each $j$, in the figure], where $(\bar{n}_{T}^{\prime})_{j}$ and
$(\bar{n}_{s}^{\prime})_{j}$ indicate the average number of particles in total
and in each squeezed-vacuum state, respectively, as input of the $j$th MZI. We
then optimize the separable strategy over each $|\alpha_{j}^{\prime}|^{2}$.
The gain factor is

$\mathcal{G}_{3}(\bm{v})=\frac{\min_{|\alpha_{1}^{\prime}|^{2},...,|\alpha_{d}^{\prime}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}}{\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}}.$ (33)

Results $d=2$ and $d=3$ are shown in Fig. 5(e) and (f), respectively. It is
particularly interesting that the entangled strategy outperforms the separable
one in this case, $\mathcal{G}_{3}(\bm{v})\geq 1$ for all $\bm{v}$. Indeed, it
clearly shows that the main responsible for the gain in sensitivity is the use
of an optimal entangled multi-mode state in place of a product of squeezed-
vacuum states.

When considering $\bm{v}=\bm{v}_{\rm ave}$, optimal strategies are obtained
taking $|\alpha_{j}^{\prime}|^{2}=|\alpha_{j}|^{2}=\bar{n}_{c}$ for all $j$
and $\tilde{\bm{u}}=1/\sqrt{d}$. In this case, $\mathcal{G}_{3}(\bm{v}_{\rm
ave})$ equals Eq. (29), for $\bar{n}_{c}\gg\bar{n}_{s}$, and the same
considerations as above can be obtained. In particular,
$\mathcal{G}_{3}(\bm{v}_{\rm ave})=d$ for $\bar{n}_{s}\gg d$, while
$\mathcal{G}_{3}(\bm{v}_{\rm ave})=e^{2r}$ for $\bar{n}_{s}\ll d$ (in
particular, in the limit $d\to\infty$).

####  IV.2.4 $\mathcal{C}_{4}$: same total average number of particles and
same squeezed-vacuum strength

Here, we compare the entangled strategy using a single squeezed-vacuum state
$|\xi\rangle$ with the separable strategy using $d$ copies of the same
squeezed-vacuum states $|\xi\rangle$. In other words, here, $r_{j}^{\prime}=r$
for all $j$. For simplicity, we set
$|\alpha_{j}^{\prime}|^{2}=\bar{n}_{c}^{\prime}$ and
$|\alpha_{j}|^{2}=\bar{n}_{c}$ for all $j$. We also consider the same total
average number of particles, $\bar{n}_{T}=\bar{n}_{T}^{\prime}$, where
$\bar{n}_{T}=d\bar{n}_{c}+\bar{n}_{s}$ and
$\bar{n}_{T}^{\prime}=d(\bar{n}_{c}^{\prime}+\bar{n}_{s})$, and recall that
$\bar{n}_{s}=\sinh^{2}r$. In this case, using Eqs. (2) and (7), the
sensitivity achieved with the separable strategy is

$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}=\frac{\bar{n}_{c}^{\prime}e^{-2r}+\bar{n}_{s}}{d(\bar{n}_{c}^{\prime}-\bar{n}_{s})^{2}}$
(34)

for all $\bm{v}$. One of the characteristic features of Eq. (34) is the
divergence at $\bar{n}_{c}^{\prime}=\bar{n}_{s}$ (namely
$\bar{n}_{T}=2d\bar{n}_{s}$) and the saturation of the Cramer-Rao bound,
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
sQCR}=1/[d(\bar{n}_{c}^{\prime}e^{2r}+\bar{n}_{s})]$, for $\bar{n}_{T}\gg
e^{2r}\bar{n}_{s}$ PezzePRL2008 . Instead, the entangled strategy, should be
optimized over the QC transformation, thus giving

$\mathcal{G}_{4}(\bm{v})=\frac{\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}}{\min_{\bm{U}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}}.$ (35)

Let us focus on the case $\bm{v}=\bm{v}_{\rm ave}$ and take the limit
$\bar{n}_{c},\bar{n}_{c}^{\prime}\gg\bar{n}_{s}$ [that, in particular imply
$\bar{n}_{T}\gg(d+1)\bar{n}_{s}$, $\bar{n}_{T}\approx d\bar{n}_{c}$ and
$\bar{n}_{T}^{\prime}\approx d\bar{n}_{c}^{\prime}$]. In this regime,
$\min_{\bm{U}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$ is given by Eq.
(20) which, taking into account Eq. (34), provides

$\mathcal{G}_{4}(\bm{v_{\rm
ave}})=\frac{\bar{n}_{T}e^{-2r}+d\bar{n}_{s}}{\bar{n}_{T}e^{-2r}+\bar{n}_{s}}.$
(36)

If $\bar{n}_{T}\gg de^{2r}\bar{n}_{s}$, Eq. (36) simplifies to

$\mathcal{G}_{4}(\bm{v_{\rm ave}})=1.$ (37)

Although there is no gain in this case, it is still interesting that the
parallel strategy using a single squeezed-vacuum state achieves the same
performance as the sequential strategy using $d$ squeezed states with the same
squeeze parameter.

The opposite regime, $\bar{n}_{T}\ll\bar{n}_{s}e^{2r}$, should be considered
with care. In particular, for
$(d+1)\bar{n}_{s}\ll\bar{n}_{T}\ll\bar{n}_{s}e^{2r}$, the second term in both
the numerator and denominator of Eq. (36) dominates, giving
$\mathcal{G}_{4}(\bm{v_{\rm ave}})=d$. Achieving $\mathcal{G}_{4}(\bm{v_{\rm
ave}})>1$ in this regime is due to the divergence of Eq. (34), as further
discussed in Appendix G. In particular, a calculation of the gain factor
following a numerical optimization of
$\min_{\bm{U}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$ shows that
$\mathcal{G}_{4}(\bm{v_{\rm ave}})$ diverges when $\bar{n}_{T}=2d\bar{n}_{s}$.
This is an artifact due to the use of Eq. (34) outside the regime where it
saturates the Cramer-Rao bound.

##  V Optimal linear combination of phases

In the previous section, we have discussed the optimal configuration of the
sensor network of Fig. 1(a) that maximizes the sensitivity for the estimation
of a fixed a linear combination of phases $\bm{v}\cdot\bm{\theta}$. Here we
consider the opposite problem. Given a specific configuration of the sensor
network, namely a specific QC transformation $U$, coherent state intensities
and squeezed parameter, we discuss the optimal linear combination of phases
$\bm{v}\cdot\bm{\theta}$ that can estimated with the smallest possible
uncertainty.

In the following, we fist provide a general framework for the minimization
problem considered here, namely we introduce the notion of Fisher and
squeezing spectra. We then apply this formalism to the Mach-Zehnder sensor
network of Fig. 1(a) for random choices of the QC and discuss different
sensitivity limits and regimes.

###  V.1 Fisher spectrum and squeezing spectrum

Finding the optimal vector $\bm{v}\in\mathbb{R}^{d}$ that minimizes Eq. (1)
and/or Eq. (3) is solved by calculating the spectrum of the matrices
$\bm{\mathcal{F}_{\rm Q}}$ and $\bm{\mathcal{M}}$, that we indicate as Fisher
and squeezing spectrum, respectively. These spectra contain useful information
regarding the multiparameter problem, in general. We have

$\min_{\bm{v}\in\mathbb{R}^{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}=\frac{1}{\mu_{\rm max}d},$ (38)

where $\mu_{\rm max}$ is the largest eigenvalue of $\bm{\mathcal{M}}$. The
corresponding optimal eigenvector $\bm{v}_{\mu_{\rm max}}$ gives the linear
combination of parameters $\bm{v}_{\mu_{\rm max}}\cdot\bm{\theta}$ that can be
estimated with the smallest possible uncertainty when using the specific
method of moments considered (namely, based on the chosen measurement
observables $\hat{X}_{j}$, probe state and phase encoding transformation).
Following the inequality $\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\leq\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$, we have $f_{\rm
max}\geq\mu_{\rm max}$, where $f_{\rm max}$ is the largest eigenvalue of the
QFIM and satisfies

$\min_{\bm{v}\in\mathbb{R}^{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}=\frac{1}{f_{\rm max}d}.$ (39)

The corresponding optimal eigenvector $\bm{v}_{f_{\rm max}}$ (in general,
$\bm{v}_{f_{\rm max}}\neq\bm{v}_{\mu_{\rm max}}$) gives the linear
combinations of parameters, $\bm{v}_{f_{\rm max}}\cdot\bm{\theta}$, that can
be estimated with the highest possible sensitivity (when optimized over all
generalized output measurements and all possible estimation strategies) for
the given probe state and phase encoding transformation. The demonstration of
Eqs. (38) and (39) is reported in Appendix H. Furthermore, a degeneracy (e.g.
in the squeezing spectrum) reveals independent linear combinations of
parameters that can be estimated with the same sensitivity. Specifically, if
$d_{\mu}$ is the degeneracy of the eigenvalue $\mu$ of $\bm{\mathcal{M}}$,
then the sensitivity $\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}=1/(\mu d)$
is the same for any $\bm{v}$ given by a linear combination of the $d_{\mu}$
orthonormal eigenvectors $\bm{v}_{\mu}^{(1)},...,\bm{v}_{\mu}^{(d_{\mu})}$. In
particular, $\bm{\mathcal{M}}^{-1}$ is defined on the subspace of
$\mathbb{R}^{d}$ generated by a basis of eigenvectors of $\bm{\mathcal{M}}$
corresponding to finite eigenvalues (and similarly for $\bm{\mathcal{F}_{\rm
Q}}^{-1}$).

###  V.2 Random choice of quantum circuit

Here, we consider random choices of the QC and find the corresponding optimal
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$. To be more explicit, we
generate random unitary QC matrices $\bm{U}$ (with uniform deHaar measure) and
calculate the largest eigenvalue $\mu_{\rm max}$ of $\bm{\mathcal{M}}$, Eq.
(5). Furthermore, without loss of generality, we take the same number of
photons in each coherent state, namely $|\alpha_{j}|^{2}=\bar{n}_{c}$ for all
$j$.

Figure 3 summarizes our findings, while different analytical limits are
discussed below. The figure shows $\mathcal{E}_{\rm QC}[1/(\mu_{\rm max}d)]$
(green dots), where $\mathcal{E}_{\rm QC}[...]$ indicates statistical
averaging. For comparison, we also consider $\mathcal{E}_{\rm QC}[1/(f_{\rm
max}d)]$ (red triangles). An analytical upper bound to Eq. (38) can be derived
by taking $\bm{v}=\bm{\tilde{u}}/\sqrt{d}$, giving
$\min_{\bm{v}\in\mathbb{R}^{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\leq\bm{\tilde{u}}^{T}\mathcal{M}(\bm{\tilde{u}})^{-1}\bm{\tilde{u}}/d$.
The inequality is valid for every QC nota3  and numerical calculations reveal
that it is tight in a wide regimes of parameters. In particular, for
$\bar{n}_{T}\gg\bar{n}_{s}$, and taking the statistical average, we find the
simplified expression nota3

$\frac{\mathcal{E}_{\rm
QC}\big{[}\tilde{\bm{u}}^{T}\mathcal{M}(\tilde{\bm{u}})^{-1}\tilde{\bm{u}}\big{]}}{d}=\frac{e^{-2r}}{\bar{n}_{T}}+\frac{\bar{n}_{s}\mathcal{S}}{\bar{n}_{T}^{2}},$
(40)

where $\mathcal{S}\equiv\mathcal{E}_{\rm
QC}[d\sum_{j=1}^{d}\tilde{u}_{j}^{4}]$. Equation (40) is plot as solid black
line in Fig. 5.

Figure 5: Optimized phase uncertainties as a function of $\bar{n}_{T}$.
Symbols show $\mathcal{E}_{\rm QC}[1/(\mu_{\rm max}d)]$ (green dots) and
$\mathcal{E}_{\rm QC}[1/(f_{\rm max}d)]$ (red triangles), where
$\mathcal{E}_{\rm QC}[...]$ indicates statistical averaging over random
choices of the QC. Error bars are root mean square fluctuations. The solid
line is Eq. (40). The dot-dashed line is $e^{-2r}/\bar{n}_{T}$, Eq. (42), the
dashed line is $\sqrt{\mathcal{S}}/\bar{n}_{T}^{3/2}$, Eq. (43), while the
dotted line is $\bar{n}_{s}\mathcal{S}/\bar{n}_{T}^{2}$, Eq. (44). The grey
regions are defined by $1/\bar{n}_{T}$ (SN) and $1/\bar{n}_{T}^{2}$ (HL).
Here, $d=10$, $\bar{n}_{s}=100$ and statistical averaging is obtained over
$10^{4}$ random choices of unitary transformation $\bm{U}$.

Regime $\bar{n}_{T}\gg\bar{n}_{s}e^{2r}$. In this regime, Eq. (5) becomes
nota2

$\bm{\mathcal{M}}^{-1}=\frac{e^{-2r}-1}{\bar{n}_{c}}\bm{\tilde{u}}\bm{\tilde{u}}^{T}+\frac{1}{\bar{n}_{c}}\bm{I}_{d},$
(41)

where $\bm{I}_{d}$ is the $d\times d$ identity matrix. Equation (41) can be
diagonalized straightforwardly: we find $\mu_{\max}=\bar{n}_{c}e^{2r}$, the
corresponding eigenvector being $\bm{v}_{\mu_{\rm
max}}=\bm{\tilde{u}}/\sqrt{d}$. In this case, the upper bound
$\bm{\tilde{u}}^{T}\mathcal{M}(\tilde{u})^{-1}\bm{\tilde{u}})/\sqrt{d}$ is
tight, with the first term in Eq. (40) dominating over the second one. The
optimal sensitivity is

$\min_{\bm{\bm{v}}\in\mathbb{R}^{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}=\frac{e^{-2r}}{\bar{n}_{T}},$ (42)

shown as dot-dashed line in Fig. 5. Equation. (42) holds for any QC. It is
worth noticing that such sensitivity is independent of $\bm{\tilde{u}}$, while
$\bm{v}_{\mu_{\rm max}}$ is independent of the numbers of particles
$\bar{n}_{c}$ and $\bar{n}_{s}$ used. Equation (42) agrees with the numerical
calculations shown in Fig. 5. Below we show that
$f_{\max}=\mu_{\max}=\bar{n}_{c}e^{2r}$ in this regime, with corresponding
eigenvectors $\bm{v}_{f_{\rm max}}=\bm{v}_{\mu_{\rm
max}}=\bm{\tilde{u}}/\sqrt{d}$. The optimal sensitivity predicted by the QFIM
is thus saturated by the practical estimation method given by the method of
moments: in the present limit, $\mathcal{E}_{\rm QC}[1/(\mu_{\rm
max}d)]=\mathcal{E}_{\rm QC}[1/(f_{\rm max}d)]=e^{-r}/\bar{n}_{T}$ with
negligible fluctuations due to random choices of the QC.

Figure 6: (a) $\mathcal{E}_{\rm QC}[\min_{r}1/(\mu_{\rm max}d)]$ as a function
of $\bar{n}_{T}$ and for $d=10$ (dots). The solid line is Eq. (43). The
corresponding optimal values of $\bar{n}_{s}$ are shown in the inset (dots).
There, the solid line is $\bar{n}_{s}=\sqrt{\bar{n}_{T}}/\sqrt{4\mathcal{S}}$.
Panel (b) shows $\mathcal{E}_{\rm QC}[\min_{\bar{n}_{s}}1/(\mu_{\rm
max}d)]\times\bar{n}_{T}^{3/2}$ as a function of $d$ and for
$\bar{n}_{T}=10^{6}$ (dots). The solid line is Eq. (43). The corresponding
optimal values of $\bar{n}_{s}$ are shown in the inset, where the solid line
corresponds to the theoretical prediction. $\mathcal{E}_{\rm QC}[]$ indicates
statistical average, here over $10^{4}$ random QC transformations, where error
bars are root mean square fluctuations.

Optimal squeezing for $\bar{n}_{T}\approx\bar{n}_{s}e^{2r}$. We now optimize
the average number of particles in the squeezed-vacuum state in order to
maximize $\mu_{\rm max}$, for a given total average number of particles
$\bar{n}_{T}$ and QC transformation $\hat{U}$. Such optimization cannot be
performed analytically and we rely on a numerical diagonalization of Eq. (5).
For each QC, we evaluate numerically the maximum eigenvalue $\mu_{\rm max}$ of
the corresponding $\bm{\mathcal{M}}$ and optimize it with respect to
$\bar{n}_{s}$. Numerical results are compared to the analytical optimization
of Eq. (40). For $\bar{n}_{s}\gg 1$ (such that $e^{2r}\approx 4\bar{n}_{s}$)
this predicts

$\min_{r}\frac{\mathcal{E}_{\rm
QC}\big{[}\tilde{u}^{T}\mathcal{M}(\tilde{u})^{-1}\tilde{u}\big{]}}{d}\approx\frac{\sqrt{\mathcal{S}}}{\bar{n}_{T}^{3/2}}$
(43)

for $\bar{n}_{s}\approx\sqrt{\bar{n}_{T}/(4\mathcal{S})}$. In Fig. 6(a) we
plot $\mathcal{E}_{\rm QC}[\min_{\bar{n}_{s}}1/(\mu_{\rm max}d)]$ as a
function of $\bar{n}_{T}$ and for fixed $d$ (dots). The solid line is Eq.
(43). The inset shows the corresponding optimal values of $\bar{n}_{s}$
(dots), the solid line being
$\bar{n}_{s}=\sqrt{\bar{n}_{T}}/\sqrt{4\mathcal{S}}$. In Fig. 6(b) we plot
$\mathcal{E}_{\rm QC}[\min_{\bar{n}_{s}}1/\mu_{\rm
max}]\times\bar{n}_{T}^{3/2}$ as a function of $d$, where the corresponding
optimal values of $\bar{n}_{s}$ are shown in the inset. The numerical results
(dots) are in excellent agreement with Eq. (43) (solid line). Equation (43) is
further shown as dashed line in Fig. 5.

Transient Heisenberg scaling for
$\bar{n}_{s}\ll\bar{n}_{T}\ll\bar{n}_{s}e^{2r}$. In this regime, the first
term in Eq. (40) can be neglected and we obtain

$\frac{\mathcal{E}_{\rm
QC}\big{[}\tilde{u}^{T}\mathcal{M}(\tilde{u})^{-1}\tilde{u}\big{]}}{d}\approx\frac{\bar{n}_{s}\mathcal{S}}{\bar{n}_{T}^{2}}.$
(44)

This predicts a transient Heisenberg scaling, for fixed $\bar{n}_{s}$, with
prefactor approximately given by $\bar{n}_{s}$. This prediction is confirmed
in Fig. 5 where Eq. (44) is shown as the dotted line.

Figure 7: Plot of $\min_{\bar{n}_{s}}1/(f_{\rm max}d)$, averaged over $10^{4}$
random choices of the QC (dots). The solid line is the analytical prediction
$1/\bar{n}_{T}^{2}$.

Quantum Cramer-Rao bound. We now study the QFIM, Eq. (6), for
$|\alpha_{j}|^{2}=\bar{n}_{c}$ for all $j$. In the limit
$\bar{n}_{c}e^{2r}\gg\bar{n}_{s}$ (which is $\bar{n}_{c}\gg 1$ for
sufficiently large values of $\bar{n}_{s}$), Eq. (6) assumes the simple form

$\bm{\mathcal{F}_{\rm
Q}}=\bar{n}_{c}(e^{2r}-1)\bm{\tilde{u}}\bm{\tilde{u}}^{T}+\bar{n}_{c}\bm{I}_{d}.$
(45)

The maximum eigenvalue is $f_{\rm max}=e^{2r}\bar{n}_{c}$ and the
corresponding eigenvector is $\bm{v}_{f_{\rm max}}=\bm{\tilde{u}}/\sqrt{d}$.
For $d\bar{n}_{c}\gg\bar{n}_{s}$, we have $\bar{n}_{T}\approx d\bar{n}_{c}$
and we find
$\min_{\bm{\bm{v}}\in\mathbb{R}^{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}=e^{-2r}/\bar{n}_{T}$. This behaviour holds also for
$\bar{n}_{T}\gg\bar{n}_{s}e^{2r}$, where
$\min_{\bm{\bm{v}}\in\mathbb{R}^{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}=\min_{\bm{\bm{v}}\in\mathbb{R}^{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}$ as given in Eq. (42), see Fig. 5. Furthermore, taking $\bar{n}_{s}\gg
1$ (so that $\bar{n}_{s}\approx e^{2r}/4$), we can optimize $f_{\rm
max}=4\bar{n}_{c}\bar{n}_{s}$ with respect to $\bar{n}_{s}$, for a fixed
$\bar{n}_{T}$: replacing $\bar{n}_{s}=\bar{n}_{T}-\bar{n}_{c}d$ and taking the
derivative with respect to $\bar{n}_{c}$, we find the optimal condition
$d\bar{n}_{c}=\bar{n}_{s}=\bar{n}_{T}/2$. This predicts the saturation of the
Heisenberg limit

$\min_{\bm{\bm{v}}\in\mathbb{R}^{d}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}=\Delta^{2}(\bm{\tilde{u}}\cdot\bm{\theta})_{\rm
eQCR}=\frac{1}{\bar{n}_{T}^{2}}$ (46)

with respect to the total number of particles $\bar{n}_{T}$. In Fig. 7 we show
the the statistical average of $\min_{\bar{n}_{s}}1/(f_{\rm max}d)$ (dots) as
a function of $\bar{n}_{T}$. The solid line is $1/\bar{n}_{T}^{2}$. In Fig. 5
we plot $1/(f_{\rm max}d)$, averaged on random choices of the QC (triangles).
Numerical simulations agree well with analytical predictions in the different
limits. In particular, we see that the QFIM tends to a sub-shot noise scaling
faster than the moment matrix, i.e. for smaller values of $\bar{n}_{T}$.

##  VI Conclusions and Discussion

This work generalizes one of the most important archetype of quantum
interferometry – namely, the single MZI with coherent$\otimes$squeezed-vacuum
light CavesPRD1981 ; ParisPLA1995 ; BarnettEPJD2003 ; PezzePRL2008 ;
LangPRL2013 ; Ruo-BercheraPRA2015 ; SparaciariPRA2016  – to a distributed
sensor network composed by $d$ MZIs, see Fig. 1(a). The multiphase estimation
analysis is based on a method of moments requiring local and independent
photocounting at the output of each MZI. This avoids to recombine the phase-
shifted modes using a second multimode beam splitter TriggianiARXIV ;
OhPRR2020 ; GePRL2018 . The scheme is thus optimal to realize a highly
spatially-separated sensor using a multimode entangled state of a large number
of particles. In particular, our MZI sensor network is characterized by
different regimes reaching $\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}=O(\bar{n}_{T}^{-3/2})$ and the Heisenberg scaling
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}=O(\bar{n}_{T}^{-2})$.

The manuscript is focused on two different problems:

i) We have optimized the full Mach-Zehnder sensor network in order to minimize
the uncertainty, $\Delta^{2}(\bm{v}\cdot\bm{\theta})$, for the estimation of
an arbitrary fixed linear combination of phase shifts
$\bm{v}\cdot\bm{\theta}$. Studying the cases $d=2$ and $d=3$, we have shown
that the optimized parallel strategy, exploiting the mode entanglement
generated by the linear multimode splitting is never surpassed by an optimized
sequential strategy using mode-separable squeezed-vacuum states, when
considering different constraints. For any number of parameters, the maximum
gain of the parallel over the sequential strategy is a factor $d$, for the
estimation of the generalized average phase $\bm{v}_{\rm
ave}\cdot\bm{\theta}=\sum_{j=1}^{d}\pm\theta_{j}/d$. While the literature on
distributed quantum sensing has mainly focused on the estimation of specific
linear combinations of different parameters, the possibility to optimize the
sensor network for the estimation of any desired $\bm{v}\cdot\bm{\theta}$, as
shown here, is generally highly desirable. This possibility has only been
shown in the configurable multimode displacement sensor of Ref. XiaPRL2020
and for multipass phase sensing GebhartARXIV  using a photonic Bell state
ZhaoPRX2021 . Further configurable approaches include the splitting and
multimode recombination of squeezed-vacuum light TriggianiARXIV  and twin-Fock
states GePRL2018 . In our case, the optimized sensing scheme uses local
measurements and avoids the recombination of the states in a second quantum
circuit.

ii) For arbitrary unitary splitting the squeezed vacuum state, we have
identified optimal linear combinations of the $d$ phases that minimize the
uncertainty $\Delta^{2}(\bm{v}\cdot\bm{\theta})$. Typical results and scalings
with the total average number of particles $\bar{n}_{T}$ hold regardless the
random choice of unitary QC transformation.

A further interesting problems raised in the context of multiparameter
estimation is whether a single sensor network allows the estimation of
multiple linear combination of phases at the same time RubioJPA2020 . This
problem is solved here by changing the mode $D$ of the QC where the squeezed
vacuum state is injected, see Fig. 1(a). Indeed, a single $d$-mode QC can be
optimized to estimate $d$ different (e.g. orthogonal) linear combinations
$\bm{v}_{D}\cdot\bm{\theta}$ (with $D=1,...,d$ and
$\bm{v}_{i}\cdot\bm{v}_{j}=\delta_{ij}$) with the same sensitivity: each input
mode $D$ of the QC corresponds to a specific optimal
$\bm{v}_{D}\cdot\bm{\theta}$.

It is also worth comparing here directly with the results of Ref. XiaPRL2020
that considered a sensor network based on the linear splitting of a squeezed-
vacuum state, displacement operations and homodyne measurements on each output
mode. It should be noticed that displacement sensing and phase sensing using a
network of MZIs (as considered here) corresponds, in general, different
parameter-encoding transformations. A precise mapping is obtained within a
(mean-field) Holstein-Primakof approximation, where the mode operator
$\hat{a}_{j}$, see Fig. 1, is replaced by the classical number $|\alpha_{j}|$.
In this case, the $j$th Mach-Zehnder transformation
$e^{-\theta_{j}(\hat{a}^{\dagger}_{j}\hat{b}_{j}-\hat{b}_{j}^{\dagger}\hat{a}_{j})/2}$
reduces to the single-mode displacement operator $\hat{D}(q_{j})$ with
$q_{j}=|\alpha_{j}|/2$. Our results predict
$\Delta^{2}(\bm{v}\cdot\bm{q})_{\rm
eQCR}=\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\times\bar{n}_{c}/4=e^{-2r}/4$ for $|\alpha_{j}|^{2}=\bar{n}_{c}$, in
agreement with Refs. XiaPRL2020 ; ZhuangPRA2018 . Yet, while the phase
uncertainty is characterized by a convenient scaling with $\bar{n}_{T}$,
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm eQCR}=e^{-2r}/\sqrt{\bar{n}_{T}}$,
the displacement uncertainty does not scale with $\bar{n}_{T}$. In this
regime, the quantum Cramer-Rao bound is saturated by the estimation strategy
based on the multiparameter method of moments GessnerNATCOMM2020  based on
local photodetection. Furthermore, our formalism allows to go beyond the
homodyne limit of Refs. XiaPRL2020 ; ZhuangPRA2018  and discuss the phase
sensitivity of a network of MZIs with respect to the total number of particles
$\bar{n}_{T}$ used: this is crucial to discuss scalings of phase variance that
are faster than $O(1/\bar{n}_{T})$.

The results of this work are relevant in current experiments realizing
squeezed vacuum-light and multimode linear splitting transformations
GuoNATPHYS2020 ; XiaPRL2020 ; NokkalaNJP2018 . They pave the way to sensor
networks using multiple MZIs – in both optical and atomic systems – with a
large variety of applications ranging from field and biological sensing,
gravitational wave detection, quantum clocks and inertial measurements.

##  VII Appendix

###  VII.1 Detailed derivation of Eqs. (5) and (6)

In this Appendix, we provide details on the derivation of Eqs. (5) and (6).
Our methods are based on a technique to calculate the QFIM that was outlined
in Ref. GagatsosPRA2016 . It should be noticed, however, that Ref.
GagatsosPRA2016  considered a different sensor network configuration: the
generalization to a network of MZIs is not straightforward and requires
additional algebraic work. Furthermore, our derivation corrects some flaws
which are present in Ref. GagatsosPRA2016  and that led to an incorrect final
expression (see discussion below). Alternative approaches to calculate the
QFIM of Gaussian states have been also considered, see Refs. JiangPRA2014 ;
BanchiPRL2015 ; NicholsPRA2018 ; OhPRR2020 .

####  VII.1.1 Preliminary definitions

We consider the general case where a product of $d^{\prime}$ single-mode
squeezed states is sent to a passive linear network $\hat{A}^{\dagger}$ and
transformed according to

$|\Psi\rangle=\hat{A}^{\dagger}\bigotimes_{k=1}^{d^{\prime}}|\beta_{k},\xi_{k}\rangle.$
(47)

Here, $|\beta_{k},\xi_{k}\rangle$ is the single-mode displaced-squeezed state
in the mode $k$: $\beta_{k}$ is the coherent amplitude of the state and
$\xi_{k}=r_{k}e^{i\varphi_{k}}$ its squeeze parameter. In the following, we
will assume that $\hat{A}^{\dagger}$ is a Gaussian unitary, that is, a unitary
operator which transforms Gaussian states into Gaussian states. Because
$\hat{A}^{\dagger}$ is also a passive, i.e. particle-number preserving,
transformation, if we set
$\hat{\bm{c}}=(\hat{c}_{1},\dots,\hat{c}_{d^{\prime}})^{T}$, a relation
$\hat{A}\hat{\bm{c}}\hat{A}^{\dagger}=\hat{\bm{c}}^{\prime}=\bm{\mathcal{A}}\hat{\bm{c}}$
must hold, with $\bm{\mathcal{A}}$ a unitary matrix (similarly,
$\hat{A}\hat{\bm{c}}^{\dagger}\hat{A}^{\dagger}=\bm{\mathcal{A}}^{*}\hat{\bm{c}}^{\dagger}$).
Here, $\hat{c}^{\dagger}_{k}$ and $\hat{c}_{k}$ are bosonic creation and
annihilation operators, respectively,

We recall the definition of the Q-function for the state $|\Phi\rangle$ of a
generic $d^{\prime}$-mode system:

$Q(\bm{\alpha})=\frac{|\langle\bm{\alpha}|\Phi\rangle|^{2}}{\pi},$ (48)

where $|\bm{\alpha}\rangle=\bigotimes_{k=1}^{d^{\prime}}|\alpha_{k}\rangle$,
$|\alpha_{k}\rangle$ being an arbitrary single-mode coherent state in mode
$k$. From the point of view of the Q-function, a transformation
$\hat{A}^{\dagger}|\Phi\rangle$ is equivalent to a transformation
$\hat{A}|\bm{\alpha}\rangle$ of the coherent states; moreover, it is a well
know property that a Gaussian passive transformation $\hat{A}$ sends a product
of coherent states into another product of coherent states, in particular:
$\hat{A}|\bm{\alpha}\rangle=|\bm{\alpha}^{\prime}\rangle=|\bm{A}\bm{\alpha}\rangle$.
As the notation just used suggests, and according to Ref. GagatsosPRA2016 ,
$\bm{A}$ is the matrix that implements the transformation
$\bm{\alpha}^{\prime}=\bm{A}\bm{\alpha}$ of the amplitudes of the coherent
states associated with the Q-function. It is possible to show that this matrix
is the same as the one which describes the transformation of the annihilation
operators implemented by $\hat{A}$, that is $\hat{A}^{\dagger}\hskip
2.84544pt\bm{\hat{c}}\hskip
2.84544pt\hat{A}=\bm{\hat{c}}^{\prime}=\bm{A}\bm{\hat{c}}$, and,
correspondingly, the hermitian conjugate of the matrix which describes the
transformation of the annihilation operators implemented by
$\hat{A}^{\dagger}$, which we have denoted as $\bm{\mathcal{A}}$ above. The
relation $\bm{\mathcal{A}}^{\dagger}=\bm{A}$ will be frequently used in what
follows.

Using the Q-function representation of the states, (47), Ref. GagatsosPRA2016
showed that

$\displaystyle\langle\hat{n}_{i}\rangle$
$\displaystyle=-1+\partial_{i}\partial_{i}^{*}G(\bm{\mu})\Bigr{\rvert}_{\bm{\mu}=0},$
(49) $\displaystyle\langle\hat{n}_{i}\hat{n}_{j}\rangle$
$\displaystyle=\left[\partial_{i}\partial_{i}^{*}\partial_{j}\partial_{j}^{*}-(1+\delta_{ij})\partial_{i}\partial_{i}^{*}-\partial_{j}\partial_{j}^{*}\right]G(\bm{\mu})\Bigr{\rvert}_{\bm{\mu}=0}+1,$
(50)

where
$\bm{\mu}=(\lambda_{1},\dots,\lambda_{d^{\prime}},\lambda_{1}^{*},\dots,\lambda_{d^{\prime}}^{*})^{T}$
is an arbitrary $2d^{\prime}$-dimensional complex vector, $\partial_{i}$,
$\partial_{i}^{*}$ are shorthand notation for $\partial/\partial\lambda_{i}$,
$\partial/\partial\lambda_{i}^{*}$,
$\hat{n}_{i}=\hat{c}^{\dagger}_{i}\hat{c}_{i}$ and the expectation values are
evaluated in state $|\Psi\rangle$. We have $G(\bm{\mu})=e^{\Delta}$, where

$\Delta\equiv\displaystyle\frac{1}{4}\left(\bm{\nu_{b}}^{\dagger}\bm{M}^{-1}\bm{\mu}+\bm{\mu}^{\dagger}\bm{M}^{-1}\bm{\nu_{b}}+\bm{\mu}^{\dagger}\bm{M}^{-1}\bm{\mu}\right),$
(51)
$\bm{\nu_{b}}=(b_{1},\dots,b_{d^{\prime}},b_{1}^{*},\dots,b_{d^{\prime}}^{*})^{T},$
(52)
$b_{j}=\sum_{k}\bm{A}^{\dagger}_{jk}(\beta_{k}+\beta_{k}^{*}e^{i\varphi_{k}}\tanh{r_{k}}),$
(53) $\bm{M}^{-1}=2\begin{pmatrix}\bm{E}&-\bm{NE}^{T}\\\
-\bm{N}^{\dagger}\bm{E}&\bm{E}^{T}\end{pmatrix},$ (54)
$\bm{N}=\bm{A}^{\dagger}\bm{D}\bm{A}^{*},$ (55)
$\bm{E}=\bm{A}^{\dagger}\bm{C}\bm{A},$ (56)
$\bm{C}_{jk}=\delta_{jk}\cosh^{2}r_{k},$ (57)

and

$\bm{D}_{jk}=\delta_{jk}e^{i\varphi_{k}}\tanh{r_{k}}.$ (58)

Using Eq. (54), we find

$\displaystyle\Delta=\frac{1}{4}\left[\begin{pmatrix}\bm{b}^{*}&\bm{b}\end{pmatrix}\cdot
2\begin{pmatrix}\bm{E}&-\bm{NE}^{T}\\\
-\bm{N}^{\dagger}\bm{E}&\bm{E}^{T}\end{pmatrix}\begin{pmatrix}\bm{\lambda}\\\
\bm{\lambda}^{*}\end{pmatrix}+\begin{pmatrix}\bm{\lambda}^{*}&\bm{\lambda}\end{pmatrix}\cdot
2\begin{pmatrix}\bm{E}&-\bm{NE}^{T}\\\
-\bm{N}^{\dagger}\bm{E}&\bm{E}^{T}\end{pmatrix}\begin{pmatrix}\bm{b}\\\
\bm{b}^{*}\end{pmatrix}+\begin{pmatrix}\bm{\lambda}^{*}&\bm{\lambda}\end{pmatrix}\cdot
2\begin{pmatrix}\bm{E}&-\bm{NE}^{T}\\\
-\bm{N}^{\dagger}\bm{E}&\bm{E}^{T}\end{pmatrix}\begin{pmatrix}\bm{\lambda}\\\
\bm{\lambda}^{*}\end{pmatrix}\right],$

which we can write more explicitly as

$\displaystyle\Delta$
$\displaystyle=\frac{1}{2}\left[\bm{b}^{*}\cdot\left(\bm{E}\bm{\lambda}-\bm{NE}^{T}\bm{\lambda}^{*}\right)+\bm{b}\cdot\left(-\bm{N}^{\dagger}\bm{E}\bm{\lambda}+\bm{E}^{T}\bm{\lambda}^{*}\right)+\bm{\lambda}^{*}\cdot\left(\bm{E}\bm{b}-\bm{NE}^{T}\bm{b}^{*}\right)+\bm{\lambda}\cdot\left(-\bm{N}^{\dagger}\bm{E}\bm{b}+\bm{E}^{T}\bm{b}^{*}\right)\right.$
$\displaystyle\left.\qquad+\bm{\lambda}^{*}\cdot\left(\bm{E}\bm{\lambda}-\bm{NE}^{T}\bm{\lambda}^{*}\right)+\bm{\lambda}\cdot\left(-\bm{N}^{\dagger}\bm{E}\bm{\lambda}+\bm{E}^{T}\bm{\lambda}^{*}\right)\right].$

We then calculate the first and second partial derivatives of $\Delta$ with
respect to $\lambda_{i}$ and $\lambda_{i}^{*}$:

$\displaystyle\begin{aligned}
\partial_{i}\Delta&=\frac{1}{2}\sum_{k}\left[-\left(\left(\bm{N}^{\dagger}\bm{E}\right)_{ik}+\left(\bm{N}^{\dagger}\bm{E}\right)_{ki}\right)(\bm{\lambda}_{k}+\bm{b}_{k})+2\bm{E}_{ki}(\bm{\lambda}_{k}^{*}+\bm{b}_{k}^{*})\right]=\sum_{k}\bm{E}^{*}_{ik}(\bm{\lambda}_{k}^{*}+\bm{b}_{k}^{*})-(\bm{EN})^{*}(\bm{\lambda}_{k}+\bm{b}_{k}),\end{aligned}$
(59) $\displaystyle\begin{aligned}
\partial_{i}^{*}\Delta&=\frac{1}{2}\sum_{k}\left[2\bm{E}_{ik}(\bm{\lambda}_{k}+\bm{b}_{k})-\left(\left(\bm{NE}^{T}\right)_{ik}+\left(\bm{NE}^{T}\right)_{ki}\right)(\bm{\lambda}_{k}^{*}+\bm{b}_{k}^{*})\right]=\sum_{k}\bm{E}_{ik}(\bm{\lambda}_{k}+\bm{b}_{k})-(\bm{EN})(\bm{\lambda}_{k}^{*}+\bm{b}_{k}^{*}),\end{aligned}$
(60) $\displaystyle\partial_{i}^{*}\partial_{j}\Delta=\bm{E}_{ij},$ (61)
$\displaystyle\partial_{i}\partial_{j}^{*}\Delta=\bm{E}_{ji}=\bm{E}^{*}_{ij},$
(62)
$\displaystyle\partial_{i}^{*}\partial_{j}^{*}\Delta=-\frac{1}{2}\left(\left(\bm{NE}^{T}\right)_{ij}+\left(\bm{NE}^{T}\right)_{ji}\right)=-(\bm{EN})_{ij},$
(63)
$\displaystyle\partial_{i}\partial_{j}\Delta=-\frac{1}{2}\left(\left(\bm{N}^{\dagger}\bm{E}\right)_{ij}+\left(\bm{N}^{\dagger}\bm{E}\right)_{ji}\right)=-(\bm{EN})^{*}_{ij}.$
(64)

To derive the above equations, we have used the following relations:

$\displaystyle\bm{E}=\bm{E}^{\dagger}\implies\bm{E}^{T}=\bm{E}^{*},$
$\displaystyle\bm{N}=\bm{N}^{T}\implies\bm{N}^{\dagger}=\bm{N}^{*},$
$\displaystyle(\bm{EN})^{T}=\bm{EN},$
$\displaystyle\bm{NE}^{T}=(\bm{EN})^{T}=\bm{EN},$
$\displaystyle\bm{N}^{\dagger}\bm{E}=(\bm{EN})^{\dagger}=\left((\bm{EN})^{T}\right)^{*}=(\bm{EN})^{*}.$

We are now ready to work out the partial derivatives of $G(\bm{\mu})$ that
appear in Eqs. (49) and (50):

$\displaystyle\partial_{i}\partial_{i}^{*}G(\bm{\mu})=(\partial_{i}\partial_{i}^{*}\Delta)e^{\Delta}+(\partial_{i}\Delta)(\partial_{i}^{*}\Delta)e^{\Delta}=\left(\bm{E}_{ii}+(\partial_{i}\Delta)(\partial_{i}^{*}\Delta)\right)e^{\Delta},$

and

$\displaystyle\partial_{i}\partial_{i}^{*}\partial_{j}\partial_{j}^{*}G(\bm{\mu})=$
$\displaystyle\partial_{i}\left[\left(\bm{E}_{jj}+(\partial_{j}\Delta)(\partial_{j}^{*}\Delta)\right)(\partial_{i}^{*}\Delta)e^{\Delta}+\left((\partial_{i}^{*}\partial_{j}\Delta)(\partial_{j}^{*}\Delta)+(\partial_{j}\Delta)(\partial_{i}^{*}\partial_{j}^{*}\Delta)\right)e^{\Delta}\right]$
$\displaystyle=$
$\displaystyle\left[(\bm{EN})_{ij}(\bm{EN})^{*}_{ij}-(\bm{EN})^{*}_{ij}(\partial_{i}^{*}\Delta)(\partial_{j}^{*}\Delta)-(\bm{EN})_{ij}(\partial_{i}\Delta)(\partial_{j}\Delta)+\bm{E}_{ij}\bm{E}^{*}_{ij}+\bm{E}_{ij}(\partial_{i}\Delta)(\partial_{j}^{*}\Delta)+\bm{E}^{*}_{ij}(\partial_{i}^{*}\Delta)(\partial_{j}\Delta)\right.$
$\displaystyle\left.+\left(\bm{E}_{ii}+(\partial_{i}\Delta)(\partial_{i}^{*}\Delta)\right)\left(\bm{E}_{jj}+(\partial_{j}\Delta)(\partial_{j}^{*}\Delta)\right)\right]e^{\Delta}.$

By evaluating the derivatives at $\bm{\mu}=0$, we get

$\partial_{i}\partial_{i}^{*}G(\bm{\mu})\Bigr{\rvert}_{\bm{\mu}=0}=\bm{E}_{ii}+(\partial_{i}\Delta)_{0}(\partial_{i}^{*}\Delta)_{0},$

where $(\partial_{i}\Delta)_{0}$ is shorthand notation for
$(\partial_{i}\Delta)\bigr{\rvert}_{\bm{\mu}=0}$, and

$\displaystyle\partial_{i}\partial_{i}^{*}\partial_{j}\partial_{j}^{*}G(\bm{\mu})\Bigr{\rvert}_{\bm{\mu}=0}=$
$\displaystyle(\bm{EN})_{ij}(\bm{EN})^{*}_{ij}-(\bm{EN})^{*}_{ij}(\partial_{i}^{*}\Delta)_{0}(\partial_{j}^{*}\Delta)_{0}-(\bm{EN})_{ij}(\partial_{i}\Delta)_{0}(\partial_{j}\Delta)_{0}+\bm{E}_{ij}\bm{E}^{*}_{ij}+\bm{E}_{ij}(\partial_{i}\Delta)_{0}(\partial_{j}^{*}\Delta)_{0}+\bm{E}^{*}_{ij}(\partial_{i}^{*}\Delta)_{0}(\partial_{j}\Delta)_{0}$
$\displaystyle+\left(\bm{E}_{ii}+(\partial_{i}\Delta)_{0}(\partial_{i}^{*}\Delta)_{0}\right)\left(\bm{E}_{jj}+(\partial_{j}\Delta)_{0}(\partial_{j}^{*}\Delta)_{0}\right).$

Finally, going back to Eqs. (49) and (50), we get

$\displaystyle\langle\hat{n}_{i}\rangle=-1+\partial_{i}\partial_{i}^{*}G(\bm{\mu})\Bigr{\rvert}_{\bm{\mu}=0}=-1+\bm{E}_{ii}+(\partial_{i}\Delta)_{0}(\partial_{i}^{*}\Delta)_{0}$
(65)

and

$\displaystyle\langle\hat{n}_{i}\hat{n}_{j}\rangle=$
$\displaystyle\left[\partial_{i}\partial_{i}^{*}\partial_{j}\partial_{j}^{*}-(1+\delta_{ij})\partial_{i}\partial_{i}^{*}-\partial_{j}\partial_{j}^{*}\right]G(\bm{\mu})\Bigr{\rvert}_{\bm{\mu}=0}+1$
$\displaystyle=$
$\displaystyle\partial_{i}\partial_{i}^{*}\partial_{j}\partial_{j}^{*}G(\bm{\mu})\Bigr{\rvert}_{\bm{\mu}=0}-(1+\delta_{ij})(\langle\hat{n}_{i}\rangle+1)-(\langle\hat{n}_{j}\rangle+1)+1$
$\displaystyle=$
$\displaystyle(\bm{EN})_{ij}(\bm{EN})^{*}_{ij}-(\bm{EN})^{*}_{ij}(\partial_{i}^{*}\Delta)_{0}(\partial_{j}^{*}\Delta)_{0}-(\bm{EN})_{ij}(\partial_{i}\Delta)_{0}(\partial_{j}\Delta)_{0}+\bm{E}_{ij}\bm{E}^{*}_{ij}+\bm{E}_{ij}(\partial_{i}\Delta)_{0}(\partial_{j}^{*}\Delta)_{0}+\bm{E}^{*}_{ij}(\partial_{i}^{*}\Delta)_{0}(\partial_{j}\Delta)_{0}$
$\displaystyle+\langle\hat{n}_{i}\rangle\langle\hat{n}_{j}\rangle-\delta_{ij}\left(\bm{E}_{ii}+(\partial_{i}\Delta)_{0}(\partial_{i}^{*}\Delta)_{0}\right).$
(66)

At this point, we introduce the $d^{\prime}\times d^{\prime}$ matrix $\bm{h}$,
with elements

$\bm{h}_{ij}=\langle\hat{n}_{i}\hat{n}_{j}\rangle-\langle\hat{n}_{i}\rangle\langle\hat{n}_{j}\rangle,$
(67)

whose expression, taking into account the above equations for
$\langle\hat{n}_{i}\rangle$ and $\langle\hat{n}_{i}\hat{n}_{j}\rangle$, can be
immediately derived:

$\displaystyle\bm{h}_{ij}=$
$\displaystyle(\bm{EN})_{ij}(\bm{EN})^{*}_{ij}-(\bm{EN})^{*}_{ij}(\partial_{i}^{*}\Delta)_{0}(\partial_{j}^{*}\Delta)_{0}-(\bm{EN})_{ij}(\partial_{i}\Delta)_{0}(\partial_{j}\Delta)_{0}+\bm{E}_{ij}\bm{E}^{*}_{ij}+\bm{E}_{ij}(\partial_{i}\Delta)_{0}(\partial_{j}^{*}\Delta)_{0}+\bm{E}^{*}_{ij}(\partial_{i}^{*}\Delta)_{0}(\partial_{j}\Delta)_{0}$
$\displaystyle-\delta_{ij}\left(\bm{E}_{ii}+(\partial_{i}\Delta)_{0}(\partial_{i}^{*}\Delta)_{0}\right).$

This equation can be rewritten in a compact form by introducing the vector

$\bm{\gamma}_{i}\equiv(\partial_{i}\Delta)_{0}.\\\ $

Notice that $(\gamma_{i})^{*}=(\partial_{i}^{*}\Delta)_{0}$. From Eqs. (59)
and (60), we get
$\bm{\gamma}_{i}=\sum_{k}\bm{E}^{*}_{ik}b_{k}^{*}-(\bm{EN})^{*}_{ik}b_{k}$,
namely $\bm{\gamma}=\bm{E}^{*}\bm{b}^{*}-(\bm{EN})^{*}\bm{b}$. Finally, making
use of vector $\bm{\gamma}$, and of the Hadamard entrywise product $\circ$, we
can rewrite $\bm{h}$ in the compact form:

$\displaystyle\bm{h}=$
$\displaystyle\bm{EN}\circ(\bm{EN})^{*}-\bm{EN}\circ\bm{\gamma}\bm{\gamma}^{T}-(\bm{EN})^{*}\circ\left(\bm{\gamma}\bm{\gamma}^{T}\right)^{*}+\bm{E}\circ\bm{E}^{*}+\bm{E}\circ\bm{\gamma}\bm{\gamma}^{\dagger}+\bm{E}^{*}\circ\left(\bm{\gamma}\bm{\gamma}^{\dagger}\right)^{*}-\left(\bm{E}+\bm{\gamma}\bm{\gamma}^{\dagger}\right)\circ\bm{I}.$
(68)

It is evident that matrix $\bm{h}$ is real and symmetric. An expression
similar to Eq. (68) was derived in Ref. GagatsosPRA2016 , see Eq. (14) in that
reference. There are however important differences with respect to Eq. (68)
due to flaws in the derivation reported in Ref. GagatsosPRA2016 .

####  VII.1.2 Quantum Fisher information matrix, Eq. (6)

Using Eq. (68), we now show how to express the QFIM in terms of the matrices
$\bm{E}$, $\bm{EN}$, etc., which were introduced in the previous paragraph and
will here be evaluated for the specific case of the Mach-Zehnder sensor
network of Fig. 1. In our sensing scheme, the initial state is given by

$|\Psi_{\rm
in}\rangle=\left(|\alpha_{1}\rangle\otimes\dots\otimes|\alpha_{d}\rangle\right)\otimes\left(|0\rangle\otimes\dots\otimes|\xi\rangle\otimes\dots\otimes|0\rangle\right).$

It is a product state of coherent states in modes $a_{1},...,a_{d}$, a
squeezed-vacuum state $|\xi\rangle$ in mode $(b_{\rm in})_{D}$ and the vacuum
$|0\rangle$ in modes $(b_{\rm in})_{j}$ for $j=1,...,d$ and $j\neq D$. This
initial state should be compared with the product state in Eq. (47). In order
to facilitate the identification of the two cases, we can set
$\hat{a}_{j}\equiv\hat{c}_{j}$ and $(\hat{b}_{\rm in})_{j}\equiv\hat{c}_{j+d}$
$(j=1,\dots,d)$, thus introducing a more homogeneous notation valid for all of
the $2d=d^{\prime}$ input modes of the sensing apparatus. We then identify the
$2d\times 2d$ unitary matrix corresponding to the mode transformation
performed by the QC as

$\bm{U}_{\rm QC}^{\dagger}=\begin{pmatrix}\bm{I}_{d}&0\\\
0&\bm{U}^{\dagger}\end{pmatrix}.$

The $d\times d$ identity matrix $\bm{I}_{d}$ describes the action of the QC on
the coherent states, while
$\hat{b}_{j}=\sum_{k}(\bm{U}^{\dagger})_{jk}(\hat{b}_{\rm in})_{k}$
$(j=1,\dots,d)$, $\bm{U}^{\dagger}$ being a unitary $d\times d$ matrix. We
denote as $|\Psi_{0}\rangle$ the output state of the QC:
$|\Psi_{0}\rangle=\hat{U}_{\rm QC}^{\dagger}|\Psi_{\rm in}\rangle$.

The phases $\theta_{1},\dots,\theta_{d}$ to be estimated are encoded in
$|\Psi_{0}\rangle$ through the unitary transformation
$\otimes_{j=1}^{d}e^{-i\theta_{j}(\hat{J}_{y})_{j}}$, where the transformation
$e^{-i\theta_{j}(\hat{J}_{y})_{j}}$ identifies the $j$th MZI in the network,
with
$\hat{H}_{j}=(\hat{J}_{y})_{j}=(\hat{a}_{j}^{\dagger}\hat{b}_{j}-\hat{b}_{j}^{\dagger}\hat{a}_{j})/2i$.
This is equivalent to the phases being encoded in the state
$|\Psi\rangle=\otimes_{j=1}^{d}e^{-i\frac{\pi}{2}(\hat{J}_{x})_{j}}|\Psi_{0}\rangle$,
where
$(\hat{J}_{x})_{j}=(\hat{a}_{j}^{\dagger}\hat{b}_{j}+\hat{b}_{j}^{\dagger}\hat{a}_{j})/2$,
through the unitary transformation
$\otimes_{j=1}^{d}e^{-i\theta_{j}(\hat{J}_{z})_{j}}$, with
$(\hat{J}_{z})_{j}=(\hat{a}_{j}^{\dagger}\hat{a}_{j}-\hat{b}_{j}^{\dagger}\hat{b}_{j})/2$.
This alternative formulation is more convenient here. Overall, $|\Psi\rangle$
can be expressed as

$|\Psi\rangle=\left(\otimes_{j=1}^{d}e^{-i\frac{\pi}{2}(\hat{J}_{x})_{j}}\right)\hat{U}^{\dagger}_{\rm
QC}|\Psi_{\rm in}\rangle.$

This equation should be compared with Eq. (47): the identification
$\hat{A}^{\dagger}=(\otimes_{j=1}^{d}e^{-i\frac{\pi}{2}(\hat{J}_{x})_{j}})\hat{U}^{\dagger}_{\rm
QC}$ is straightforward. The action of $\hat{A}^{\dagger}$ on the annihilation
operators of the input modes can be represented by the relation
$\hat{c}_{j}^{\prime}=\sum_{k}\bm{\mathcal{A}}_{jk}\hat{c}_{k}$
($j=1,\dots,2d$) with

$\displaystyle\bm{\mathcal{A}}=$
$\displaystyle\frac{1}{\sqrt{2}}\begin{pmatrix}\bm{I}_{d}&-i\bm{I}_{d}\\\
-i\bm{I}_{d}&\bm{I}_{d}\end{pmatrix}\begin{pmatrix}\bm{I}_{d}&0\\\
0&\bm{U}^{\dagger}\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}\bm{I}_{d}&-i\bm{U}^{\dagger}\\\
-i\bm{I}_{d}&\bm{U}^{\dagger}\end{pmatrix}.$ (69)

Notice that the matrix
$\frac{1}{\sqrt{2}}\begin{pmatrix}\bm{I}_{d}&-i\bm{I}_{d}\\\
-i\bm{I}_{d}&\bm{I}_{d}\end{pmatrix}$ in the product above describes an array
of balanced beam splitters working in parallel, corresponding to
$\otimes_{j=1}^{d}e^{-i\frac{\pi}{2}(\hat{J}_{x})_{j}}$. Under the hypotheses
of a pure state $|\Psi\rangle$ for the system and of a phase-imprinting
transformation of the form
$\otimes_{j=1}^{d}e^{-i\theta_{j}(\hat{J}_{z})_{j}}$ (see above), the QFIM is
given by

$\left(\bm{\mathcal{F}_{\rm
Q}}\right)_{ij}=4\left(\langle\Psi|(\hat{J}_{z})_{i}(\hat{J}_{z})_{j}|\Psi\rangle-\langle\Psi|(\hat{J}_{z})_{i}|\Psi\rangle\langle\Psi|(\hat{J}_{z})_{j}|\Psi\rangle\right)\hskip
14.22636pt(i,j=1,\dots,d).$

Since
$(\hat{J}_{z})_{j}=(\hat{a}_{j}^{\dagger}\hat{a}_{j}-\hat{b}_{j}^{\dagger}\hat{b}_{j})/2=(\hat{c}_{j}^{\dagger}\hat{c}_{j}-\hat{c}_{j+d}^{\dagger}\hat{c}_{j+d})/2=(\hat{n}_{j}-\hat{n}_{j+d})/2$,
after simple calculations we get to express $\left(\bm{\mathcal{F}_{\rm
Q}}\right)_{ij}$ as

$\left(\bm{\mathcal{F}_{\rm
Q}}\right)_{ij}=\bm{h}_{i,j}+\bm{h}_{i+d,j+d}-\bm{h}_{i,j+d}-\bm{h}_{i+d,j},\hskip
14.22636pt(i,j=1,\dots,d)$ (70)

where the expression of $\bm{h}_{ij}$,
$\bm{h}_{ij}=\langle\Psi|\hat{n}_{i}\hat{n}_{j}|\Psi\rangle-\langle\Psi|\hat{n}_{i}|\Psi\rangle\langle\Psi|\hat{n}_{j}|\Psi\rangle$,
is exactly the one already given in Eq. (67). Notice that
$\bm{\mathcal{F}_{\rm Q}}$ is a $d\times d$ matrix, whose elements, according
to Eq. (70), can be obtained as combinations of the elements of $\bm{h}$, a
$2d\times 2d$ matrix.

Equation (68) from the previous paragraph expresses $\bm{h}$ in terms of the
two matrices $\bm{E}$ and $\bm{EN}$ and the vector $\bm{\gamma}$. The two
matrices are derived by referring to Eqs. from (55) to (58). In particular,
from Eqs. (57) and (58) we get

$\displaystyle\bm{C}=\begin{pmatrix}\bm{I}_{d}&0\\\
0&\bm{C}_{1}\end{pmatrix},$ (71)

and

$\displaystyle\bm{CD}=\begin{pmatrix}0&0\\\
0&\bm{C}_{1}\bm{D}_{1}\end{pmatrix}.$ (72)

In the above equations, $\bm{C_{1}}$ and $\bm{C_{1}D_{1}}$ are $d\times d$
matrices with elements

$(\bm{C_{1}})_{ij}=\delta_{ij}\left[\left(1-\delta_{Dj}\right)+\delta_{Dj}c^{2}\right],$
(73)

and

$(\bm{C_{1}D_{1}})_{ij}=\delta_{ij}\delta_{Dj}e^{i\varphi}sc,$ (74)

respectively, where $D$ is the index of the input port into which
$|\xi\rangle$ is injected, and $s\equiv\sinh{r}$, $c\equiv\cosh{r}$ (these
shortcuts will be repeatedly used below). On account of Eqs. (55), (56) and
(69) and recalling the fundamental relation
$\bm{A}=\bm{\mathcal{A}}^{\dagger}$ discussed above, we have
($\bm{U}\equiv\bm{A_{1}}$ in the following)

$\displaystyle\bm{E}=$
$\displaystyle\bm{A}^{\dagger}\bm{C}\bm{A}=\frac{1}{2}\begin{pmatrix}\bm{I}_{d}+\bm{E_{1}}&i(\bm{I}_{d}-\bm{E_{1}})\\\
-i(\bm{I}_{d}-\bm{E_{1}})&\bm{I}_{d}+\bm{E_{1}}\end{pmatrix},$

with $\bm{E_{1}}=\bm{A_{1}}^{\dagger}\bm{C_{1}}\bm{A_{1}}$, and

$\displaystyle\bm{EN}=$
$\displaystyle\bm{A}^{\dagger}\bm{CD}\bm{A}^{*}=\frac{1}{2}\begin{pmatrix}-\bm{E_{1}N_{1}}&-i\bm{E_{1}N_{1}}\\\
-i\bm{E_{1}N_{1}}&\bm{E_{1}N_{1}}\end{pmatrix},$

with $\bm{E_{1}N_{1}}=\bm{A_{1}}^{\dagger}\bm{C_{1}D_{1}}\bm{A_{1}}^{*}$. It
is not difficult to see that, when none of the input modes is in a state
$|\beta,\xi\rangle$ with both $\beta\neq 0$ and $\xi\neq 0$ – that is, each
mode is either in a coherent or in a squeezed vacuum state – the expression of
$\bm{\gamma}$ can be simplified as $\bm{\gamma}=\bm{b}^{*}$. Indeed, if that
is the case, then $\bm{E}^{*}\bm{b}^{*}=\bm{b}^{*}$ and
$(\bm{EN})^{*}\bm{b}=0$. In such case one also has
$\bm{b}=\bm{A}^{\dagger}\bm{\beta}$ and $\bm{\beta}=(\bm{\beta_{0}}\hskip
2.84544pt0)^{T}$, with

$\displaystyle\bm{\beta_{0}}=\begin{pmatrix}\alpha_{1}\\\ \vdots\\\
\alpha_{d}\end{pmatrix}=\begin{pmatrix}|\alpha_{1}|e^{i\phi_{1}}\\\ \vdots\\\
|\alpha_{d}|e^{i\phi_{d}}\end{pmatrix}.$

Thus, one finds

$\displaystyle\bm{\gamma}=\frac{1}{\sqrt{2}}\begin{pmatrix}\bm{\beta_{0}}^{*}\\\
i\bm{\beta_{0}}^{*}\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}\bm{\gamma_{0}}\\\
i\bm{\gamma_{0}}\end{pmatrix},$

where the symbol $\bm{\gamma_{0}}\equiv\bm{\beta_{0}}^{*}$ was introduced. We
are interested in the following combinations of $\bm{\gamma}$:

$\bm{\gamma}\bm{\gamma}^{T}=\frac{1}{2}\begin{pmatrix}\bm{\gamma_{0}}\bm{\gamma_{0}}^{T}&i\bm{\gamma_{0}}\bm{\gamma_{0}}^{T}\\\
i\bm{\gamma_{0}}\bm{\gamma_{0}}^{T}&-\bm{\gamma_{0}}\bm{\gamma_{0}}^{T}\end{pmatrix},$

and

$\bm{\gamma}\bm{\gamma}^{\dagger}=\frac{1}{2}\begin{pmatrix}\bm{\gamma_{0}}\bm{\gamma_{0}}^{\dagger}&-i\bm{\gamma_{0}}\bm{\gamma_{0}}^{\dagger}\\\
i\bm{\gamma_{0}}\bm{\gamma_{0}}^{\dagger}&\bm{\gamma_{0}}\bm{\gamma_{0}}^{\dagger}\end{pmatrix},$

where

$\displaystyle\bm{\gamma_{0}}\bm{\gamma_{0}}^{T}=$
$\displaystyle\begin{pmatrix}|\alpha_{1}|^{2}e^{-2i\phi_{1}}&|\alpha_{1}||\alpha_{2}|e^{-i(\phi_{1}+\phi_{2})}&\ldots&|\alpha_{1}||\alpha_{d}|e^{-i(\phi_{1}+\phi_{d})}\\\
|\alpha_{2}||\alpha_{1}|e^{-i(\phi_{2}+\phi_{1})}&|\alpha_{2}|^{2}e^{-2i\phi_{2}}&\ldots&|\alpha_{2}||\alpha_{d}|e^{-i(\phi_{2}+\phi_{d})}\\\
\vdots&\vdots&&\vdots\\\
|\alpha_{d}||\alpha_{1}|e^{-i(\phi_{d}+\phi_{1})}&|\alpha_{d}||\alpha_{2}|e^{-i(\phi_{d}+\phi_{2})}&\ldots&|\alpha_{d}|^{2}e^{-2i\phi_{d}}\end{pmatrix}\equiv\bm{\phi_{+}},$

and

$\displaystyle\bm{\gamma_{0}}\bm{\gamma_{0}}^{\dagger}=$
$\displaystyle\begin{pmatrix}|\alpha_{1}|^{2}&|\alpha_{1}||\alpha_{2}|e^{-i(\phi_{1}-\phi_{2})}&\ldots&|\alpha_{1}||\alpha_{d}|e^{-i(\phi_{1}-\phi_{d})}\\\
|\alpha_{2}||\alpha_{1}|e^{-i(\phi_{2}-\phi_{1})}&|\alpha_{2}|^{2}&\ldots&|\alpha_{2}||\alpha_{d}|e^{-i(\phi_{2}-\phi_{d})}\\\
\vdots&\vdots&&\vdots\\\
|\alpha_{d}||\alpha_{1}|e^{-i(\phi_{d}-\phi_{1})}&|\alpha_{d}||\alpha_{2}|e^{-i(\phi_{d}-\phi_{2})}&\ldots&|\alpha_{d}|^{2}\end{pmatrix}\equiv\bm{\phi_{-}}.$

According to Eq. (68), which reads

$\displaystyle\bm{h}=$
$\displaystyle\bm{EN}\circ(\bm{EN})^{*}-\bm{EN}\circ\bm{\gamma}\bm{\gamma}^{T}-(\bm{EN})^{*}\circ\left(\bm{\gamma}\bm{\gamma}^{T}\right)^{*}+\bm{E}\circ\bm{E}^{*}+\bm{E}\circ\bm{\gamma}\bm{\gamma}^{\dagger}+\bm{E}^{*}\circ\left(\bm{\gamma}\bm{\gamma}^{\dagger}\right)^{*}-\left(\bm{E}+\bm{\gamma}\bm{\gamma}^{\dagger}\right)\circ\bm{I},$

$\bm{h}$ is obtained by referring to the following matrices:

$\displaystyle\bm{EN}\circ(\bm{EN})^{*}=\frac{1}{4}\begin{pmatrix}\bm{E_{1}N_{1}}\circ(\bm{E_{1}N_{1}})^{*}&\bm{E_{1}N_{1}}\circ(\bm{E_{1}N_{1}})^{*}\\\
\bm{E_{1}N_{1}}\circ(\bm{E_{1}N_{1}})^{*}&\bm{E_{1}N_{1}}\circ(\bm{E_{1}N_{1}})^{*}\end{pmatrix},$
(75)
$\displaystyle-\bm{EN}\circ\bm{\gamma}\bm{\gamma}^{T}=\frac{1}{4}\begin{pmatrix}\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}&-\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}\\\
-\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}&\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}\end{pmatrix},$
(76)
$\displaystyle\bm{E}\circ\bm{E}^{*}=\frac{1}{4}\begin{pmatrix}\bm{I}_{d}+\left(\bm{E_{1}}+\bm{E_{1}}^{*}\right)\circ\bm{I}_{d}+\bm{E_{1}}\circ\bm{E_{1}}^{*}&\bm{I}_{d}-\left(\bm{E_{1}}+\bm{E_{1}}^{*}\right)\circ\bm{I}_{d}+\bm{E_{1}}\circ\bm{E_{1}}^{*}\\\
\bm{I}_{d}-\left(\bm{E_{1}}+\bm{E_{1}}^{*}\right)\circ\bm{I}_{d}+\bm{E_{1}}\circ\bm{E_{1}}^{*}&\bm{I}_{d}+\left(\bm{E_{1}}+\bm{E_{1}}^{*}\right)\circ\bm{I}_{d}+\bm{E_{1}}\circ\bm{E_{1}}^{*}\end{pmatrix},$
(77)
$\displaystyle\bm{E}\circ\bm{\gamma}\bm{\gamma}^{\dagger}=\frac{1}{4}\begin{pmatrix}\bm{I}_{d}\circ\bm{\phi_{-}}+\bm{E_{1}}\circ\bm{\phi_{-}}&\bm{I}_{d}\circ\bm{\phi_{-}}-\bm{E_{1}}\circ\bm{\phi_{-}}\\\
\bm{I}_{d}\circ\bm{\phi_{-}}-\bm{E_{1}}\circ\bm{\phi_{-}}&\bm{I}_{d}\circ\bm{\phi_{-}}+\bm{E_{1}}\circ\bm{\phi_{-}}\end{pmatrix},$
(78)
$\displaystyle-\left(\bm{E}+\bm{\gamma}\bm{\gamma}^{\dagger}\right)\circ\bm{I}_{d}=\frac{1}{2}\begin{pmatrix}\bm{I}_{d}+\bm{E_{1}}\circ\bm{I}_{d}+\bm{I}_{d}\circ\bm{\phi_{-}}&0\\\
0&\bm{I}_{d}+\bm{E_{1}}\circ\bm{I}_{d}+\bm{I}_{d}\circ\bm{\phi_{-}}\end{pmatrix}.$
(79)

The above set of equations shows that $\bm{h}$ is a sum of $2d\times 2d$
matrices of the form $\begin{pmatrix}\bm{X}^{(k)}&\bm{Y}^{(k)}\\\
\bm{Y}^{(k)}&\bm{X}^{(k)}\end{pmatrix}$, $\bm{X}^{(k)}$, $\bm{Y}^{(k)}$ being
$d\times d$ matrices. Consequently, computing the two sums
$\bm{X}=\sum_{k}\bm{X}^{(k)}$ and $\bm{Y}=\sum_{k}\bm{Y}^{(k)}$ is enough to
determine $\bm{h}$ completely. We find

$\displaystyle\bm{X}=\frac{1}{4}\bm{E_{1}N_{1}}\circ(\bm{E_{1}N_{1}})^{*}+\frac{1}{4}\bm{E_{1}}\circ\bm{E_{1}}^{*}+\frac{1}{4}\left[\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}+(\bm{E_{1}N_{1}})^{*}\circ\bm{\phi_{+}}^{*}\right]+\frac{1}{4}(\bm{E_{1}}\circ\bm{\phi_{-}}+\bm{E_{1}}^{*}\circ\bm{\phi_{-}}^{*})+\frac{1}{4}(\bm{E_{1}}+\bm{E_{1}}^{*})\circ\bm{I}_{d}-\frac{1}{2}\bm{E_{1}}\circ\bm{I}_{d}-\frac{1}{4}\bm{I}_{d}$

and

$\displaystyle\bm{Y}=\frac{1}{4}\bm{E_{1}N_{1}}\circ(\bm{E_{1}N_{1}})^{*}+\frac{1}{4}\bm{E_{1}}\circ\bm{E_{1}}^{*}-\frac{1}{4}\left[\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}+(\bm{E_{1}N_{1}})^{*}\circ\bm{\phi_{+}}^{*}\right]-\frac{1}{4}(\bm{E_{1}}\circ\bm{\phi_{-}}+\bm{E_{1}}^{*}\circ\bm{\phi_{-}}^{*})-\frac{1}{4}(\bm{E_{1}}+\bm{E_{1}}^{*})\circ\bm{I}_{d}+\frac{1}{2}\bm{I}_{d}\circ\bm{\phi_{-}}+\frac{1}{4}\bm{I}_{d}.$

Finally, from Eq. (70), and using $\bm{h}_{i,j}=\bm{h}_{i+d,j+d}=\bm{X}_{ij}$
and $\bm{h}_{i,j+d}=\bm{h}_{i+d,j}=\bm{Y}_{ij}$, we find

$\displaystyle\bm{\mathcal{F}_{\rm Q}}=$ $\displaystyle
2\bm{X}-2\bm{Y}=\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}^{*}+(\bm{E_{1}N_{1}})^{*}\circ\bm{\phi_{+}}^{*}+\bm{E_{1}}\circ\bm{\phi_{-}}+\bm{E_{1}}^{*}\circ\bm{\phi_{-}}^{*}-\bm{\phi_{-}}\circ\bm{I}_{d}+\bm{E_{1}}^{*}\circ\bm{I}_{d}-\bm{I}_{d}.$
(80)

Furthermore, we use

$\displaystyle(\bm{E_{1}N_{1}})_{ij}$ $\displaystyle=$
$\displaystyle\left(\bm{A_{1}}^{\dagger}\bm{C_{1}D_{1}}\bm{A_{1}}^{*}\right)_{ij}=\sum_{kl}\left(\bm{U}^{\dagger}\right)_{ik}\left(e^{i\varphi}sc\right)\delta_{kl}\delta_{Dl}\left(\bm{U}^{*}\right)_{lj}=e^{i\varphi}sc\hskip
2.84544pt\left(\bm{U}^{\dagger}\right)_{iD}\left(\bm{U}^{*}\right)_{Dj}$
$\displaystyle=$ $\displaystyle e^{i\varphi}sc\hskip
2.84544pt\left(\bm{U}^{*}\right)_{Di}\left(\bm{U}^{*}\right)_{Dj}$

and

$\displaystyle(\bm{E_{1}})_{ij}$ $\displaystyle=$
$\displaystyle\left(\bm{A_{1}}^{\dagger}\bm{C_{1}}\bm{A_{1}}\right)_{ij}=\sum_{kl}\left(\bm{U}^{\dagger}\right)_{ik}\delta_{kl}\left[\left(1-\delta_{Dl}\right)+\delta_{Dl}c^{2}\right]\left(\bm{U}\right)_{lj}=\sum_{kl}\left(\bm{U}^{\dagger}\right)_{ik}\delta_{kl}\left(\bm{U}\right)_{lj}+\sum_{kl}\left(\bm{U}^{\dagger}\right)_{ik}(c^{2}-1)\delta_{kl}\delta_{Dl}\left(\bm{U}\right)_{lj}$
$\displaystyle=$
$\displaystyle\delta_{ij}+s^{2}\left(\bm{U}^{*}\right)_{Di}\left(\bm{U}\right)_{Dj},$

derived on the basis of Eq. (74) and (73), respectively. Also, since
$(\bm{\phi_{+}})_{ij}=|\alpha_{i}||\alpha_{j}|e^{-i\phi_{i}}e^{-i\phi_{j}}$,
and
$(\bm{\phi_{-}})_{ij}=|\alpha_{i}||\alpha_{j}|e^{-i\phi_{i}}e^{i\phi_{j}}$,
one has

$(\bm{E_{1}N_{1}}\circ\bm{\phi_{+}})_{ij}=e^{i\varphi}sc\hskip
2.84544pt\left(\bm{U}^{*}\right)_{Di}\left(\bm{U}^{*}\right)_{Dj}|\alpha_{i}||\alpha_{j}|e^{-i\phi_{i}}e^{-i\phi_{j}},$
(81)

and

$\displaystyle(\bm{E_{1}}\circ\bm{\phi_{-}})_{ij}=$
$\displaystyle\left[\delta_{ij}+s^{2}\left(\bm{U}^{*}\right)_{Di}\left(\bm{U}\right)_{Dj}\right]|\alpha_{i}||\alpha_{j}|e^{-i\phi_{i}}e^{i\phi_{j}}=|\alpha_{i}|^{2}\delta_{ij}+s^{2}\left(\bm{U}^{*}\right)_{Di}\left(\bm{U}\right)_{Dj}|\alpha_{i}||\alpha_{j}|e^{-i\phi_{i}}e^{i\phi_{j}}.$
(82)

For the sake of a lighter notation, we set $(\bm{U})_{Di}\equiv u_{i}$; the
same notation was used in the main text. Plugging Eqs. (81) and (82) into Eq.
(80) and by some algebraic manipulation, we find

$\displaystyle(\bm{\mathcal{F}_{\rm Q}})_{ij}=$
$\displaystyle|\alpha_{i}||\alpha_{j}|\left[\left(e^{-i\varphi}u_{i}e^{i\phi_{i}}u_{j}e^{i\phi_{j}}+e^{i\varphi}u^{*}_{i}e^{-i\phi_{i}}u^{*}_{j}e^{-i\phi_{j}}\right)sc+2\delta_{ij}+\left(u^{*}_{i}e^{-i\phi_{i}}u_{j}e^{i\phi_{j}}+u_{i}e^{i\phi_{i}}u^{*}_{j}e^{-i\phi_{j}}\right)s^{2}-\delta_{ij}\right]+|u_{i}|^{2}\delta_{ij}s^{2}$
$\displaystyle=$
$\displaystyle|\alpha_{i}||\alpha_{j}|\left[\Re\left(e^{i\chi_{i}}u_{i}\right)\Re\left(e^{i\chi_{j}}u_{j}\right)\left(e^{2r}-1\right)+\Im\left(e^{i\chi_{i}}u_{i}\right)\Im\left(e^{i\chi_{j}}u_{j}\right)\left(e^{-2r}-1\right)\right]+\left(|\alpha_{i}|^{2}+|u_{i}|^{2}\bar{n}_{s}\right)\delta_{ij},$
(83)

where $\chi_{i}\equiv\phi_{i}-\varphi/2$ and
$s^{2}\equiv\sinh^{2}r=\bar{n}_{s}$. By imposing the condition
$\Im(e^{i\chi_{j}}u_{j})=0$, for $j=1,...,d$, one finally recovers Eq. (6) of
the main text.

####  VII.1.3 Inverse moment matrix, Eq. (5)

The so called moment matrix corresponds to the covariance matrix of a
particular set of estimators of the unknown parameters $\theta_{i}$,
$i=1,\dots,d$. It is defined as follows GessnerNATCOMM2020 :

$\bm{\mathcal{M}}=\bm{G}^{T}\bm{\Gamma}^{-1}\bm{G},$ (84)

where

$\bm{G}_{ij}=\frac{\partial\langle\hat{X}_{i}\rangle_{\rho(\bm{\theta})}}{\partial\theta_{j}}\hskip
14.22636pt(i=1,\dots,K;\hskip 2.84544ptj=1,\dots,d),$ (85)

and

$\bm{\Gamma}_{ij}=\langle\hat{X}_{i}\hat{X}_{j}\rangle_{\rho(\bm{\theta})}-\langle\hat{X}_{i}\rangle_{\rho(\bm{\theta})}\langle\hat{X}_{j}\rangle_{\rho(\bm{\theta})}\hskip
14.22636pt(i,j=1,\dots,K),$ (86)

$\rho(\bm{\theta})$ being the output state of the whole sensor network. The
$\hat{X}_{i}$, $i=1,\dots,K$, are $K$ hermitian operators which correspond to
measurements performed on the output state $\rho(\bm{\theta})$. We choose
$\hat{X}_{i}\equiv(\hat{J}_{z})_{i}=\left(\hat{a}_{i}^{\dagger}\hat{a}_{i}-\hat{b}_{i}^{\dagger}\hat{b}_{i}\right)/2$,
with $i=1,\dots,d$. In such case, the two matrices $\bm{G}$ and $\bm{\Gamma}$
can be evaluated through the Q function-based technique already illustrated in
the previous section. This time, we need to take into account the complete
transformation performed by the network on the input state, including the
phase–encoding stage. Such transformation is described by the matrix

$\displaystyle\bm{A}^{\dagger}=$
$\displaystyle\begin{pmatrix}\bm{\tilde{C}}&\bm{\tilde{S}}\\\
-\bm{\tilde{S}}&\bm{\tilde{C}}\end{pmatrix}\begin{pmatrix}\bm{I}_{d}&0\\\
0&\bm{A_{1}}^{\dagger}\end{pmatrix}=\begin{pmatrix}\bm{\tilde{C}}&\bm{\tilde{S}}\bm{A_{1}}^{\dagger}\\\
-\bm{\tilde{S}}&\bm{\tilde{C}}\bm{A_{1}}^{\dagger}\end{pmatrix},$ (87)

where $\bm{\tilde{C}}_{ij}=\cos{\left(\theta_{i}/2\right)}\delta_{ij}$ and
$\bm{\tilde{S}}_{ij}=\sin{\left(\theta_{i}/2\right)}\delta_{ij}$. Let
$|\Psi(\bm{\theta})\rangle$ be the output state of our MZI sensor network. On
the basis of Eq. (86) and having set $\hat{X}_{i}\equiv(\hat{J}_{z})_{i}$, we
have

$\bm{\Gamma}_{ij}=\langle\Psi(\bm{\theta})|(\hat{J}_{z})_{i}(\hat{J}_{z})_{j}|\Psi(\bm{\theta})\rangle-\langle\Psi(\bm{\theta})|(\hat{J}_{z})_{i}|\Psi(\bm{\theta}\rangle)\langle\Psi(\bm{\theta})|(\hat{J}_{z})_{j}|\Psi(\bm{\theta})\rangle\hskip
14.22636pt(i,j=1,\dots,d),$ (88)

which can be rewritten as

$\bm{\Gamma}_{ij}=\frac{1}{4}\left(\bm{h}_{i,j}+\bm{h}_{i+d,j+d}-\bm{h}_{i,j+d}-\bm{h}_{i+d,j}\right)\hskip
14.22636pt(i,j=1,\dots,d),$ (89)

where
$\bm{h}_{ij}=\langle\Psi(\bm{\theta})|\hat{n}_{i}\hat{n}_{j}|\Psi(\bm{\theta})\rangle-\langle\Psi(\bm{\theta})|\hat{n}_{i}|\Psi(\bm{\theta})\rangle\langle\Psi(\bm{\theta})|\hat{n}_{j}|\Psi(\bm{\theta})\rangle$.
Notice that $\bm{h}_{ij}$ can still be computed by means of Eq. (68), provided
that the expression of matrix $\bm{A}^{\dagger}$ is that given in Eq. (87).
Matrices $\bm{C}$ and $\bm{CD}$, on the other hand, being only related to the
form of the input state, have the same expressions as in Eqs. from (71) to
(74). So we get

$\displaystyle\bm{E}=\bm{A}^{\dagger}\bm{C}\bm{A}=\begin{pmatrix}\bm{\tilde{C}}^{2}+\bm{\tilde{S}}\bm{E_{1}}\bm{\tilde{S}}&-\left(\bm{\tilde{S}}\bm{\tilde{C}}-\bm{\tilde{S}}\bm{E_{1}}\bm{\tilde{C}}\right)\\\
-\left(\bm{\tilde{S}}\bm{\tilde{C}}-\bm{\tilde{C}}\bm{E_{1}}\bm{\tilde{S}}\right)&\bm{\tilde{S}}^{2}+\bm{\tilde{C}}\bm{E_{1}}\bm{\tilde{C}}\end{pmatrix}$
(90)

and

$\displaystyle\bm{EN}=$
$\displaystyle\bm{A}^{\dagger}\bm{CD}\bm{A}^{*}=\begin{pmatrix}\bm{\tilde{S}}\bm{E_{1}N_{1}}\bm{\tilde{S}}&\bm{\tilde{S}}\bm{E_{1}N_{1}}\bm{\tilde{C}}\\\
\bm{\tilde{C}}\bm{E_{1}N_{1}}\bm{\tilde{S}}&\bm{\tilde{C}}\bm{E_{1}N_{1}}\bm{\tilde{C}}\end{pmatrix}.$

We also have

$\displaystyle\bm{\gamma}=$
$\displaystyle\bm{b}^{*}=\begin{pmatrix}\bm{\tilde{C}}\bm{\beta_{0}}^{*}\\\
-\bm{\tilde{S}}\bm{\beta_{0}}^{*}\end{pmatrix}=\begin{pmatrix}\bm{\tilde{C}}\bm{\gamma_{0}}\\\
-\bm{\tilde{S}}\bm{\gamma_{0}}\end{pmatrix},$

so that

$\displaystyle\bm{\gamma}\bm{\gamma}^{\dagger}=\begin{pmatrix}\bm{\tilde{C}}\bm{\phi_{-}}\bm{\tilde{C}}&-\bm{\tilde{C}}\bm{\phi_{-}}\bm{\tilde{S}}\\\
-\bm{\tilde{S}}\bm{\phi_{-}}\bm{\tilde{C}}&\bm{\tilde{S}}\bm{\phi_{-}}\bm{\tilde{S}}\end{pmatrix},$
(91)

and

$\displaystyle\bm{\gamma}\bm{\gamma}^{T}=\begin{pmatrix}\bm{\tilde{C}}\bm{\phi_{+}}\bm{\tilde{C}}&-\bm{\tilde{C}}\bm{\phi_{+}}\bm{\tilde{S}}\\\
-\bm{\tilde{S}}\bm{\phi_{+}}\bm{\tilde{C}}&\bm{\tilde{S}}\bm{\phi_{+}}\bm{\tilde{S}}\end{pmatrix}.$

At this point, we set $\theta_{i}=\pi/2$, a choice which is expected to lead
to minimum estimation uncertainty [maximum slope condition, see Eq. (93)
below]. We thus obtain the same matrices as in Eqs. (75)-(79), while, for what
concerns Eq. (76), we find

$\displaystyle-\bm{EN}\circ\bm{\gamma}\bm{\gamma}^{T}=\frac{1}{4}\begin{pmatrix}-\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}&\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}\\\
\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}&-\bm{E_{1}N_{1}}\circ\bm{\phi_{+}}\end{pmatrix},$

which has opposite sign with respect to the corresponding matrix used in the
QFIM case. As a consequence, we eventually obtain an expression for
$\bm{\Gamma}$ which is almost identical to the QFIM, the most notable
difference being the swap of the real and imaginary part of the term
$e^{i\chi_{i}}u_{i}$ (compare with Eq. (83)):

$\displaystyle\bm{\Gamma}_{ij}=\frac{1}{4}\left\\{|\alpha_{i}||\alpha_{j}|\left[\Im\left(e^{i\chi_{i}}u_{i}\right)\Im\left(e^{i\chi_{j}}u_{j}\right)\left(e^{2r}-1\right)+\Re\left(e^{i\chi_{i}}u_{i}\right)\Re\left(e^{i\chi_{j}}u_{j}\right)\left(e^{-2r}-1\right)\right]+\left(|\alpha_{i}|^{2}+|u_{i}|^{2}\bar{n}_{s}\right)\delta_{ij}\right\\}.$
(92)

Consider now Eq. (85), rewritten for the specific choice
$\hat{X}_{i}\equiv(\hat{J}_{z})_{i}$:

$\displaystyle\bm{G}_{ij}=$
$\displaystyle\frac{\partial\langle\Psi(\bm{\theta})|(\hat{J}_{z})_{i}|\Psi(\bm{\theta})\rangle}{\partial\theta_{j}}=\frac{1}{2}\frac{\partial\langle\Psi(\bm{\theta})|(\hat{n}_{i}-\hat{n}_{i+d})|\Psi(\bm{\theta})\rangle}{\partial\theta_{i}}\delta_{ij},\hskip
14.22636pt(i=1,\dots,d).$ (93)

A formula for the evaluation of
$\langle\Psi(\bm{\theta})|\hat{n}_{i}|\Psi(\bm{\theta})\rangle$ was derived in
Eq. (65):

$\displaystyle\langle\Psi(\bm{\theta})|\hat{n}_{i}|\Psi(\bm{\theta})\rangle=-1+\bm{E}_{ii}+\bm{\gamma}_{i}\bm{\gamma}_{i}^{*}\hskip
14.22636pt(i=1,\dots,2d).$

The same information is conveniently condensed in the (diagonal) matrix

$\displaystyle-\bm{I}_{d}+\left(\bm{E}+\bm{\gamma}\bm{\gamma}^{\dagger}\right)\circ\bm{I}_{d},$

whose diagonal elements are the
$\langle\Psi(\bm{\theta})|\hat{n}_{i}|\Psi(\bm{\theta})\rangle$,
$i=1,\dots,2d$. Making use of Eqs. (90) and (91), one finds

$\displaystyle-\bm{I}_{d}+\left(\bm{E}+\bm{\gamma}\bm{\gamma}^{\dagger}\right)\circ\bm{I}_{d}=$
$\displaystyle\begin{pmatrix}-\bm{I}_{d}+\bm{\tilde{C}}^{2}+\left(\bm{\tilde{S}}\bm{E_{1}}\bm{\tilde{S}}\right)\circ\bm{I}_{d}+\bm{\tilde{C}}\bm{\phi_{-}}\bm{\tilde{C}}&0\\\
0&-\bm{I}_{d}+\bm{\tilde{S}}^{2}+\left(\bm{\tilde{C}}\bm{E_{1}}\bm{\tilde{C}}\right)\circ\bm{I}_{d}+\bm{\tilde{S}}\bm{\phi_{-}}\bm{\tilde{S}}\end{pmatrix}$
$\displaystyle\begin{pmatrix}-\bm{I}_{d}+\bm{\tilde{C}}^{2}+\bm{\tilde{S}}^{2}\left(\bm{E_{1}}\circ\bm{I}_{d}\right)+\bm{\tilde{C}}\bm{\phi_{-}}\bm{\tilde{C}}&0\\\
0&-\bm{I}_{d}+\bm{\tilde{S}}^{2}+\bm{\tilde{C}}^{2}\left(\bm{E_{1}}\circ\bm{I}_{d}\right)+\bm{\tilde{S}}\bm{\phi_{-}}\bm{\tilde{S}}\end{pmatrix},$

where the general result
$(\bm{Z}\bm{Y}\bm{Z})\circ\bm{I}_{d}=\bm{Z}^{2}(\bm{Y}\circ\bm{I}_{d})$ has
been used, $\bm{Z}$ being a diagonal matrix. Note that, for $i=1,\dots,d$, the
expectation value
$\langle\Psi(\bm{\theta})|(\hat{n}_{i}-\hat{n}_{i+d})|\Psi(\bm{\theta})\rangle$
corresponds to the difference between the two diagonal blocks of this matrix,
so that, after simple algebraic manipulation, we get

$\displaystyle\bm{G}_{ij}=$
$\displaystyle\frac{1}{2}\frac{\partial\langle\Psi(\bm{\theta})|(\hat{n}_{i}-\hat{n}_{i+d})|\Psi(\bm{\theta})\rangle}{\partial\theta_{i}}\delta_{ij}=\frac{1}{2}\sin\theta_{i}\left[\bm{E_{1}}\circ\bm{I}_{d}-\left(1+|\alpha_{i}|^{2}\right)\right]_{ii}\delta_{ij}.$

We set again $\theta_{i}=\pi/2$, ending up with:

$\displaystyle\bm{G}_{ij}=\frac{1}{2}\left(|u_{i}|^{2}s^{2}-|\alpha_{i}|^{2}\right)\delta_{ij}.$
(94)

At this point, if we aimed at obtaining the explicit expression of the moment
matrix, we would still have to compute the inverse of matrix $\bm{\Gamma}$. On
the other hand, we only need $\bm{G}^{-1}$ to derive
$\bm{M}^{-1}=(\bm{G}^{-1})^{T}\bm{\Gamma}\bm{G}^{-1}$, the advantage being
that $\bm{G}$ is diagonal and thus easily inverted. From Eqs. (92) and (94),
one readily gets

$\displaystyle\left(\bm{M}^{-1}\right)_{ij}=$
$\displaystyle|\alpha_{i}||\alpha_{j}|\frac{\Im{(e^{i\chi_{i}}u_{i})}}{|\alpha_{i}|^{2}-|u_{i}|^{2}\bar{n}_{s}}\frac{\Im{(e^{i\chi_{j}}u_{j})}}{|\alpha|_{j}^{2}-|u_{j}|^{2}\bar{n}_{s}}(e^{2r}-1)+|\alpha_{i}||\alpha_{j}|\frac{\Re{(e^{i\chi_{i}}u_{i})}}{|\alpha_{i}|^{2}-|u_{i}|^{2}\bar{n}_{s}}\frac{\Re{(e^{i\chi_{j}}u_{j})}}{|\alpha|_{j}^{2}-|u_{j}|^{2}\bar{n}_{s}}(e^{-2r}-1)+\frac{|\alpha_{i}|^{2}+|u_{i}|^{2}\bar{n}_{s}}{(|\alpha_{i}|^{2}-|u_{i}|^{2}\bar{n}_{s})^{2}}\delta_{ij},$
(95)

where $\chi_{i}\equiv\phi_{i}-\varphi/2$ and
$s^{2}\equiv\sinh^{2}r=\bar{n}_{s}$. Under the condition
$\Im(e^{i\chi_{j}}u_{j})=0$, for $j=1,...,d$, one finally recovers Eq. (5).

###  VII.2 Demonstration of Eqs. (11) and (12)

First, let us recall the general expression of
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$, obtained using Eq. (5):

$\displaystyle\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}=\bm{v}^{T}\bm{\mathcal{M}}(\tilde{\bm{u}},|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},\bar{n}_{s})^{-1}\bm{v}$
$\displaystyle\quad=(e^{-2r}-1)\left(\sum_{j=1}^{d}\frac{|\alpha_{j}|\tilde{u}_{j}v_{j}}{|\alpha_{j}|^{2}-\tilde{u}_{j}^{2}\bar{n}_{s}}\right)^{2}+\sum_{j=1}^{d}\frac{|\alpha_{j}|^{2}+\tilde{u}_{j}^{2}\bar{n}_{s}}{\left(|\alpha_{j}|^{2}-\tilde{u}_{j}^{2}\bar{n}_{s}\right)^{2}}v_{j}^{2}.$

Also, as a general result, we notice that the minimization problem
$\min_{\bm{U}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$ is characterized
by a symmetry of a simple type: the value of the minimum is invariant under
the transformation $v_{i}\to-v_{i}$ performed on an arbitrary number of
components of $\bm{v}$. This can be proved as follows. Let
$\bm{\tilde{u}}_{\bm{v}}$ be the specific vector $\bm{\tilde{u}}$ that
minimizes $\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$. Inverting the sign
of $v_{i}$ for a certain $i\in\\{1,\dots,d\\}$ [$v_{i}^{\prime}=-v_{i}$] will
leave the second term in the second line of Eq. (VII.2) unaltered, only the
first term will change. However, it is easy to realize that the only effect of
this transformation is to modify the minimum point of
$\Delta^{2}(\bm{v}^{\prime}\cdot\bm{\theta})_{\rm emom}$, which is now
$\bm{\tilde{u}}_{\bm{v^{\prime}}}$ with
$(\bm{\tilde{u}}_{\bm{v^{\prime}}})_{i}=-(\bm{\tilde{u}}_{\bm{v}})_{i}$, while
the minimum itself will remain the same as for
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm emom}$. Also, this argument can be
immediately extended to a transformation that inverts the sign of an arbitrary
number of components of $\bm{v}$. The same symmetry is easily seen to apply
also to
$\min_{r_{1}^{\prime},\dots,r_{d}^{\prime}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
smom}$, the sequential strategy case.

The lower bound (11) is obtained by noticing that

$\min_{\bm{U},|\alpha_{1}|^{2},...|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\geq\min_{\bm{U},|\alpha_{1}|^{2},...|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}_{j}\cdot\bm{\theta})_{\rm
emom}$

where we have indicated $\bm{v}_{j}$ as the vector with elements
$v_{j}=1/\sqrt{d}$ and $v_{i\neq j}=0$ and corresponds to the estimation of a
single phase shift $\bm{v}_{j}\cdot\bm{\theta}=\theta_{j}/\sqrt{d}$. This
inequality is supported by numerical evidence, see Fig. 3. Specifically, from
Eq. (VII.2), we have

$\displaystyle\Delta^{2}(\bm{v}_{j}\cdot\bm{\theta})_{\rm emom}$
$\displaystyle=$
$\displaystyle\frac{(e^{-2r}-1)|\alpha_{j}|^{2}\tilde{u}_{j}^{2}}{d(|\alpha_{j}|^{2}-\tilde{u}_{j}^{2}\bar{n}_{s})^{2}}+\frac{|\alpha_{j}|^{2}+\tilde{u}_{j}^{2}\bar{n}_{s}}{d(|\alpha_{j}|^{2}-\tilde{u}_{j}^{2}\bar{n}_{s})^{2}}$
(97)
$\displaystyle\approx\frac{(e^{-2r}-1)\tilde{u}_{j}^{2}+1}{d|\alpha_{j}|^{2}}+\frac{\tilde{u}_{j}^{2}\bar{n}_{s}}{d|\alpha_{j}|^{4}},$

where the approximate expression in the second line is obtained for
$|\alpha_{j}|^{2}\gg\bar{n}_{s}$. Rewriting the numerator of Eq. (97) as
$|\alpha_{j}|^{2}+[\bar{n}_{s}-(1-e^{-2r})|\alpha_{j}|^{2}]\tilde{u}_{j}^{2}$,
we see that the minimum in the interval $0\leqslant\tilde{u}^{2}_{j}\leqslant
1$ is clearly achieved for $\tilde{u}^{2}_{j}=1$ if the term between square
brackets is negative, namely $|\alpha_{j}|^{2}>\bar{n}_{s}/(1-e^{-2r})$, or
simply $|\alpha_{j}|^{2}>\bar{n}_{s}$ for $r\gg 1$. Moreover, the minimum of
Eq. (97) with respect to $|\alpha_{j}|^{2}$ is clearly obtained by setting
$|\alpha_{j}|^{2}$ to its maximum value, which corresponds to
$|\alpha_{j}|^{2}=d\bar{n}_{c}$. Thus, we have

$\min_{\bm{U},|\alpha_{1}|^{2},...|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}_{j}\cdot\bm{\theta})_{\rm
emom}=\frac{e^{-2r}}{d^{2}\bar{n}_{c}}+\frac{\bar{n}_{s}}{d^{3}\bar{n}_{c}^{2}}.$

We recall that the above equation has been derived under the condition
$|\alpha_{j}|^{2}=d\bar{n}_{c}\gg\bar{n}_{s}$ that implies
$\bar{n}_{T}=d\bar{n}_{c}+\bar{n}_{s}\approx d\bar{n}_{c}$ and this
$\bar{n}_{T}\gg\bar{n}_{s}$. We thus recover Eq. (11).

The upper bound (12) is obtained for the sub-optimal conditions i)
$|\alpha_{j}|^{2}=\bar{n}_{c}$ for all $j$, where
$\bar{n}_{c}=(\bar{n}_{T}-\bar{n}_{s})/d$, and ii)
$\bm{\tilde{u}}=\bm{v}\sqrt{d}$. The bound reads

$\displaystyle\min_{\bm{U},|\alpha_{1}|^{2},...|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\leq\bm{v}^{T}\bm{\mathcal{M}}(\bm{v}\sqrt{d},\bar{n}_{c},\bar{n}_{s})^{-1}\bm{v}$
$\displaystyle=d\bar{n}_{c}(e^{-2r}-1)\left(\sum_{j=1}^{d}\frac{v_{j}^{2}}{\bar{n}_{c}-dv_{j}^{2}\bar{n}_{s}}\right)^{2}+\sum_{j=1}^{d}\frac{\bar{n}_{c}+dv_{j}^{2}\bar{n}_{s}}{(\bar{n}_{c}-dv_{j}^{2}\bar{n}_{s})^{2}}v_{j}^{2}.$

Taking into account the normalization $\sum_{j=1}^{d}v_{j}^{2}=1/d$, we have
$dv_{j}^{2}\leq 1$ for all $j$. Therefore,

$\left(\sum_{j=1}^{d}\frac{v_{j}^{2}}{\bar{n}_{c}-dv_{j}^{2}\bar{n}_{s}}\right)^{2}\leq\frac{1}{d^{2}(\bar{n}_{c}-\bar{n}_{s})^{2}}$

and

$\sum_{j=1}^{d}\frac{\bar{n}_{c}+dv_{j}^{2}\bar{n}_{s}}{(\bar{n}_{c}-dv_{j}^{2}\bar{n}_{s})^{2}}v_{j}^{2}\leq\frac{d\bar{n}_{c}+\mathcal{W}\bar{n}_{s}}{d^{2}(\bar{n}_{c}-\bar{n}_{s})^{2}},$

where $\mathcal{W}=d^{3}\sum_{j=1}^{d}v_{j}^{4}$. Notice that $\mathcal{W}\geq
1$, with $\mathcal{W}=1$ for $v_{j}=\pm 1/d$ for all $j$. Combining the above
equations and using $d\bar{n}_{c}=\bar{n}_{T}-\bar{n}_{s}$, we find

$\min_{\bm{U},|\alpha_{1}|^{2},...|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}\leq\frac{(\bar{n}_{T}-\bar{n}_{s})e^{-2r}+\bar{n}_{s}\mathcal{W}}{[\bar{n}_{T}-(d+1)\bar{n}_{s}]^{2}}.$

Finally, we recover Eq. (12) under the condition
$\bar{n}_{T}\gg(d+1)\bar{n}_{s}$.

###  VII.3 Demonstration of Eqs. (16) and (17)

The lower bound to
$\min_{\bm{U},|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}$ is obtained by first using the Cauchy-Schwarz inequality

$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm eQCR}=\bm{v}^{T}\bm{\mathcal{F}_{\rm
Q}}^{-1}\bm{v}\geq\frac{|\bm{v}|^{4}}{(\bm{v}^{T}\bm{\mathcal{F}_{\rm
Q}}\bm{v})}.$ (98)

From Eq. (6) we then have

$\bm{v}^{T}\bm{\mathcal{F}_{\rm
Q}}\bm{v}=(e^{2r}-1)\bigg{(}\sum_{j=1}^{d}|\alpha_{j}|\tilde{u}_{j}v_{j}\bigg{)}^{2}+\sum_{j=1}^{d}|\alpha_{j}|^{2}v_{j}^{2}+\bar{n}_{s}\sum_{j=1}^{d}\tilde{u}_{j}^{2}v_{j}^{2}.$
(99)

We have

$\bigg{(}\sum_{j=1}^{d}|\alpha_{j}|\tilde{u}_{j}v_{j}\bigg{)}^{2}\leq\bigg{(}\sum_{j=1}^{d}v_{j}^{2}\bigg{)}\bigg{(}\sum_{j=1}^{d}|\alpha_{j}|^{2}\tilde{u}_{j}^{2}\bigg{)}\leq\frac{\sum_{j=1}^{d}|\alpha_{j}|^{2}}{d},$

where the first inequality is due to Cauchy-Schwarz and the second is a
consequence of $\tilde{u}_{j}^{2}\leq 1$ and $\sum_{j=1}^{d}v_{j}^{2}=1/d$.
Using $v_{j}^{2}\leq 1/d$, we also have

$\sum_{j=1}^{d}|\alpha_{j}|^{2}v_{j}^{2}\leq\frac{\sum_{j=1}^{d}|\alpha_{j}|^{2}}{d},$

and

$\sum_{j=1}^{d}\tilde{u}_{j}^{2}v_{j}^{2}\leq\frac{1}{d}.$

Combining Eq. (99) with the above inequalities gives

$\bm{v}^{T}\bm{\mathcal{F}_{\rm
Q}}\bm{v}\leq\frac{e^{2r}\sum_{j=1}^{d}|\alpha_{j}|^{2}+\bar{n}_{s}}{d}.$

Taking into account that
$\sum_{j=1}^{d}|\alpha_{j}|^{2}=\bar{n}_{T}-\bar{n}_{s}$ and
$|\bm{v}|^{4}=1/d^{2}$, from Eq. (98), we obtain

$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\geq\frac{1}{d[\bar{n}_{T}e^{2r}-\bar{n}_{s}(e^{2r}-1)]}.$

The lower bound is valid for all $\bm{U}$ and all
$|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}$ and thus also for the optimal
configuration,

$\min_{\bm{U},\,|\alpha_{1}|^{2},...,|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\geq\frac{1}{d[\bar{n}_{T}e^{2r}-\bar{n}_{s}(e^{2r}-1)]},$ (100)

which corresponds to the lower bound (16).

To derive Eq. (17), we first use the Sherman-Morrison formula to invert Eq.
(6):

$\displaystyle\bm{\mathcal{F}_{\rm Q}}^{-1}=$
$\displaystyle\frac{1}{|\alpha_{i}|^{2}+\tilde{u}_{i}^{2}\bar{n}_{s}}\delta_{ij}$
$\displaystyle-\frac{\left(e^{2r}-1\right)}{1+\mathcal{K}\left(e^{2r}-1\right)}\frac{|\alpha_{i}||\alpha_{j}|\tilde{u}_{i}\tilde{u}_{j}}{\left(|\alpha_{i}|^{2}+\tilde{u}_{i}^{2}\bar{n}_{s}\right)\left(|\alpha_{j}|^{2}+\tilde{u}_{j}^{2}\bar{n}_{s}\right)},$
(101)

where
$\mathcal{K}=\sum_{j=1}^{d}\frac{|\alpha_{j}|^{2}\tilde{u}_{j}^{2}}{|\alpha_{j}|^{2}+\tilde{u}_{j}^{2}\bar{n}_{s}}$.
From Eq. (VII.3), it is possible to derive the general expression for
$\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm eQCR}$:

$\displaystyle\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}=\bm{v}^{T}\bm{\mathcal{F}_{\rm
Q}}(\bm{\tilde{u}},|\alpha_{1}|^{2},...,|\alpha_{d}|^{2},\bar{n}_{s})^{-1}\bm{v}$
$\displaystyle\quad=\sum_{j=1}^{d}\frac{v_{j}^{2}}{|\alpha_{j}|^{2}+\tilde{u}_{j}^{2}\bar{n}_{s}}-\frac{\left(e^{2r}-1\right)}{1+\mathcal{K}\left(e^{2r}-1\right)}\left(\sum_{j=1}^{d}\frac{|\alpha_{j}|\tilde{u}_{j}v_{j}}{|\alpha_{j}|^{2}+\tilde{u}_{j}^{2}\bar{n}_{s}}\right)^{2}.$

Similar to the derivation of Eq. (12), the upper bound (17) is obtained by a
specific configuration of the sensor network that simplifies the above
equation: i) $|\alpha_{j}|^{2}=\bar{n}_{c}$ for all $j$, where
$\bar{n}_{c}=(\bar{n}_{T}-\bar{n}_{s})/d$, and ii)
$\bm{\tilde{u}}=\bm{v}\sqrt{d}$. In this case, we obtain

$\min_{\bm{U},\,|\alpha_{1}|^{2},...|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\leq\frac{\mathcal{K}^{\prime}}{1+d\bar{n}_{c}(e^{2r}-1)\mathcal{K}^{\prime}}.$
(103)

with
$\mathcal{K}^{\prime}=\sum_{j=1}^{d}\frac{v_{j}^{2}}{\bar{n}_{c}+dv_{j}^{2}\bar{n}_{s}}$.
To derive an upper bound that is independent from $\bm{v}$, we notice that
$\mathcal{K}^{\prime}/(1+x\mathcal{K}^{\prime})$ is a monotonic growing
function of $\mathcal{K}^{\prime}$, implying that
$f(\mathcal{K}_{1}^{\prime})<f(\mathcal{K}_{2}^{\prime})$ if
$\mathcal{K}_{1}^{\prime}<\mathcal{K}_{2}^{\prime}$. We have
$\bar{n}_{c}+dv_{j}^{2}\bar{n}_{s}\geq\bar{n}_{c}$ and thus
$\mathcal{K}^{\prime}\leq 1/(d\bar{n}_{c})$. Therefore,

$\min_{\bm{U},\,|\alpha_{1}|^{2},...|\alpha_{d}|^{2}}\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
eQCR}\leq\frac{1}{d\bar{n}_{c}e^{2r}}.$ (104)

We recover the upper bound (17) when noticing that
$d\bar{n}_{c}=\bar{n}_{T}-\bar{n}_{s}$.

###  VII.4 Demonstration of Eqs. (20) and (21)

Let us consider Eq. (VII.2) for $|\alpha_{j}|^{2}=\bar{n}_{c}$, for all $j$.
In the regime $\bar{n}_{c}\gg\bar{n}_{s}$ [or, equivalently,
$\bar{n}_{T}=d\bar{n}_{c}+\bar{n}_{s}\gg(d+1)\bar{n}_{s}$], we find

$\displaystyle\Delta^{2}(\bm{v}\cdot\bm{\theta})_{\rm
emom}=\frac{e^{-2r}-1}{\bar{n}_{c}}\left(\sum_{i=1}^{d}\tilde{u}_{i}v_{i}\right)^{2}+\frac{1}{d\bar{n}_{c}}+\frac{\bar{n}_{s}}{\bar{n}_{c}^{2}}\sum_{i=1}^{d}\tilde{u}_{i}^{2}v_{i}^{2}.$
(105)

We now show that the choice $\bm{\tilde{u}}=\sqrt{d}\bm{v}$ is optimal for the
estimation of the generalized average phase $\bm{v}_{\rm
ave}\cdot\bm{\theta}$, thus proving Eq. (20).

Let us denote Eq. (105) with $f(\bm{\tilde{u}})$, as a function of
$\bm{\tilde{u}}$. The goal here is to minimize $f(\bm{\tilde{u}})$ with
respect to $\bm{\tilde{u}}$ under the normalization condition
$g(\bm{\tilde{u}})=\sum_{i=1}^{d}\tilde{u}_{i}^{2}-1=0$. This problem can be
solved through the method of Lagrange multipliers. Let $\lambda$ denote the
Lagrange multiplier, the computation of the partial derivative with respect to
$\tilde{u}_{i}$ of the Lagrangian function
$\mathcal{L}(\bm{\tilde{u}},\lambda)=f(\bm{\tilde{u}})-\lambda
g(\bm{\tilde{u}})$ gives

$\frac{\partial}{\partial\tilde{u}_{i}}\mathcal{L}(\bm{\tilde{u}},\lambda)=\frac{e^{-2r}-1}{\bar{n}_{c}}2\left(\sum_{j=1}^{d}\tilde{u}_{j}v_{j}\right)v_{i}+2\frac{\bar{n}_{s}}{\bar{n}_{c}^{2}}\tilde{u}_{i}v_{i}^{2}-2\lambda\tilde{u}_{i}.$
(106)

The constrained minimum problem is solved by the pair
$(\bm{\tilde{u}},\lambda)$ which satisfies
$(\partial/\partial\tilde{u}_{i})\mathcal{L}(\bm{\tilde{u}},\lambda)=0$ and
$g(\bm{\tilde{u}})=0$ at the same time. Note that
$\bm{\tilde{u}}=\sqrt{d}\bm{v}$ satisfies the constraint
$g(\bm{\tilde{u}})=0$. Therefore, the choice $\bm{\tilde{u}}=\sqrt{d}\bm{v}$
is optimal if and only if a value of $\lambda$ can be found such that
$(\partial/\partial\tilde{u}_{i})\mathcal{L}\left(\bm{\tilde{u}},\lambda\right)|_{\bm{\tilde{u}}=\sqrt{d}\bm{v}}=0$.
We thus search for the unique value of $\lambda$ that simultaneously solves
all the equations in the following set:

$\frac{e^{-2r}-1}{\bar{n}_{c}}+2d\frac{\bar{n}_{s}}{\bar{n}_{c}^{2}}v_{i}^{2}-2d\lambda=0,$
(107)

with $i$ such that $v_{i}\neq 0$, since we have divided by $v_{i}\neq 0$ – if
$v_{i}=0$, the $i$-th equation is identically satisfied. If such $\lambda$
actually exists, by subtracting consecutive pairs of equations contained in
(107), we will end up with

$2d\frac{\bar{n}_{s}}{\bar{n}_{c}^{2}}\left(v_{i}^{2}-v_{i+1}^{2}\right)=0,$
(108)

where $v_{i},v_{i+1}\neq 0$. On the other hand, if the above conditions are
met, we can immediately determine the sought value of $\lambda$ by solving any
of the [linear] equations in (107). Therefore, Eq. (108) contains necessary
and sufficient conditions for $\bm{\tilde{u}}=\sqrt{d}\bm{v}$ to give the
optimal QC, which can be fulfilled in two different ways:

$\frac{{\bar{n}}_{s}}{\bar{n}_{c}^{2}}=0,\qquad{\rm or}\qquad
v_{i}^{2}-v_{i+1}^{2}=0,$

with $v_{i},v_{i+1}\neq 0$. The first alternative may be interpreted as
referring to the situation when the last term in Eq. (105) is so small with
respect to the rest that it can be neglected. If that is the case, the system
in (107) admits one solution $\lambda$ for every value of $\bm{v}$, implying
that $\bm{\tilde{u}}=\sqrt{d}\bm{v}$ must be optimal in the estimation of any
linear combination of phases. On the other hand, if such term cannot be
neglected, we need that $v_{i}^{2}-v_{i+1}^{2}=0$ for all $v_{i},v_{i+1}\neq
0$, namely, that all the non-vanishing components of $\bm{v}$ are equal in