Tensor-network approach for quantum metrology in many-body quantum systems.

Identification of the optimal quantum metrological protocols in realistic many particle quantum models is in general a challenge that cannot be efficiently addressed by the state-of-the-art numerical and analytical methods. Here we provide a comprehensive framework exploiting matrix product operators (MPO) type tensor networks for quantum metrological problems. The maximal achievable estimation precision as well as the optimal probe states in previously inaccessible regimes can be identified including models with short-range noise correlations. Moreover, the application of infinite MPO (iMPO) techniques allows for a direct and efficient determination of the asymptotic precision in the limit of infinite particle numbers. We illustrate the potential of our framework in terms of an atomic clock stabilization (temporal noise correlation) example as well as magnetic field sensing (spatial noise correlations). As a byproduct, the developed methods may be used to calculate the fidelity susceptibility-a parameter widely used to study phase transitions.

Introduction

Quantum metrology[1-6] is plagued by the same computational difficulties afflicting all quantum information processing technologies, namely, the exponential growth of the dimension of many particle Hilbert space[7,8]. Apart from idealized noiseless models[1] as well as models operating within the fully symmetric subspace[9] (where the Hilbert space dimension grows linearly with the number of particles) only small-scale problems are feasible via direct numerical study, and even a slight increase in the number of elementary objects makes such an approach intractable. Interestingly, in case of uncorrelated noise models, easily computable fundamental precision bounds are available[10-16] and hence a deep-physical insight may be obtained even if direct numerical optimization is infeasible. However, in cases when one deals with metrological models involving correlated noise, or whenever states outside the fully symmetric subspace are involved, there are no efficient methods that can be applied.

Correlated noise appears naturally in a number of highly relevant metrological problems. Temporal noise correlations are present in the atomic clock stabilization problem[17], making identification of the optimal quantum clock stabilization strategies a highly non-trivial task[18-20]. An even more challenging case involves models where time-correlated noise cannot be effectively described via some classical stochastic process[21] and as such manifests non-Markovian features of quantum dynamics[22], as e.g., in NV-center sensing models[23]. A second natural setting exhibiting non-trivial noise correlations is that of many-body systems such as, e.g., spin chains. Here, typically, spatially correlated noise emerges, which is of crucial relevance for any models where the effective signal is obtained from spatially distributed probes[24-26].

Temporal noise correlations usually decay rapidly. Similarly, in the spatially correlated case one expects on dimensional and energetic grounds that noise will be short-range correlated. The most successful approach for classical simulating short-range correlated many-body systems is via the variational tensor-network state (TNS) ansatz[27]. Among many ansatz classes, that led to unparalleled insights into the physics of quantum many-body systems, the most relevant for the present work is the matrix product operator (MPO) ansatz for density operators[28] and also its infinite particle limit known as infinite MPO (iMPO)[29].

In this paper, building upon experience obtained from uncorrelated noise metrological studies[30], where the optimal input states where shown to be efficiently described as matrix product states (MPS), we develop a comprehensive tensor-network-based framework allowing to (i) calculate relevant metrological quantities (such as the Quantum Fisher Information (QFI) or a Bayesian-type cost), (ii) optimize input probe states and as a result (iii) identify the optimal metrological protocol. All this is accomplished while remaining fully within the tensor-network formalism and thus avoiding the curse of dimensionality along the way.

Results

Efficient identification of the optimal quantum metrological protocol

A paradigmatic task in quantum metrology is schematically depicted in Fig. 1. The central goal is to find the best input probe state , the best measurement and estimator so as to minimize the average uncertainty , where the expectation is over all measurement results x.

Fig. 1

A scheme of a standard quantum metrological task.

A probe state is subject to a parameter-dependent quantum evolution, mathematically represented by a parameter-dependent quantum channel . A POVM type measurement is then carried out yielding a conditional probability distribution , where . Given the conditional probability distribution the objective is to estimate the value of the unknown parameter . To this end one employs an estimator function to produce a given estimate for given the measurement outcome x.

This task is somewhat facilitated by exploiting a fundamental result in quantum metrology, namely, the Cramér-Rao inequality[31,32], which lower-bounds the average uncertainty of the best possible estimator where is the QFI of the output state.

In this paper, we will use the following formula for the QFI[9,33]:where is the derivative of with respect to . This form is equivalent to the standard definition of the QFI[31,32], as can be seen by solving the above maximization problem with respect to L—this is formally a quadratic function in matrix L and the resulting extremum condition yields the standard linear equation for the symmetric logarithmic derivative (SLD) L, . When the solution for L is plugged into the above formula it yields in agreement with the standard definition of the QFI.

