Optimizing quantum-enhanced Bayesian multiparameter estimation of phase and noise in practical sensors

Achieving quantum-enhanced performances when measuring unknown quantities requires developing suitable methodologies for practical scenarios, that include noise and the availability of a limited amount of resources. Here, we report on the optimization of sub-standard quantum limit Bayesian multiparameter estimation in a scenario where a subset of the parameters describes unavoidable noise processes in an experimental photonic sensor. We explore how the optimization of the estimation changes depending on which parameters are either of interest or are treated as nuisance ones. Our results show that optimizing the multiparameter approach in noisy apparata represents a significant tool to fully exploit the potential of practical sensors operating beyond the standard quantum limit for broad resources range.


INTRODUCTION
The goal of quantum metrology is to estimate a set of physical parameters, exploiting quantum resources to achieve improved performances beyond those achievable by classical methods.The use of quantum probes discloses the capability to reach the Heisenberg limit (HL), gaining a quadratic scaling advantage over the standard quantum limit (SQL) corresponding to the use of N independent probes [1][2][3][4][5].Often in a real scenario, even if the interest relies on a single parameter, the process is unavoidably affected by the presence of unknown noises.For these reasons, it is usually more effective to treat these estimations using a multiparameter approach [6][7][8][9][10].Despite their importance, experimental demonstrations of quantum enhanced estimation in the multiparameter case are still few and limited to modest amounts of coherent quantum resources [5,[11][12][13][14][15].The following two extremal scenarios are prototypical for multiparameter metrology.In the first case, all the unknown parameters are treated on the same level, and thus one needs to optimize the overall amount of information extracted.Here, the adoption of quantum probes can provide improved performances with respect to strategies where each parameter is estimated separately [16][17][18][19].In the second extremal case, only one parameter is considered to be of interest, but the dynamics of the metrological evolution intrinsically involves other nuisance parameters, of which an approximate knowledge is although necessary to retrieve a good estimator for the desired one.For instance, we have to deal with this scenario when different sources of noise affect the evolution: phase and visibility [20][21][22][23], phase and phase diffusion [24,25], magnetic field and decoherence [26] for example.The optimal strategy in this case is very different from the optimal one in the former scenario, since now the interest is to maximize the information extracted on one parameter at the expense of all the others [27].In the general case, intermediate configurations between these two extremal scenarios can be defined, corresponding to different choices of the cost function.For example, a couple of parameters could be considered of interest while the others are treated as nuisance.For each specific scenario, different strategies may thus turn out to be optimal.In general, the importance of the different parameters can be weighted arbitrarily.
Another crucial aspect of quantum metrology in a practical scenario regards the availability of a finite amount of resources N in the estimation process.The standard approach is based on a theoretical framework dedicated to defining bounds and strategies in the asymptotic limit of large N .However, when only a finite amount of resources is available, any estimation strategy needs to be tailored to optimize the convergence for low values of N [28][29][30].A powerful tool here is represented by adaptive protocols, which enable faster convergences to the ultimate limits [31,32].These have been implemented both through online [33,34] and offline [35,36] approaches also resorting to the use of different machine learning algorithms [37][38][39][40].These techniques demonstrated two relevant characteristics, namely fast convergence to the ultimate bounds and performances independent of the particular value of the parameter of interest.Different experimental applications of adaptive techniques have been reported, first in single-parameter estimation problems [38,39,[41][42][43][44] and then in a multiparameter setting [14,45].In this scenario, the versatility of the multiparameter approach allows to choose the optimal allocation of resources, depending on which are the parameters of interest and which are the ones associated to noise processes, treated instead as nuisances.Furthermore, the application of adaptive strategies in the quantum metrology context is beneficial also when the quantum probe state is in a continuous variable (CV) Gaussian state.In particular, squeezed light has emerged as a valuable resource as demonstrated by its use in various domains, from gravitational wave detection [46] to recent research on distributed quantum sensing, where CV states enable enhanced precision and sensitiv-ity in multiparameter estimation tasks [47,48].However, also in the CV framework, real-world measurement are often dynamic and subject to various sources of noise and fluctuation.Adaptive strategies, such as adaptive Bayesian estimation, allow systems to dynamically respond to changing conditions, optimizing measurement protocols in real-time ensuring enhanced measurement performances.
In this work, we investigate a multiparameter estimation scenario, where the parameters of interest are a physical rotation angle [49,50] together with the noise values involved in the interferometric measurements.To this end, we employ orbital angular momentum (OAM) of single photons, carrying tunable OAM values up to 50, able to show N00N-like sensitivities for rotation estimations [51][52][53][54][55][56].Importantly, we extend the single parameter study [52] to a multiparameter approach within a Bayesian framework [30,[57][58][59] for all the aforementioned scenarios by employing an adaptive strategy ensuring the optimal allocation of resources [60].Such an approach allows us to extend the multiparameter estimation problems to the regime of O(30, 000) number of quantumlike resources, experimentally demonstrating sub-SQL precision in the estimation of the rotation angle for wide resources ranges even with nuisance parameters.In order to quantify the quality of the reached performances, we define non-tight Bayesian bounds on the estimations in each considered scenario.This work is applicable to all those cases where the visibilities must be estimated on the fly, for example because they fluctuate from one rotation estimation to the other.However, also if a pre-calibration of the visibilities before the rotation estimation is possible, our approach will be optimal with respect to the figures of merit and resources considered, if the resources consumed in the pre-calibration are also counted.Indeed, if a pre-calibration would be more efficient, our protocol would choose it, as a first stage of the estimation.

PRECISION BOUNDS FOR THE MULTIPARAMETER ESTIMATION
The goal of multiparameter quantum metrology is to identify regimes where the estimation precision outperforms the one achievable by classical probes.A crucial aspect is to keep such enhanced performances as the employed resources increase.It becomes key to develop a platform able to investigate such a regime with quantum scaling of the precision.Here, we study the simultaneous estimation, in the large resource regime, of a rotation angle θ ∈ [0, π) and of a collection of parameters (the fringe visibilities V s1 , • • • , V s4 defined below) that affect the efficiency of the detection process [51].More specifically, in our scheme, before each measurement step, we select a control parameter s out of a set of four possible values {s 1 = 1, s 2 = 2, s 3 = 11, s 4 = 51}, each corresponding to a device which can be switched on at will and produces the quantum-like resource.In the ideal (noiseless) scenario such choice is meant to force the interferometer to produce a single-photon, N00N-like output state analogous to those employed in [42] to achieve quantum limited precision, i.e. the vector |ψ s (θ being orthogonal circular polarization states.Unfortunately, the selection of high values of s also has the indirect effect to add noise into the model which ultimately deteriorates the associated visibility V s ∈ [0, 1] of the measurements we perform on |ψ s (θ)⟩ to recover θ.Our scheme relays on two different types of detections, the first corresponding to the projection of |ψ s (θ)⟩ on the basis {(|0⟩ ± |1⟩)/ √ 2}, while the second uses {(|0⟩ ± i |1⟩)/ √ 2} as reference basis.Our interest lies in determining how the estimation precision changes depending on the different perspective from which the multiparameter problem is addressed.To this end, we apply the Bayesian algorithm in Ref. [60].This procedure identifies the most effective adaptive strategy depending on the different roles assigned to each parameter (θ, V s1 , V s2 , V s3 , and V s4 ), whether they are of interest or are treated as nuisance parameters.In this protocol, the Bayesian posterior probability distribution for all the parameters is represented by a particle filter.For more details see Appendix B and [60].From the information accumulated in the reconstructed posterior distribution, at each step, a greedy strategy selects the experimental settings that minimize the future expected estimator variance [60].This is computed through a brute force procedure that consists in simulating the evolution of the posterior distribution for all the possible measurements settings, and picking the one with the smaller expected variance.We optimize both the value s, that determines the successive probe state and the basis of the polarization measurement we perform on the output state |ψ s (θ)⟩.The information flow in the Bayesian algorithm is represented in Fig. 1, here the measurement outcomes are used to update the posterior by the computer unit that also calculates the next optimal measurement.The parameters tuned by the optimization algorithm are the used quantum resource and the appropriate polarization measurement basis.To quantify the efficiency of the estimation, we define a suitable figure of merit and derive bounds limiting its minimal achievable value for the investigated metrological task.Focusing on a generic experimental run r N composed by a series of individual measurements, where, given i ∈ {1, • • • , 4}, the control value s i is used ν i times, we define the total number of resources N = 4 i=1 ν i s i [2].Hence, indicating with θ(r N ) , the estimated values of θ, V s1 , • • • , V s4 , we get from the Bayesian procedure, we gauge the associated experimental error via the quantity where G is a weight matrix that controls which parameters are to be treated as nuisance and which are the parameters of interest [27].We repeat the whole estimation for a collection θ 1 , θ 2 , • • • , θ J of J different values of θ uniformly distributed in [0, π), and for each θ j , we run the whole procedure M > 200 times.Indeed, in this protocol, we can not fix the total number N of resources, because the stochastic nature of the measurement outcomes propagates to the online choice of the next q-plate and then to N , which is, therefore, a stochastic variable whose exact value is only known at runtime.However, we can repeat M ≫ 1 times the simulation and collet the precision after each used photon in the set of tuple {(∆ 2 rn,G (θ j ), n)}, where n is the total resource number used up to that point to get to precision ∆ 2 rn,G (θ j ).We will then look at all the points with total resources falling in the interval [n, n + ∆n] with ∆n small, and the median error of this cluster is the error we associate with n.We choose to cluster only the point having n > 100, with the chosen interval being ∆n = 50.This finally leads us to the identification of the error figure of merit: with the median taken on the M executions.This is the median square error over the different estimates, each obtained, as usual in quantum metrology, as the mean of the posterior distribution.The choice to compute the median instead of the average is due to the fact, that the former is much less sensitive to outliers which are always present in such estimation protocols (see Appendix A).It is possible to impose a theoretical lower limit on (1) which can be used to evaluate the quality of our experiment.Specifically, we can write: where the constant C G is defined and computed in Appendix A via non-linear programming on the Cramér-Rao bound and connecting the mean square error to the median square error.This bound is legit if we assume the local asymptotic normality (LAN) of the estimators.

EXPERIMENTAL APPARATUS
In order to investigate a multiparameter metrological problem, we employ a state-of-the-art experimental apparatus, that allows to generate single-photon high-OAM value states showing quantum-like sensitivities for the estimation of a rotation angle.Importantly, we extend the result from the single parameter case [52] to a multiparameter one, exploring unprecedented regimes for multiparameter quantum metrology.The experimental setup is composed of a series of q-plates [52,61] with an increasing topological charge q arranged in a cascade configuration as reported in Fig. 1.We test the protocol through single-photon states generated via a spontaneous parametric down-conversion (SPDC) source in a Sagnac configuration.Photon pairs are emitted at 808 nm pumping a periodically poled titanyl phosphate (ppKTP) crystal with a continuous laser with a wavelength equal to 404 nm.Once generated, one photon is sent through the apparatus while the other acts like a trigger allowing to work with coincidences events when the two photons of the pairs are detected by avalanche photodiodes within a time window of few nanoseconds.Starting from single photons prepared in the horizontal polarization state, we can generate a N00N-like state in the OAM degree of freedom of the form: where |0, 1⟩ refers to the circular polarization, while |±m⟩ with m = 2q represents the OAM state and depends on the FIG. 1. Sketch of the experimental setup.Single photons are sent through the apparatus consisting of a generation stage (blue rectangle) and a measurement stage (orange rectangle).The first is composed of a polarizing beam splitter (PBS) and three q-plates with different topological charges q = 1/2; 5; 25, respectively, followed by a half-waveplate (HWP).The measurement stage is composed of the same elements of the preparation mounted in reverse order in a compact and motorized cage which can be rotated by an angle θ.After the final PBS, the photons are measured through avalanche photodiodes (APDs).The measurement HWP allows to set the polarization measurement basis.The outcomes of the APDs are used by a computer to update the posterior distribution and optimize, for the next measurement, the tunable parameters consisting of the activated qplate and the polarization basis.activated q-plate.The prepared state is then sent to a measurement stage composed of the same set of q-plates in reverse order, allowing to reconvert the OAM state into a polarization state.Such measurement stage can be rotated by an angle θ ∈ [0, π) by means of a fully motorized rotation cage, as realized in Ref. [52].In this way, the photon state passing through the full setup becomes the vector |ψ s (θ)⟩ defined earlier, where s = 2q + 1 is the total angular momentum of the photon.By appropriately choosing the active q-plate, determining m, the frequency of the oscillation interference fringes can be tuned.Finally, through a half-waveplate after the q-plate, it is possible to select also the measurement polarization basis.Given the employed devices, we have access to states with s = 1, 2, 11, 51.The first is obtained when no q-plate is activated, while the others are achieved activating in turn one and only one of the three mounted q-plates in the preparation stage, and the corresponding one in the receiver.These states generate oscillation patterns, retrieved from measurements in polarization basis and characterized by visibilities V s which depend upon the selected s, this is schematically represented by the insert "Model" in Fig. 1.

EXPERIMENTAL RESULTS
Here, we investigate the two scenarios mentioned in the introduction: in the first one the interest is focused only on the rotation angle while the visibilities involved in the process are only considered as unavoidable noises.In order to achieve a reliable estimation of the rotation, it is indeed useful to keep track of all the visibilities values which are therefore treated in such a scenario as nuisance parameters.In the second scenario inspected, we consider also the visibilities as interesting parameters.This happens for instance if the user is interested in the complete characterization of an interferometer, and therefore needs an estimation of the noise levels too.In this context, depending on the selected configuration, the optimal multiparameter estimation protocol can be found by assigning positive weights to both the rotation angle and the visibilities.Whatever the scenario, the optimization is performed in a Bayesian framework that follows the strategy developed in [60].Our procedure is however very different when it comes to counting the resources.A single experiment consumes a single photon, however, we do not count this as the consumption of a single resource but as the use of a number of resources equal to the total angular momentum of the photon in that experiment.This is done to set a correspondence between the photon states and the entangled N00N states with sizes equal to the available total angular momentum of the photons.Our N00N-like states can achieve quantum-limited precision, similarly as multipass protocols do [42].It is worth noting that, similar to multipass protocols, which quantify the resources invested in estimation based on the number of interactions between the probe and the sample, it is reasonable to consider the total angular momentum as a valuable resource in the estimation protocol.There are several reasons for this choice.Generating and measuring higher-order orbital angular momentum requires more complex devices such as qplates with higher topological charges or spatial light modulators with larger effective areas.Moreover, the propagation of these states requires also to address their divergence.
We start measuring the angle θ while treating the visibilities as nuisance.To this end, we set the weight matrix G in Eq. ( 1) to have G 11 = 1 as the only non-null entry.The online execution of the Bayesian algorithm is only simulated.This means we first run a stage where we perform a large number of measurements for all the possible q-plates, polarization basis, and rotation angles θ j , and then we queue the outcomes.In the second stage, in running the Bayesian analysis offline, we evoke the needed measurement at each step from the appropriate queue.We collect such outcomes for J = 8 different rotation angles in [0, π) (reported in the Table of Appendix D) , which is the periodicity interval of the system, reporting in panel a) -b) and c) of Fig. 2 the experimental median squared error for some of the different investigated scenarios.The optimized median square error when the rotation angle is the only parameter of interest is reported in panel a) of Fig. 2 where the experimental data are compared with the Cramér-Rao bound on the median of Eq. ( 2), the SQL and the HL [62] (both corrected to refer to the median error).Along with the precision for the optimized Bayesian strategy we report, as a comparison, also the experimental precision for a non-optimized strategy (violet line in Fig. 2a)) in which the activated q-plate and polarization measurement basis are randomly chosen.While it shows for large resources an advantage with respect to the SQL, it is in the transient much worse than the optimized strategy.
In the Bayesian estimation the prior distributions on the visibilities is uniform in [0, 1] and the prior on the angle is likewise uniform in [0, π).With the goal of minimizing the posterior variance on the estimation of the rotation angle, the adaptive algorithm selects at each step the most appropriate quantum-like (q-plates) resource among the available ones.The resources that the Bayesian algorithm selects for the phase estimation are reported in panel d) -e) and f) of Fig. 2. For N > 10 4 only the q-plate with larger topological charge q is used, while for N very small only s = 1 is selected.In between, we have a transient regime where multiple q-plates are used together.This is the region that benefits the most from the optimization.
The experimental results show that the obtained error approaches the computed bound, which proves to be a valid reference even if non-tight.Notably, for the evaluation of the phase, even if the visibility values are completely unknown, the implemented multiparameter protocol shows an enhanced estimation precision compared to the SQL for a large resources range, previously unexplored by multiparameter estimation experiments.Importantly, differently from [52] (see Appendix C for more details) where the measurement strategy has been precalibrated according to the visibility values, here we show that it is still possible to obtain a large region showing sub-SQL performances in the medium-large resource range (1 < N < 10 4 ), treating the 4 visibilities as nuisance parameters.
We now consider the scenario where the visibilities are not treated anymore as nuisances but count in the evaluation of the error.This happens for instance if the user is interested in the full characterization of the apparatus and its configuration, therefore needing an estimation of the noise levels as well as of the rotation angle.It is interesting to see how the optimization of the available resources changes in these new configurations.In particular, we shall focus on the scenario where one is interested in the estimation of the angle and one of the visibilities (the i 0 -th one), while the three remaining ones are still treated as nuisance parameters.Under this assumption, the parameters of interest are collected in the pair (θ, V si 0 ) and the median error of Eq. ( 1) is computed with the weight matrix having G 1,1 = G i0,i0 = 1 as the only nonzero elements.As shown in Fig. 2, under these circumstances the Bayesian protocol, although switching to higher dimensional OAM states to decrease the error on the rotation angle, continues to use also the q-plate related to the visibility chosen: this is necessary to obtain a good precision on the joint error.The results of the estimation of each of the four possible couples are reported in Fig. 2b, c, d and e.In particular, the plateaus in Fig. 2d and Fig. 2e highlighted in green appear since, for few resources, the q-plates corresponding to s = 11 and s = 51 are not significantly used, therefore, the estimator of the corresponding visibilities remains the mean FIG. 2. a) Median square error of Eq.( 1) for the phase only (blue line), with all the other parameters treated as nuisances.The green line is the SQL= 1/N , and the red one is the HL= π 2 /N 2 , both converted to bound on the median square error.The violet line is the precision for non-adaptive measurements.In all the plots the black line is the bound of Eq. ( 2) computed for the appropriate weight matrix G. b)-c) Median square error for the phase and one of the visibilities, respectively Vs 1 , Vs 2 .For each plot, all the other three visibilities are treated as nuisance parameters.In light red we plot the 99% confidence interval for the experimentally evaluated median error.d) -e) -f) Frequency of use of each q-plates as a function of the number of resources, computed for a batch of 8 • 10 3 experiments with the measurement optimization described in the main text.
value of the uniform prior distribution in [0, 1], i.e. 0.5, while the error on the phase decreases, thereby reaching a plateau in the error determined by the value of the visibility V i0 itself.This changes only when the algorithm starts to use the high charge q-plates, and the error finally decreases.In Fig. 3 other scenarios where the visibilities are not treated as nuisance parameters but become themselves parameters of study are reported.It is interesting to see how the optimization of the available resources changes in these new configurations.When we are interested in the estimation of the rotation angle and one of the visibilities, the Bayesian protocol, although switching to higher dimensional OAM states to decrease the error on the rotation angle, will continue to use the q-plate related to the visibility chosen (see Fig. 3 panel d)-e)).In particular, the initial plateau appearing in Fig. 3 a)-b) are due to the fact that for few resources the q-plates corresponding to s = 11 and s = 51 are not significantly used and therefore the estimator of the corresponding visibilities remains the mean value of the uniform distribution in [0, 1] (the prior), i.e. 0.5, until when the algorithm starts to use the high charge q-plates, and the error finally decreases.Lastly, we try to estimate all the five parameters (θ, V 1 , V 2 , V 3 , V 4 ), that means setting all the diagonal elements of the matrix G in equal to 1, obtaining the results reported in panel c) and f) in Fig. 3.
As a final observation, we notice that in the case where one focuses only on the phase θ there is still a pretty large region in which the slope of the optimized estimation error scales with sub-SQL performances (specifically looking at Fig. 2a, this happens for values of N between 2000 and 5000).The region where sub-SQL scaling (in N ) can be observed is defined by the maximum amount of the topological charge q employed in the apparatus.This region can be thus extended by increasing the maximum value of q.Such behaviour on the contrary is damped (but still present) when considering the estimation of the pairs (θ, V si 0 ).The reason for this is that the visibility is an inherently classical parameter and cannot benefit from the high angular momenta being our quantum resources.

