General expressions for the quantum Fisher information matrix with applications to discrete quantum imaging

The quantum Fisher information matrix is a central object in multiparameter quantum estimation theory. It is usually challenging to obtain analytical expressions for it because most calculation methods rely on the diagonalization of the density matrix. In this paper, we derive general expressions for the quantum Fisher information matrix which bypass matrix diagonalization and do not require the expansion of operators on an orthonormal set of states. Additionally, we can tackle density matrices of arbitrary rank. The methods presented here simplify analytical calculations considerably when, for example, the density matrix is more naturally expressed in terms of non-orthogonal states, such as coherent states. Our derivation relies on two matrix inverses which, in principle, can be evaluated analytically even when the density matrix is not diagonalizable in closed form. We demonstrate the power of our approach by deriving novel results in the timely field of discrete quantum imaging: the estimation of positions and intensities of incoherent point sources. We find analytical expressions for the full estimation problem of two point sources with different intensities, and for specific examples with three point sources. We expect that our method will become standard in quantum metrology.

Quantum metrology deals with the estimation of unknown parameters from the measurement outcomes of quantum experiments. Naturally, the goal is to design experiments and data-analysis strategies such that the uncertainty in the estimated parameters is minimized. One of the appeals of quantum metrology is that it promises reduced uncertainties compared to what is possible with comparable classical resources [1][2][3][4][5]. In order to understand and quantify the advantage offered by quantum metrology, we typically calculate the quantum Fisher information (QFI) or, in the general case of simultaneous estimation of multiple parameters, the QFI matrix (QFIM). The inverse of the QFIM yields the quantum Cramér-Rao bound (QCRB), a lower bound on the uncertainty of any unbiased estimator of the parameters. In particular, analytical solutions for the QFIM are desirable as they provide valuable insights into how the estimation error depends on and scales with tunable parameters. Since the QFIM is a function of the density matrix, existing general approaches to compute the QFIM conventionally assume that the density matrix is expressed with respect to an orthogonal basis [6], or even start with a density matrix in its diagonalized form [7]. However, analytical matrix diagonalization is typically limited to low dimensions. On top of that, the density matrix often has a natural and elegant representation with respect to some non-orthogonal basis. In such cases it is not only cumbersome to choose a suitable orthogonal basis and to expand the density matrix in such a basis, but it also complicates the further computation of the QFIM.
In this work, we address these problems and provide a general and efficient method to analytically compute the QFIM. Our approach is based on a new formal solution for the QFIM. Compared to previous methods, our solution relies on a non-orthogonal basis approach [8,9] which allows us to express all matrices with respect to an arbitrary, possibly non-orthogonal basis. Further, our solution does not rely on matrix diagonalization but on matrix inversion. Our results improve over the analysis of Refs. [8,9] by providing an expression for the QFIM which relies on the general solution of the associated Lyapunov equations and, thus, avoids solving the Lyapunov equations for each parameter separately. In comparison with Ref. [10], our expressions for the QFIM are general and do not depend on particular properties of the problem under consideration.
There are many situations in quantum metrology where the natural representation of ρ involves a non-orthogonal set of states, e.g., when coherent states are involved (such as for noisy Schrödinger cat states or entangled coherent states [9]). In this work, we demonstrate the power of our approach by deriving novel analytical results in the field of discrete quantum imaging, a rapidly developing branch of quantum metrology [11]. For the problem of imaging two incoherent point sources of monochromatic light in the paraxial regime, we provide a complete solution for the QFIM when all parameters, i.e., the spatial locations and relative intensity of the sources, are estimated simultaneously. Our method further allows us to derive analytical results for some special arrangements of three point sources, providing new insights into the imaging of multiple sources.
The paper is organized as follows. In Section 1, we provide a precise formulation of the problem and discuss existing approaches to compute the QFIM. In Section 2, we present our formal solution for multiparameter quantum estimation with respect to non-orthogonal bases. In Section 3, we introduce quantum imaging and present our results for the QFIM for two and three point sources. In Section 4, we discuss the computational requirements of our method and close the paper with some concluding remarks.