The above QFI formula has an advantage over the standard one, when one wants to additionally perform optimization of the QFI over the input states in order to find the optimal quantum metrological protocol. This problem can be written as a double maximization problem:With fixed is effectively a function of , so in what follows for the sake of simplicity of notation we will write , and instead of , . The Figure of Merit (FoM) is linear in and quadratic in L. This formulation leads to an extremely efficient iterative numerical procedure for determining the optimal input probe state: start with some random (or an educated guess for an) input state, determine the corresponding optimal L by performing the relevant optimization. Then, for the L just found, reverse the procedure and look for the optimal input state. This procedure converges very quickly and yields the optimal input probe state, as well as the corresponding QFI. This approach was first proposed in refs. [9,34] in the Bayesian estimation context, and then applied to the QFI FoM in ref. [33] (recently the method has been rediscovered in a slightly modified incarnation in ref. [35] and proved useful in studying metrological properties of states with a positive partial transpose). The above considerations are valid whenever the quantity to be optimized is given in the form of Eq. (2): need not necessarily be the derivative of the state with respect to the estimated parameter. As such, this procedure is applicable in the Bayesian approach, and also in case of less trivial FoMs (see the atomic clock stabilization example).

We consider metrological models involving N distinguishable d-dimensional probes, so that the total Hilbert space is . We assume that the parameter is unitarily encoded in the output state according to a product of unitaries given by the exponential of local generators (or Hamiltonians):where is the generator acting on the nth particle. Most importantly for this paper, the noise, represented above by the operator , is not assumed to be local, but admits the possibility of noise correlations. In order to be able to make use of MPO approach efficiently we assume that the noise correlations are short-range, i.e., may be effectively approximated as a product of single, two-, three-, etc. particle maps up to some cutoff point after which to is assumed that noise correlations do not extend beyond r neighboring particles—see Supplementary Note 1 for a more detailed and more physical discussion of this approximation.

Expressing the optimization problem using the MPO formalism

The optimization scheme described above is completely general and, provided the systems considered are small enough, it can be implemented numerically using the standard quantum state representation. For larger systems, however, this will not be possible in general, and one way to circumvent this difficulty is to implement the above algorithm using the MPO formalism.

In order to do so, we first write the initial N particle state as an MPO (see Supplementary Note 4 for a short introduction to tensor networks):where denotes the standard product basis for N-particle states, with being the physical indices, while are matrices (for which entries are identified by virtual indices ) and is referred to as the bond-dimension of the MPO. The above state is depicted diagrammatically belowwhere the contracted lines represent virtual indices and uncontracted lines physical indices. We may vectorize and obtain its MPS representation by formally bending the vertical legs upward (), which results inwhere a thick vertical line ranges over a doubled physical index .

Our next step is to find the tensor-network representation of . For definiteness, we focus on the situation when the noise operator can be described by a subsequent action of singe-particle and two-particle terms —which in the following are denoted by Y and X, respectively. Physically this is the case when single and two-particle evolution terms commute and no particle is distinguished—generalizations to more complex situations are tedious but straightforward. In the tensor-network representation the action of the operator takes the following form:where we place X operators in a skewed orientation in order to maintain a manifestly translation invariant model.

The operator X is defined with respect to a product basis for the doubled legs and it acts on two neighboring subsystems and , i.e., it is the tensor . We may now perform the singular value decomposition (SVD) of tensor X, which has the following graphical representation:Here, we regard the two vertical legs of X on the left as a doubled leg acting on a single virtual system and the two vertical legs on the right as acting on a second virtual system: in this way one can think of X as a simple matrix acting on a virtual system and we can then apply the SVD.

We apply the SVD to each X operator and absorb into the U tensor from the right (respectively, into the tensor from the left). Let (the upper index indicates two-particle nature of the noise) be the number of non-zero (or more practically non-negligible) singular values of matrix S. Introduction of , finalizes the decomposition of the X-layer of the tensor-network:

The final step to obtain an MPO representation for is to combine the tensors , T, W, Y, and a tensor —which represents the unitary phase encoding process—into a single new MPS tensor :The doubled horizontal legs in the scheme above can be now combined into thicker horizontal legs to obtain the MPO representation of (note that we have split back the vertical legs so we have a proper density matrix, and not its vectorized form):

