Typicality of Heisenberg scaling precision in multi-mode quantum metrology
Giovanni Gramegna Danilo Triggiani Paolo Facchi Frank A. Narducci Vincenzo Tamma
quant-ph








giovanni.gramegna@ba.infn.it
Dipartimento di Fisica and MECENAS, Università di Bari, I-70126 Bari, Italy 
INFN, Sezione di Bari, I-70126 Bari, Italy



danilo.triggiani@port.ac.uk
School of Mathematics and Physics, University of Portsmouth, Portsmouth PO1 3QL, UK



Dipartimento di Fisica and MECENAS, Università di Bari, I-70126 Bari, Italy 
INFN, Sezione di Bari, I-70126 Bari, Italy


Department of Physics, Naval Postgraduate School, Monterey, CA, United States 



vincenzo.tamma@port.ac.uk
School of Mathematics and Physics, University of Portsmouth, Portsmouth PO1 3QL, UK
Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK





We propose a measurement setup reaching Heisenberg scaling precision for the estimation of any distributed parameter φ (not necessarily a phase) encoded into a generic M-port linear network composed only of passive elements. The scheme proposed can be easily implemented from an experimental point of view since it employs only Gaussian states and Gaussian measurements. Due to the complete generality of the estimation problem considered, it was predicted that one would need to carry out an adaptive procedure which involves both the input states employed and the measurement performed at the output; we show that this is not necessary: Heisenberg scaling precision is still achievable by only adapting a single stage. The non-adapted stage only affects the value of a pre-factor multiplying the Heisenberg scaling precision: we show that, for large values of M and a random (unbiased) choice of the non-adapted stage, this pre-factor takes a typical value which can be controlled through the encoding of the parameter φ into the linear network. 


Valid PACS appear here
Typicality of Heisenberg scaling precision in multi-mode quantum metrology
    Vincenzo Tamma
    December 30, 2023
==========================================================================







§ INTRODUCTION

The precision achievable in a measurement when all experimental noise sources are minimized is ultimately determined by the discreteness of all physical phenomena: electronic devices will suffer the discreteness of the electric charge, whereas the quantum nature of light will affect optical devices. Due to this quantum noise, 
the error in the estimation of a physical parameter φ through a measurement employing N probes (e.g. photons, electrons) is strongly limited by the so-called “shot noise” factor of 1/√(N). However, it has been proven that quantum features such as entanglement and squeezing can be exploited to go beyond the shot-noise limit  and reach a precision of order 1/N, the so-called Heisenberg limit (HL) <cit.>. 

Several quantum metrological problems have been largely studied and a few approaches have been proposed to reach a HL sensitivity. However, these protocols are usually difficult to implement experimentally due to the convoluted and challenging measurement procedures <cit.> and the fragile quantum coherence needed in the input states <cit.>. 
Gaussian states, on the other hand, provide a promising avenue for quantum optical technologies <cit.>, since they are easier to create and manipulate experimentally compared to non-Gaussian ones, such as Fock states. Moreover, they allow a complete analytical treatment from a theoretical point of view <cit.>. In particular, the squeezing of a Gaussian state, which allows for highly reduced-noise signals, appears to be a valuable tool to reach quantum super-sensitive precision <cit.>. From a metrological perspective, squeezed states are often used along with Gaussian measurements <cit.>, defined as measurement schemes producing a Gaussian probability distribution of the outcomes for any Gaussian state <cit.>. Homodyne and heterodyne detection represent paradigmatic examples of Gaussian measurements. It has been shown, both theoretically <cit.> and experimentally <cit.>, that an adaptive homodyne phase estimation performs better than heterodyne detection, and approaches closer to the intrinsic quantum uncertainty than any previous technique when no prior knowledge of the phase is given. The importance of feedback and adaptivity in quantum estimation protocols has been underlined also in subsequent works <cit.>. Adaptiveness can be avoided in an optimal protocol (or near optimal) only if some constraint in the range of variation of the parameter is given <cit.>.


Within the domain of quantum optics, photons are sent as probes through an interferometer where a parameter φ to be estimated is encoded. The information about the parameter is imprinted then in the output state of the photons, and it can be extracted by a suitable measurement. The situation which has been often considered is the case where φ is an optical phase <cit.> or a phase-like parameter <cit.>. These results clearly apply also to situations in which other quantities of interest (e.g. a distance) can be converted into an optical phase <cit.>, but they fail to cover more general situations, e.g. the unknown parameter is distributed among several components of the interferometer. Recently, some progress has been made along this direction concerning the estimation of particular functions of multiple parameters distributed in a specific manner within a particular network. <cit.>. It has been also shown in a recent work <cit.> that the presence of a single unknown parameter distributed in multiple nodes of an arbitrary network introduces non-trivial complications if no constraints are given on the range of values the parameter is allowed to assume: in fact, it appears that a simultaneous adaptive procedure both in the input probe and in the measurement is needed in order to reach the HL, making the whole scheme quite unfeasible from a practical point of view. Furthermore, the proposed scheme requires an unquantified precision and number of resources  in the adaptive procedure.