DISCUSSIONS AND CONCLUSIONS
In summary, quantum sensing promises to be one of the first quantum technologies exploited to enhance tasks with respect to what is achievable with classical resources.Most of the realistic metrological problems involve more than one unknown parameter, which led to the birth of multiparameter quantum metrology.In this context, a fundamental problem is to optimally allocate the finite available resources, depending on which parameters are treated as nuisance noises and which are the parameters of interest, in order to unlock wider regions of sub-SQL scaling.In this work, we accomplished both these tasks considering the scenario where the parameter of interest is either only the rotation angle or the angle and the fringes visibility.We experimentally showed that this approach is able to reach sub-SQL performances on the estimation precision even when unknown nuisance noises are present, for a FIG. 3. a)-c) Median square error for the phase and one of the visibilities, respectively Vs 3 and Vs 4 (blue lines).For each plot, all the other three visibilities are instead treated as nuisance parameters.c) Median square error on the joint estimation of the phase and all four visibilities.In light red, we plot the 99% confidence interval for the experimentally evaluated median error.d)-f) The frequency of use of each q-plates as a function of the number of resources, computed for a batch of 10 3 experiments with the measurement optimization described in the main text.d)-e) refer to the joint estimation of the rotation angle and one visibility, respectively Vs 3 and Vs 4 .Plot f) refers to the joint estimation of the rotation angle and all the visibilities.resources range O(30, 000).The obtained results have shown the possibility of extending the advantages of multiparameter quantum metrology in the large resource domain.On the other, the methodology here described can find application in large varieties of experimental platforms for quantum sensing with analogous noise models, thus representing a tool for future generations of quantum sensors, that can be employed to boost sensing scenarios like the dynamical tracking of biological reactions [23], but also for estimating the relative rotation of communicating stations [63].Finally, another key advantage of this approach is that it enables the system to be adaptable and resource-efficient, tailoring the estimation process to the specific characteristics and parameters of interest in various experimental setups, including more complex scenarios involving multiple entangled sensors [64].This adaptability allows for improved accuracy and scalability, making it a valuable tool in the pursuit of ultra-precise sensing across diverse applications.

