Estimation with Heisenberg-Scaling Sensitivity of a Single Parameter Distributed in an Arbitrary Linear Optical Network

Let us consider a M-channel passive linear network, whose action on the state of the probe is described by the unitary operator ${\widehat{U}}_{\phi }$, $\phi$ being a single, generally distributed, unknown parameter we are interested to estimate. Due to its passive and linear nature, this network can be represented by a unitary matrix. Let then, $U_{\phi }$ be the $M\times M$ unitary matrix representing the action of ${\widehat{U}}_{\phi }$ on the annihilation operators ${\widehat{a}}_i$, $i=1,\dots ,M$, associated with each channel of the network, satisfying the commutation relations $\left[{\widehat{a}}_i{\widehat{a}}_j\right]=0$, $\left[{\widehat{a}}_i{\widehat{a}}_j^{†}\right]={\delta }_{ij}$, where we denote with ${\delta }_{ij}$ the Kronecker delta. The matrix $U_{\phi }$ is thus defined by the transformation

We will consider a single-mode squeezed state as a probe, with $\widehat{S}\left(\mathit{r}\right)={\mathrm{e}}^{\frac{1}{2}r({\widehat{a}}_1^{†2}-{\widehat{a}}_1^2)}$ squeezing operator and $N={sinh}^{2}r$ average number of photons, all injected in a single channel of the apparatus, i.e., the first one with the choice of squeezing parameters $\mathit{r}$ in Equation (2). In other words, the state $|{\psi }_0\rangle$ presents a non-vanishing number of photons only in the first mode. As discussed in Section 1, the approach of squeezing-based estimation strategies is to infer the value of $\phi$ from the transformation of the covariance matrix $\mathsf{\Gamma }=\mathrm{diag}({\mathrm{e}}^{2r},1,\dots ,{\mathrm{e}}^{-2r},1,\dots )/2$ of the state $|\psi \rangle$ after the interferometric evolution ${\widehat{U}}_{\phi }$. To do so, we will consider the model where a single output channel, say the first, is measured through homodyne detection. We will denote with $\theta$ the phase of the local oscillator, which coincides with the phase of the measured quadrature ${\widehat{x}}_{\theta }$. We will assume that $r>0$ without loss of generality. We notice that, with this assumption, the squeezed quadrature of the state $|{\psi }_0\rangle$ is ${\widehat{p}}_1\equiv -\mathrm{i}({\widehat{a}}_1-{\widehat{a}}_1^{†})/\sqrt{2}$, with $\mathrm{Var}\left[{\widehat{p}}_1\right]={\mathrm{e}}^{-2r}/2$, while the anti-squeezed quadrature field is ${\widehat{x}}_1\equiv ({\widehat{a}}_1+{\widehat{a}}_1^{†})/\sqrt{2}$, with $\mathrm{Var}\left[{\widehat{x}}_1\right]={\mathrm{e}}^{2r}/2$. In terms of the creation and annihilation operators, the measured quadrature field can be expressed as ${\widehat{x}}_{\theta }=({\mathrm{e}}^{-\mathrm{i}\theta }{\widehat{a}}_1+{\mathrm{e}}^{\mathrm{i}\theta }{\widehat{a}}_1^{†})/\sqrt{2}$.

Since the linear network ${\widehat{U}}_{\phi }$ is arbitrary, the average number of photons that can be actually detected after the interferometric evolution ranges between 0 and N. Naturally, if this number were to be small, or far from N, we would expect a sub-optimal performance of the estimation scheme, since most of the photons would come out of the network ${\widehat{U}}_{\phi }$ from channels that are not observed, and information on $\phi$ would in this way be lost. Moreover, it may happen that the transformation that it imposes on the probe in the transition to the first output port is trivial, namely that the element ${\left(U_{\phi }\right)}_{11}$ of the transition amplitude matrix $U_{\phi }$ does not depend on $\phi$. This occurrence would preclude the probe from acquiring any observable information on the parameter. In order to prevent these conditions from happening, this model includes the presence of two auxiliary linear and passive networks acting on the probe, ${\widehat{V}}_{\mathrm{in}}$ before and ${\widehat{V}}_{\mathrm{out}}$ after the linear network ${\widehat{U}}_{\phi }$ (see Figure 2). The first auxiliary stage ${\widehat{V}}_{\mathrm{in}}$ can be understood as network scattering, which distributes the probe through multiple input channels of the parameter-dependent network ${\widehat{U}}_{\phi }$. The purpose of the second stage ${\widehat{V}}_{\mathrm{out}}$ is instead to refocus the probe, after the interaction with ${\widehat{U}}_{\phi }$, into the only observed channel. The unitary matrix describing the overall network is thus given by the matrix product of the three single matrix representations $V_{\mathrm{out}}$, $U_{\phi }$ and $V_{\mathrm{in}}$.

Since in this model all the photons are injected in the first input channel, which is also the only channel observed at the output, the only relevant transition amplitude is the element ${\left(u_{\phi }\right)}_{11}$ of the overall unitary matrix in Equation (3). We can then rewrite where $P_{\phi }={\left|{\left(V_{\mathrm{out}}{U}_{\phi }{V}_{\mathrm{in}}\right)}_{11}\right|}^2$ is the probability that a single photon injected in the first port is detected at the first output port of the overall network, and ${\gamma }_{\phi }=arg\left({\left(V_{\mathrm{out}}{U}_{\phi }{V}_{\mathrm{in}}\right)}_{11}\right)$ is the phase acquired by the probe during the interferometric evolution. The Gaussian nature of the probe and of the homodyne measurements yield a Gaussian probability density function which governs the outcomes of the homodyne detection [18,19,20,21]. The univariate Gaussian probability density function in Equation (5) is centred at zero due to the absence of a displacement in the probe, while its variance is given by (see Appendix A) where $\theta$ is the phase of the local oscillator.