As a result we end up with an MPO with the bond-dimension . The generalization of this derivation beyond the case of nearest-neighbor correlations will lead to an MPO representation of the density matrix with bond-dimension . Here, represents the contribution to the effective bond-dimension of the output state resulting from the action of the correlated noise, where is the number of non-zero singular values that will appear when considering the s-particle noise term.

Operator can be efficiently written as an MPO as well thanks to the fact that it is as a commutator of with H, where H is a sum of local HamiltoniansIn what follows we assume that the basis associated with the physical indices j is chosen to be the eigenbasis of the local Hamiltonian h, , where are the corresponding eigenvalues. The MPO representation of can be written at the cost of doubling the bond-dimension ():where is equal toThe matrices that appear in the above construction are responsible for the increase of the bond-dimension but guarantee that the effect of trace in Eq. (7) is equivalent to that resulting from the sum of MPO acted upon consecutively by the commutator of the local Hamiltonians corresponding to different particles.

Finally, the optimization algorithm needs to be implemented using the MPO structures introduced above. Here, we just sketch how the two-step iterative process of identifying the optimal L and is implemented within the MPO framework in order to arrive at the solution of Eq. (3)—see the Methods section for a detailed description. In the first iteration step we seek the maximum of the over a Hermitian operator L with a fixed initial density matrix (QFI is a quadratic function of L). We use the MPO representation of L and, instead of trying to optimize the whole operator L at once, we iteratively optimize each tensor site by site in a loop up to the desired convergence of the FoM. During the one-site optimization step we express our FoM as a quadratic function of the vectorized tensor, which we want to optimize. This gives us a condition for the extremum in a form of a linear equation for the components of this tensor, which we solve using the Moore–Penrose pseudo-inverse. The second iteration step, where the optimization over the initial state is performed, is similar in concept but varies in details. One can show that the optimal input state is pure and then using the Heisenberg picture of the evolution process express the optimization problem as a minimization of a Rayleight quotient involving a Hermitian operator and the respective input state. When the Hermitian operator is interpreted as a Hamiltonian then such a problem reduces to the well known variational method in quantum mechanics, which when the MPS formalism is incorporated can be efficiently solved using the DMRG algorithm[27].

Note that in the above process we have an outer iteration loop related with the iteration process where we subsequently optimize L and in the FoM quantity and at the same time we have an inner iteration loop in each of these steps where we iterate over sites of the relevant tensor-network in order to find the optimal solution within a given class of MPO. Furthermore, while running the algorithm, one has to choose the bond-dimensions for the input state, as well as for the L operator. We start with low bond-dimensions (typically 1 or 2) and then increase them subsequently only when the increase yields a noticable change () in the optimized FoM. In this way we achieve the desired accuracy while keeping the bond-dimensions as low as possible thus guaranteeing the efficiency of the algorithm. This procedure may be implemented both in the finite particle number regime as well in the asymptotic limit of infinite number of particles by utilizing the iMPO—see the Methods section. This latter approach provides us with a unique insight into the asymptotic efficiency of the quantum enhanced metrological protocols.

Below we present three applications of our framework chosen in a way so as to highlight the possibility of applying the framework to a variety of qualitatively different physical problems and therefore demonstrate versatility of the approach.

Magnetic field sensing with locally correlated noise

Consider N particles each with spin-, which interact for a fixed time t with an external magnetic field B (assumed to be in the z direction) whose strength fluctuates. The fluctuations induce an effective dephasing process on the particles and apart from the uncorrelated dephasing contribution[10,36] we will also take into account correlations between field fluctuations at the nearest-neighbor particle sites. The magnetic field at site n is given as , where B is the mean value of the field to be estimated. Here, we assume that fluctuations are Gaussian and have no relevant temporal correlations. The corresponding variance as well as the nearest-neighbor correlation functions of the fluctuating field read: , , where represents the strength of correlations and may be both positive and negative (anti-correlation)—for simplicity, we assume periodic boundary conditions, so in fact also the particles N and 1 are correlated. This model corresponds to the choice , (where g is the gyromagnetic ratio for the particle), whereas using the standard cumulant expansion techniques the field fluctuations lead to:where are column vectors (), and C is the correlation matrixwhere , represent, respectively, local and correlated dephasing strength—see Supplementary Note 2 for extended discussion.