PRELIMINARIES
Let us consider the general case of simultaneously estimating n parameters, θ = (θ 1 , . . . , θ n ). Further, letθ = (θ 1 , . . . ,θ n ) be an estimator of θ, withθ µ the estimator of θ µ . The uncertainty in the estimatorθ can be characterized by the covariance matrix Cov θ and is lower bounded by the QCRB [12,13] with M the number of repetitions of the experiment (i.e. the statistical ensemble size), and H the quantum Fisher information matrix (QFIM). We use the notation that vectors and matrices are printed in bold while scalars and operators are not. Further, θ j denotes the jth coefficient of vector θ and H µ,ν the (µ, ν)th coefficient of matrix H. The QCRB represents a matrix inequality in the sense that Cov θ − 1 M H −1 has to be a positive semi-definite matrix. Instead of the matrix inequality (1), one often considers a lower bound on the summed variances of eachθ j which is obtained by taking the trace on both sides of Eq. (1).
For a parameter-dependent density operator ρ, the coefficients of the QFIM H are defined as dependence in our notation for the sake of brevity. The SLDs are defined via Lyapunov equations and by inserting Eq. (3) for the derivative in Eq. (2) it becomes apparent that H is a symmetric matrix. Analytical solutions for the QFIM are usually obtained based on the following classical result: for any density operator ρ with nonzero eigenvalues, a formal solution for Eq. (3) is given by [14] Similarly, the QFIM can be written as [2] The most common method to solve the integral in Eq. (4) relies on expanding ρ in its eigendecomposition ρ = ∑ d j=1 λ j λ j ⟩ ⟨λ j , with nonzero eigenvalues λ j and eigenvectors λ j ⟩. Inserting the solution of the integral (4) into Eq. (2) yields the QFIM see Ref. [7] for an overview of analytical results for the QFIM based on the eigendecomposition of ρ, including formulas which hold if ρ does not have full rank. In particular, Eq. (6) holds also for non-full-rank ρ if one restricts the summation to indices for which λ l + λ m > 0 holds. However, the density operator ρ is usually not given in its eigendecomposition, and the diagonalization of ρ is known to be a hard problem; it requires solving the characteristic equation, i.e., a polynomial equation of order d. Abel's impossibility theorem [15] states that algebraic solutions for polynomial equations with arbitrary coefficients are impossible for dimension d ≥ 5. Therefore, finding analytical solutions for the QFIM by expanding ρ in its eigenbasis usually works only for low dimensions or some special cases.
An alternative to the formal solution (4) is to expand the Lyapunov equations (3) in an orthonormal basis and to vectorize the corresponding matrix equation, which is a well-known approach in the mathematical literature [16] but has received little attention in quantum metrology [6,10]. Using vectorization, the Lyapunov equations are transformed into a linear system which can be solved without diagonalizing ρ. In particular, solving a d-dimensional linear system can in principle be done analytically for any finite d and does not suffer from the same limitations as solving polynomial equations of order d, i.e., computing eigenvalues of ρ. If there exists an orthonormal basis, i.e, an orthonormal set of states supporting ρ and its derivatives (by the parameters of interest) such that the density matrix ρ (formed by the coefficients of ρ in that basis) has full rank,Šafránek gives the following formal solution for the QFIM [6] where ∂ µ ρ (∂ ν ρ) is the matrix formed by the coefficients of ∂ µ ρ (∂ ν ρ) expanded in the orthonormal basis, 1 denotes the identity matrix of the same dimension as ρ, A denotes the complex conjugate of A, A † denotes the conjugate transpose of A, and vec(A) is defined as a column vector obtained by concatenating the columns a j of A, In the derivation ofŠafránek's formula, Eq. (7), the inverse matrix in Eq. (7) is part of the formal solution of the Lyapunov equations. Compared to Eq. (6),Šafránek's formula has the advantage that it does not rely on matrix diagonalization but instead uses the inverse of a d 2 × d 2 matrix; matrix inversion is equivalent to solving a linear system and thus does not share the limitations of analytical matrix diagonalization. A drawback of Eq. (7) is that it requires ρ to be invertible, i.e., ρ needs to be given with respect to an orthogonal basis such that its matrix has full rank. If ρ does not have full rank and, thus, is not invertible, one can replace ρ with an invertible matrix ρ s = (1 − s)ρ + s d 1 such that the QFIM is then given by [6] However, Eq. (9) involves additional analytical computations, and if the dimension d of ρ is much larger than its rank, the matrix to be inverted is much larger than necessary. Additionally, ρ is often given with respect to non-orthogonal states which form a basis spanning the support of ρ. Then, the matrix of ρ with respect to this non-orthogonal basis has a compact form and full rank. In such cases, it is nontrivial to find a suitable orthogonal basis which spans only the support of ρ. Instead, one often has to rely on bases for the whole Hilbert space which leads to an inefficient representation of ρ if the rank of ρ is smaller than the dimension of the Hilbert space.
In the following, we address these problems by deriving a general formal solution for the QFIM using a nonorthogonal basis approach [8,9]. Similarly to Eq. (7), our solution does not rely on matrix diagonalization but on matrix inversion, and we will show that it can be seen as a generalization ofŠafránek's formula to non-orthogonal bases.

THE QUANTUM FISHER INFORMATION MATRIX FOR GENERAL BASES
In this section we will present our general expressions for the QFIM and discuss some special cases. The derivation (see Appendix A) uses block-vectorization, a variation of standard vectorization, which allows us to separate a basis into different parts, corresponding to different subspaces. In this way, we define all relevant matrices on their support such that they are invertible.