Once the probability density function $p_{\phi }\left(x\right)$ is known, it is possible to evaluate the Fisher information [44,45] of the estimation scheme, which in turn fixes the ultimate precision $\delta {\phi }_{\min }$ achievable in the estimation of $\phi$ through $\nu$ iteration of the measurement, given by the Cramer-Rao bound [44,45]

For a Gaussian distribution centred on zero, the Fisher information reads (see Appendix B) where ${\partial }_{\phi }:=\mathrm{d}/\mathrm{d}\phi$. We notice from Equation (6) that all the information on the parameter $\phi$ is encoded in the variance ${\sigma }_{\phi }^2$ of the measured quadrature ${\widehat{x}}_{\theta }$ through the two quantities $P_{\phi }$ and ${\gamma }_{\phi }$. Thus, we can split ${\partial }_{\phi }{\sigma }_{\phi }^2$ into two contributions, one containing the derivative of $P_{\phi }$, the other the derivative of ${\gamma }_{\phi }$, namely where ${\partial }_P$ and ${\partial }_{\gamma }$ are derivatives with respect to $P_{\phi }$ and ${\gamma }_{\phi }$, respectively, so that

Generally, without imposing any condition on the setup, this model does not achieve the Heisenberg scaling in the precision for the estimation of $\phi$. In fact, we can explicitly rewrite the variance ${\sigma }_{\phi }^2$ in Equation (6) in terms of the average number of photons $N={sinh}^{2}r$ in the probe where $\mathcal{O}\left(1\right)$ is a term of order equal to or smaller than 1, negligible in the asymptotic regime of N large (We will say that given two functions $f\left(N\right)$ and $g\left(N\right)$, $f\left(N\right)=\mathcal{O}\left(g\right(N\left)\right)$ when ${lim}_{N\to \infty }\left|f\left(N\right)/g\left(N\right)\right|<+\infty$). We can also rewrite the derivatives ${\partial }_P{\sigma }_{\phi }^2$ and ${\partial }_{\gamma }{\sigma }_{\phi }^2$ in terms of N

Plugging the asymptotics shown in Equations (12) and (13) into the expression of the Fisher information in Equation (9), we notice that the numerator of $\mathcal{F}\left(\phi \right)$ can be of the order $N^2$ at most. Since the denominator is in general of the order $N^2$ as well, it yields an overall general scaling of the Fisher information of $\mathcal{O}\left(1\right)$—i.e., even lower than the SQL.

In order for this setup to reach the Heisenberg scaling, we thus need to impose some constraints which prevent the denominator of the Fisher information, i.e., the variance ${\sigma }_{\phi }^2$, to grow with N. We show in Appendix C that the asymptotic conditions (Given a function $f\left(N\right)$ and a finite sum $p\left(N\right)$ of powers of N, we will say that $f\left(N\right)\sim p\left(N\right)$ when they show the same asymptotic behaviour. In formulas, $f\left(N\right)\sim p\left(N\right)$ when $f\left(N\right)=p\left(N\right)+\mathcal{O}\left(N^{s-\epsilon }\right)$, $\forall \epsilon >0$, with s exponent of the smallest power of N appearing in the sum $p\left(N\right)$) need to be satisfied for large N, with $0\le \ell <N$ and $k\ne 0$ arbitrary constant independent of N. We will discuss more in detail the physical meaning of these conditions in Section 2.3, and we will see that Equation (14) is a minimum-resolution requirement on the tuning of the local oscillator, while Equation (15) is the condition on the refocusing of the probe. Intuitively, Equations (14) and (15) assure that ${\sigma }_{\phi }^2$ at the denominator of $\mathcal{F}\left(\phi \right)$ does not grow as fast as ${\partial }_{\phi }{\sigma }_{\phi }^2$ at the numerator: instead, we can see in Appendix C that, when these conditions hold, the variance ${\sigma }_{\phi }^2$ becomes of the order $1/N$, while its derivative ${\partial }_{\phi }{\sigma }_{\phi }^2$ remains constant for large N. In particular, we show in Appendix C that the Fisher information in Equation (9) asymptotically reads proving the achievement of the Heisenberg scaling, with positive factor reaching its maximum value for $k=\pm 1/4$ and $\ell =0$, namely $\varrho (1/4,0)=1$.

Compared with the conditions found in the literature for single-parameter Gaussian estimation schemes based on squeezed-vacuum probes which, adjusted to the notation employed so far, can be translated into ${\gamma }_{\phi }-\theta ={tan}^{-1}{\mathrm{e}}^{2r}\sim \pi /2+{\left(4N\right)}^{-1}$ and $P_{\phi }=1$ [26,27], we see that Equations (14) and (15) achieve two important further results

Lastly, we recall that it is always possible to asymptotically saturate the Cramér-Rao bound in Equation (8) in the limit $\nu \to +\infty$ of samples with a large number $\nu$ of observations. In particular, the maximum-likelihood estimator ${\tilde{\phi }}_{\mathrm{MLE}}$ is an asymptotically efficient and Gaussian estimator which can be obtained through the maximisation of the Likelihood function associated with the set ${\left\{x_i\right\}}_{i=1,\dots ,\nu }$ of the $\nu$ measurement outcomes of the quadrature field ${\widehat{x}}_{\theta }$ [44,45]. In Appendix D we see that the non-trivial solution which maximises the Likelihood function $\mathcal{L}(\phi ;\mathit{x})$ in Equation (18) is simply given by the estimator ${\tilde{\phi }}_{\mathrm{MLE}}$ satisfying where ${\sigma }^2\left(\phi \right)$ is the variance ${\sigma }_{\phi }^2$ in Equation (6) as a function of $\phi$, and $S{\left(\mathit{x}\right)}^2$ is the usual sample variance