In order to write the evolution manifestly in the MPO formalism, we will replace , which forms a basis for the vectorized input density matrix . The action of on is identical to the action of the operator on , i.e., , wherewithandNote that the and mutually commute with each other, so thatDenoting and we finally obtainwhich is the form of evolution the same as discussed in Result section guaranteeing efficient MPO description.

After evolution through quantum channel , the phase is imprinted in our state through unitary evolution according to Eq. (4) with local Hamiltonians —in the tensor-network picture this is represented by the action of , where . Written in the basis , reads:

Figure 2(a) presents a comparison of results of the QFI optimization procedure for exemplary noise parameters obtained using the finite number of particles N MPO approach and the asymptotic value of the QFI per particle obtained using the iMPO approach. The results obtained via the two approaches are in very good agreement. This is a numerical confirmation that indeed the iMPO approach, which, as described in the Methods section is much more conceptually involved, yields correct results. In Fig. 2b we present the contour plot depicting the asymptotic value of the QFI per particle as a function of noise parameters, where bottom left inset also demonstrates that state-of-the-art methods developed with uncorrelated noise models in mind yield bounds that are far from satisfactory in case of correlated noise models. The main qualitative feature that clearly emerges is the decrease of the optimal QFI with the increase of correlated noise part parameter . At the same time going into the anti-correlation regime (negative ) allows for a significant increase in the achievable QFI—this is to be expected based on intuition obtained from quantum error-correction-based metrology where purely anti-correlated noise may be even completely removed recovering the Heisenberg scaling[26].

Fig. 2

Quantum Fisher Information for the magnetic field sensing model.

a Comparison of the QFI per particle for a magnetic field sensing problem in presence of locally correlated dephasing as a function of the number of spins in a chain N (for dephasing noise model with local noise parameter and correlation parameter ) calculated using the finite MPO approach (black dots) with asymptotic value obtained using the iMPO approach (black solid line). Gray crosses indicate results obtained via the standard full-Hilbert space description. Gray lines show state-of-the-art bounds on obtained by decomposing the dynamics into effectively independent channels of increasing complexity (dotted, dash-dotted, dashed)—see Supplementary Note 2 for details. For comparison, the solid gray line corresponds to the bound obtained when all correlations are neglected and only local dephasing noise is taken into account. b Asymptotic value (obtained using iMPOs) of QFI per particle for dephasing type noise in function of local and between nearest neighbors noise parameters. Black equipotential lines are in logarithmic scale. Left inset shows a slice of the main plot along and presents the results obtained using the iMPO approach (black dots) compared with the exact asymptotic result for a weakly squeezed state strategy (light gray line), and the gray lines show state-of-the-art bounds on QFI/N (dotted, dash-dotted, dashed). Right inset shows a slice of the main plot along and presents result obtained using iMPO approach (black dots) compared with the known exact result for strictly local noise with (light gray line).

In case of purely local dephasing it is known that in the limit of large number of particles the fundamental bound, , can be saturated by protocols involving weakly spin-squeezed states[10,37]. We have performed analogous analysis in case of correlated noise, see Supplementary Note 2, and found out that asymptotically the strategy involving optimally squeezed states and standard Ramsey measurement lead to the asymptotic precision , which can be related with the corresponding Fisher information per particle equal to:We have checked that this formula agrees with our numerical results up to the desired accuracy (), and the representative comparison of the numerical data and this formula is provided in the left inset of Fig. 2b. This implies that similarly as in the uncorrelated dephasing models, weakly spin-squeezed states are asymptotically optimal.

Using the above formula, we may also go back to the original problem of magnetic field sensing. Utilizing the relation we get the corresponding magnetic field sensing precision:The above formula assumes a fixed interrogation time t. We may generalize the considerations, and fix the total interrogation time T, which we allow to split into T/t independent interrogation steps. The corresponding estimation uncertainty reads:which when optimized over t reaches the minimal value when and yields:

Based on the above results we can expect that weakly spin-squeezed states should also be optimal in case of a more general dephasing noise, provided the range of correlations r is finite and we consider the asymptotic limit . In this case, following analogous calculations, we would arrive at the optimal magnetic field sensing precision of the formwhere represent magnetic field correlations for particles at distance : . Comparing this result with the performance of the GHZ states for the same model, see ref. [38], we notice that there is the factor improvement in performance of the optimally spin-squeezed states over the GHZ states familiar from uncorrelated dephasing considerations[36,37].

Atomic clock stabilization

