Simultaneous estimation of multiple phases in generalised Mach–Zehnder interferometer

In this work we investigate the problem of simultaneous estimation of phases using generalised three- and four-mode Mach–Zehnder interferometer. In our setup, we assume that the phases are placed in each of the modes in the interferometer, which introduces correlations between estimators of the phases. These correlations prevent simultaneous estimation of all these phases, however we show that we can still obtain the Heisenberg-like scaling of precision of joint estimation of any subset of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d-1$$\end{document}d-1 phases, d being the number of modes, within completely fixed experimental setup, namely with the same initial state and set of measurements. Our estimation scheme can be applied to the task of quantum-enhanced sensing in three-dimensional interferometric configurations.

Multiparameter estimation is a rapidly growing field of quantum metrology [1][2][3][4][5][6][7][8] , in which an initial quantum state undergoes an evolution which depends on several parameters, and the final task of the process relies on estimating these parameters based on measurement statistics of the final state with as low variance as possible. Multiparameter estimation encounters specific difficulties, which do not occur in the single-parameter case 9,10 . One of the most significant difficulties is the possibility of an occurrence of correlations between the estimators corresponding to different parameters, which could decrease the overall precision of joint estimation 3 . Also, in this regard, it is worth mentioning that in Ref. 11 it has been shown that the number of simultaneously estimatable parameters reduces when an external reference mode is absent. The main tool used to evaluate the precision of multiparameter estimation is the generalisation of the Quantum Fisher Information (QFI) 12 to the multiparameter case, known as the Quantum Fisher Information Matrix (QFIM) 5 . For sufficiently uncorrelated parameters the QFIM is invertible and the covariance matrix of the estimated parameters is bounded from below by the inverse of the QFIM. This bound is a multiparameter version of the Quantum Cramer-Rao bound (QCRB) 3 . The elements of the QFIM depend on the size of the initial probe state. If the state is N-partite the elements of QFIM can depend at most quadratically on N, which is known as the Heisenberg-like (HL) scaling 13,14 . Several works were showing that the Heisenberg-like scaling of precision of estimation of all the parameters is possible for an entangled input state and some measurement strategy, which in principle demands the use of arbitrary multiports 1,4,[15][16][17][18] .
In this work we state the problem of simultaneous estimation of multiple phases using 3-and 4-port generalised Mach-Zehnder interferometer 19 . We consider a scenario in which the phases are placed within each of the internal ports of the d-mode interferometer (see Fig. 1a), and the task is to simultaneously estimate any (d − 1) -element subset of them, whereas the remaining one is known, and serves as a phase reference.
Note that such configuration implies that the phases are strongly correlated, and the QFIM for all the phases in the interferometer is singular. The singularity of the QFIM reflects the impossibility of simultaneous estimation of all the phases without an external reference mode 11,20,21 .
We show that with the use of a fixed initial entangled probe state and a fixed interferometer one can obtain the Heisenberg-like scaling of precision of simultaneous estimation of any (d − 1)-element subset of the phases, without any change in the setup. This means that we have both the same initial state as well as the same set of local measurements when estimating each subset.
In a typical approach to quantum phase estimation one treats all the unitary part before the phaseshifts as a preparation of the initial state, whereas all the part after them is treated as an implementation of the measurement. In this work we apply a different point of view, treating the entire interferometer as a fixed single unitary operation, with the aim of investigating the metrological properties of a generalised Mach-Zehnder interferometer as a whole.

