Experimental investigation of Bayesian bounds in multiparameter estimation

Quantum parameter estimation offers solid conceptual grounds for the design of sensors enjoying quantum advantage. This is realised not only by means of hardware supporting and exploiting quantum properties, but data analysis has its impact and relevance, too. In this respect, Bayesian methods have emerged as an effective and elegant solution, with the perk of incorporating naturally the availability of a priori information. In this article we present an evaluation of Bayesian methods for multiple phase estimation, assessed based on bounds that work beyond the usual limit of large samples assumed in parameter estimation. Importantly, such methods are applied to experimental data generated from the output statistics of a three-arm interferometer seeded by single photons. Our studies provide a blueprint for a more comprehensive data analysis in quantum metrology.


I. INTRODUCTION
The aim of quantum metrology is to study how the use of quantum phenomena can be beneficial in parameter estimation [1][2][3][4][5][6][7][8][9]. As the technological level of quantum sensors increases, ensuring reliable and efficient data processing becomes a crucial aspect in the operation of devices. While in proof-of-principle demonstrations the enhancement in precision can be shown with simple data analysis, for real applications optimal use of the available data is required [10]: this is instrumental to reach the necessary accuracy, as well as the ultimate limits on precision.
The estimation process relies on two sources of information. On the one hand, we may have some a priori information available on our parameter from modelling or preliminary measurements. On the other hand, the collected data allow to refine our previous knowledge. Bayesian data analysis thus provides the most natural setting to combine both contributions to deliver a high-quality estimator [11,12]. The standard approach to assess its performance in terms of precision is to compare the experimental variance with that of the Cramér-Rao bound [11]. This is determined by the Fisher information, which indicates the amount of information embedded in the output probability of the quantum sensor. However, this approach is valid only for local estimation, i.e. in a small range around a known value of the parameter. This may not be the case in general, and suitable methods to achieve local conditions are needed [13][14][15][16][17][18][19][20][21][22][23][24]. More exhaustive data processing thus has to incorporate cases with relatively broad a priori information. This is particularly relevant when limited metrological resources can be invested for the measurement [25][26][27][28][29].
The interplay between a priori information and collected data becomes even more complex in the multiparameter case. Here, the estimation procedure is affected not only by the initial uncertainties on the individual parameters, but also on their correlations. Indeed, the multiparameter Fisher information sets a given degree of statistical correlation between parameters [2,30], which may contrast with the one obtained from a priori considerations, especially in its orientation in the parameter space to be estimated.
Adequate bounds on the variance of Bayes estimators have been introduced by Ziv and Zakai (ZZ) [31], and by Van Trees (VT) [32], but have never been tested on data from a quantum sensing experiment. In this article, we show the applicability of the ZZ and VT bounds in multiphase integrated photonic quantum sensing. This is emerging as one of the most solid technologies for the realisation of quantum sensors [33][34][35][36]. Our results show that both bounds capture the limits in multiphase estimation, especially for limited resources. The VT approach, besides being less computationally demanding, provides a tighter bound on the experimental variance. The ZZ approach still provides a useful reference, and captures the dependence on the variance on the correlations present in the a priori information. Our findings contribute to shaping actual data analysis for future real-life quantum sensors.