QFIM for general bases
Our general solution for the QFIM relies only on one assumption, which is that the density operator ρ is given with respect to a d-dimensional basis B = { ψ j ⟩} d j=1 , where B is a set of linearly independent states ψ j ⟩ spanning the support of ρ. Note that we do not assume that the basis B is orthogonal. Then, ρ = ∑ j,k=1 ρ B j,k ψ j ⟩ ⟨ψ k can be represented by a full-rank hermitian matrix ρ B with coefficients ρ B j,k . In the following, we write A B for the matrix representation of A with respect to the basis B.
Matrix equations can be rewritten for matrices which are given with respect to a general, non-orthogonal basis B using the Gramian G B j,k = ⟨ψ j ψ k ⟩, defined with respect to the basis B. Then, for matrices A and B, defined with respect to an orthogonal basis, we have the following replacements for the matrix product and the trace operation [9] AB Clearly, if B is orthonormal, then G B = 1 and we retrieve the standard matrix operations. With this, the Lyapunov equation for the parameter θ µ becomes where B µ denotes a basis which spans the support of ρ and ∂ µ ρ, and G Bµ is the Gramian of B µ . Let us denote the (i, j)th matrix block of a matrix A as A ij . In comparison, A i,j denotes a matrix coefficient. In the following, all matrices are divided in 4 blocks as follows: the A 11 block is always of size B × B where • denotes cardinality and B the given basis spanning the support of ρ. The remaining blocks complete the total matrix. For example, we have where the zeros denote blocks of zeros such that ρ Bµ is of size B µ × B µ . If the derivative of ∂ µ ρ lies within the support of ρ, all but the ρ Bµ 11 = ρ B block vanish. In Appendix A, we use block-vectorization to derive the following general solution for the SLD where we defined the matrices and where mat(⋅) is defined to take a vector of length n 2 and rearrange its coefficients to a matrix of size n × n by inserting the first n coefficients of the vector to the first column of the matrix, the next n coefficients to the second column, and so forth; for example for n = 2 All that remains to be done is to use Eq. (14) to calculate the QFIM which, using our replacement rules (10-11), is given by where B µ,ν is a basis spanning the support of ρ, ∂ µ ρ, and ∂ ν ρ. The matrix L Bµ,ν µ is obtained by padding the matrix L Bµ µ in Eq. (14) with zeros, i.e., (L Bµ,ν µ ) i,j = 0 if one or both of the states ψ i ⟩ and ψ j ⟩ lie outside B µ . Eqs. (14) and (19) represent our general solution for the QFIM. The only nontrivial calculations are to find the inverses of C and D, as given by Eqs. (15) and (16). Similar methods have been used in Ref. [10], although our result is more general because the result of Ref. [10] depends on particular properties of the problem considered in that work.
We proceed by expressing the compatibility conditions [17] with respect to a non-orthogonal basis. This facilitates the evaluation of the compatibility conditions if the QFIM is computed using Eqs. (14) and (19), and in Section 3 we will discuss the compatibility conditions in the context of discrete quantum imaging.

Compatibility conditions
When considering a "multiparameter scenario", i.e., the simultaneous estimation of n > 1 parameters, an interesting question is how it compares to an (overly optimistic) "separate scenario" where one assumes that each parameter can be estimated independently, i.e., that there is an estimation scheme for each parameter, and for each individual estimation scheme one assumes that (i) all other parameters are known and (ii) the same resources are available as for the multiparameter scenario. While such separate scenario is an idealization which usually cannot be implemented (see Refs. [18,19] for a rigorous treatment of quantum parameter estimation with nuisance parameters), it represents a useful benchmark and performs at least as well as the multiparameter scenario. If the so-called compatibility conditions are fulfilled, the multiparameter scenario matches the performance of the separate scenario (while using only the resources of one of the estimation schemes) [17]. The compatibility conditions consist of (i) the commutation condition, (ii) the independence condition, H µ,ν = 0 for all µ ≠ ν, and (iii) the initial-state condition, i.e., that there exists a single initial probe state which is optimal for every estimation scheme in the separate scenario. Note that the QFIM H and Γ can be seen as real and imaginary part of the same quantity since If the commutation condition is fulfilled, there exist optimal measurements [20] such that the QCRB can be saturated (although the estimators of the parameters might not be independent). It is possible to find independent estimators for the parameters if the independence condition is fulfilled, i.e., if the off-diagonal coefficients of the QFIM vanish. If, in addition to the commutation and independence conditions, the initial-state condition is fulfilled, an optimal multiparameter scenario can be constructed which matches the performance of the corresponding separate scenario.
The remainder of this section is devoted to showing how to recover special cases of relevance from our general results.

RetrievingŠafránek's formula as a special case
In order to retrieveŠafránek's formula, as given in Eq. (7), from our general solution, we have to make the additional assumption that the basis B, which spans the support of ρ, is orthogonal, and that the derivatives ∂ µ ρ are supported by B. Then, all matrices can be expressed with respect to B, and G B = 1. This means that C = ρ B , D = 1 B ⊗ρ B +ρ B ⊗1 B , and E vanishes because B µ = B for all µ. The SLD is then given by L B µ = 2 mat D −1 vec ∂ µ ρ B and the QFIM is found to be where we used tr [AB] = vec (A) † vec (B) to get to Eq. (24), and we used that the inverse of an invertible hermitian matrix is hermitian to find Eq. (25) which is indeed identical with Eq. (7). Also, since Eq. (6) can be seen as a special case of Eq. (7) [6], it is a special case of our general solution as well.