In this work we demonstrate the typicality of the Heisenberg limited sensitivity with a simple metrological technique which overcomes at the same time all these serious drawbacks. In particular, we consider a general scenario in which φ can be any parameter embedded into an arbitrary linear passive M-modes interferometer: it can be a parameter characterizing any specific component of the interferometer, or arbitrarily distributed among different components of the circuit. We will show that an experimentally feasible scheme achieving Heisenberg scaling is typically possible in such a general scenario. 
In our scheme (see <ref>), a single-mode squeezed vacuum state is sent through a linear, passive preliminary stage which scatters the input photons among all the M channels of the interferometer, in order to extract the information on φ which is distributed among all the modes. A second auxiliary stage at the output of the interferometer refocuses the photons in the only observed output port. By employing a single-mode homodyne detection, we present two broad conditions which together suffice to reach the HL: the first being the requirement that most of the injected photons are successfully refocused on the observed output mode; the second simply being a minimal-resolution requirement on the homodyne measurement. Remarkably, these conditions allow for imperfections both in the refocusing and in the measurement. Heisenberg scaling is thus achievable by choosing two additional passive and linear stages, whose roles are to conveniently scatter the input probe to all the M modes, and then to refocus the photons. Despite the fact that the choice of these unitary stages can be in general φ-dependent, we show that it is always sufficient to adapt only one of the two stages, which we will thus call optimized, while the other stage can be chosen arbitrarily and independent of the parameter. Moreover we show that the optimized stage can be prepared
with a precision which is achievable using only classical resources, or by means of a preliminary classical estimation. 
This is also consistent with the result obtained in <cit.>; namely, that a preliminary classical  estimation of φ yields enough information to correctly prepare the optimized stage and thus to achieve Heisenberg scaling in the estimation protocol. Finally we show that the non-optimized stage affects the precision simply by a constant pre-factor. Using typicality and results of measure concentration in high-dimensional vector spaces, we show that distributing the unknown parameter among an high number of modes M allows this pre-factor to typically take non-vanishing values.

The rest of the paper is organized as follows. In Section <ref> we describe the proposed optical interferometer and the relative Fisher Information. In Section <ref> we use the Fisher Information to prove that the Heisenberg scaling can be achieved under suitable physical conditions; we then show how, even in the most general case, all the adaptivity can be confined within one of the auxiliary stages. In Section <ref> we discuss the typicality of our results for interferometers with a large number of channels. Finally, in Section <ref> we draw some conclusions and discuss the outlook.





§ THE PROPOSED SETUP


Let us consider a metrological scheme where the parameter to be estimated is encoded into an M-ports passive linear network described by the unitary Û_φ acting on M bosonic modes â_j  (j=1,…,M) obeying the canonical commutation relations [â_j,â_k^†]=δ_jk, and [â_j,â_k]=[â_j^†,â_k^†]=0. For a passive linear network, the action of Û_φ on the annihilation operators is associated with an M× M unitary matrix via:

    Û^†_φâ_jÛ_φ = ∑_k=1^M (U_φ)_jkâ_k.

The unitarity of the matrix U_φ is strictly related to the conservation of the number of photons injected. By  definition, U_φ is the matrix of the single-photon transition amplitudes, i.e. (U_φ)_jk^2 is the probability that a single photon injected into the k-th input channel ends up in the j-th output channel due to the action of the network.



We now propose an estimation scheme reaching Heisenberg scaling if suitable conditions are satisfied. As shown in  <ref>, the preparation of the input probe consists in two steps: first, we inject a single-mode squeezed vacuum in the first port of .
Then, the unitary stage  is used to scatter the photons injected among all the modes.
The input state of the network  in our protocol is therefore given by |ψ_0⟩ = Ŝ_1(r)|vac⟩, 
where Ŝ_1(r)=^r/2(â_1^2-â_1^† 2) is a single-mode squeezing operator with squeezing parameter r>0, and |vac⟩=|0⟩^⊗ M is the M-mode vacuum state. The average number of photons injected in the apparatus is thus N=sinh^2 r. The state |ψ_0⟩ at the output of the network  undergoes the unitary  which refocuses all the photons into a single mode, namely the first one, where a homodyne measurement of the field quadrature x̂_θ is performed. If the refocusing procedure is not perfect there will be some photons scattered into other channels with probability 1 - P_φ, where 

    P_φ = 
    ( U_φ)_11^2

is defined by the the probability amplitude ()_11 = ()_11  for the transition from the first input to the first output port in the overall interferometer =, with  and  being the single-photon unitary matrix representatives of  and  respectively, obtained analogously to (<ref>).

The homodyne measurement is described by a Positive Operator Valued Measure (POVM) ℳ={Π̂_x}, whose elements are defined by

    Π̂_x =  ^θâ_1^†â_1|x⟩_11⟨x|^-θâ_1^†â_1.


The probability of obtaining a value x  from a measurement of the quadrature x̂_θ = ^θâ^†_1â_1x̂_1^-θâ^†_1â_1 is then given by Born's rule

    p(x|φ)=(Π̂_x Ŝ_1(r)|vac⟩⟨ ̂|S^†_1(r)^†),
 
over the output state Ŝ_1(r)|vac⟩ after the overall interferometric evolution =, which yields (see Appendix <ref>)

    p(x|φ) = 1√(2πΔ_φ)exp(-x^2/2Δ_φ),

where the variance of the Gaussian distribution

    Δ_φ = 12(1+()_11^2(cosh2r -1)+
    
    +[^-2iθ()_11^2]sinh 2r ),

encodes the parameter φ through the interferometric transition amplitude ()_11 = ()_11.

It is known from classical estimation theory that the maximum precision attainable when inferring the value of the unknown parameter φ, is given by the  so-called Cramér-Rao bound <cit.>

    δφ⩾1√(ν F(φ))=δφ_min,

where ν is the number of measurments performed, while 
F(φ) is the Fisher Information (FI):

    F(φ)=∫ p(x|φ) (∂log p(x|φ)/∂φ)^2  x.

The bound (<ref>) can be asymptotically saturated through post-processing Bayesian data analysis <cit.>.

The Fisher Information related to a Gaussian probability distribution with variance Δ_φ can be evaluated by inserting (<ref>) into (<ref>), and it reads

    F(φ)=12(∂_φΔ_φΔ_φ)^2.

 
Plugging (<ref>) into (<ref>), it is possible to explicitly evaluate the FI (see Appendix <ref>) obtaining:


    F(φ) = 2((∂_φ P_φ)f(N) -2P_φ(∂_φγ_φ) h(N)1+2P_φ f(N))^2,

where

    f(N):=N(1+cos[2(γ_φ-θ)]√(1+1/N)),


    h(N):=Nsin[2(γ_φ-θ)]√(1+1/N),

with 

    γ_φ= = ( U_φ)_11
 
being the accumulated phase through the interferometric evolution. 



§ HEISENBERG SCALING