Generally, Equation (19) cannot be solved analytically, so that numerical methods need to be employed to find non-trivial solutions. Nevertheless, it is possible to find some exceptions, particularly for elementary functional dependencies of $P_{\phi }$ and ${\gamma }_{\phi }$ on the unknown parameter $\phi$. For example, in the case for which $P_{\phi }\equiv P$ is independent of $\phi$, and the functional dependence of $\gamma \left(\phi \right):={\gamma }_{\phi }$ on $\phi$ of the phase acquired by the probe is invertible, the function $\sigma {\left(\phi \right)}^2$ in Equation (19) can be easily inverted as well, and the maximum-likelihood estimator reads

We can notice how, due to the presence of the cosine in ${\sigma }_{\phi }^2$ in Equation (6), some prior knowledge on the parameter $\phi$ is required in order to correctly choose the invertibility interval for $cos\left(2\gamma \right(\phi )-2\theta )$—i.e., to choose the correct value of $n\in \mathbb{N}$ and the sign of the arccos function in Equation (21). In the next section, we will see how a classical prior knowledge of the parameter $\phi$ is required to satisfy condition (15), achievable with a prior coarse estimation reaching an uncertainty of the order of $1/\sqrt{N}$. Such prior knowledge on the parameter, for a large enough N, can be employed to choose the correct invertibility interval.

We can see that both conditions in Equations (14) and (15) are $\phi$-dependent, suggesting that an adaptive procedure must take place in order to employ the estimation scheme described in Section 2.1, as it is customary for ab-initio Gaussian estimation strategies [24,27,35,36,37]. However, some considerations can be made in this regard.

Condition (14) fixes the phase of the quadrature ${\widehat{x}}_{\theta }$ which needs to be measured. The quantity ${\gamma }_{\phi }$ is in fact the phase acquired by the squeezed vacuum during the interferometric evolution from the first input port to the first output port, and for ${\gamma }_{\phi }\equiv \phi$ Equation (14) resembles the condition $\theta =\phi +{tan}^{-1}{\mathrm{e}}^{2r}$ found in the literature for single-phase estimation [24]. On the other hand, condition (14) is a looser condition to reach the Heisenberg scaling, and it puts in relation the precision with which we are able to choose the phase $\theta$ of the local oscillator—given by the resolution of the homodyne detection apparatus—with the precision achievable in the estimation of $\phi$. In particular, it is evident how the minimum resolution for the homodyne detector required to reach an uncertainty $\delta \phi$ of order $1/N$ must be, in turn, of order $1/N$. This is in agreement with the common notion in metrology for which a sensor cannot detect changes in the quantity that is being measured which are smaller than its resolution.

Interestingly, we notice from Equation (14) that the constant k cannot be equal to zero. Counterintuitively, the value $k=0$ coincides with the choice of measuring the quadrature ${\widehat{x}}_{{\gamma }_{\phi }+\pi /2}$, namely the minimum-variance quadrature after the squeezed vacuum undergoes a phase-shift of magnitude ${\gamma }_{\phi }$, i.e., after the interferometric evolution given by ${\widehat{u}}_{\phi }$ in Equation (3). This apparent incongruity can be explained by observing the expression of ${\partial }_{\phi }{\sigma }_{\phi }^2$ in Equation (10). Differently from displacement-encoding approaches—in which the information on the parameter is obtained from the transformation of the displacement of the probe, and thus minimizing the noise of the signal is always the optimal choice—here the value of $\phi$ is encoded in the variance of the quadrature itself. For us to be able to extract information on the parameter, the variance of the signal needs to be sensitive to small variations in $\phi$—i.e., the derivative ${\partial }_{\phi }{\sigma }_{\phi }^2$ must be non-vanishing. Of the two contributions of ${\partial }_{\phi }{\sigma }_{\phi }^2$ in Equation (10), the one originating from the variations in $P_{\phi }$ is identically vanishing when condition (15) is satisfied, since ${\partial }_{\phi }{P}_{\phi }\simeq 0$ for $P_{\phi }$ close to its maximum. The remaining contribution derives from the variations in the overall phase ${\gamma }_{\phi }$, and the variance of the maximally squeezed quadrature ${\widehat{x}}_{{\gamma }_{\phi }+\pi /2}$ is a stationary point with respect to variations in the phase ${\gamma }_{\phi }$ and is thus insensitive to $\phi$, namely ${\partial }_{\gamma }{\sigma }_{\phi }^2=0$ for ${\gamma }_{\phi }-\theta =\pi /2$ in Equation (11) (See Figure 3).

Condition (15) is the requirement that most of the photons injected into the network end up in the observed output channel. In fact, this condition can be rewritten in terms of the average number of photons that are not correctly refocused $N(1-P_{\phi })\sim \ell$. Thus, condition (15) tells us that the number of photons which are not observed must be a constant ℓ, not growing with N. In other words, this condition assures that most of the information on $\phi$ encoded in the probe is not lost in channels that are not observed. As a matter of fact, we can see from Equation (6) that the variance of the observed quadrature ${\widehat{x}}_{\theta }$ after the interferometric evolution is the convex combination of the variances of a squeezed vacuum and the pure vacuum, with coefficients $P_{\phi }$ and $1-P_{\phi }$, respectively. In order for this variance to be ‘squeezed’, in the sense that it is of order $1/N$, the contribution from the pure vacuum must be of order $1/N$, namely $1-P_{\phi }=\mathcal{O}(1/N)$.