ACKNOWLEDGMENTS
The Bayesian data analysis has been programmed with the Python framework PyTorch and ran on a GPU.The code can be found on github [65].We gratefully acknowledge computational resources of the Center for High Performance Computing (CHPC) at SNS.This work is supported by the ERC Advanced grant PHOSPhOR (Photonics of Spin-Orbit Optical Phenomena; Grant Agreement No. 828978), by the Amaldi Research Center funded by the Ministero dell'Istruzione dell'Università e della Ricerca (Ministry of Education, University and Research) program "Dipartimento di Eccellenza" (CUP:B81I18001170001) and by MIUR (Ministero dell'Istruzione, dell'Università e della Ricerca) via project PRIN 2017 "Taming complexity via QUantum Strategies a Hybrid Integrated Photonic approach" (QUSHIP) Id. 2017SRNBRK.

Appendix A: Derivation of multiparameter precision bounds
In this section, we derive the multiparameter precision bound of Eq. ( 2) in the main text.In order to do so, we start from the quantity that we simulate, that are the measurements outcomes.The greedy algorithm selects for each photon consumed in the experiment the best value of s among the available ones and the polarization basis b that give us the maximum information gain.We called r N the list of tuples containing these choices, together with the relative measurement outcome o i.e.
Notice that the number of total resources used N and the number of measurements K, i.e. the number of photons are different.We also call this string the trajectory of the estimation run.The most widely employed figure of merit for the precision of an estimator is the squared deviation from the true values of the parameters.Given the estimators θ(r N ) and V (r N ) si for θ and V si respectively, and a weight matrix G that codifies which parameters are of interests, we have introduced in the main text the error quantity We then take the expectation value of this precision on the prior distribution for θ and define ∆ 2 We can approximate these expressions by means of M simulations for each of the J angle θ j .So that we can write where r m,j N is the trajectory of the m-th experimental run for the j-th angle.The figure of merit in Eq. (1) of the main text is computed by taking the median of the M quantities J j=1 ∆ 2 r m,j N ,G (θ j ) instead of the mean.We now see how the Cramèr-Rao (CR) bound sets a limit to ∆ 2 G and how this can become a bound for the median error M 2 G .Depending on the measurement basis chosen, the results o 1 , o 2 ∈ {−1, +1} of the two polarization measurements for the i−th q-plate are distributed respectively according to We have ν i measurements for the q-plate s i in total, that we assume being evenly split between the two polarization basis.From these probabilities we can write the 5 × 5 Fisher information matrix I (FI matrix) for the five parame-ters (θ, V s1 , V s2 , V s3 , V s4 ), whose non-zero elements are , for i = 1, 2, 3, 4, and I i+1,1 = I 1,i+1 for symmetry.The Cramér-Rao bound, holding true for asymptotically unbiased estimators, is then expressed by the following inequality involving ∆ 2 G (θ): by taking the expectation value on the prior on θ we have (A6) We now want to renormalize the uses of each q-plate ν i in such a way to highlight the dependence on the total number of resources N , i.e. ν i := x i N .The FI matrix I becomes I = N I, where the entries of I are similar to that of I, only that ν i is substituted with x i .The CR bound reads now where the expectation value of the matrix I is diagonal with entries and C G is the solution of the following minimization problem: In order to get a reference value for the median square error, we suppose that the estimators θ(r N ) and V (r N ) si are asymptotically normal and unbiased, we give some empirical evidence of this later in this section.Because the square of the deviations | θ(r N ) − θ j | 2 and | V (r N ) si − V si | 2 are left-skewed and independent variables, we observe that The variable θ(r N ) − θ j is normally distributed and centred in zero.Under this hypothesis it is easy to show that the median of its square is proportional to the variance with a factor ξ ≃ 0.4549 that can be estimated numerically.Therefore, the bound on the median error of the estimation is with C G = ξC G .We will now briefly prove the relation between the median and the mean square error.Under the aforementioned assumptions we have θ(r N ) − θ ∼ N 0, σ 2 .From the formula for the transformation of the probability distributions we compute the probability density function for y := | θ(r N ) − θ| 2 , which is The mean value of y is by definition σ 2 , i.e., By introducing k = y σ 2 and M = ξσ 2 and changing variable in the above integral we get In this integral the only unknown is ξ, and it can be solved numerically to get Eq.(A12), since by definition xi is the proportionality factor between the median and the variance.In Fig. 4 we report the comparison of the estimator for the rotation angle produced by the Bayesian procedure with a reference Gaussian.This justifies empirically the affirmation we made in the paper that the distribution for the estimator is asymptotically a bias-less Gaussian affected by outliers.What we measure with the median square error is the width of the central Gaussian-like part of the blue distribution, while ignoring the outliers.