We now demonstrate that Heisenberg scaling sensitivity can be achieved in the proposed metrological setup shown in  <ref>, if two conditions are met. The first condition is the constraint that the average number of photons scattered into channels which are not measured is a finite quantity ℓ_φ, independent of N, which translates into the condition

    P_φ= 1 - ℓ_φN +  O (1N^2),  ℓ_φ≥0

on the probability P_φ in (<ref>). Here, ℓ_φ  depends in general on the linear network  in which the parameter is embedded,  and on the auxiliary stages  and  in <ref>. In subsection <ref> we will show how, for any given arbitrary , it is possible to optimize only a single stage so that the probability distribution in (<ref>) can be expressed as (<ref>). The second condition relates the accumulated phase γ_φ = through the whole setup and the phase θ=θ_φ of the measured quadrature field x̂_θ, according to

    θ_φ = γ_φ±π2 + k_φN +  O (1N^2),  𝑘_φ≠ 0,

where k_φ can depend on φ, but is assumed to be independent of N. In practice, one can even fix  k_φ to a constant value without using additional resources. 

A heuristic explanation behind this condition can be found in  <ref>: in order to maximise the ratio in (<ref>) while keeping constant N=sinh^2 r, the choice of the quadrature  x̂_θ to be measured is a  trade-off between two opposite behaviours. One consists in minimizing the variance Δ_φ in the denominator of (<ref>), while the other consists in maximizing the sensitivity of  the variance with respect to the variations of φ, namely choosing θ such that ∂_φΔ_φ in the numerator is maximal. The former is met for θ as close as possible to γ±π/2, since x̂_γ±π/2 are the squeezed quadratures after the rotation in phase space by the phase γ_φ  accumulated through the interferometer; the latter instead requires a choice of θ-γ  far enough from the stationary points of the variance at γ±π/2, where Δ_φ and therefore the overall probability distribution p(x|φ) are insensitive to variation in the parameter φ. Noticeably, the larger is N, and thus the squeezing parameter, the closer to the squeezed direction the quadrature field should be measured, as can be seen in (<ref>)


In order to prove the claim of HL scaling, we will evaluate the asymptotics of the Fisher information (<ref>) as N→∞. Substituting the value θ=θ_φ in (<ref>) into (<ref>) and (<ref>), we get

    f(N)= -12+2k_φ^2N+18N+  O (1N^2),


    h(N) = 2k_φ(1+12N)+ O (1N^2).

Hence, substituting (<ref>) and (<ref>) in (<ref>), and neglecting higher order terms, the asymptotic behavior of the Fisher information reads

    F(φ)
    	∼ 8ϱ(k_φ,ℓ_φ)(∂_φγ_φ)^2N^2,

with

    ϱ(k,ℓ)=(8k1+16k^2+4ℓ)^2.

The quadratic scaling in the mean number photons N in (<ref>) finally proves that conditions (<ref>) and (<ref>) suffice to reach the Heisenberg scaling.

The asymptotics for the Fisher Information carries two pre-factors, ϱ(k_φ,ℓ_φ) and (∂_φγ_φ)^2. We easily notice that the pre-factor ϱ(k,ℓ) vanishes only at k=0, and attains its maximum  at ℓ=0, k=± 1/4:

    ϱ(k,ℓ)≤ϱ(± 1/4,0)=1,

so that, with this choice of the constants k and ℓ, the Fisher Information asymptotically reads

    F(φ)|_k=±1/4
    ℓ=0∼ 8(∂_φγ_φ)^2N^2

Moreover, ϱ(k,ℓ) is a decreasing function of ℓ independent of k, so that ℓ=0 is always the best case, meaning that the less photons that are scattered in different channels, the higher the sensitivity in the estimation. Instead, for a fixed arbitrary positive value of ℓ, the maximum of ϱ(k,ℓ) is reached for k=±√(4ℓ+1)/4. 



 §.§ One-sided adaptivity


Since P_φ=( U_φ)_11^2, condition (<ref>) may appear to require a simultaneous optimization of the input  and the output  in a parameter-dependent way.
This two-sided adaptation can be quite difficult to realize in practice.

However, we are going to show that in fact conditions (<ref>) and (<ref>) can always be satisfied with just a one-sided parameter-dependent adaptation, which can be performed either at the input or at the output of the network equivalently. And remarkably, this adaptation can be accomplished by performing a preliminary classical, shot-noise limited, estimation of φ. 

In particular, one can choose to adaptively optimize only  and fix  to an arbitrary parameter-independent unitary stage: in this case, one can set the parameter-dependent condition

    ()_1i=(^† U_φ^†)_1i+ O (1/√(N)).

Alternatively, it is possible to adaptively optimize only  with the condition

    ()_i1 = (U_φ^†^†)_i1+ O (1/√(N)),

while 
 can be  arbitrarily chosen. 

Remarkably both equations (<ref>) and (<ref>) imply that an error of the order of O (1/√(N)) is allowed to prepare the optimized stage to reach condition (<ref>). To show that both equations (<ref>) and (<ref>) satisfy condition (<ref>), we notice that P_φ can be expressed as a transition probability

    P_φ=( U_φ)_11^2=⟨v_out|v_in⟩^2,

between the two normalized vectors |v_in⟩ = U_φ|e_1⟩ and |v_out⟩ = ^†|e_1⟩,
with |e_1⟩=(1,0,…,0)^T. Then, equation (<ref>) translates into

    |v_out⟩ =  |v_in⟩ + |δ v⟩ =  ^iε H|v_in⟩,   ε= O (1/√(N)),
 
with some H=H^†, by unitarity.
Therefore, we can see that

    P   =⟨v_out|v_in⟩^2=⟨v_in|^-iε Hv_in⟩^2
       = 1-i ε⟨v_in|H v_in⟩ - O(ε^2) ^2 = 1-O(ε^2)
       =1- O (1/N)