This condition can also be seen as a requirement of the performance of the refocusing stage ${\widehat{V}}_{\mathrm{out}}$. In fact, in order to satisfy condition (15) for a given choice of ${\widehat{V}}_{\mathrm{in}}$, the auxiliary stage ${\widehat{V}}_{\mathrm{out}}$ must be chosen so that $|V_{\mathrm{out}}{U}_{\phi }{V}_{\mathrm{in}}{|}^2\sim 1-\ell /N$. As discussed earlier, this implies that, in general, the auxiliary stage ${\widehat{V}}_{\mathrm{out}}$ which satisfies this condition depends on the value of the parameter itself, requiring an adaptive approach to find an optimal refocusing stage to reach the Heisenberg scaling. We show now that the information on $\phi$ required to engineer an adequate refocusing stage to reach the Heisenberg scaling can be obtained through a classical estimation strategy, namely that which achieves the scaling $1/\sqrt{N}$ typical of the shot-noise limit. This result is due to the structure of $P_{\phi }={\left|{\left(V_{\mathrm{out}}{U}_{\phi }{V}_{\mathrm{in}}\right)}_{11}\right|}^2$, which is essentially a transition probability $P=|{\mathit{v}}_{\mathrm{out}}·{\mathit{v}}_{\mathrm{in}}{|}^2$ between the unitary vectors ${\mathit{v}}_{\mathrm{in}}=U_{\phi }{V}_{\mathrm{in}}{\mathit{e}}_1$ and ${\mathit{v}}_{\mathrm{out}}=V_{\mathrm{out}}^{†}{\mathit{e}}_1$, with ${\mathit{e}}_1={(1,0,\dots ,0)}^T$. Hence, a small tilt of order $O(1/\sqrt{N})$ between the unit vectors ${\mathit{v}}_{\mathrm{in}}$ and ${\mathit{v}}_{\mathrm{out}}$ yields a quadratic reduction in their transition probability. To prove this, we will call ${\phi }_{\mathrm{cl}}$ the rough guess of the value of $\phi$ that is sufficiently precise to engineer a refocusing stage ${\widehat{V}}_{\mathrm{out}}$ which satisfy condition (15), and we will show that the estimation strategy to obtain this rough estimate of $\phi$ is classical, namely that the error $\delta \phi =\phi -{\phi }_{\mathrm{cl}}$ associated with the prior rough estimation is allowed to be of order $1/\sqrt{N}$. For a given choice of ${\widehat{V}}_{\mathrm{in}}$ and ℓ, we will call ${\widehat{V}}_{\mathrm{out}}\left({\phi }_{\mathrm{cl}}\right)$ a solution of Equation (15). The single-photon transition probability $P_{\phi }$ appearing in this condition can be written as the squared complex modulus of the scalar product of two M-dimensional complex vectors $U_{\phi }{V}_{\mathrm{in}}{\mathit{e}}_1$ and $V_{\mathrm{out}}{\left({\phi }_{\mathrm{cl}}\right)}^{†}{\mathit{e}}_1$, with ${\mathit{e}}_1={(1,0,\dots ,0)}^{\mathrm{T}}$. We can then write where the transition probability $\eta (\phi ,{\phi }_{\mathrm{cl}})$ is a smooth function of $\phi$ and ${\phi }_{\mathrm{cl}}$, with a locus of points of maxima along the condition ${\phi }_{\mathrm{cl}}=\phi$, since, for a perfect knowledge of the parameter the auxiliary stage, ${\widehat{V}}_{\mathrm{out}}({\phi }_{\mathrm{cl}}=\phi )$ would satisfy $\eta (\phi ,\phi )=1$. If the prior estimation ${\phi }_{\mathrm{cl}}$ slightly deviates from the real value of the parameter $\phi ={\phi }_{\mathrm{cl}}+\delta \phi$, we can write the expansion where the derivative of $\eta (\phi ,{\phi }_{\mathrm{cl}})$ is zero along the condition $\phi ={\phi }_{\mathrm{cl}}$. We can see, comparing Equations (23) and (15), that an error $\delta \phi$ of order $1/\sqrt{N}$ suffices to correctly engineer a refocusing stage ${\widehat{V}}_{\mathrm{out}}$ that allows for the Heisenberg scaling. It is then possible to conceive two-step ab initio protocols exploiting the model presented in this section: a first, coarse, classical estimation of the parameter $\phi$ is performed and the rough estimate ${\phi }_{\mathrm{cl}}$ is obtained, with an error $\delta \phi =\phi -{\phi }_{\mathrm{cl}}=\mathcal{O}(1/\sqrt{N})$ of the same order of the shot-noise limit. Then, the classical information obtained on $\phi$ can be employed to engineer the refocusing stage ${\widehat{V}}_{\mathrm{out}}$, once ${\widehat{V}}_{\mathrm{in}}$ is fixed, so that the overall network satisfies condition (15), allowing us to reach the Heisenberg scaling through the quantum strategy described in Section 2.1.

Lastly, we notice that, in order to satisfy condition (15), it is also possible to optimize the input auxiliary stage ${\widehat{V}}_{\mathrm{in}}$ while arbitrarily fixing the refocusing stage ${\widehat{V}}_{\mathrm{out}}$. In such a case, identical considerations can be made regarding the possibility of a two-step protocol, since the optimization ${\widehat{V}}_{\mathrm{in}}$ still requires only a classical coarse estimation ${\phi }_{\mathrm{cl}}$ of the parameter. Interestingly, only one of the auxiliary stages needs to be optimized, and thus depends on ${\phi }_{\mathrm{cl}}$, whether it is ${\widehat{V}}_{\mathrm{in}}$ or ${\widehat{V}}_{\mathrm{out}}$. This leaves the choice of the second auxiliary stage completely arbitrary, notwithstanding that the pre-factor ${\left({\partial }_{\phi }{\gamma }_{\phi }\right)}^2$ appearing in the Fisher information in Equation (16) is not vanishing. Indeed, the condition ${\left({\partial }_{\phi }{\gamma }_{\phi }\right)}^2=0$ corresponds to the situation in which the optimized network ${\widehat{u}}_{\phi }={\widehat{V}}_{\mathrm{out}}\left({\phi }_{\mathrm{cl}}\right){\widehat{U}}_{\phi }{\widehat{V}}_{\mathrm{in}}$ acts trivially, namely without imprinting any information about $\phi$, on the probe. Remarkably, it has been shown that, for a random choice of the non-optimized auxiliary network, sampled uniformly within the set of all the possible linear networks, the pre-factor ${\left({\partial }_{\phi }{\gamma }_{\phi }\right)}^2$ multiplying the scaling $N^2$ in the Fisher information is typically different from zero [34]. In other words, within certain non-restrictive regularity conditions and for linear networks with a large enough number of channels, the value of the pre-factor becomes essentially unaffected by the choice of the non-optimised auxiliary network. This important feature can be exploited for experimental applications, for example, employing the arbitrary non-optimised stage to manipulate the information encoded into the probe regarding the structure of a linear network with multiple unknown parameters, ultimately allowing the choice of functions of such parameters to be estimated at the Heisenberg scaling sensitivity [46].