A typical atomic clock operates in a feedback loop where the local oscillator (LO, e.g., laser) is stabilized to atomic reference frequency by periodically interrogating atoms (using radiation from the LO) and based on the measured response, the frequency of the LO is corrected[17]. One of the main goals in the design of the clock interrogation scheme is to achieve the lowest instability typically quantified by the Allan variance (AVAR)where represents averaging over frequency fluctuations of the LO described by some stochastic process, denotes averaging time, atomic reference angular frequency and time-dependent angular frequency of the LO. The key feature from our perspective is the fact that the LO frequency fluctuations are temporally correlated and lead effectively to a temporarily correlated dephasing of atoms.

Fixing the physical properties of the atoms the goal is to optimize their initial states, interrogations times, measurements and feedback corrections in order to minimize the AVAR. Performing such a comprehensive optimization is not feasible. In ref. [20] a lower bound on the achievable AVAR was introduced the quantum Allan Variance (QAVAR):where is the AVAR of free running LO, represents a correction to it from the feedback loop and T is the interrogation time. We do not provide here explicit forms of the operators and , and refer the interested reader to ref. [20], but just note that the structure of the formulas are similar to that in Eq. (2). The important information is that, if the atoms with which the atomic clock interacts are described via states on some d dimensional Hilbert space , then the and objects act on a tensor space , where is the number of atomic cycles that need to be considered in order to calculate QAVAR.

In Fig. 3a we present the exemplary results for an atomic clock operating on one two-level atom, where the LO fluctuations are characterized by the autocorrelation function , which is a combination of Ornstein-Uhlenbeck (OU) process and white Gaussian frequency noise (, , ) (see more detailed discussion of the model in Supplementary Note 3). In the long averaging time limit QAVAR takes the form with some constant c, which we will refer to as asymptotic coefficient. The results indicate that completely neglecting noise correlations in the analysis of the clock performance is unjustified. These calculations have been attempted in ref. [20] using the full-Hilbert space description, but were not capable of approaching the regime where the character of the scaling of the QAVAR and the coefficient could be unambiguously read out.

Fig. 3

Quantum Allan Variance in the atomic clock stabilization model.

a QAVAR (times ) as a function of the averaging time for the atomic clock (based on one atom) with the local oscillator (LO) noise, which is strictly local (yellow line) or also includes the nearest neighbors correlations (orange line), plotted against the AVAR of uncorrected LO (black line). b QAVAR asymptotic coefficient as a function of the number of atoms in the atomic clock with LO noise, which is strictly local (yellow dots connected by solid line/light gray dots connected by dotted line) or also includes the nearest neighbors correlations (orange dots connected by solid line/gray dots connected by dotted line) for the optimal/NOON state, plotted against the AVAR asymptotic coefficient of uncorrected LO (black dots connected by solid line). The sphere depicts Husimi Q distribution on the Bloch sphere for the optimal state of ten atoms with marked equipotential lines where quasiprobability is equal to 0.1 for coherent spin state (black dashed line) and for this optimal state (black solid line), which shows that the state is squeezed.

Following this approach we calculate the QAVAR asymptotic coefficient c and the corresponding optimal interrogation time T as a function of the number of atoms in the clock, see Fig. 3b. From this figure we see that the differences in QAVAR between cases with strictly local noise and when nearest-neighbor noise correlations are included only grow with the increasing number of atoms. This implies that noise correlations play an important role in the accurate analysis of clock performance. We confront the results (which are optimized over the input state) with values obtained using a NOON/GHZ states as an input, . The NOON states are highly prone to dephasing noise, and hence the optimal interrogation times will be necessary reduced compared to the optimal, more robust states. This is a manifestation of a generic poor performance of the NOON/GHZ states in realistic (noisy) scenarios with increasing particle number N[10,11,36].

Fidelity susceptibility for many-body thermal states

In condensed matter context, fidelity between many-body states and , that differ by a small variation of a parameter in a Hamiltonian, is a mean to identify the location of a phase transition[41,42]. This is where the fidelity susceptibility , defined by , has a maximum indicating a fundamental change in the state of the system. This concept was employed in ref. [43] to evaluate the usefulness of a quantum phase transition, that happens at zero-temperature, for precise sensing of the parameter in a realistic system at a finite temperature. QFI defines a metric in the space of quantum states (the Bures metric)[32] and is directly related with the fidelity susceptibility, namely .