We can notice that in both equations (<ref>) and (<ref>), no assumption on the non-optimized stage is made, so that its choice is completely arbitrary. This freedom affects the precision of the estimation of φ through the N-independent pre-factor (∂_φγ_φ)^2 which appears in the Fisher Information (<ref>). At this point, one may argue that this pre-factor may be vanishing if a poor choice for the non-adapted unitary is made. Remarkably, in the next Section we will show that the pre-factor is typically non-vanishing for random choices of the non-adapted stage and suitably well-behaved given linear networks .




§ TYPICAL SENSITIVITY

In this section we will address in more detail the study of the pre-factor (∂_φγ_φ)^2 in the Fisher information (<ref>), clarifying under what circumstances it can be safely considered non-vanishing and characterizing its magnitude for random choices of the non-optimized stage. First of all, we can link (∂_φγ_φ)^2 to the derivative of the matrix element (u_φ)_11=( U_φ)_11=√(P_φ)^γ_φ:

    (∂_φ u_φ)_11^2   =(∂_φ√(P_φ)+i(∂_φγ_φ)√(P_φ))^γ_φ^2
       =(∂_φ√(P_φ))^2+(∂_φγ_φ)^2 P_φ.

If condition (<ref>) is satisfied, equation (<ref>) simplifies to

    (∂_φγ_φ)^2 = (∂_φ u_φ)_11^2+  O (1/N),

so that the two quantities are equal up to order 1/N.

Now, if the adaptation is performed in the output, i.e. we choose an arbitrary  and adapt  according to equation (<ref>), we see that

    (∂_φ u_φ)_11^2=(^† G_φ)_11^2+ O (1/N),

where the Hermitian operator

    G_φ:= U_φ^†∂_φ U_φ/∂φ

is the (φ-dependent) generator of U_φ.
If, on the other hand, condition (<ref>) is realized through an adaptation on the input while taking an arbitrary , then equation (<ref>) implies that

    (∂_φ u_φ)_11^2=( U_φ G_φ U_φ^†^†)_11^2+ O (1/N).


Using equations (<ref>)-(<ref>), we can finally rewrite the asymptotic expression of the Fisher information (<ref>) as

    F(φ)∼ϱ(k_φ,ℓ_φ) f(U,G_φ)N^2,

as N→∞, where

    f(U,G_φ)=(U^† G_φ U)_11^2,

with U= if the optimization is performed on the output, while U=U_φ^†^† if the optimization is carried out on the input.  We emphasize that the pre-factor  is completely independent of the choice of the optimized stage.

The maximization of the prefactor (<ref>) can be realized, for example, if U=V_φ is some unitary diagonalizing G_φ i.e. satisfying equation V_φ^† G_φ V_φ=D_φ with D_φ=(g_1,g_2,…,g_M) being the diagonal matrix of the eigenvalues of G_φ, ordered in such a way that g_1=G_φ is the maximum eigenvalue in absolute value <cit.>.  Actually, it is not necessary to take a diagonalizing unitary to maximize (<ref>), since only the first column of U enters in the definition of f(U,G_φ); hence, to maximize f(U,G_φ) it is sufficient to require this column to be the eigenvector of G_φ corresponding to the maximum eigenvalue G_φ. 
However, even that requirement would necessitate the complete knowledge of G_φ, which in general depends on the unknown parameter φ.
Therefore,  it is more relevant to consider arbitrary choices of the non-adapted network (the unitary U) independently of φ in order to determine the practical advantages of the obtained Heisenberg scaling precision for finite values of N and only one (classically) adapted stage. 

For this reason, we will perform now a statistical analysis on the typical values which can be assumed by the prefactor f(U,G_φ) for random choices of the unitary U.
Assuming no prior knowledge of the unitary U, we sample it from the unitary group U(M) according to the  unbiased uniform distribution probability, i.e. the unitarily invariant Haar measure 𝒫.

For a random unitary  U, sampled according to this distribution, the average value of the prefactor f(U,G_φ) can be computed using techniques from random matrix theory (see Appendix <ref>):

    𝔼[f(U,G_φ)] 
    =(G_φ^2)+(G_φ)^2/M(M+1),

where [·] denotes the expectation value over U(M)  with respect to the Haar measure. 

In the trivial case of a generator proportional to the identity, G_φ=G_φ1, which corresponds to the case of a network U_φ=^G_φ1 acting as a φ-dependent global phase shifter, we have (G_φ)^2=M^2G_φ^2 and (G_φ^2)=MG_φ^2, so that the average value of the pre-factor  equals the maximum one, f_max=G_φ^2, in accordance to the fact that in this particular case every unitary in U(M) diagonalizes G_φ. 

In general, we are interested in determining the conditions which make this average value in (<ref>) as large as possible. First of all, we can note directly from expression (<ref>) that eigenvalues of opposite signs can have a detrimental effect on this average, since they lower the value of (G_φ). In general, we can find a lower bound on the average value (<ref>) using Jensen's inequality [X^2]⩾[X]^2 to obtain

    [f(U,G_φ)]⩾[(U^† G_φ U)_11]^2=[  (G_φ)/M]^2,

where again, the average [(U^† G_φ U)_11] has been computed using standard techniques (see Appendix <ref>). 


 


Notice that the right-hand side of this inequality is nothing but the square of the average between G_φ's eigenvalues. Hence, if we have some degree of control on the eigenvalues, we can achieve a result which is a certain fraction α of the maximum value f_max if the average of G_φ's eigenvalues is at least a fraction √(α) of the maximum eigenvalue, namely:

    [  (G_φ)/M]^2⩾αG_φ^2   ⇒  [f(U,G_φ)]⩾α f_max.


However, this may be not sufficient for our purposes, since the average of a random variable alone does not determine its typical behaviour: a paradigmatic elementary example is that of a real random variable taking only the values 0 or 1 with equal probabilities, thus having an average of 1/2 even if it never takes values close to 1/2. 

We will show now that this is not the case for the pre-factor f(U,G_φ), thanks to the fact that it is a sufficiently well behaved function with respect to the random unitary U. In fact, by using  results on concentration of measure in high-dimensional probability spaces, we prove in Appendix <ref> that for an network with a large number M of ports the pre-factor f(U,G_φ) becomes typical, meaning that it becomes almost constant with respect to random choices of U∈U(M) (according to the unitarily invariant measure), hence concentrating around its average value (<ref>),  bounded below by (<ref>).