In this section, we will apply our model for the estimation of distributed parameters to a particular example of a 2-channel network, in which the unknown parameter $\phi$ influences the reflectivity ${\eta }_{\phi }$ of a beam-splitter and the magnitudes ${\lambda }_{\phi }$ and ${\lambda }_{\phi }^{\prime }$ of two phase-shifts (see Figure 4). We can think of the global parameter $\phi$ as an external physical property, such as the temperature or the magnitude of the electromagnetic field, affecting the components of the network ${\widehat{U}}_{\phi }$. We will suppose that the functional dependence of the phase-shifts ${\lambda }_{\phi }$, ${\lambda }_{\phi }^{\prime }$ and of the reflectivity ${\eta }_{\phi }$ on the true value of the parameter $\phi$ are known and smooth, whether given by some law of nature or opportunely engineered. With reference to Figure 4, we can write the matrices representing the action of the beam-splitter and the phase-shifts as respectively, where ${\sigma }_i$, $i=1,2,3$, is the i-th Pauli matrix and ${\mathbb{I}}_2$ is the $2\times 2$ identity matrix, so that the network ${\widehat{U}}_{\phi }$ is represented by the matrix

We easily notice that $|{\left(U_{\phi }\right)}_{11}{|}^2={cos}^2{\eta }_{\phi }$, which, in general, is different from one and thus does not satisfy the condition (15), with the exception of the two values ${\eta }_{\phi }=0,\pi$ which correspond to the absence of the mixing between the two modes. As described in the model earlier, we then add two auxiliary stages $V_{\mathrm{in}}$ and $V_{\mathrm{out}}\equiv V_{\mathrm{out}}\left({\phi }_{\mathrm{cl}}\right)$, of which only one depends on a prior coarse estimation ${\phi }_{\mathrm{cl}}$ of $\phi$ realised with a classical strategy, so that $\phi -{\phi }_{\mathrm{cl}}=\delta \phi =\mathcal{O}(1/\sqrt{N})$. In particular we choose as input stage and as output stage where ${\alpha }_{{\phi }_{\mathrm{cl}}}=({\lambda }_{{\phi }_{\mathrm{cl}}}-{\lambda }_{{\phi }_{\mathrm{cl}}}^{\prime })/2-\pi /4$ is a quantity which can be obtained through a classical estimation ${\phi }_{\mathrm{cl}}$ of $\phi$. A straightforward calculation of $P_{\phi }={\left|{\left(V_{\mathrm{out}}\left({\phi }_{\mathrm{cl}}\right)U_{\phi }{V}_{\mathrm{in}}\right)}_{11}\right|}^2$ yields where $\delta \lambda ={\lambda }_{\phi }-{\lambda }_{{\phi }_{\mathrm{cl}}}$ and $\delta {\lambda }^{\prime }={\lambda }_{\phi }^{\prime }-{\lambda }_{{\phi }_{\mathrm{cl}}}^{\prime }$ are the error in the estimates of ${\lambda }_{\phi }$ and ${\lambda }_{\phi }^{\prime }$ due to the imprecision of the classical estimation ${\phi }_{\mathrm{cl}}$. We can then easily see that $P_{\phi }$ in Equation (28) satisfies condition (15), since both the errors $\delta \lambda$ and $\delta {\lambda }^{\prime }$ are of order $1/\sqrt{N}$,—i.e., $\delta \lambda =\left({\partial }_{\phi }{\lambda }_{\phi }\right)\delta \phi =\mathcal{O}(1/\sqrt{N})$, and similarly for $\delta {\lambda }^{\prime }$—and thus we obtain from Equation (28)

In order to evaluate the Fisher information in Equation (16), we need to calculate both the phase acquired by the probe throughout the whole interferometric evolution ${\gamma }_{\phi }$, and the coefficient ℓ. The phase ${\gamma }_{\phi }$ is easily obtained as the complex phase of ${\left(u_{\phi }\right)}_{11}\equiv {\left(V_{\mathrm{out}}\left({\phi }_{\mathrm{cl}}\right)U_{\phi }{V}_{\mathrm{in}}\right)}_{11}$

Since $\delta \phi =\mathcal{O}(1/\sqrt{N})$, we call h the finite N-independent constant such that $\delta \phi \sim h/\sqrt{N}$. The transition probability $P_{\phi }$ can then be written as so that the factor $\ell =h^2{({\partial }_{\phi }{\lambda }_{\phi }-{\partial }_{\phi }{\lambda }_{\phi }^{\prime })}^2/4$ appearing in the Fisher information in Equation (16) is easily evaluated comparing Equations (15) and (31). The Fisher information can be obtained from Equation (16), with ${\partial }_{\phi }{\gamma }_{\phi }$ given by Equation (30), and $\varrho (k,\ell )$ given by Equation (17), with k given by the condition on the local oscillator phase and $\ell =h^2{({\partial }_{\phi }{\lambda }_{\phi }-{\partial }_{\phi }{\lambda }_{\phi }^{\prime })}^2/4$.