II. BAYESIAN MULTIPARAMETER ESTIMATION WITH A PRIORI KNOWLEDGE
Bayesian estimation theory [11,12,37] represents a paradigm for parameter estimation protocols and has been extensively and fruitfully employed [8,17,18,[36][37][38][39][40][41][42][43][44][45][46]. In general, in multiparameter problems one has to deal with a set of physical parameters φ = (φ 1 , . . . , φ n ), whose (unknown) values have to be estimated by the user. Typically, information on such set φ is obtained by preparing a set of N probe states ρ 0 of an auxiliary system. The evolution of the auxiliary system depends on the set φ via a unitary operation U (φ), or more generally via a physical map L φ . Measurement Π x of the N probe states ρ φ after evolution provides information on the parameters φ.
Such process can be embedded in a Bayesian framework, where the parameters φ to be estimated are treated as random variables distributed according to some distribution p(φ), which encodes the actual knowledge on the parameters values (see Fig. 1). Bayesian estimation permits to naturally encode in such framework the availability of an arbitrary a priori information on the parameters, that may occur in several scenarios. This initial knowledge is encoded in the so-called prior distribution A(φ), which acts as the starting point for the estimation process. After preparing N probe states, the measurement outcomes x = (x 1 , . . . , x N ) are collected. Here, this vector represents the set of the N measurements obtained by sending the sequence of inputs, where the outcome on a single probe x k = 1, . . . , s can take one of s possible values according to the system behavior. Prior knowledge is updated according to the Bayes' rule: p(φ|x) ∝ p(x|φ)A(φ).
Here p(x|φ) is the conditional probability of detecting the sequence x of outcomes provided a given set of values φ for the parameters. For N independent probes, and individual measurements, such conditional probability is expressed as p(x|φ) = N k=1 p(x k |φ), where p(x k |φ) is the likelihood of the given outcome x k for a single copy of the probe.
The output distribution p(φ|x) fully encodes all the final knowledge on the parameters at the end of the estimation process. A common choice for an estimator of φ is provided by its expectation valueφ = φp(φ|x) dφ, that can be proven to be an asymptotically efficient and unbiased estimator [47][48][49]. In single-parameter estimation, the precision can be quantified by the mean square error MSE(φ) = x p(x|φ)(φ − φ) 2 . This generalizes in the multiparameter scenario to a matrix C(φ), with elements C ij (φ) = . Such quantities substantially encode the precision in the estimate with respect to the true values of the parameters. Different theoretical studies focused on the extension of the Bayesian framework to the multiparameter case, also paying attention to the limited data regime [27,[50][51][52][53].
In general, the mean square error (or its corresponding matrix form for a multiparameter scenario), is not accessible by the user, given the inherent absence of knowledge of the true parameter value. Bayesian approaches provide an additional framework to describe the estimation process. In particular, for a given experimental instance with N probe states and the string x of measurement outcomes, the posterior distribution p(φ|x) encodes the degree of confidence resulting from the estimation process via its covariance matrix Σ, whose elements are given by The diagonal terms of Σ represent the variances of the individual parameters. Conversely, off-diagonal terms contain information about the correlation between the parameters. In the limit of large N , given that Bayesian estimation is unbiased, it can be shown that the covariance matrix Σ and matrix C asymptotically coincide.
Different quantities can be used to bound the estimation error, as a function of the resources N employed in the process. The most common bound applying to the covariance matrix is the celebrated Cramér-Rao bound [11] where N is the number of repetitions of the experiment, and the Fisher information matrix F (φ) has entries F i,j (φ) = x ∂ φi p(x|φ)∂ φj p(x|φ)/p(x|φ). This limit holds for local estimation: this means that the measurement assumes previous knowledge of the parameters φ, up to some uncertainty we aim at improving. Such an improvement may come from performing adaptive measurements, or, as it is often the case, by increasing the number of events N . Referring to the discussion above, in this large N limit, the width of the posterior distribution has little to do with the a priori distribution A(φ), and is governed uniquely by the Bayesian update.
When a limited number of repetitions are available, the Cramér-Rao limit may indeed not provide relevant information. However, if the prior distribution is regular and derivable, there exists a simple result due to Van Trees, which states that the matrix sets a limit to the covariance as which is called Van Trees (VT) bound [54]. This considers how both the average Fisher information available following the measurement and the prior enter in the final information; when N is large, the usual Cramér-Rao bound is recovered, as the contribution from the a priori distribution becomes negligible.
For the VT bound to exist, A(φ) must be regular; in order to obtain a valid bound also for generalised distributions, we adopt the Ziv-Zakai (ZZ) bound [31], in the multiparameter form proposed by Bell et al. [55]. Differently from the VT bound, the ZZ bound is written explicitly in a scalar form, by introducing a unit vector u in the parameter space, and considering the error u Σu: We can define the error probability of distinguishing φ from φ as: based on our measurements and the probabilities π 0 = 1 − ). These expressions hint at how the ZZ method adopts binary hypothesis testing as a way to bound the error, and avoids the need of a regular a priori distribution.

III. PLATFORM AND METHODOLOGY FOR BAYESIAN MULTIPHASE ESTIMATION
We have tested Bayesian estimation with a priori knowledge in an experimental platform provided by an integrated multiarm interferometer realized via the femtosecond laser writing technique [56] on a glass chip, where the unknown parameters are provided by two relative phase shifts inside the structure (see Fig. 2). In this way we apply the Bayesian framework and study the corresponding bounds for the paradigmatic problem of multiphase estimation, having several applications to sensing and microscopy [33,35,[57][58][59][60][61][62]. More specifically, the interferometer is the three-mode generalization of a Mach-Zehnder one [34,35]. This is realized by replacing the input and output beam-splitters with two cascaded tunable tritters [63]. The output state of the interferometer will then depend on two different relative phase shifts between the internal modes. Within the circuit, phases can be actively reconfigured via thermo-optic phase shifters, obtained by fabricating resistive microheaters on top of the interferometer [35,36,64]. In the present device, resistors R A and R B are used to fine-tune the input and output tritter transformations. The phase shifts for the estimation process are the two relative phases (φ 1 , φ 2 ) referred to the central internal arm of the interferometer. These phases can be modified via the set of thermo-optic shifters corresponding to resistors R 1 − R 6 and whose response functions can be also calibrated using machine learning techniques [65]. Redundancy in such set of resistors can be used to implement adaptive estimation protocols, as shown in Ref. [36]. Single-photon probes are conditionally generated from a two-photon source based on the spontaneous parametric down-conversion process. Namely, one of the two generated photon is directly measured, acting as a trigger for the experiments, while the other photon is injected in the integrated device via a single-mode fiber array in the first input port. After evolution through the interferometer, the output modes are coupled to an array of multimode fibers and measured via single-photon detectors.  ). b, Scheme of the reconfigurable interferometer employed for multiphase estimation. The first transformation UA is a balanced tritter, and includes a thermo-optic phase shifter placed on resistor RA used to adjust the transmittivities for balanced operations. The phases φ1 and φ2 to be estimated, which can be changed via the set of resistors R1 − R6, are the two relative phases with respect to the central mode, acting as a reference. Output transformation UB is a second balanced tritter, which can be fine-tuned via resistor RB, and it is used to recombine the output modes before the measurement.
As a first step, prior to the estimation experiments, singlephoton states are employed to characterize the response of the interferometer. In particular, this amounts to experimentally reconstruct the correspondence between the current propagating on each resistor and the applied optical phase on the waveguide below, including potential cross-talks between the different tunable phase shifts. Furthermore, the same measurements are employed to reconstruct the transmittivity values for each directional coupler. This characterization is performed via the procedure reported in [35,36]. More in details, single-photon input states are sent in the interferometer, and the output statistics for different values of the current are recorded. These data are then employed to perform a fitting procedure, using a mathematical model for the interferometer. Such a model includes both static parameters such as directional coupler transmittivities (T j ) and static phases obtained at zero current in the resistors (φ 0j ), and the dynamical response function φ l = φ l ({i j , R j }), being R j the resistances and i j the applied currents [35,36]. This procedure finally leads to knowledge of the likelihood function P M L (x|φ 1 , φ 2 ).
Having access to such response function has a fundamental role for Bayesian estimation experiments, since this approach is based on progressively updating the posterior distribution according to the measurement outcomes and the likelihood function. Indeed, an accurate knowledge of the system likelihood is required to avoid biases or additional errors in phase estimation processes, and is a common need shared by quantum sensor platforms based on different technologies.
Then, phase estimation experiments are performed on the same system, characterizing the process for different values of the prior distribution A(φ 1 , φ 2 ). The posterior distribution after each experiment with N resources is updated from Bayesian update as described above, and carries all the available information on the unknown parameters. Relevant quantities, such as the estimate as well as the confidence interval related to the covariance matrix Σ, are extracted from this function. Bayesian protocols require evaluation of multidimensional integral, whose dimension is given by the number of unknown parameters. Such integrals become progressively more difficult to handle from a computational point of view. To address such issue, a possible solution is provided by performing an appropriate discretization of the parameters space, which can be performed by using the particle approximation [66]. Such a technique corresponds to replacing the continuous parameter space with a discrete set of M p points y i = (φ 1,i , φ 2,i ), named particles. Each particle has its own weight w i , related to the conditional probability associated with the corresponding values y i . More specifically, the initial particle set is obtained by randomly generating M p values y i according to the prior distribution. The corresponding initial weights are uniformly set to w i = 1/M p . Bayesian update is then performed by changing at each step the weights of the particle according to the Bayes' rule. Thus, at each step k the weights of the particles are updated as w i → w i P M L (x k |φ 1,i , φ i,2 ), followed by proper renormalization. After the N steps of the estimation have been performed, the estimates are obtained as: while the covariance matrix can be estimated as: In our case, given the intermediate N regime and the intrinsic periodicity of the phase parameters, we employed the corresponding circular counterparts for the estimate and the covariance matrix [67]. For a two-parameter space we employed a value of M p = 1600, which is sufficient to provide a good approximation of the estimation experiments in the investigated regimes. Furthermore, we applied the resampling technique [68] throughout the data analysis. Such technique is employed to avoid a common issue that may arise in particle approximations. In particular, for progressively larger N and corresponding lower uncertainty on the parameters, most of the weights w i will be zero-valued. In this case, the set of parti-cles is no more informative on the experiment, and the quantities evaluated from Eqs. (6-7) fail to provide a good approximation of the exact values. To avoid such issue, a new set of particles is periodically resampled [68] throughout the estimation experiments. This method consists in checking, at each step of the protocol, the concentration of the weight values. If 1/ i w 2 i < M th (in our case, the choice is M th = M p /2), a new set of particle positions and weights ({y i }, w i ) is generated according to the current knowledge on the parameters. More specifically, M p new particle positions are generated by random choice from the current ones {y i } according to the weights w i . Then, the new set of particle positions is perturbed according to a multivariate normal distribution having mean µ i = ay i + (1 − a)φ and covariance matrix Σ = (1 − a 2 )Σ. Here,φ and Σ are the mean and covariance of particles before resampling ({y i }, w i ) as defined in Eqs. (6)(7), and a is a parameter included in the interval a ∈ [0, 1]. In our case, we used a = 0.98 according to Refs. [36,66]. Once the new particle positions are determined, the new weights are reset to uniform w i = 1/M p .

IV. RESULTS
Two-phase estimation is performed by Bayesian updating of the a priori distribution A(φ 1 , φ 2 ) based on the reconstructed probabilities P M L (x|φ 1 , φ 2 ). In the asymptotic limit of large samples N 1, the only relevant figure to assess the uncertainty of the measurement is the local Fisher information F (φ 1 , φ 2 ): usually, in this regime the output posterior p(φ 1 , φ 2 |x) is sufficiently narrow that details on the prior distribution become irrelevant -except for ruling out possible ambiguities occurring however at a large scale. With a smaller sample, N 200, both P M L (x|φ 1 , φ 2 ) and A(φ 1 , φ 2 ) are relevant: this regime is the one we will be exploring.
We first consider an a priori distribution in the Gaussian form: , where µ 1 and µ 2 are the mean values, and Γ captures the initial uncertainties and the correlations on the parameters. For equal uncertainties Γ 1,1 = Γ 2,2 , the width of the prior distribution is quantified by a single parameter σ = Γ 1,1 . In particular, we quantify the correlations by means of the Pearson coefficient [69] ρ = Γ 1,2 /(Γ 1,1 Γ 2,2 ) 1/2 . Similarly, we can introduce a quantity for the Fisher information ν = −F 1,2 /(F 1,1 F 2,2 ) 1/2 , where the additional minus sign here originates from the inversion of F in order to obtain a covariance matrix.
The relevant figure in our case is the total variance V = Tr [Σ]. Fig. 3 shows the scaling of this variance on the phases as a function of N . The points, corresponding to our estimates based on the experimental probabilities, are efficiently captured by both the ZZ and the VT bounds, with the former becoming marginally looser for increasing N . Since we have to consider two parameters, the minimal variance at the ZZ bound V ZZ writes V ZZ = Z(u 1 ) + Z(u 2 ), where u 1 and u 2 are orthogonal unit vectors. We found that the tighter bound on V corresponds to choosing u 1 and u 2 as the eigen- vectors of the covariance matrix Σ, since these are associated to independent parameters.
A first case in this plot considers the scenario of a symmetric A(φ 1 , φ 2 ), Γ 1,1 = Γ 2,2 , presenting no correlation between the two parameters, ρ = 0. Other possibilities are illustrated in Fig. 3, for the same widths, but different values of ρ: the bounds are still able to capture the variance efficiently, but it must be remarked that the uncertainties themselves are varied. In particular, it can be verified that the minimal uncertainties are achieved when the correlations in the prior distribution match those of the Fisher information: ρ = ν, i.e. in the absence of competing symmetries in the updated probability. The particular case of Fig. 3 exemplifies a general behaviour. The grids in Fig. 4 show that the relative variation between the assessed variance at N = 1 and the corresponding ZZ bound remains below 1%, regardless the initial correlation ρ. These results also provide upper limits for the behaviour of the VT bound, which appears to be tighter (see Fig. 3).
The second key parameter in the a priori is its width, which is varied in the plots of Fig. 5a. Here we show the variance as a function of N , for different values of Γ 1,1 . As expected, the narrower the prior distribution, the more precise the estimate is. Crucially, both the ZZ and the VT bounds become less strict for wider prior distributions. The origin of this behaviour can be traced in the form of the posterior distributions, shown in Fig. 5b-c. When A(φ 1 , φ 2 ) is narrower, it eventually produces a well-behaved a posteriori distribution, with a single peak. In the opposite limit of a wide a priori distribution, the final distribution is affected by ambiguities in the outcome probabilities, on which A(φ 1 , φ 2 ) does not provide enough information for a discrimination. This emphasises the impact of ambiguities in the output distributions. In fact, different phase settings may lead to identical outcomes. The settings can only be discriminated if the a priori distribution is sufficiently narrow, regardless of the number of copies. The bounds, on the other hand, do not account for this increased complexity in the posterior.
In general, the VT bound, which is considerably simpler to calculate, performs better than the more sophisticated ZZ bound. However, this latter presents the advantage of working also for non-derivable a priori distributions, as for the case presented in Fig. 6. We show the bounds as a function of N when the a priori distribution is now in the form of a 2D Heaviside rectangle of width ∆. The bounds are not as accurate as in the analytical cases of regular a priori functions, providing only an estimate for the magnitude of the variance.

V. DISCUSSION
Parameter estimation using quantum resources promises the development of a novel generation of sensors with enhanced precision capabilities. Besides the technological advances, this development must be inevitably accompanied by appropriate data analysis techniques to fully exploit the sensors precision. Notable frameworks to be addressed require going beyond local estimation and asymptotic scenarios, thus dealing with limited resources and different prior knowledge. In this general framework, suitable bounds need to be identified to properly address the performance of quantum estimation experiments.
Here, we have investigated Bayesian multiphase estimation with different prior knowledge in integrated multiarm interferometers. The latter represents a promising platform for quantum sensing, and at the same time can be used as a testbench to develop proper techniques for multiparameter estimation. Within this system, we have investigated the interplay between the natural correlation on the unknown parameters, related to the Fisher information matrix, and the correlation encoded in the prior distribution. Such scenario has been assessed by using two different bounds for Bayesian estimation, that can naturally take into account the availability of different prior knowledge. We observe that both Van Trees and Ziv-Zakai bounds provide a strict estimate of the estimation error for regular and narrow prior distribution. This is observed independently of the relative orientation of the prior distribution correlation with respect to the system Fisher matrix. For larger priors, the bounds become less accurate in the presence of ambiguities in the system output probabilities, which affect the estimation precision. Finally, when the prior distribution is not regular the Ziv-Zakai bound still provides a reasonable, albeit not strict, estimate of the process features.
The present analysis provides a detailed insight on the role of the different quantities underlying multiparameter estimation, in the general regime with limited resources and different prior knowledge. Such results are expected to be of relevance for the development of quantum sensors, and can stimulate further investigation in defining a general and comprehensive framework for multiparameter quantum estimation.