In formulas, we have that

    𝒫(f-𝔼[f]⩾ε)⩽ 2 exp(-A M/G_φ^4ε^2).

where A=(72π^3)^-1. This result tells us that for large interferometers it is extremely unlikely to obtain a pre-factor sensibly different from its average, since for large values of M the probability of f(U,G_φ) being different from its average is exponentially suppressed. 


This result can be also seen from the exact distribution of the pre-factor computed for some particular cases, which is shown in  <ref>. It can be seen from these figures that as M is increased the distribution of the pre-factor concentrates around its average: in particular, for the chosen configuration, a value of M=20 is already sufficient to get this concentration. Thus, for any well-behaved linear network  such that the expectation value in (<ref>) is far enough from zero, any random choice of the non-adapted stage in the proposed interferometric setup typically yields an Heisenberg-scaling precision for the estimation of φ if the number M of interferometric channels is large enough.





§ CONCLUSIONS

We demonstrated by using a simple metrological technique the typicality of  Heisenberg scaling precision for the estimation of a generic parameter φ encoded into an arbitrary M-mode network. Our scheme can be applied regardless of the nature of the parameter, which can even  be distributed among several components of the network. In particular, the proposed scheme  makes use of a single-mode squeezed state as a probe, scattered throughout all the modes by means of an auxiliary passive linear stage. Once the information on the parameter is gathered by the probe, this gets refocused on a single output channel by a second auxiliary stage, and then detected with homodyne measurement. The analysis of the Fisher information associated with such scheme reveals that, if a constant average number of photons (not scaling with the total number of photons injected) is scattered into channels different from the one measured, due to an imperfect refocusing procedure, the Heisenberg limit can be asymptotically reached, provided that the homodyne detection is performed with a sufficient resolution. For a distributed parameter, the refocusing is generally parameter-dependent, implying some sort of adaptive procedure in order to correctly refocus the probe. However, we have shown that all the dependence on the parameter can be entirely bounded to only one of the two auxiliary stages, while the other only affects the estimation through a multiplicative pre-factor. Moreover, we have also discussed how all the information on the parameter needed to sufficiently refocus the probe can be obtained with a classical shot-noise precision, meaning that the number of resources required to adaptively optimize the auxiliary stages is not detrimental for the Heisenberg scaling precision. 
Finally, we have shown that, for a large number of modes, Heisenberg scaling is typically obtained by an arbitrary non-adapted stage, with an overwhelming probability, i.e. an exponentially suppressed probability of failure.



§ ACKNOWLEDGEMENTS

This work was supported by the Office of
Naval Research Global (N62909-18-1-2153). PF and GG are partially supported by Istituto Nazionale di Fisica Nucleare (INFN) through the project “QUANTUM”, and by the Italian National Group of Mathematical Physics (GNFM-INdAM).

GG and DT contributed equally to the drafting of this work.






§ DERIVATION OF THE PROBABILITY DISTRIBUTION IN (<REF>)


In order to evaluate the probability distribution associated to the measure of the quadrature field x_θ, it is natural to proceed with the phase-space formalism for Gaussian quantum optics. For a deep and systematic overview in this topic, several works and reviews can be found in literature <cit.>. In this appendix, we will briefly introduce only the concepts and tools needed to obtain the expression (<ref>) for the probability distribution p(x|φ) of measuring the quadrature value x at the output of the proposed interferometer in  <ref>. We recall first the definition (<ref>)

    p(x|φ)=⟨vac|Ŝ^†_1(r )û^†_φ^θâ_1^†â_1|x⟩_11⟨ x|^-θâ_1^†â_1û_φŜ_1(r )|vac⟩,

where  describes the overall interferometric evolution of the single mode squeezed state Ŝ_1(r)|vac⟩. To evaluate (<ref>), it is useful to firstly recover its Fourier transform,

    χ(ξ|φ)   =∫x̣  p(x|φ)^- xξ
       =⟨vac|Ŝ^†_1(r )û^†_φ^θâ_1^†â_1^-x̂_1ξ^-θâ_1^†â_1û_φŜ_1(r )|vac⟩.

It is possible to write this characteristic function in a more canonical way. Indeed, first notice that we can write (the derivation is given in Appendix <ref> for completeness)

    ^θâ_1^†â_1^-x̂_1ξ^-θâ_1^†â_1=
    D̂(ξ⃗_θ),

where 

    D̂(ξ⃗_θ)=^-ξ⃗_θ·z⃗
 
is the displacement operator,  with

    ξ⃗_θ=
    [ ξcosθ; ξsinθ ],  ẑ⃗=
    [ x̂_1; p̂_1 ].

Then, using equation (<ref>) we can write the characteristic function (<ref>) as

    χ(ξ|φ)=[D̂(ξ⃗_θ)û_φŜ_1(r )|vac⟩⟨ ̂|S^†_1(r )û^†_φ],

Due to the Gaussian nature of the squeezed vacuum state and the linearity of the interferometric setup, the characteristic function (<ref>) is a Gaussian bivariate function centred in zero, of the form <cit.>

    χ(ξ|φ)=^-1/2ξ⃗^T_θσ_φξ⃗^T_θ,

where σ_φ is the 2× 2 covariance matrix of the whole interferometer output state Ŝ_1(r)|vac⟩, reduced to the first mode. 

In order to evaluate this matrix, we firstly recover the covariance matrix Γ_0 of the input state Ŝ_1(r )|vac⟩, which reads

    Γ_0 =12[  ^2ℛ    0;    0 ^-2ℛ ],

where ℛ is the M× M diagonal matrix with a single non-zero entry ℛ_11≡ r. After the action of the interferometer, the covariance matrix transforms into

    Γ_φ = R_φΓ_0 R_φ^T,