Appendix B: The Bayesian algorithm
In this section we complement the presentation of the Bayesian algorithm that we gave in the main text, with the application to our q-plates setup in mind.With respect to the original formulation, we made a few corrections necessary because of the circular nature of the angular variable that we are FIG.4. Probability density of the error for the rotation angle estimator (blue bars), compared with a reference Gaussian distributions (orange bars).This plot is for a estimation with N = 3000 resource and a batchsize b = 3000 (total number of runs).The outliers of the estimator distribution (which appear even further away from the central peak than the ones showed) are very impactful on the mean square error, while not affecting at all the median square error.going to measure.In every experiment, n p = 5000 particles have been used.The parameters to estimate are collected in the vector x := (θ, V s1 , V s2 , V s3 , V s4 ), that contains the phase in the first entry and the four visibilities in the other ones.The Granade's approach [60] is based on a particle filter, which represents on a computer the posterior probability distribution with the ensemble E := {(x k , w k )}, of n p pairs composed of a particle k in the position x k and with a weight w i .
The weights are updated after each measurements with the Bayes rule, i.e.
and depend therefore on the measurement step t.The probability p(o t |s i , x k ) is either p 1 or p 2 according to the basis of the polarization measurement chosen.The weights w k t+1 are renormalized after each update.In the following, we will drop the dependence on the measurement step t except when needed.The j-th component of the k-th particle of the ensemble will be represented as x k j , and could correspond to the phase if j = 0, that is x k 0 = θ k or to one of the visibilities if j = 1, 2, 3, 4, that is x k j = V k sj .The mean of the angular values is computed as while the mean values of the visibilities are Together they form the vectorial mean of the posterior distribution μ = (μ 0 , μ1 , μ2 , μ3 , μ4 ).Its components will also be our estimators for the parameters i.e. θ := μ0 and Vsi := μi .The covariance matrix is defined as where for i, j = 1 we use the circular distance, i.e.
while for i, j > 4 we use the regular distance The scalar variance of the posterior distribution is by simulating each one of the possible 8 experiments (4 qplates and 2 polarization basis among which we have to choose) and selecting the one with the lowest expected variance we realize the greedy optimization described in the main text.In particular, we compute for each combination of qplate and polarization basis the expected scalar variance, i.e.
where σ 2 (o) is the variance of the posterior distribution where the weigh update of Eq. (B1) is performed assuming the outcome o for the measurement.This variance is weighted with an estimation for the probability of getting the outcome o.In this way the algorithm attempts to concentrate as much as possible the distribution around its mean, without however planning for more than one step ahead in the future.The simplest possible non-greedy extension of this would be to simulate two measurement steps, this would mean computing 64 possible expected variances.Simulating many more steps becomes quickly unfeasible, and a more refined approach, possibly involving machine learning is required.The resampling strategy of the Granade procedure in [60] has also undergone minor changes to adapt it to the phase estimation problem.
based on a particle filter, which represents on a computer the posterior probability distribution with the ensemble E := {(x k , w k )}, of n p pairs composed of a particle k in the position x k and with a weight w i .
The weights are updated after each measurements with the Bayes rule, i.e.
and depend therefore on the measurement step t.The probability p(o t |s i , x k ) is either p 1 or p 2 according to the basis of the polarization measurement chosen.The weights w k t+1 are renormalized after each update.In the following, we will drop the dependence on the measurement step t except when needed.The j-th component of the k-th particle of the ensemble will be represented as x k j , and could correspond to the phase if j = 0, that is x k 0 = θ k or to one of the visibilities if j = 1, 2, 3, 4, that is x k j = V k sj .The mean of the angular values is computed as μ0 := arg np k=1 while the mean values of the visibilities are μj = np k=1 Together they form the vectorial mean of the posterior distribution μ = (μ 0 , μ1 , μ2 , μ3 , μ4 ).Its components will also be our estimators for the parameters i.e. θ := μ0 and Vsi := μi .The covariance matrix is defined as where for i, j = 1 we use the circular distance, i.e.
while i, j > we use the regular distance The scalar variance of the posterior distribution is by simulating each one of the possible 8 experiments (4 q-plates and 2 polarization basis among which we have to choose) and selecting the one with the lowest expected variance we realize the greedy optimization described in the main text.In particular, we compute for each combination of q-plate and polarization basis the expected scalar variance, i.e.
E σ 2 := o=±1 p 1/2 (o|s i , Vsi , θj )σ 2 (o) , where σ 2 (o) is the variance of the posterior distribution where the weigh update of Eq. ( 17) is performed assuming the outcome o for the measurement.This variance is weighted with an estimation for the probability of getting the outcome o.In this way the algorithm attempts to concentrate as much as possible the distribution around its mean, without however planning for more than one step ahead in the future.The simplest possible non-greedy extension of this would be to simulate two measurement steps, this would mean computing 64 possible expected variances.Simulating many more steps becomes quickly unfeasible, and a more refined approach, possibly involving machine learning is required.The resampling strategy of the Granade procedure in [1] has also undergone minor changes to adapt it to the phase estimation problem.

Details on the data analysis
In this protocol, we can not fix the total number N of resources, because the stochastic nature of the measurement outcomes propagates to the online choice of the next q-plate and then to N , which is therefore a stochastic variable whose exact value is only known at runtime.However, we can repeat M ≫ 1 times the simulation and collet the precision after each used photon in the set of tuple {(∆ 2 rn,G (θ j ), n)}, where n is the total resource number used up to that point to get to precision ∆ 2 rn,G (θ j ).We will then look at all the points with total resources falling in the interval [n, n + ∆n] with ∆n small, and the median error of this cluster is the error we associate to n.We choose to cluster only the point having n > 100, with the chosen interval being ∆n = 50.