Results
Analytical description of generalised Mach-Zehnder interferometer. In this section, we describe our setup for estimation of multiple phases. For the convenience of the presentation, we will use the standard Hilbert space description, which demands the introduction of N interferometers, however, as we discuss in details in "Single-interferometer optical implementation" section, the entire scenario can be as well described using the second-quantised framework, in which only one interferometer suffices. The setup (see Fig. 1b) consists of an N-partite GHZ source, which sends single photons to N measurement stations, such that a d-mode bundle goes to each of the stations, and the basis states are encoded by a path of the photon within the bundle. Therefore the initial state reads: in which the preparation symmetric multiport U d will be specified later. Further each of the N measurement stations apply d-mode interferometers, which consist of generalised Mach-Zehnder interferometers involving d phases, the (d − 1)-element subsets of which will be estimated. A generalised Mach-Zehnder interferometer 19 consists of two symmetric multiports intertwinned with a series of phaseshifts on each of the d modes linking the multiports (see Fig. 1a). By a symmetric multiport we mean a d-mode multiport with the property that a single photon entering by any of the d input modes has a uniform probability 1 d to be detected in any of the d output modes (see e.g. 19 ).
For metrological considerations one usually needs a description of an interferometer in the Hamiltonianlike form U = e −ih i α i , with an explicit form of generators corresponding to phaseshifts. The typical choice of the operator basis for finding the generators h i is the set of Gell-Mann matrices, however we found that much more convenient choice for describing interferometers based on symmetric multiports is the Heisenberg-Weyl operator basis, defined as: where ω = exp(2iπ/d) . In the following we will present explicit form of the generators for d = 2, 3, 4.
Two-mode case. In the two-mode case the Heisenberg-Weyl operators (2) are equivalent to the standard Hermitian Pauli matrix basis. The symmetric two-port can be expressed in the following form: whereas the phase-imprinting part of the interferometer has the following representation: The evolution of the entire interferometer can be expressed in a concise way using only Y operators: The U 2 evolution is generated by two Hamiltonians: the eigenvalues of which read: {0, 1}.
Three-mode case. The description of the 3-mode case in the Heisenberg-Weyl basis turns out to be completely analogous to the 2-mode one. Firstly, the symmetric multiport has the following presentation in terms of generalised X operators: www.nature.com/scientificreports/ whereas the phase part depends solely on generalised Z operators: The entire evolution is, in full analogy to the two-mode case, generated solely by the generalised Y matrices: The U 3 evolution is generated by three Hamiltonians: The above Hamiltonians fulfill h 1 + h 2 + h 3 = 0 , and their eigenvalues read: {2/3, −1/3, −1/3} . Note that despite the fact that in the three-mode case the Heisenberg-Weyl operators (2) are no longer Hermitian, their appropriate combinations give rise to a proper Hermitian Hamiltonians (11).
Four-mode case. The description of a 4-mode case is a bit more complicated. The symmetric multiport still has analogous simple form in terms of generalised X operator: whereas the phase-imprinting part can be presented in the following concise way: where entries of the vector z fulfill relations: and K mn = ω (m−1)(n−1) . The entire evolution is generated by 4 Hamiltonians: however their exact form, although still depending only on generalised Y matrices, is much more complicated than in previous cases: where the star * denotes complex conjugation and for clarity we introduced constants: iπ(X + X 2 ) , www.nature.com/scientificreports/ The eigenvalues of h i are: {3/4, −1/4, −1/4, −1/4} . In this case the same remark applies as in the previous one: properly defined functions of the Heisenberg-Weyl non-Hermitian operators give rise to Hermitian Hamiltonians.