where R_φ is the orthogonal and symplectic matrix associated with the interferometer unitary matrix 

    R_φ = W^†[   u_φ     0;     0 u^*_φ ]
    W,

and

    W = 1√(2)[   1  i1;   1 -i1 ],

where 1 is the M× M identity matrix. R_φ can be easily evaluated to be

    R_φ =
    [  u_φ -u_φ;  u_φ  u_φ ],

so that Γ_φ in (<ref>) reads

    Γ_φ =
    [    Δ X^2_φ     Δ XP_φ; (Δ XP_φ)^T    Δ P^2_φ ].

where we have defined the M× M matrices

    Δ X^2_φ   ≡12[^2ℛ^† - ^-2ℛ^†]
       =12[[cosh(2ℛ)^†]+[sinh(2ℛ)^T]],
    Δ P^2_φ   ≡12[ -^2ℛ^† + ^-2ℛ^†]
       = 12[[cosh(2ℛ)^†]-[sinh(2ℛ)^T]],
    Δ XP_φ   ≡12[ -^2ℛ^† - ^-2ℛ^†]
       =12[-[cosh(2ℛ)^†]+[sinh(2ℛ)^T]].

In the second lines of each of the previous expression, we have exploited the fact that ℛ is real. We are interested to evaluate σ_φ, the covariance matrix reduced to the first mode, which we can now readily write

    σ_φ =
    [ (Δ X^2_φ)_11  (Δ XP_φ)_11;  (Δ XP_φ)_11 (Δ P^2_φ)_11 ],

and insert in (<ref>). Our final step is to invert the Fourier transform to finally get the expression of the probability distribution p(x|φ) given in (<ref>). In order to do that, we introduce the 2×2 orthogonal matrix

    O_θ =
    [  cosθ -sinθ;  sinθ  cosθ ],

such that ξ⃗_θ = O_θξ⃗_0, with ξ⃗_0 = (ξ,0)^T. Then, the characteristic function (<ref>) can be written in a more convenient way, namely

    χ(ξ|φ) = ^-1/2ξ⃗^T_0 O_θ^Tσ_φ O_θξ⃗^T_0=^-1/2(O_θ^Tσ_φ O_θ)_11ξ^2.

Exploiting, by the definition of ℛ, the identities

    (cosh(2ℛ)^†)_11   = cosh(2r )()_11^2+ ∑_i=2^M ()_1i^2 
    
    = cosh(2r )   ()_11^2+(1-()_11^2),
    (sinh(2ℛ)^T)_11   =sinh(2r )()_11^2,

and some elementary trigonometry, the term (O_θ^Tσ_φ O_θ)_11 can be further manipulated to match the expression of Δ_φ given in (<ref>). In fact

    (   O_θ^Tσ_φ O_θ)_11 = ∑_i,j=1,2(O_θ)_i1(O_θ)_j1(σ_φ)_ij
       =cos^2θ(Δ X^2_φ)_11 + sin^2θ(Δ P^2_φ)_11+
         +2cosθsinθ(Δ XP_φ)_11
       =12cos^2θ(cosh(2r )()_11^2+sinh(2r )(()_11^2))+
         +12sin^2θ(cosh(2r )()_11^2-sinh(2r )(()_11^2))+
         +12(cos^2θ+sin^2θ)(1-()_11^2)+
         +cosθsinθsinh(2r )[()_11^2]
       =12(1+()_11^2(cosh(2r )-1)+
         +sinh(2r )(cos(2θ)(()_11^2)+sin(2θ)[()_11^2]))
       =12(1+()_11^2(cosh(2r )-1)+
         +[^-2iθ()_11^2]sinh(2r ))≡Δ_φ.

After applying the inverse Fourier transformation on the Gaussian characteristic function (<ref>), the probability distribution reads

    p(x|φ) = 12π∫dξ χ(ξ|φ)^ xξ =1√(2πΔ_φ)exp(-x^2/2Δ_φ)

as displayed in (<ref>).

 


 §.§ Derivation of (<ref>) 
 

Exploiting the unitarity of ^θâ_1^†â_1^, we can write

    ^θâ_1^†â_1^-ξx̂_1^-θâ_1^†â_1   =exp (-ξ^θâ_1^†â_1x̂_1^-θâ_1^†â_1).

By using the definition x̂_1 = (â_1 + â_1^†)/√(2), the first-mode quadrature operator along θ reads

    ^θâ_1^†â_1x̂_1^-θâ_1^†â_1 = 
    1/√(2) (â_1(θ) + â_1(θ)^†),

where â_1(θ) = ^θâ_1^†â_1â_1^-θâ_1^†â_1 is the first-mode annihilation operator at time θ in the Heisenberg picture. The Heisenberg equation is obtained by taking the derivative of â_1(θ) with respect to θ, and reads

    ̣̂a_1(θ) /dθ =  [â_1(θ)^†â_1(θ), â_1(θ)] = - â_1(θ),

since [â_1(θ),â_1(θ)^†] = [â_1,â_1^†] =1 and [â_1(θ),â_1(θ)] = [â_1,â_1] =0.
Therefore,

    â_1(θ) = ^-θâ_1,     â_1(θ)^† = ^-θâ_1^†,

the second equality being obtained by taking the adjoint of the first. By plugging (<ref>) into (<ref>) one gets

    ξ^θâ_1^†â_1x̂_1^-θâ_1^†â_1 =  ξcosθ x̂_1 + ξ sinθp̂_1 = ξ⃗_θ·z⃗,

where the vectors ξ⃗_⃗θ⃗ and z⃗ are given in (<ref>).
Inserting this expression into (<ref>), we finally obtain (<ref>).


§ DERIVATION OF THE FISHER INFORMATION IN (<REF>)

In this appendix we will evaluate the FI in (<ref>) from the expression in (<ref>). Let us recall that the variance Δ_φ of the Gaussian probability density function (<ref>) reads

    Δ_φ = 12(1+()_11^2(cosh2r -1)+
    
    +[^-2iθ()_11^2]sinh 2r ).