We notice from the expression of ${\alpha }_{{\phi }_{\mathrm{cl}}}=({\lambda }_{{\phi }_{\mathrm{cl}}}-{\lambda }_{{\phi }_{\mathrm{cl}}}^{\prime })/2-\pi /4$ that the unknown reflectivity of the beam splitter ${\eta }_{\phi }$ does not influence the refocusing stage ${\widehat{V}}_{\mathrm{out}}\left({\phi }_{\mathrm{cl}}\right)$, but it appears in the phase ${\gamma }_{\phi }$ acquired by the probe in Equation (30). In particular, if the two phases ${\lambda }_{\phi }$ and ${\lambda }_{\phi }^{\prime }$ are vanishing, the dependence of ${\alpha }_{{\phi }_{\mathrm{cl}}}$, and thus of the refocusing stage ${\widehat{V}}_{\mathrm{out}}\left({\phi }_{\mathrm{cl}}\right)$, on the classical estimation ${\phi }_{\mathrm{cl}}$ of the parameter disappears completely. In other words, this network for ${\lambda }_{\phi }={\lambda }_{\phi }^{\prime }=0$ transforms the reflectivity ${\eta }_{\phi }$ of a beam splitter into the magnitude of a phase shift, independently from ${\eta }_{\phi }$.

We will consider a linear and passive network ${\widehat{U}}_{\phi }$, which depends on a single and generally distributed parameter, for example affecting several components of the network, as shown in Figure 1. Once again, ${\widehat{U}}_{\phi }$ admits a unitary matrix representation $U_{\phi }$ obtained through Equation (1). Differently from Section 2 though, we consider as a probe the single-mode squeezed state with $\widehat{D}\left(\mathit{\alpha }\right)={\mathrm{e}}^{\alpha ({\widehat{a}}_1^{†}-{\widehat{a}}_1)}$, i.e., a squeezed coherent state with a mean number of photons $N=N_{\mathrm{D}}+N_{\mathrm{S}}={\alpha }^2+{sinh}^{2}r$, $\alpha ,r>0$, where we can introduce the displacement $\mathit{d}$, given by ${\alpha }_i=(d_i+\mathrm{i}{d}_{i+M})/\sqrt{2}$, $i=1,\dots ,M$. We remark that the choice $\alpha ,r>0$ is a specific ($\phi$-independent) condition which is required in displacement-encoding approaches: squeezing and displacement can, in general, have different complex phases, and the condition $\alpha ,r>0$ assures that the squeezed quadrature (${\widehat{x}}_{\pi /2}$ for $r>0$) is orthogonal—i.e., conjugated—to the displaced quadrature (${\widehat{x}}_0$ for $\alpha >0$), and hence it is the most sensible to the presence of phases. Moreover, in our model, we will consider homodyne detection in all M output channels of the linear network, so that M quadrature fields ${\widehat{x}}_{i,{\theta }_i}$ are measured, where ${\theta }_i$ is the phase of the i-th local oscillator, $i=1,\dots ,M$.

Since we are observing all the output ports of the network, but a non-vanishing number of photons is injected only in the first channel, only the first column of the unitary matrix $U_{\phi }$ is relevant in our model, consisting of the transition amplitude of single photons from the first to every channel of the network (see Equation (1)). We can thus employ the parametrisation where $P_i$ is the probability that a single photon is transmitted through the linear network from the first to the i-th channel, and ${\gamma }_i$ is the phase that it would acquire (In order to keep the notation lighter, we are dropping the subscript $\phi$ in $P_i$, ${\gamma }_i$, $\mathit{\mu }$ and $\Sigma$. Nonetheless, it rests assured that these quantities, in general, depend on the unknown parameter). Once again, due to the Gaussian nature of the model, the (joint) probability distribution associated with the outcome $\mathit{x}$ of the M homodyne detectors is also Gaussian and reads [18,19,20,21]

In Equation (34), both the covariance matrix $\Sigma$ and the mean $\mathit{\mu }$ depend on the parameter $\phi$. The elements of the covariance matrix $\Sigma$ are evaluated in Appendix A, and read where ${\delta }_{ij}$ is the Kronecker delta, ${\overline{\gamma }}_i={\gamma }_i-{\theta }_i$ is the phase acquired at the output of the i-th channel relative to the correspondent local oscillator, and $\mathit{\mu }$ is the mean vector (see Appendix A)

The determinant $det[\Sigma ]$ can also be written in compact form (see Appendix A)

For a multivariate Gaussian distribution of the form shown in Equation (34), the Fisher information can be easily evaluated from its definition in Equation (7) and employ the expression of the probability distribution $p_{\phi }\left(\mathit{x}\right)$ in Equation (34) as (see Appendix B) where $C=det[\Sigma ]{\Sigma }^{-1}$ is the cofactor matrix of $\Sigma$ and $\left[A\right]={\sum }_i{A}_{ii}$ denotes the trace of the matrix A. Compared with the Fisher information shown in Equation (9) for the model described in Section 2, the Fisher information for this setup includes an additional term ${\mathcal{F}}_{\mathrm{D}}\left(\phi \right)$ which depends on the derivative of the average $\mathit{\mu }$ with respect to the parameter $\phi$. Moreover, the contribution ${\mathcal{F}}_{\mathrm{S}}\left(\phi \right)$, representing the information on $\phi$ encoded in the covariance matrix $\Sigma$, can be split into two terms, of which the first resembles the Fisher information in Equation (9) once we perform the substitution ${\sigma }_{\phi }^2\to det[\Sigma ]$: we will show in the following that, in the asymptotic regime of large $N_{\mathrm{S}}$, this is the only contribution of ${\mathcal{F}}_{\mathrm{S}}\left(\phi \right)$ which, besides ${\mathcal{F}}_{\mathrm{D}}\left(\phi \right)$, reaches the Heisenberg scaling.

Similarly to what happens to the setup described in Section 2, the Fisher information in Equation (38) does not generally reach the Heisenberg scaling unless certain conditions are met. To show this, it is necessary to evaluate all the contributions of $\mathcal{F}\left(\phi \right)$ in terms of the number of photons $N_{\mathrm{D}}$ and $N_{\mathrm{S}}$ in the asymptotic regime. First, it is convenient to express the cofactor matrix C explicitly in terms of the squeezing parameter r. In Appendix B, we see that a closed form for C in such terms exists where and ${\Sigma }_{ij}$ are shown in Equation (35). We then notice that, in order to analyse the generic asymptotic behaviour of the Fisher information in Equation (38), it suffices to separately examine the asymptotics of the terms $\mathit{\mu }$, ${\Sigma }_{ij}$, $S_{sti}$, $det[\Sigma ]$, and of their derivatives.