Precision of estimation of multiple phases with generalised Mach-Zehnder interferometer.
In this section, we analyse the estimation precision of our proposed estimation scheme for d = 3, 4 and an arbitrary number N of photons in the initial state. There are two kinds of phases in our experiment: the (d − 1)-element subset of unknown phases to be jointly estimated, and one remaining phase, assumed to be fixed and known, serving as a phase reference. Such a situation, in which the set of all parameters determining the final probability distribution is divided into parameters of interest (the ones we estimate) and additional parameters, is well known in estimation theory 30 (see "Methods" section for more detailed presentation). In our case the additional parameter is the reference phase. Since we assume entirely fixed experimental setup with no optimisation of measurements, we estimate the precision of estimation using classical Fisher Information Matrix techniques. The precision of joint estimation of several parameters in the presence of fixed and known additional parameters is specified by the Cramer-Rao bound based on the inverese of the Fisher Information Submatrix corresponding to the parameters of interest (55). In order to describe the precision of estimation by a single quantity we take the trace of both sides of the Cramer-Rao bound (55) 27 : where in our case the matrix F ( α) I ,I is a Fisher Information submatrix corresponding to the subset of jointly estimated phases (denoted here by a symbol I meaning parameters of interest, see "Methods" section). Therefore our task would be to analyse the behaviour of the quantity Tr F (� α) I ,I −1 as a function of the number of photons N in order to find asymptotic scaling of precision. In our setup we allow for an optimisation of the initial state, which depends solely on the dimension of the local multiport. As an optimisation strategy we take the maximisation of the mean QFI per parameter. We utilise the following inequality: which follows directly from (50) (see "Methods" section). Note that the LHS of the above inequality is just the mean QFI per parameter, which can be treated as an approximate measure of average estimation performance per parameter. The total collective Hamiltonian corresponding to the action of N local interferometers reads: where h i denotes any of the local Hamiltonians from formulas (11) and (16). The inequality (19) implies that the optimal state should be an eigenstate of the operator i H 2 i .
Three-mode case. Optimal state. The N-qutrit optimal state that maximizes the trace of the QFIM F Q is given by: where One can easily prove that fact by noticing that the operator U 3 simultaneously diagonalises the local Hamiltonians (11). Indeed, the operators (11) are expressible solely by the operators Y and Y 2 , and we have the relations: i is also diagonal with the same set of eigenstates.
Achievable precision. As described in "Methods" section, in order to calculate the estimation precision via classical Fisher Information matrix (46) we have to determine the parameter-dependent probability distribution for measurement outcomes p(k|� α) . In our setup the outcomes are labeled by the numbers i k ∈ {0, 1, 2} which denote detector clicks in local modes i k in measurement stations k, therefore the distribution has the form p(i 1 , . . . , i N |� α) . Further, as we show in "Methods" section, the final probability distribution (59) depends only on the total number of clicks in local modes {0, 1, 2} specified by integers z, j, d respectively, therefore p(i 1 , . . . , i N |� α) = p(z, j, d|� α) . Using the final form of the probability distribution (59) we can directly calculate the classical 2 × 2 Fisher Information submatrix (46) corresponding to joint estimation of two of the three phases {α 1 , α 2 } = I , whereas the third phase is set to zero as a reference mode: www.nature.com/scientificreports/ where the multinomial coefficient counts the number of separate detection situations giving rise to the total of z, j, d clicks in modes {0, 1, 2} . In the above formula we took the third mode as a reference mode, however the form of the above Fisher information submatrix does not depend on this choice due to symmetry of the final probability distribution (59) with respect to the parameters α i . Therefore the following analysis holds for estimation of any 2-mode subset of modes of the interferometer. The exact analytical expression for the above defined Fisher information submatrix for arbitrary values of α 's is complicated, however we were able to calculate its inverse and find the optimal scaling of the quantity Tr F (� α) I ,I −1 as a function of the number of photons N. It reads: for the optimal values of estimated phases: Assuming that the estimated phases are equal, α 1 = α 2 = α , the trace of the inverse Fisher information submatrix has the following form: In order to visualise how robust our strategy is for estimation of arbitrary values of the phases we plot the value of Tr F (α 1 , α 2 ) I ,I −1 as a function of the phases α 1 , α 2 for N = 8 , see Fig. 2.
All the above analysis indicates the local character of the Fisher Information-based approach to precision of estimation: the precision strongly depends on the values of the estimated phases. Therefore in realistic applications one needs to obtain some prior knowledge of the phases in order to tune the interferometer in a way that the unknown phases are close to the optimal values for which the error is the lowest.
Notice that even though we do not allow for optimization of final measurements, we still obtain the Heisenberg-like scaling of precision of joint estimation for each of the parameters around its optimal values.
Four-mode case. Optimal state. In analogy to the previous case the N-ququart state which maximises the trace of the QFIM reads: The proof of optimality of (27) follows the same logic as (21). U 4 simultaneosly diagonalizes local Hamiltonians h i (16) with the eigenstates being the standard basis |0�, |1�, |2� , and |3� . Consequently, following the action of the unitary operation U ⊗N 4 , the total collective Hamiltonians H i (20) are diagonal with the eigenstates |0...0�, |1...1�, |2...2� , and |3...3� , which implies that the operator i H 2 i is also diagonal with the same set of eigenstates.
Achievable precision. In analogy to the previous case we have to determine the final probability distribution p(k|� α) . In the current case the outcomes are labeled by the numbers i k ∈ {0, 1, 2, 3} which denote detector clicks in local modes i k in measurement stations k, therefore the distribution has the form p(i 1 , . . . , i N |� α) . As shown in "Methods" section, the final probability distribution (62) depends only on the total number of clicks in local modes {0, 1, 2, 3} specified by integers z, j, d, t respectively, therefore p(i 1 , . . . , i N |� α) = p(z, j, d, t|� α) . Using the final form of the probability distribution (62) we can directly calculate the classical 3 × 3 Fisher Information submatrix corresponding to joint estimation of the three phases {α 1 , α 2 , α 3 } = I , whereas the fourth phase is set to zero as a reference mode: where the multinomial coefficient counts the number of separate detection situations giving rise to the total of z, j, d, t clicks in modes {0, 1, 2, 3} . As in the previous case the above Fisher information submatrix has the same form for any choice of the reference mode, therefore the following analysis of precision of estimation holds for estimating any triple of phases chosen from all the four ones.
Analogously to the previous case we calculated the optimal scaling of the quantity Tr F (� α) I ,I −1 as a function of the number of photons N. It reads: for the optimal values of estimated phases: Assuming that all the estimated phases are equal, α 1 = α 2 = α 3 = α , the trace of the inverse Fisher information submatrix scales with the number of photons N as follows: In this case, we also obtain the Heisenberg-like scaling of precision of joint estimation of any triple of the phases around their optimal values. The same remark on the local character of precision of estimation applies here: one needs to gain a prior knowledge of the unknown phases in order to estimate them around the optimal working point specified by the optimal phases (31).
Symmetric 5-mode multiport with symmetric Hamiltonian. We found that symmetric multiports for d = 5, 6 can be also generated by the powers of generalised X operators in analogy to formulas (3), (8) and (12): As pointed out in seminal paper 29 it is sometimes advisable to analyse the optical multiports from the Hamiltonian perspective. Such an analysis can be necessary in implementations of optical multiports with active optical devices. Let us move to the second quantisation description, in which we define the symmetric Hamiltonian for d-mode optical instrument as:   (33) www.nature.com/scientificreports/ which is a slight generalisation of the definition used in 29 which additionally includes phases. In Ref. 29 , Sect. III it is stated that Hamiltonian of the form (35) cannot generate evolution of a symmetric multiport for d > 4 . However, it can be easily seen that the evolution (33) of a symmetric 5-port is in fact generated by such a Hamiltonian. To see this let us notice that the Hamiltonian of the symmetric multiport (33) reads up to a constant factor: Using the Jordan-Schwinger map: which maps matrix operators on C d into second-quantised operators on d-mode Fock space, we obtain that H 5 has the symmetric representation (35) with the following phases: On the other hand the Hamiltonian: generating the evolution of the 6-mode symmetric multiport (34) does not have symmetric representation (35).
Instead it can be represented as: with the following amplitudes: Single-interferometer optical implementation. Although our scheme in the version presented in the Fig. 1 can be directly implemented using optical interferometry, the main difficulty of such an implementation lies in preparing the multiphoton GHZ state. Despite the progress in realising multiphoton entanglement in recent years (cf. e.g. 31 ) it is still chalenging to prepare such states for higher number of subsystems, which may suggest that the proposed scheme is unfeasible. However our scheme can be transformed into simpler one which is completely feasible within current optical technology by performing second quantisation of the scheme. Note that the final states in our setup are symmetric states of photons which are distinguishable by path degree of freedom: The second-quantised version of the above states, which assumes that a state of N indistinguishable photons is send to a single d-mode interferometer can be expressed as: where the |NOON d � N states are defined by the formula: The detection probabilities for the second-quantised version of the scheme are identical to the ones for the original scheme (59) and (62), only the meaning of the detection events changes: now the numbers z, j, d (and t for the 4-mode case) denote photon counts in modes {0, 1, 2} (and respectively in mode 3) in a single interferometer consisting of the initial-state-correcting multiports U 3 (22) (or U 4 (28)) and the generalised Mach-Zehnder interferometer U 3 (10) (or U 4 (15)). Notice that the equivalence of the original scheme with the second-quantised (39) (42) www.nature.com/scientificreports/ one can be already seen at the level of probabilities (59) and (62), since they do not distinguish in which of the N stations there was a click in a given mode, but depend solely on the total number of clicks in given modes across all the labs. The single-interferometer version of the scheme is experimentally feasible at the current stage since, in contrast to the multiphoton GHZ sources, there are experimental methods to produce N-photon NOON states (44) for higher values of N 23 . Note that a similar idea of using fixed multimode Mach-Zehnder interferometer for multiphase estimation already appeared in several works [24][25][26][27][28] , however in all these works input states with small definite number of photons are considered, therefore they lack discussion about scaling of precision with the size of the initial state.