The derivative of Δ_φ is written as a sum of two contributions

    ∂_φΔ_φ = 12(cosh 2r  -1)∂_φ()_11^2+
    
    +12sinh 2r [^-2iθ∂φ()_11^2].

The derivative in the first contribution is thus evaluated

    ∂_φ^2 
       =^* + ^*
       = 2[ ^*]

while the derivative in the second contribution reads

    [^-2iθ∂_φ^2] = 2[^-2iθ]

Then, defining γ_φ as the phase of , and recalling that sinh^2 r  = N,  (<ref>) reads

    ∂_φΔ_φ =
       = [^*((cosh 2r  -1)+^2i(γ_φ-θ)sinh 2r )]
       =2[^*(N+^2i(γ_φ-θ)√(N^2+N))],

while (<ref>) can be written as

    Δ_φ = 12(1+2^2×
    ×(N+cos(2γ_φ-2θ)√(N^2+N))).

Inserting the expressions (<ref>) and (<ref>) into (<ref>), we get the FI 

    F(φ)=
       =8([^*(N + ^2i(γ_φ - θ)√(N^2+N))]1+2(u_φ)_11^2(N+cos(2γ_φ-2θ)√(N^2+N)))^2.

Moreover we can write = ^γ_φ√(P_φ), so that 

    ^* = ∂_φ P_φ/2 +  P_φ∂_φγ_φ.

Since P_φ, γ_φ and their derivatives are real, once we define the quantities

    f(N):=N(1+cos[2(γ_φ-θ)]√(1+1/N)),


    h(N):=Nsin[2(γ_φ-θ)]√(1+1/N),

we easily obtain 

    F(φ) = 2((∂_φ P_φ)f(N) -2P_φ(∂_φγ_φ) h(N)1+2P_φ f(N))^2

as displayed in (<ref>).


§ ANALYTIC DISTRIBUTION OF THE PRE-FACTOR (<REF>) IN THE FISHER INFORMATION (<REF>) FOR GENERATORS WITH ONLY TWO DISTINCT EIGENVALUES

We will derive here the explicit form of the probability density function for the pre-factor = (U^† G_φ U)_11^2 in (<ref>) for a fixed generator G_φ as U is sampled from U(M) with the Haar measure. First of all, note that this distribution depends only on the eigenvalues of G_φ, which we denote with g=(g_1,…,g_M), dropping the φ subscript for notation simplicity. This can be seen using the spectral decomposition of G_φ=V_φ^† D_φ V_φ, where D_φ=(g),  yielding:

    U^† G_φ U=(V_φ U)^† D_φ (V_φ U) U^† D_φ U,
 
where in the last step we used the invariance property of the Haar measure and we used the notation  to say that the two random variables have the same distribution. In light of this remark, we have that

    f(U,D_φ)
       =(U^† D_φ U)_11^2
       =(∑_j=1^M u_j^2 g_j)^2,

having defined the random vector u=Ue_1 obtained by the application of the random matrix U∈U(M) to the fixed basis vector e_1=(1,…,0)^T∈^M, where ^M denotes the set of M-tuples of complex numbers. We see that f(U,D_φ) can be interpreted as a weighted average of the eigenvalues of G_φ with random weights; these weights are given by the square modulus of the components of a random vector drawn from the unit sphere in ^M with the Haar measure. The distribution of this random variable can be quite complicated for a generic choice of the G_φ's eigenvalues g=(g_1,…,g_M). We will consider here the situation in which there are at most two distinct eigenvalues g_1⩾ g_2⩾ 0, i.e:

    g=(g_1,…,g_1_k, g_2,…,g_2_M-k)

so that 

    f(U,D_φ)   =(g_1 ∑_j=1^k u_j^2 +g_2 ∑_j=k+1^M u_j^2)^2
       =[(g_1-g_2) ∑_j=1^k u_j^2 +g_2]^2

where we used the normalization constraint

    ∑_j=1^M u_j^2=1.

In order to get the distribution of (<ref>), let us first consider the random quantity τ(U) defined by the sum inside the brackets, namely:

    τ(U)=∑_j=1^k u_j^2

We start from the distribution q(t) of τ(U), defined in such a way that q(t)ṭ is the probability to have

    t ⩽τ(U)= ∑_j=1^k u_j^2 ⩽ t+ṭ

or, defining x_2j-1:= u_j and x_2j:= u_j, the probability to have:

    t ⩽∑_j=1^k (x_2j-1^2+ x_2j^2) ⩽ t+ṭ.

This probability can be interpreted as the geometrical surface of a 2k-dimensional hyperspherical cap of a (2M-1)-dimensional hypersphere sitting in ^2M. Using this interpretation, one then finds that <cit.>:

    q(t)=(M-1)!/(k-1)!(M-k-1)!t^k-1(1-t)^M-k-1χ_[0,1](t),

where 

    χ_[0,1](t)=
    	1     0⩽ t⩽ 1,
    
    	0    otherwise.

Starting from the distribution (<ref>) of τ(U), the probability density function p(x) of 

    f(U,(g))=[(g_1-g_2)τ(U)+g_2]^2,

can be found with a change of variables to be:

    p(x)=1/2√(x)Δ gq(√(x)-g_2/Δ g),

where Δ g:=g_1-g_2. We then have explicitly:

    p(x)=C(g_1-√(x))^M-k-1(√(x)-g_2)^k-1/√(x)χ_[g_2,g_1](√(x))

where C is a normalization constant given by:

    C=1/2(Δ g)^M-1(M-1)!/(k-1)!(M-k-1)!.

This distribution is valid whenever G_φ has only two distinct positive eigenvalues g_1⩾ g_2⩾ 0. Numerical results are compared with the probability density function (<ref>) in  <ref>.



§ DERIVATION OF THE AVERAGE VALUE (<REF>) OF THE PRE-FACTOR IN (<REF>) FOR RANDOM UNITARY MATRICES U

We collect here some results needed for the computation of averages over the unitary group which have been used in the main text. Denoting with 𝒫 the Haar probability measure,  the average of a function f:U(M)→ is defined as 

    𝔼[f(U)]=∫ f(U)𝒫(U)