Unlike at zero-temperature, the fidelity between a thermal many-body states represented by MPOs is not tractable in general. This is why a quasi-fidelity was employed defining a quasi-susceptibility, , that provides bounds for the exact fidelity susceptibility[40]The Hamiltonian considered in ref. [43] was the spin- XX modelwith a quantum critical point at . Taking the MPOs studied in ref. [43], we bypass the tractability problem employing the part of our scheme with to calculate the QFI . This exact susceptibility for the chain with 64 spins is shown in Fig. 4, together with the upper and lower bounds (Eq. 25). The accuracy of the fidelity susceptibility is limited by the finite bond-dimension as well the finite parameter difference . Nevertheless, we obtain satisfying results with relative error around for and .

Fig. 4

Fidelity susceptibility for the thermal state.

Exact fidelity susceptibility for a thermal many-body state (dots connected by the line) at the critical point in function of dimensionless inverse temperature in the XX model (Eq. 26) with 64 spins. The shaded band shows the bounds (Eq. 25). As predicted in refs. [39,40], the exact value tends to the upper/lower bound for high/low temperatures.

This example demonstrates that the scheme for calculating fidelity susceptibility is a useful byproduct of our general algorithm. Beyond the present metrological context, it paves a way to generalize the zero-temperature fidelity approach to detecting quantum phase transitions[41,42]—by now standard in condensed matter physics—to phase transitions in quantum many-body systems at finite temperature.

Discussion

We have provided a comprehensive framework for optimization of quantum metrological protocols using the tensor-network formalism. The potential to deal effectively with correlated noise models, as well as directly access the asymptotic is what makes this framework unique. We also expect that this framework may be adapted to deal with even more challenging metrological problems including noisy multiparameter estimation[44,45], waveform estimation[46,47], or the study of the effectiveness of adaptive metrological protocols, including quantum error-correction-based schemes[16,26,48,49]. We also expect that this numerical framework may be crucial for understanding better metrological models with temporally correlated noise especially of non-Markovian nature[22], where effective tools to find the optimal metrological protocols in such cases are yet to be developed.

Methods

We have implemented all of our algorithms in MATLAB with the help of ncon() function[50] for tensor contraction.

Optimization in the finite number of particles regime

As discussed in the main text maximization of our FoM , see Eq. (3), leading to the maximal possible QFI for a given metrological model is a two-step iterative process. First, we show how to maximize the FoM over a Hermitian operator L with fixed , and then we focus on the maximization over the input state with fixed L.

We search for the optimal L in the form of an MPOwith a finite bond-dimension . Without loss of generality, each is assumed Hermitian in its physical indices,to ensure that L is Hermitian. We call this condition a Hermitian gauge.

The bond-dimension for L is expected to be small for weakly correlated noise models. Indeed, in the limit of an uncorrelated product state , L is a sum of local operators, . Here, is the SLD for the single-particle problem applied to particle n. In a similar way as for presented in the main text, the sum can be represented by an MPO with a bond-dimension 2. Therefore, is the limiting value for uncorrelated noise models.

Graphical representation of the FoM reads:Maximization of the FoM over L is equivalent to a joint maximization over each tensor . We relax this optimization problem by iterating an optimization loop. In the loop we first find the optimal , then , and so on up to , after which we go back to . The optimization over each is performed with all other tensors fixed. The loop is repeated until the FoM converges.

After fixing the other tensors, the FoM becomes quadratic in and the optimal is found as a solution to a linear equation. For definiteness, to explain the procedure, we focus on the generic example of . We start by vectorizing :This allows us to represent the -FoM as:which can be also written in a compact way asHere, are the elements of the vector , and are the elements of the matrix A. Both and A describe the entire tensor-network complementing the distinguished vector in the two respective terms of the -FoM. After taking a derivative with respect to , we obtain a linear equation for the extremum:The matrix typically has a non-zero kernel and the linear equation does not have a unique solution. We use the Moore–Penrose pseudo-inverse, , to obtain a solution that does not contain any zero modes of .

If the linear equation was non-singular, then its exact solution would satisfy the Hermitian gauge (Eq. 28). For the typical singular case, using an SVD of to construct its pseudo-inverse, we have to truncate singular values falling below a small but finite cutoff, set by multiplied by the highest singular value. As the cutoff solution need not satisfy the Hermitian gauge condition exactly, we filter out its small anti-Hermitian part with the substitution:From experience, this substitution can improve numerical stability but is not necessary when all initial are in the Hermitian gauge (Eq. 28) and is large enough to suppress the anti-Hermitian part of the solution. However, with too large a cutoff the final optimized L does not reach the maximal possible value of the QFI. Therefore, we adjust to obtain the highest QFI achievable without compromising the stability.