Unitary parameterization with commuting generators
Let us introduce a simple unitary parameterization model using operator equations. The parameter-dependent density operator is given by ρ = U (θ)ρ 0 U † (θ), with ρ 0 an parameter-independent initial state, U (θ) = exp −i ∑ n µ=1 K µ θ µ a unitary operator which encodes the parameter dependence into ρ 0 , and we assume that the generators K µ commute.
It is then easy to show that the QFIM can be written as where ρ ′ µ = [K µ , ρ 0 ] replaces the derivative when compared with Eq. (2), and L ′ µ is given via We can solve Eq. (27) 14) which, together with Eq. (26) in matrix form, constitute a solution for the QFIM for the special case of unitary parameterization. Note that this solution does not depend on the parameters θ [since it does not depend on U † (θ)] but only on the generators K µ and the initial state ρ. This is in accordance with known special cases for orthogonal bases [6] and solutions based on matrix diagonalization [7].

Single-parameter estimation
It is worth noticing that the presented method of analytically computing the QFIM also applies to the evaluation of the QFI for single-parameter estimation. The QFI H for estimation of a parameter θ µ is the scalar special case of the QFIM, In particular, avoiding matrix diagonalization and making use of a non-orthogonal basis can be advantageous for computing the QFI as well as the QFIM. In some cases of single-parameter estimation (or even for a few parameters), the analytical computation can be simpler if one solves the linear system [see Eq. (A10) in Appendix A] obtained from the Lyapunov equations explicitly for the particular parameter of interest instead of using the general solution (14) which holds for any parameter. On the other hand, if one wants to estimate multiple parameters, it will often be more convenient to compute the general solution, i.e., to find the inverses of C and D [Eqs. (15)(16)] instead of solving the linear system (A10) for each parameter.

DISCRETE QUANTUM IMAGING
In the previous section we derived a formal solution for the QFIM without diagonalizing the density operator ρ and without expressing ρ with respect to an orthogonal set of states. We continue by showing through the technologically relevant example of discrete quantum imaging that our formal solution can be very useful for finding novel analytical expressions for the QFIM.