whenever this integral is defined. In order to derive the results of this work we are interested only in the moments of the matrix elements U_ij, i.e. the averages of some powers of the matrix elements and their complex conjugates, which we give here in the following lemma.


	Denoting with U_ij the matrix elements of a Haar distributed random unitary matrix U∈U(M), the only non-vanishing moments up to the fourth order are given by:
	
    𝔼[U_ij^2]   =1/M   (1⩽ i,j ⩽ M),
    𝔼[U_ij^4]   =2/M(M+1)   (1⩽ i,j ⩽ M),
    𝔼[U_ij^2 U_kj^2]   =1/M(M+1)    (i≠ k),
    𝔼[U_ij^2 U_il^2]   =1/M(M+1)    (j≠ l),
    𝔼[U_ij^2U_kl^2]   =1/M^2-1    (i≠ k, j≠ l),
    𝔼[U_ij U_kl U_il^* U_kj^*]   =-1/M(M^2-1)    (i≠ k, j≠ l).



The results gathered in the previous lemma can be expressed with two compact formulas <cit.> as:

    𝔼[U_ijU_kl^*]=δ_ikδ_jl/M
       𝔼[U_ijU_klU_mn^* U_pq^*]=δ_imδ_jnδ_kpδ_lq+δ_ipδ_jqδ_kmδ_ln/M^2-1+
               -δ_imδ_jqδ_kpδ_ln+δ_ipδ_jnδ_kmδ_lq/M(M^2-1)

These formulas allow us to prove the following lemma, which has been used to derive the averages (<ref>) and (<ref>) of the main text.

	Given a generic M× M complex matrix A, the following results hold:
	
    𝔼[(U^† A U)_ij]=(A)/Mδ_ij

	
    𝔼[(U^† A U)_ij^2]=(A)^2+(A^2)/M(M+1)δ_ij



	Equation (<ref>) can be proved using (<ref>):
	
    𝔼[(U^† A U)_ij]   =∑_k,l𝔼[U^†_ikA_klU_lj]
       =∑_k,lA_kl𝔼[U^*_kiU_lj]
       =∑_k,lA_klδ_klδ_ij/M
       =(A)/Mδ_ij,

	while equation (<ref>) can be proved using (<ref>):
	
    𝔼[(U^† A U)_ij^2]=
           =∑_k,l,m,n𝔼[U^†_ikA_klU_lj U^†_imA_mnU_nj]
           =∑_k,l,m,nA_klA_mn𝔼[U_lj U_njU_mi^*U_ki^*]
           =[(A^2)+(A)^2](1/M^2-1-1/M(M^2-1))δ_ij
           =(A^2)+(A)^2/M(M+1)δ_ij


For A=G_φ and i=j=1, the expressions (<ref>) and (<ref>) reduces to the equalities in (<ref>) and to (<ref>), respectively.


§ DERIVATION OF THE TYPICALITY RESULT IN  (<REF>)

In this appendix we will show how to derive equation (<ref>) starting from a standard result on concentration of measure in high-dimensional spaces known as Levy's Lemma, which we report in the following theorem for the sake of completeness.


	Let f:𝕊^n-1→ be a function defined over the unit euclidean sphere
	
    𝕊^n-1={ x∈^n|∑_k=1^n x_k^2=1}

	endowed with the invariant Haar probability measure 𝒫. Denote with  L the Lipschitz constant of the function, i.e. the minimum L such that
	
    f(x)-f(y)⩽ L x-y_2,

	for all x,y∈𝕊^n-1, where x_2=√(∑_k=1^n x_k^2) is the Euclidean norm.
	Then:
	
    𝒫(f-[f]⩾ε )⩽ 2 e^-n ε^2/CL^2,

	where C is some positive constant which can be taken to be C=9π^3 <cit.>.


In order to apply Theorem <ref> to our case, we need to compute the Lipschitz constant associated with the pre-factor (<ref>). First, note that  can be interpreted as a function defined on a real unit sphere. In fact, it can be written as

    =(∑_j=1^M u_j^2 g_j )^2.

where u is a complex vector on the unit sphere, given by u=Ue with e=(1,0,…,0)^T∈^M. Since only the squared moduli u_j^2 appear in this expression, we can recast the problem in terms of a real vector x∈^2M whose components are defined by:

    x_2j-1= u_j,      x_2j= u_j,      j=1,…,M.

The normalization constraint ∑_j=1^M u_j^2=1 becomes 

    ∑_j=1^2M x_j^2=1,
 
so that x∈𝕊^2M-1, the unit sphere sitting inside ^2M. We see then that the random factor in equation (<ref>) can be envisioned as a function defined over the unit sphere 𝕊^2M-1:

    (∑_j=1^Mu_j^2g_j)^2   =(∑_j=1^M (x_2j-1^2+x_2j^2)g_j)^2
       =(x^T G x)^2=:f(x),

where we have defined the diagonal matrix G=(g) with g=(g_1,g_1,…,g_M,g_M)∈^2M.
In order to apply Theorem <ref> we need to estimate the Lipschitz constant L of the function f; to this aim, we  evaluate the gradient of f, which is given by:

    ∇ f(x)
       =4(x^TGx) Gx.

The Lipschitz constant for f can be then obtained as:

    L=max_x∈𝕊^2M-1∇ f(x)_2=4G_φ^2.

To see this, note simply that:

    ∇ f(x)_2   =√([∇ f(x)]^T [∇ f(x)])
       =4|x^T G x|√(x^T G^2 x)
       ⩽ 4 G^2
       = 4 G_φ^2

where in the inequality we used the fact that |x^TG x| ⩽G and x^T G^2 x=Gx^2⩽G^2, while in the last equality we used the fact that G=G_φ. The value ∇ f(x)_2=4G_φ^2 can be obtained with x=(1,0…,0)^T, which together with (<ref>) proves (<ref>). Applying Theorem <ref> to our case with n=2M and L=4G_φ^2 finally yields equation (<ref>).




*