Now we move on to the maximization of the FoM over the input state for a fixed L. We start by rewriting aswhere by we denote the channel, which is dual to (the evolution written in the Heisenberg picture). We can rewrite this aswhere we introduce , and . By analogy with the construction of the MPO representation for and , we can easily construct the MPO representation of and from the known MPO form of L. The tensors determining the MPO form of and are denoted by and , respectively, and their respective bond-dimensions are and .

The quantity in Eq. (33) is maximal when is a projection on the eigenvector associated with the maximal eigenvalue of the Hermitian operator . Hence, without loss of generality, we can assume a pure input state with being an MPS with bond-dimension :The input state has bond-dimension and its MPO tensors are .

The maximization of over the input state is equivalent to the variational optimization for the ground state of a many-body “Hamiltonian” , a problem widely discussed in the many-body physics MPS literature[27,51,52]. After reinterpreting our problem as a variational minimization of the “energy”we proceed iteratively in a similar way as in the case of the maximization of over L.

For example, in order to find the minimum over , we begin by vectorizing the tensor and expressing the “energy” (Eq. 35) as a Rayleigh quotientThe tensor-network form of matrices and is given byAfter taking the derivative of Eq. (36) we obtain the condition for the extremum:which is a generalized eigenvalue problem with eigenvalue . By multiplying it with a pseudo-inverse of matrix we bring it into the form of ordinary eigenvalue problem, for which we obtain the lowest eigenvalue and its corresponding eigenvector using the Lanczos algorithm.

Modification of the entire tensor-network framework to calculate the maximal QFI for systems with open boundary conditions (OBC) poses no problem and, for systems, which are not sensitive to boundary conditions, is even advisable. In OBC the MPS representing can be brought to a canonical form, where the matrix becomes an identity and Eq. (37) reduces to a standard eigenvalue problem. There is no need to pseudo-invert .

To summarize the procedure: one needs to iteratively determine the optimal L for a given and then the optimal for a given L, until one observes convergence of the final result, e.g., the FoM does not change more than, say, after a fixed number of steps. From our numerical experience this happens very rapidly, typically after five iterations of the and L optimization steps.

While running the algorithm, one has to choose the bond-dimensions for the input state, , as well as for the SLD, , over which the optimization is performed. As in all tensor-network algorithms, keeping the bond-dimension as low as possible is essential for their efficiency. In our calculations, we typically started with a product input state, , and optimized for L with a minimal non-trivial . Then we increased one of the bond-dimensions, either or , each time repeating the optimization procedure, until we found that the QFI did not change when increasing the D’s more than by, e.g., and hence assumed that relative error of our method is around .

Optimization in the asymptotic limit

The previous subsection described an algorithm that functions for a system with a finite number of particles N. For quantum metrological problems in the presence of decoherence, it is the generic situation that the optimal quantum enhancement thanks to the use of entanglement leads asymptotically (for large N) to an improvement by a constant factor over product-state strategies. Even though the finite-system MPO approach allows us to achieve values of N that are inaccessible via exact full-Hilbert space computation, it may sometimes be not enough to reach the asymptotic limit and determine the quantum enhancement coefficient with the desired precision. For this reason, we would like to have a procedure that allows us to go directly to the infinite particle limit, calculate the maximal achievable QFI per particle and, as a result, determine the maximal quantum enhancement coefficient.

For this purpose, we exploit the infinite MPO/MPS (iMPO/iMPS) approach, see e.g., ref. [29]. We assume that all tensors are translationally invariant (TI). Then we take the limit of infinite N. Technically this limit is most natural in the case of PBCs, where it is enough to notice that, for any TI transfer matrix E, the spectral decomposition of is dominated by the leading eigenvalue and eigenvector of E. This is why in the following discussion we proceed with PBCs. In the OBC case, which is arguably more natural in some metrological contexts, one has in principle to consider the boundary conditions at infinity. However, applied to a boundary vector gives the leading eigenvector of E when k becomes longer than the finite correlation range (just as in the Lanczos algorithm). Therefore, in the bulk (i.e., far from the boundaries), all equations the TI tensors have to satisfy become the same as for the PBC.