From Equations (35)–(37), (39) and (40) we see that $\Sigma$, S and $det[\Sigma ]$ are all of order $N_{\mathrm{S}}$ in general, while $\mathit{\mu }$ is of order $\sqrt{N_{\mathrm{D}}}$. The same asymptotic behaviours are also kept for their respective derivatives, since both ${\partial }_{\phi }{\overline{\gamma }}_i\equiv {\partial }_{\phi }{\gamma }_i$ and ${\partial }_{\phi }{P}_i$ are independent of the probe, and thus of $N_{\mathrm{D}}$ and $N_{\mathrm{S}}$. We can thus see from the expression of $\mathcal{F}\left(\phi \right)$ in Equation (38) that its numerator grows at most with $N^2$, while the denominator—i.e., the denominator $det[\Sigma ]$—in general grows with $N_{\mathrm{S}}$. Therefore, in order for $\mathcal{F}\left(\phi \right)$ to reach the Heisenberg scaling, namely a scaling of order $N^2$, some conditions must be imposed so that the denominators in ${\mathcal{F}}_{\mathrm{D}}\left(\phi \right)$ and ${\mathcal{F}}_{\mathrm{S}}\left(\phi \right)$ do not grow for large $N_{\mathrm{S}}$. In Appendix C, it can be seen that, to achieve the Heisenberg scaling, the determinant $det[\Sigma ]$ must be of the order $N^{-1}$, similarly to what happens to ${\sigma }_{\phi }^2$ in the single-homodyne setup in Section 2.2, and the conditions for this to occur are

When these conditions hold, we can introduce the finite (possibly vanishing) quantities $k_i={lim}_{N_{\mathrm{S}}\to \infty }{N}_{\mathrm{S}}({\gamma }_i\mp \pi /2)$, so that the determinant $det[\Sigma ]$ reduces to (see Appendix C) while ${\partial }_{\phi }det[\Sigma ]$, ${\partial }_{\phi }\Sigma$, ${\partial }_{\phi }C$ and C tend to constant values, and ${\partial }_{\phi }\mathit{\mu }$ scales as $\sqrt{N_{\mathrm{D}}}$. We can easily see that this makes only the first two terms of $\mathcal{F}\left(\phi \right)$ in Equation (38) dominant for large $N_{\mathrm{D}}$ and $N_{\mathrm{S}}$. When conditions (41) are met, we can thus neglect the last term in Equation (38), and the Fisher information asymptotically reaches the Heisenberg scaling in $N_{\mathrm{D}}{N}_{\mathrm{S}}$ and $N_{\mathrm{S}}^2$ (see Appendix C), where we introduced the quantities while are positive, even functions which reach their maxima at $x=\pm 1/4$ and $x=0$, respectively, namely $\varrho (1/4)=1$ and $\zeta \left(0\right)=1$.

It is now possible to compare the Fisher information in Equation (43) achieved with the present scheme, with the Fisher information in Equation (16) obtained with the setup for distributed metrology employing a squeezed vacuum state and homodyne detection on a single channel. Since this setup employs a squeezed probe with a non-vanishing displacement, it lends itself to both displacement-based and squeezing-based estimation approaches. This is reflected by the presence of two separate contributions to the Fisher information in Equation (43): the first is given by the variations in the displacement $\mathit{\mu }$, the second by the variations in the determinant $det[\Sigma ]$, with respect to changes in the value of $\phi$. Both terms present a pre-factor ${(\partial \gamma )}_{\mathrm{avg}}^2$ shown in Equation (44) which resembles the pre-factor ${\left({\partial }_{\phi }{\gamma }_{\phi }\right)}^2$ in Equation (16) for the single-homodyne counterpart: in particular, ${(\partial \gamma )}_{\mathrm{avg}}^2$ is a weighted average of the sensitivities of the complex phases ${\gamma }_i$ to changes in the parameter $\phi$, each one weighted by the ‘fraction’ $P_i$ of the probe undergoing the phase shift ${\gamma }_i$. Indeed, we can see that ${(\partial \gamma )}_{\mathrm{avg}}^2$ reduces to the single-homodyne counterpart when the whole probe is refocused into a single channel—say the first—so that $P_1=1$ and $P_i=0$ for $i=2,\dots ,M$, i.e., ${(\partial \gamma )}_{\mathrm{avg}}^2={\left({\partial }_{\phi }{\gamma }_1\right)}^2$.

In the term of the Fisher information associated with the squeezing-encoding in Equation (43), the factor $\varrho \left(k_{\mathrm{avg}}\right)$ coincides in turn with the function $\varrho (k,\ell )$ in Equation (17) for $\ell =0$—i.e., the ideal case of perfect refocusing of the probe, and all photons being observed—and after the substitution $k\to k_{\mathrm{avg}}$, namely the weighted average of the coefficients $k_i$ with the same weights $P_i$. In other words, this term can be thought of as a generalisation of the Fisher information for a single homodyne due to the presence of multiple observed channels. Noticeably, this term can be set as equal to zero only for $N_{\mathrm{S}}=0$, in which case, the whole expression in Equation (43) vanishes, ruining the Heisenberg scaling, or for $k_{\mathrm{avg}}=0$. In fact, we can see from Equation (37) that the condition $k_{\mathrm{avg}}=0$—i.e., ${\overline{\gamma }}_i=\pm \pi /2$, namely when the quadratures with minimum variances are measured at each channel minimises $det[\Sigma ]$, which becomes a stationary point with respect to the variations in $\phi$, and ${\partial }_{\phi }det[\Sigma ]$ in Equation (43) vanishes. On the other hand, the first contribution to Equation (43) is instead a new term not present in the Fisher information for the single homodyne, deriving from the information on the parameter encoded in the displacement $\mathit{\mu }$. This term achieves the Heisenberg scaling in $N_{\mathrm{D}}{N}_{\mathrm{S}}$, in the sense that it reaches the Heisenberg scaling in $N=N_{\mathrm{D}}+N_{\mathrm{S}}$ if both $N_{\mathrm{D}}$ and $N_{\mathrm{S}}$ grow with N, i.e., if $N_{\mathrm{D}}=\beta N$ and $N_{\mathrm{S}}=(1-\beta )N$ with $0\le \beta <1$ independently of N. The function $\zeta \left(k_{\mathrm{avg}}\right)$ in Equation (45) does not have roots, hence the first contribution to $\mathcal{F}$ vanishes only for $N_{\mathrm{D}}=0$, i.e., for a squeezed vacuum as a probe.