Introduction to quantum imaging
Quantum imaging is concerned with using quantum-enhanced detection schemes to image point sources and objects with the highest possible resolution and finding ultimate bounds on the achievable resolution. Potential applications of quantum imaging lie in astronomy, biology, medicine, materials science, and in the semiconductor industry [11].
In 1879, Lord Rayleigh formulated a criterion for the limitations of traditional direct imaging based on classical wave optics: two incoherent point sources cannot be resolved if their separation is significantly smaller than their emission wavelength [21]. Since then, several superresolution techniques such as fluorescent microscopy [22] have been introduced to overcome Rayleigh's criterion.
In order to understand the ultimate fundamental limits of imaging, Tsang et al. developed a quantum metrological framework based on the QCRB which relies on a full quantum description of the imaging process [18,23]. A key result of Tsang et al. has been that there exist detection schemes such as spatial-mode demultiplexing [23] which resolve two incoherent point sources with an error independent of their separation and, thus, completely bypass Rayleigh's principle. This has been corroborated by several proof-of-concept experiments [24][25][26][27][28][29][30][31][32][33]. More detailed studies have shown that for any type of non-adaptive measurement, one-and two-dimensional images of multiple sources in the subdiffraction limit remain unaffected by Rayleigh's criterion only up to the second moment, while for the estimation of higher-order moments a quantum version of Rayleigh's principle reappears [34,35]. Other research has addressed the problem of optimal detection schemes and measurements [36,37] and the implications of practical imperfections such as a misalignment of the detection apparatus [38,39] or noisy detectors [40]; for a review of recent progress see Ref. [11].
Here, we address the problem of deriving analytical expressions for the QFIM for discrete quantum imaging, i.e., for imaging a discrete set of incoherent point sources. In a series of works, analytical expressions for the QFIM have been found for localizing the arbitrary three-dimensional positions of two incoherent point sources of known and possibly unequal brightness [8,[41][42][43][44]. Here we go beyond those results by considering the arguably most general detection problem of two incoherent point sources: the joint estimation of their positions and relative intensity, a total of 7 parameters. To the best of our knowledge, we are the first to deliver a fully analytical solution to this problem. Moreover we note that, when considering more than two sources, results have been of numerical nature so far [10]. Using our previously derived expressions for the QFIM, we are able to obtain analytical solutions for some classes of estimation problems involving three incoherent point sources.
Let us consider the problem of imaging N S weak, incoherent, point-like light sources with intensities {I j } N S j=1 and positions {r s = (x s , y s , z s )} N S s=1 . The light emitted by the sources is collected in the collection plane, see Fig. 1. Conventional methods of collecting light use an aperture in the collection plane, for example a circular aperture. More generally, we can imagine that light is collected at N C points in the collection plane where arbitrary apertures can be retrieved by taking a continuous limit [45]. Let {c j = (v j , w j )} N C j=1 be the collection coordinates of the N C collection points.
We follow the formulation of discrete quantum imaging as given by Lupo et al. [45], for a detailed derivation see Appendix B. We assume that at most one photon is collected per collection window, known as the limit of weak sources. Further, we consider the paraxial regime where the distance z 0 of the sources from the collection plane is much larger than the source and collection coordinates, i.e., x s , y s , z s , v j , w j ≪ z 0 , cf. Fig. 1. For multiple sources, the state of a photon impinging on the collection plane is given by [45] ρ(p, r) = N S s=1 p s ψ(r s )⟩ ⟨ψ(r s ) .
The statistical mixture in Eq. (29) takes into account that the photon must have been emitted by one of the N S sources. The probability p j that the photon has been emitted by source j is given by the relative intensity of the jth source, p j = I j I tot with the total intensity I tot = ∑ N S j=1 I j . Vectors of probabilities and source locations are defined as p = (p 1 , . . . , p N S ) and r = (r 1 , . . . , r N S ), respectively.
The photon states in the collection plane ψ(r s )⟩ emitted by the sources at r s can be parameterized as where is a reference states in the collection plane which contains the information about the location of the N C collection points, and the unitary operator is defined via the operators which generate a commutative group. Here V and W are position operators in the collection plane such that V j⟩ = v j j⟩ and W j⟩ = w j j⟩ for all j = 1, . . . , N S . After a photon has been collected at the collection plane, it is processed coherently in a general interferometer (for details see Ref. [45]) before it is measured using photodetection, cf. Fig. 1. Since the interferometer corresponds to a unitary transformation which is assumed to be independent of the source parameters, it does not change the QFIM [7] and we can proceed with calculating the QFIM from Eq. (29). Potential parameters of interest are the relative intensities of the sources, i.e., the probabilities p s , and the positions of the sources r s , which are parameters of the unitary transformation in Eq. (32).

Analytical results for discrete quantum imaging
Calculating the QFIM with our formal solution, Eqs. (14)(15)(16)(17)(18)(19), is a matter of tedious but straightforward algebra. We choose the basis B = { ψ(r s )⟩} N S s=1 , with basis states as given in Eq. (30), such that we can express the density operator (29) by a diagonal matrix ρ B (p, r) where ρ B s,s (p, r) = p s . If we want to estimate one of the probabilities p s , for example θ µ = p 1 , the corresponding basis B µ is identical with B because ∂ µ ρ is supported by B. On the other hand, if we want to estimate one of the position parameters, e.g. θ ν = x 1 , we extend B by ∂ ν ψ(r s )⟩ to obtain B ν . Note that ∂ ν ψ(r s )⟩ is linearly independent from all vectors in B only for almost all values of the parameters, for more details see Appendix C.

Compatibility conditions
In Section 2.2, we formulated the compatibility conditions with respect to non-orthogonal bases. In the case of discrete quantum imaging, the emission properties of the sources cannot be changed and therefore the initial state is fixed for all estimation schemes. Then, the compatibility conditions reduce to (i) the commutation condition, Γ µ,ν = 0 for all µ ≠ ν, where Γ is given in Eq. (20), and (ii) the independence condition, H µ,ν = 0 for all µ ≠ ν. In the following, we will refer to the commutation and independence conditions in order to interpret the results for two and three point sources.

Two sources
Here we consider the general problem of imaging two sources, i.e., the estimation of their positions and their relative intensity. The collected photon state is then given by It is convenient to reparameterize ρ using centroid and relative coordinates for the two sources: we define the relative coordinates as θ 1 = (x 1 − x 2 ) 2, θ 2 = (y 1 − y 2 ) 2, and θ 3 = (z 1 − z 2 ) 2 and the centroid coordinates as θ 4 = (x 1 + x 2 ) 2, θ 5 = (y 1 + y 2 ) 2, and θ 6 = (z 1 + z 2 ) 2. This means we want to estimate 7 parameters: the centroid and relative coordinates of the sources and θ 7 = p 1 . Note that p 2 is fixed due to normalization, p 1 + p 2 = 1, and p 1 = I 1 I tot indeed corresponds to the relative intensity. For example, if we know I tot , e.g., from an independent intensity measurement, we directly obtain I 1 from estimating p 1 . The matrices which must be inverted are the matrix C, see Eq. (15), which is of size 2 × 2, and the matrix D, see Eq. (16), which is of size 4 × 4. We find where g = (G x , G y , G z ) ⊺ and δ = (θ 1 , θ 2 , θ 3 ) ⊺ summarize the generators and relative source coordinates in vector notation, ⟨⋅⟩ = ⟨ψ(0) ⋅ ψ(0)⟩ denotes the expectation value with respect to the reference state ψ(0)⟩, the covariance of the generators is defined as Cov(g) j,k = ⟨g j g k ⟩ − ⟨g j ⟩ ⟨g k ⟩, and the variance as Var(A) = ⟨A 2 ⟩ − ⟨A⟩ 2 .
The first two columns (rows) in Eq. (35) each consist of three columns (rows), summarized using matrix and vector notation, such that the QFIM is actually a 7-dimensional matrix. The jth row and column correspond to the parameter θ j . Since we are in the paraxial regime, Eq. (35) includes only the lowest-order non-zero terms [46] in the rescaled source coordinates δ j z 0 .
By fixing certain parameters, we recover known results: if we set θ 7 = p 1 = 1 2, i.e., we consider sources of equal intensity, the upper-left 6×6 block corresponding to the centroid and relative coordinates becomes block diagonal which reproduces the well-known result that, according to the independence condition, centroid and relative coordinates can be estimated independently from each other. On the other hand, the estimation errors among different centroid coordinates are not independent (the same holds for different relative coordinates) [8,43,45]. Further, from the zero blocks in Eq. (35) it follows that the relative coordinates can be estimated independently of the relative intensity of the sources. This is in agreement with the results for two sources constrained to one dimension [41]. Finally, Eq. (35) generalizes the result for the estimation of centroid and relative coordinates of two sources of known unequal brightness [44]. We also reproduce the result of Ref. [44] that for unequal source brightness the estimation of centroid and relative coordinates is no longer independent; in particular, since the upper-left 6×6 block is proportional to Cov(g), there exists no collection pattern or aperture which makes the off-diagonal vanish while the diagonal blocks are nonzero.
We note that the QFIM is independent of the centroid coordinates, i.e., the information content of the collected light does not change by jointly moving the sources or the collection instrument as a whole. Further, the covariance Cov(g) which characterizes the QFIM for relative and centroid coordinates shows that distributions of collection points with a large variance along a x-axis (y-axis) are better suited to estimate the corresponding source coordinates along the x (y) direction, while the variance along the x-and y-axis contributes equally for the estimation of source coordinates along the z direction. Similarly, a better precision in estimating the relative intensity can be achieved if the collection points exhibit a larger variance along the x-axis (y-axis) if the sources have a larger separation along the x-axis (y-axis) compared to the y-axis (x-axis), because the variances in H 7,7 are scaled with the squared sources separation along the corresponding direction; if the sources are only separated along the z-axis, the variances of the source coordinates along the x-and y-axis contribute equally. For example, for a circular aperture of diameter D, the variance along the x-or y-axis is proportional to D 2 , i.e., the diagonal elements of the QFIM scale as D 2 .
To answer the question of the existence of optimal measurements, we compute the matrix Γ. Using the same block-partitioning as for the QFIM in Eq. (35), we find where Similarly to the lowest-order approximation of the QFIM, the coefficients of Γ are given at the lowest non-vanishing order with respect to the rescaled source coordinates δ j z 0 . We find that the diagonal blocks of Γ vanish. Thus, there exist optimal measurements for the estimation of centroid or relative coordinates (but not for the combination of both) and the QCRB can be saturated although the estimators for different centroid or relative coordinates are generally correlated since the QFIM is not diagonal. Further, we note that Γ 23 and Γ 13 are second-order with respect to δ z 0 while Γ 12 is first order. This means that in the the paraxial regime Γ 23 and Γ 13 are approximately zero and there exist optimal measurements for the relative intensity and the relative coordinates or the relative intensity and the centroid coordinates. We can infer from the zeros in Eq. (35) that one centroid coordinate and the relative intensity can be jointly estimated optimally, i.e., with independent estimators and measurements which saturate the QCRB. Important examples of apertures are the circular aperture, e.g., studied in Refs. [42][43][44], and the Gaussian-beam assumption [8,10] which corresponds to a Gaussian distribution of collection points. Since we found a general solution for any aperture or distribution of collection points, we can make more general statements. For example, for any distribution of collection points which is radially symmetric around the optical axis, i.e., such that for any collection point at position (v, w) there exists another collection point at (−v, −w), we find that Eq. (36) becomes Γ = 0, i.e., the QCRB can in principle be saturated. Clearly, this includes the aforementioned cases of circular aperture or a Gaussian distribution of collection points.

Three sources
We further consider two examples with three equidistant sources. We assume that the three sources are aligned along the x-axis, i.e., x 1 = c x − δ x , x 2 = c x , and x 3 = c x + δ x , with known centroid coordinate c x . First, let us consider the problem of estimating the distance δ x where we assume that the relative intensities are known. In lowest order with respect to δ x z 0 , the QFI is then obtained as where p 1 + p 2 + p 3 = 1. In accordance with our results for two sources, we find that a non-zero variance in the collection points in the x direction is crucial for the estimation of the relative intensity if the sources are aligned along the x-axis. For comparison, in the distance estimation of two sources of known (unequal) relative intensity along the x-axis, i.e., x 1 = c x − δ x 2, and x 2 = c x + δ x 2, with known centroid coordinate c x , one finds H = Var (g x ) for the QFI. At first sight, this is surprising because it does not seem to match Eq. (40) which yields H = 4p 3 Var (g x ) if we set p 1 to zero. However, the two estimation problems are different: in case (i), the three-sources problem with p 1 = 0 corresponds to a situation where the position of source 2 is known, while the position of source 3 has to be estimated (if δ x is varied, only source 3 moves). Accordingly, the QFI equals the QFI for the position estimation of a single source rescaled with its relative intensity p 3 ; source two effectively acts as a source of noise which reduces the relevant signal. On the other hand, in case (ii), the problem of the two sources placed symmetrically around a known centroid coordinate requires a joint estimation of the sources in order to estimate their distance (if δ x is varied, both source move). Both cases can be summarized in one formula by parameterizing the two sources as x 1 = c x + δ x (q − 1 2), and x 2 = c x + δ x (q + 1 2) where q ∈ R is an additional scaling parameter such that we obtain case (i) for q = 1 2 and case (ii) for q = 0. Then, we find for the QFI in lowest order with respect to δ x z 0 , Note that H grows quadratically with q, a typical effect when parameters are rescaled because scaling factors increase the sensitivity with respect to changes of the parameters. Finally, it is worth noticing that we can find analog results for the QFI in Eqs. (40) and (41) if the sources are aligned along the y-axis (z-axis) where Var (g x ) is replaced with Var (g y ) (Var [g z ]). For our next example, we consider again three equidistant sources of unequal brightness along the x-axis, x 1 = c x −δ x , x 2 = c x , and x 3 = c x + δ x , with known centroid coordinate. This time, we assume that their distance is known and we estimate p 1 and p 2 while p 3 is determined by p 3 = 1 − p 1 − p 2 . In lowest order with respect to δ x z 0 , the QFIM is found to be The QFIM is proportional to δ 2 x Var (g x ), which means that a nonzero source separation as well as a nonzero separation of collection points along the x-axis is necessary to estimate the relative intensities. The fact that source 2 has a special place as the middle source breaks the symmetry between source 1 and 2. This is reflected in the asymmetry between the diagonal coefficients of the QFIM which correspond to the estimation of p 1 and p 2 . For example, if we evaluate the QFIM for equal source intensities, p 1 = p 2 = 1 3 (i.e., p 3 = 1 3), the asymmetry persists: Further, note that there are statistical correlations between the estimators of p 1 and p 2 because the off-diagonal coefficients in Eq. (42) are nonzero for sources of finite brightness. We find that Γ is zero up to second order with respect to the rescaled source coordinates (while there are non-zero higher-order terms). This means that, in the paraxial regime, there exist optimal measurements for the joint estimation of the relative intensities.

DISCUSSION AND CONCLUSION
In this paper we obtained fully general solutions for the quantum Fisher information (QFI) and the QFI matrix (QFIM) [see Eq. (14) and (19)], which are figure of merits of core importance in quantum metrology. Based on these solutions, we provided a general method to analytically calculate the QFI and QFIM. Our solutions generalize previous results [9,10] and we show in particular thatŠafránek's QFIM expression [6] is a special case of our solution for orthogonal bases.
Compared to the conventional method of calculating the QFIM which relies on matrix diagonalization, our method does not share the rather strict limitations of analytical matrix diagonalization. Instead, our method relies on the computation of matrix inverses which amounts to solving linear systems. Then, finding analytical solutions is usually only limited by our ability to handle long algebraic expressions, see Appendix E for more details.
WhileŠafránek's QFIM expression [6] shares the aforementioned advantages from avoiding matrix diagonalization, our method has the additional advantage that it does not rely on expanding operators in an orthogonal basis. This can be very convenient when the density operator is given in a non-orthogonal basis. In such cases good choices for an orthogonal basis for analytical computations are often hard to find and lead to an inefficient representation which hampers analytical computation of the QFIM. In particular, switching to an orthogonal basis often leads to larger matrices which do not have full rank. Our method accepts any density matrix with respect to any given, possibly non-orthogonal basis where the only requirement is that the density matrix has full rank in this basis. Then, our method keeps the dimension of matrices as low as possible which facilitates analytical computations significantly. While our method avoids diagonalization, it requires the computation of two inverse matrices of dimension d and d 2 where d is the rank of the density matrix. Since matrix inversion is simpler than diagonalization in many respects, our method should be applied whenever it is impossible to diagonalize the density matrix. We remark once more that our method and the above considerations apply to the analytical evaluation of the QFIM. A detailed discussion about numerical computation is provided in Appendix D.
We demonstrated the usefulness of our method by deriving new analytical solutions for discrete quantum imaging which generalize existing results for two points sources and provide insights about selected problems with three point sources. Thanks to its generality and its advantages over previous methods when the density operator is given in a non-orthogonal basis, we expect that our method of computing analytical solutions for the QFI and the QFIM will find widespread application and can become a standard tool in quantum metrology.

Appendix A: Derivation of the QFIM for general bases
In this appendix, we will derive general expressions for the QFIM. The derivation involves four steps: (i) we write the Lyapunov equation for matrices defined with respect to an arbitrary (non-orthogonal) basis, (ii) we rewrite the Lyapunov equations using block-vectorization, (iii) we derive a formal solution for the SLD, and (iv) we insert the solution in the expression for the QFIM.
We would like to refer to the main text [from the beginning of Section 2.1 up to but not including Eq. (14)] for Lyapunov equations expressed in a general (non-orthogonal) basis and some notation which will be used in the following. We continue with deriving a formal solution of the Lyapunov equations. We start by introducing more notation. For a given block decomposition of matrix A, vecb(A) denotes a block-wise vectorization: Note that vecb depends on the particular partitioning of matrix A. We use the notation that A ij denotes the (i, j)th block of A in contrast to A i,j which denotes the (i, j)th coefficients of matrix A. In particular, we apply a convention for dividing matrices in blocks which has been introduced in the main text (Section 2.1).
Next we rewrite the SLD equations using an important identity for block-vectorization [47, p.49], see also Ref. [48,Lemma 4]: where C ⊺ denotes the transpose of matrix C, and ⊙ denotes the Tracy-Singh product, a generalization of the tensor product for block-partitioned matrices. For two partitioned matrices A and B, it is defined as All matrices in Eq. (A4) are give with respect to the basis B µ , however, for better readability we drop the superscripts B µ in Eq. (A4) in the following until we reach Eq. (A14). With the identity (A2), and using that ρ and G are hermitian, we find where A denotes the complex conjugate of A. By decomposing ρ and G in blocks as stated below Eq. (13), Eq. (A9) becomes (still, all matrices with respect to B µ ) where it followed from the chain rule of differentiation that vec (∂ µ ρ) 22 is zero. We continue by solving Eq. (A10) for the SLD L µ . Note that the big matrix of block matrices in Eq. (A10) is already upper triangular, and the bottom row contains only zero blocks indicating that L µ is underdetermined. We pick a solution for L µ by setting vec (L µ ) 22 = 0. This yields the following solution for the SLD: vec (L µ ) 11 = 2 D −1 vec (∂ µ ρ) 11 vec (L µ ) 21 = 2 vec (∂ µ ρ) 21 vec (L µ ) 12 = 2 vec C −1 (∂ µ ρ) 12 , vec (L µ ) 22 = 0, where we used identity (A2) for vectorization, and the matrices C, D, and E are defined in the main text, see Eqs. (15)(16)(17).
Eqs. (A11-A14) constitute a general solution to the Lyapunov equations (12) with respect to the basis B µ . In particular, using identity (A2) and the mat(⋅) operation [defined in the main text, cf. Eq. (18)], we can write the solution for L µ in matrix form where we keep in mind that all matrices starting from Eq. (A4) up to the following equation are given with respect to B µ , though, we dropped the superscript B µ for better readability. The solution of L µ in matrix form is given in Eq. (14) in the main text where we also explain how to obtain the QFIM [Eq. (19)].
An alternative representation of the QFIM to Eq. (19) can be obtained by writing the QFIM in vectorized form, where we used the cyclic property of the trace, and tr [AB] = vecb (A) † vecb (B). The block-vectorized solution Eqs. (A11-A14) can be extended to the basis B µ,ν by padding with zeros and can then be plugged into Eq. (A17).

Appendix B: A quantum formulation of discrete quantum imaging
In this appendix, we will reproduce the formulation of discrete quantum imaging as given by Lupo et al. [45]. We assume that at most one photon is collected per collection window, known as the limit of weak sources. Then, a general single-mode photon state emitted by a source at r s and collected in the collection plane can be expressed as ψ(r s )⟩ = N C j=1 γ (c j , r s ) j⟩ . (B1) The orthonormal states j⟩ contain the information about the location of the collection points where j labels the collection points, and γ (c j , r s ) are general complex amplitudes which must fulfill normalization, ∑ j γ (c j , r s ) 2 = 1. The relative phases ϕ (c j , r s ) = arg γ (c j , r s ) = kl (c j , r s ) (B2) depend on the wave number k and the path length l (c j , r s ) from the source at r s to the collection point c j , where, in the second line (B4), the distance z 0 of the sources from the collection plane, cf. Fig. 1, has been factored out and the primed variables equal the unprimed ones scaled with a factor 1 z 0 . We consider the paraxial regime where the distance of the sources from the collection plane is much larger than the source and collection coordinates, i.e., x s , y s , z s , v j , w j ≪ z 0 ; this means we can expand ϕ (c j , r s ) for all primed variables around zero. In order to get non-trivial terms for each of the source coordinates, we compute a multivariate Taylor expansion up to the third order: respect to the parameters. Only at the critical points, we find removable singularities which can be eliminated by separately calculating the QFIM at the critical points using a reduced set of basis vectors. When we try to reduce a higher order expression for the QFIM in Section 3 to first or second order, the problem of critical points becomes relevant if we try to expand around such points. In order to calculate a Taylor expansion around such removable singularities, we take the limit towards the singularity. The limit can be calculated making repeatedly use of L'Hôpital's rule. The only drawback of this method is that iterating L'Hôpital's rule requires to take higher orders of derivatives which leads to even longer expressions. In some cases, this turned out to be a computational bottleneck.