When the input state is TI then the final state is TI as well. There is a problem, however, with the operator . Its construction as an MPO in Eq. (8) is not TI. Because of this, instead of calculating the derivative exactly, we approximate it by a difference of two TI iMPOs:with infinitesimal parameter . Motivated by the defining equation for the SLDwe can consider a similar to Eq. (38) expansion of the operator L:Here, and are solutions of Eq. (39) when is replaced by, respectively, and . We search for the optimal , which is better suited for the TI formalism than the original operator L. Let us denote bythe QFI per particle that we want to maximize:A TI iMPO is defined by only one tensor, which we assign, respectively, as: , , , , , .

Optimization of a tensor-network consisting of identical tensors A is a highly nonlinear problem in A and one might think that an approach similar to the one used in the previous subsection is not applicable here. Fortunately, an efficient method for the problem was developed in ref. [53]. The main idea is to find the optimal tensor at one-site, treating all other tensors A as fixed, and then rather then replacing all tensors by perform a flexible substitution,with a mixing angle . The angle is optimized to yield the best possible FoM.

We now apply this approach to our two-step iterative procedure. First we need to find the optimal , which is equivalent to the determination of the optimal local tensor .

As explained in Supplementary Note 4, the trace of an operator represented as an iMPO is defined by its transfer matrix, so we start by introducing transfer matrices and associated with, respectively, and , which are depicted below alongside their eigendecomposition:Using this transfer matrices we can write , and with the fact that is determined by its leading eigenvalue we can write . Now we can construct iMPO representation of the FoM:where we depicted and its version with tensor distinguished (after performing vectorization).

This FoM can be also expressed in form of equation:which give us condition for an extremum:

For the powers of the eigenvalues may seem to pose a problem. Fortunately, however, this problem can be circumvented. For a given we can calculate the associated value of our FoM per particle:but going back to Eq. (2) we see that FoM per particle should have form , where and are of the same order of magnitude as the asymptotic limit of the QFI per particle. Assuming that our calculations are in the regime of , , and , and remembering binomial expansionit is to be expected that the highest eigenvalues of the transfer matrices have the form:which after inserting into Eq. (46) and using binomial expansion to the first order give us exactly . Note that, it means that we can calculate value of FoM per particle in a simple way:It is clear now that for the purpose of solving Eq. (45) we can approximate and by ones and, hence, bring the condition for the optimal to a simpler form:This equation is solved with a pseudo-inverse and its anti-Hermitian part is filtered out at every iteration step.

SLD is always traceless in any unitary parameter estimation problem, which can be seen by solving Eq. (39) for the SLD using eigenbasis of :and taking into account that . We can ensure that solution has proper normalization using the condition or, equivalently , which in the language of transfer matrices means that the highest eigenvalue of the transfer matrix associated with operator has to be equal to d.

Now we turn to the second part of our optimization procedure, namely the variational minimization over the input state. This step does not introduce any qualitatively new challenges, so we only briefly discuss it for completeness. As for the , we approximate the exact derivative of by its discrete version:After the expansion , our task becomes equivalent to minimization of the “energy density”:over . We can depict Eq. (53) in a diagrammatic form:where we have used transfer matrices , , and associated with, respectively, , , and that are represented graphically asAs previously, we expect that and for the purpose of finding the optimal tensor P, we can approximate and by ones. After taking the derivative we obtain the condition for the extremum:where is a generalized eigenvalue, whereas the matrices and have the following tensor-network representationThe matrix is a tensor product of three matrices, , the identity, and , hence its pseudo-inverse can be obtained as a tensor product of (pseudo-)inverses of the smaller matrices. Applying to Eq. (54) we bring it into the form of a standard eigenvalue problem. We solve this eigenproblem with respect to the smallest eigenvalue and its corresponding eigenvector and require that is normalized so that . Then we calculate the asymptotic value of the QFI per particle:

Just as for the finite N, iterating the and optimization steps leads to the optimal solution with the maximal QFI per particle.

While performing the numerics one should choose to be small but not too small as too small values may lead to numerical instabilities. Our general strategy in obtaining numerical results in the examples presented, was to lower the value of until we observed no noticeable change in the obtained results, while still remaining in the regime where algorithm was stable. In all the examples we studied in this paper this approach resulted in the choice of (instabilities started to appear for ). Notice that in the asymptotic iMPO approach described above we required to be small—on the order of the precision we expect from the numerical results. In other words setting the precision requirements to this implies that , and hence . What this physically means is that in our set-up the QFI per particle does not change in any noticeable way for larger N and hence the asymptotic behavior in the actual limit may be inferred from this results.