Discussion
In this work we have investigated the metrological properties of a generalised Mach-Zehnder interferometer for the number of modes equal to 3 and 4, with the emphasis put on the possibility of simultaneous estimation of (d − 1)-element subset of phases placed in arbitrary configuration across the modes. We have shown that estimation of each of the subsets can be performed with Heisenberg-like scaling of precision in an entirely fixed interferometric setup, namely with the same initial state and measurement strategy. To prove the Heisenberg-like scaling of precision we developed an analytical description of the generalised Mach-Zehnder interferometer in terms of Heisenberg-Weyl operators. This approach allows for analytical calculation of the inverse of the classical Fisher information matrix related with each of the subsets of parameters, which, in contrast to methods implicitly involving optimisation over measurement strategies based on Quantum Cramer-Rao or Holevo bounds, provides a factual limit for the efficiency of estimation within assumed concrete measurement setup.
Since our scheme allows for estimation of any (d − 1)-element subset of unknown phases placed arbitrarily across a fixed interferometer (the remaining phase is assumed to be known), the single-interferometer version of the scheme can be well suited for enhancing the performance of 3-dimensional quantum sensing tasks similar to the ones presented in Ref. 24 , in which the estimation is performed using only input states with a small definite number of photons.

Methods
General introduction to multiphase estimation. Standard approach to multiparameter estimation assumes the following estimation scheme 2-4 : an initial multipartite state ρ in undergoes an evolution � α , which depends on a vector of unknown parameters � α = (α 1 , . . . , α d ) . Finally single-particle projective measurements { } k with outcomes labeled as k's are performed, leading to final probability distribution: Having the parameter-dependent probability distribution p(k|� α) one can construct estimators {A i } of the unknown parameters {α i } . In order to estimate the joint accuracy of these estimators one has to introduce joint measure of sensitivity of the distribution p(k|� α) on the parameters {α i } . In classical estimation theory such a measure is provided by the Fisher Information Matrix F (FIM), defined as: If the estimators {A i } are unbiased, namely their mean values equal to {α i } for the entire range of α's, and the F matrix is invertible, the quality of estimation of {α i } based on distribution p(k|� α) is described by the Cramer-Rao bound: where Cov({A i }) is a covariance matrix for estimators, namely Cov({A i }) mn = Cov(A m , A n ) , and ν is the number of repetitions of the experiment. The above description of the efficiency of estimation assumes fixed measurements { } k . In quantum estimation theory one is usually interested in description of efficiency of estimation of α 's from an evolved quantum state ρ out (� α) = � � α (ρ in ) in a way which assumes optimisation over all possible measurements. This idea is encoded in the Quantum Fisher Information Matrix (QFIM), defined in an operatorbased way: where the braces denote anticommutator of operators and the operators L i are defined implicitly by the equation: In the case of pure input state |ψ� and the unitary evolution of the form U = e −iH i α i the QFIM can be expressed in an explicit form: Assuming that the QFIM is invertible the Quantum Cramer-Rao bound holds: www.nature.com/scientificreports/ It is worth to mention that there exists another version of the Quantum Cramer-Rao bound, namely the Holevo bound 3,32,33 , which does not utilise the QFIM. This bound on the precision of estimation is expressed using the notion of a cost matrix C , which is a positive matrix providing weights to the uncertainties related with different parameters and the matrix V defined elementwise by the relation V ij = Tr(X i X j ρ out (� α)) , where the Hermitian matrices X i fulfill the constraint Tr({X i , L j }ρ out (� α)) = 2δ ij and the operators L i are defined in (49). Then the Holevo bound has the following form 3 : where the norm in the RHS is the trace norm. The Holevo bound (52) is tighter than the QCRB bound (51). Despite the fact that the Holevo bound does not utilise the QFIM, it is also ill-defined in the case when the corresponding QFIM-based Cramer-Rao bound is ill-defined due to singularity of the QFIM 33 .
Multiparameter Cramer-Rao bound in the presence of additional parameters. Application of both the bounds (51) and (52) always raises the question whether found precision limitations can be saturated by experimentally accessible measurement schemes. Our main aim is to investigate the process of estimation of multiple phases within a fixed interferometer and a simple fixed measurement scheme consisting of single output mode measurements. Therefore our approach to evaluate the precision of estimation would be based on basic tools in estimation theory, namely on the classical Fisher information matrix techniques. In this way we do not need to care about optimisation of measurements.
However another issue remains to be solved concerning our setup. Namely our task is to estimate multiple phases, which constitute a subset of all the parameters on which the final probability distribution depends (being the phases to be estimated and the reference phase). Such a situation in estimation theory has been discussed extensively in Ref. 30,34 . In our case the chosen subset of phases is referred to as the set of parameters of interest, whereas the reference phase is known and fixed. In general two different cases have to be considered when additional parameters than the estimated ones appear in the setup: (i) the additional parameters are fixed and known, (ii) the additional parameters are fixed but unknown (and are called in this context the nuisance parameters). Let us denote the division of the set of all parameters into parameters of interest and additional parameters by a vector ( α I , α O ) in the case additional parameters are known and as ( α I , α N ) if they are the nuisance parameters. Then the Fisher Information Matrix and its inverse can be expressed in a block form with respect to the fixed partition of parameters ( α I , α O ) and ( α I , α N ) 30 : The Cramer-Rao bound for the parameters of interest on condition that the additional parameters are known has the following form: whereas for the case they are unknown (they are the nuisance parameters) it reads: In simple words the covariance of the estimators of parameters of interest is bounded from below either by the inverse of the submatrix of the FIM (when additional parameters are known) or by the submatrix of the inverese of the FIM (if they are unknown). It is worth to mention, that all of this holds also in the case of a single parameter of interest: if the additional parameters are unknown, still one needs to take into account entire FIM and take its inverse 34 . Probability distributions for the particular outcomes of the experimental setup. Three-mode case. Let i k ∈ {0, 1, 2} denote the measurement outcome at k-th station corresponding to detector click in i k -th local mode. Then the conditional probability distribution for the outcomes conditioned on the values of the phaseshifts reads:  (58) U 3 (α 1 , α 2 , α 3 ).U 3 = 1 √ 3   e iα 2 e iα 3 e iα 1 ω 2 e iα 2 e iα 3 ωe iα 1 e iα 2 ω 2 e iα 3 ωe iα 1   e − 1 3 i(α 1 +α 2 +α 3 ) . (59) p(0, · · · , 0 z , 1, · · · , 1 j , 2, · · · , 2 d |α 1 , α 2 , α 3 ) = 1 3 N+1 3 + 2 cos p(0, · · · , 0 z , 1, · · · , 1 j , 2, · · · , 2 d , 3, · · · , 3 t |α 1 , α 2 , α 3 , α 4 ) = 1 4 N+1 4 + 2 cos (d + j)π + N(α 2 − α 1 ) + 2 cos [(d − z)π + N(α 3 − α 1 )] + 2 cos (j − z)π + N(α 4 − α 1 ) + 2 cos (j + z)π + N(α 2 − α 3 ) + 2 cos [(d + z)π + N(α 2 − α 4 )] + 2 cos (d − j)π + N(α 3 − α 4 ) .