Finally, we will now write the Likelihood function $\mathcal{L}(\phi ;\mathit{x})$ for the setup described here, and discuss the maximum-likelihood estimator ${\tilde{\phi }}_{\mathrm{MLE}}$ saturating the Cramér-Rao bound shown in Equation (8). After performing $\nu$ measurements of the field quadratures ${\widehat{x}}_{i,{\theta }_i}$, the Likelihood of the outcomes $({\mathit{x}}_1,\dots ,{\mathit{x}}_{\nu })$ is given by with $p_{\phi }\left({\mathit{x}}_j\right)$ found in Equation (34). In Appendix D, we show that the maximum-likelihood estimator, a non-trivial solution of the maximisation of the Likelihood function in Equation (46), is implicitly given by the estimator ${\tilde{\phi }}_{\mathrm{MLE}}$ which satisfies the equation where $\Sigma$ is the covariance matrix in Equation (35), and $\mathit{\mu }$ the mean vector in Equation (36). This equation cannot, in general, be solved analytically, hence numerical methods typically need to be in place to find ${\tilde{\phi }}_{\mathrm{MLE}}$. On the other hand, Equation (47) simplifies for certain particular cases. For example, when measuring all the minimum-variance quadratures so that ${\overline{\gamma }}_i=\pm \pi /2$ in Equation (41), and the probabilities $P_i$ are independent of $\phi$, we can see from Equation (35) that ${\partial }_{\phi }\Sigma =0$. In this case, the Likelihood Equation becomes where the term in the right-hand side can be seen as a linear combination of the quantities $\mathit{\mu }-\tilde{\mathit{\mu }}$, where $\tilde{\mathit{\mu }}={\sum }_{j=1}^{\nu }{\mathit{x}}_j/\nu$ are estimators of the mean $\mathit{\mu }$. We see here how the displacement of the probe allows us to perform the estimation of $\phi$ through the sample mean $\mathit{\mu }$, i.e., the average of the outcomes of the homodyne measurements. On the other hand, for a squeezed vacuum such as a probe, $\mathit{\mu }=0$, so that Equation (47) becomes where it is possible to recognise the sample covariance matrix $\tilde{\Sigma }={\sum }_{j=1}^{\nu }{\mathit{x}}_j{\mathit{x}}_j^{\mathrm{T}}/\nu$, estimator of the covariance matrix $\Sigma$.

For the model introduced in this section, some considerations regarding the conditions in Equation (41) can also be drawn, especially in light of the feature of the single-homodyne, squeezed-vacuum scheme discussed in Section 2.3.

In fact, these conditions do nothing but fix the phases ${\theta }_i$ of the quadratures ${\widehat{x}}_{i,{\theta }_i}$ that need to be measured to achieve the Heisenberg scaling. Regarding the minimum resolution of the homodyne sensors, it appears that there is no evident advantage in this setup compared to the single-homodyne scheme, since the resolution required to tune the local oscillators of each channel is still of order $1/N$. On the other hand, introducing displacement in the probe allows, as previously discussed, the value of the parameter to be inferred from the information inscribed in the average of the measured quadratures $\mathit{\mu }$. Therefore, it is possible with this scheme to exactly measure the minimum-variance quadratures ${\widehat{x}}_{i,{\gamma }_i\pm \pi /2}$, a possibility that was prevented in the previous scheme due to the requirement $k\ne 0$ in Equation (14). Although the situation for which $k_i=0$ for $i=1,\dots ,M$ sets to zero the contribution ${\mathcal{F}}_{\mathrm{S}}\left(\phi \right)$ in Equation (43), associated with the information encoded in the covariance matrix $\Sigma$, the first contribution ${\mathcal{F}}_{\mathrm{D}}\left(\phi \right)$ of the Fisher information still reaches the Heisenberg scaling.

Another interesting feature of this protocol, which differentiates it from its single-homodyne counterpart, is that it does not require any adaptation of the network: since every output channel is observed, no condition on the transition probabilities $P_i$ is required. Therefore, not only there is no requirement for a $\phi$-dependent auxiliary stage, but there is no need for an auxiliary stage to reach the Heisenberg scaling in the first place. However, the precision in the estimation is still affected by the network through the terms $k_{\mathrm{avg}}$ and ${(\partial \gamma )}_{\mathrm{avg}}$, which appear in the constant factors multiplying the scaling in the Fisher information in Equation (38). In fact, we can see from Equation (45) that these two quantities depend on the transition probabilities $P_i$ and on the derivatives of the complex phases ${\partial }_{\phi }{\gamma }_i$. In particular, $\mathcal{F}\left(\phi \right)$ can be vanishing for exceptionally poorly conceived networks, for which ${(\partial \gamma )}_{\mathrm{avg}}=0$: for example, a network for which ${\gamma }_i$ is independent of $\phi$ for all values of i such that $P_i\ne 0$, is associated with a vanishing factor ${(\partial \gamma )}_{\mathrm{avg}}$, as we can see from its definition in Equation (44). In this case, adding a $\phi$-independent auxiliary network V, either at the input or at the output of $U_{\phi }$, would modify the transition amplitude of the overall interferometric evolution, and ultimately yield a non-vanishing value of ${(\partial \gamma )}_{\mathrm{avg}}$.