Tensor-Network Approach for Quantum Metrology in Many-Body Quantum Systems

TL;DR

This paper introduces a tensor network framework using matrix product operators to efficiently identify optimal quantum metrological protocols in many-body systems with noise, enabling analysis of asymptotic precision and probe states.

Contribution

It develops a comprehensive MPO-based approach for quantum metrology that handles noise correlations and allows direct computation of asymptotic precision in large systems.

Findings

Efficient identification of optimal quantum states under noise.

Application to atomic clock stabilization and magnetic sensing.

Methods for calculating quantum Fisher information and fidelity susceptibility.

Abstract

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. Thanks to the fact that the MPO formalism allows for an efficient description of short-range spatial and temporal noise correlations, the maximal achievable estimation precision in such models, as well as the optimal probe states in previously inaccessible regimes can be identified. Moreover, the application of infinite MPO (iMPO) techniques allows for a direct and efficient determination of the asymptotic precision of optimal protocols 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 for magnetic field sensing in the presence of locally correlated magnetic field fluctuations (spatial noise correlations). As a byproduct, the developed methods for calculating the quantum Fisher information via MPOs may be used to calculate the fidelity susceptibility - a parameter widely used in many-body physics to study phase transitions.

Equations

$$p(x|\phi )=\mathrm{Tr}({\rho }_{\phi }{\Pi }_x),$$

$${\rho }_{\phi }={\Lambda }_{\phi }\left({\rho }_0\right).$$

$${\Delta }^2\tilde{\phi }=⟨(\tilde{\phi }-\phi {)}^2⟩,$$

$${\Delta }^2\tilde{\phi }\ge \frac{1}{F({\rho }_0)},$$

$$F({\rho }_0)=\mathrm{Tr}({\rho }_{\phi }{L}^2),$$

$${\rho }_{\phi }^{\prime }=\frac{1}{2}\left(L{\rho }_{\phi }+{\rho }_{\phi }L\right),$$

$$L=\sum _{j,k}\frac{2⟨{\lambda }_j|{\rho }_{\phi }^{\prime }|{\lambda }_k⟩}{{\lambda }_j+{\lambda }_k}|{\lambda }_j⟩⟨{\lambda }_k|.$$

$$F({\rho }_0)=\underset{L}{\sup }F({\rho }_0,L),F({\rho }_0,L)=2\mathrm{Tr}\left({\rho }_{\phi }^{\prime }L\right)-\mathrm{Tr}\left({\rho }_{\phi }{L}^2\right),$$

$$F=\underset{{\rho }_0}{\sup }F({\rho }_0)=\underset{{\rho }_0,L}{\sup }F({\rho }_0,L),$$

$${\rho }_{\phi }={\Lambda }_{\phi }({\rho }_0)=e^{-\mathrm{i}H\phi }\Lambda ({\rho }_0)e^{\mathrm{i}H\phi },H=\sum _{n=1}^N{h}^{[n]},$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$${\mathcal{L}}^{(k,n)}(\rho )$$

$$\mathcal{D}[c](\rho )=c\rho c^{†}-\frac{1}{2}\left(\rho c^{†}c+c^{†}c\rho \right).$$

$$\Lambda =\exp \left(\sum _{k=1}^r\sum _{n=1}^N{\mathcal{L}}^{(k,n)}t\right).$$

$$\sum _{l=0}^{d-1}{T}_{l{j}_2\dots j_N}{W}_{k_1\dots k_{3}l{k}_5\dots k_M},$$

$$\begin{gathered} \mu =\sum _{\begin{array}{c}{j}_1,\dots ,j_N,\\ k_1,\dots ,k_N=0\end{array}}^{d-1}\mathrm{Tr}\left(A[1{]}_{k_1}^{j_1}A[2{]}_{k_2}^{j_2}\dots A[N{]}_{k_N}^{j_N}\right)\times \\ |j_1,j_2,\dots ,j_N⟩⟨k_1,k_2,\dots ,k_N|. \end{gathered}$$

$$S({\rho }_{[j,k]})\le 2{\log }_2(D),$$

$$T=US{V}^{†}.$$

$$\mathrm{Tr}(\mu )=\underset{n\to \infty }{\lim }\mathrm{Tr}(E^n).$$

$$E^n\sim {\lambda }_1^n|r_1)(l_1|,$$

$$\underset{n\to \infty }{\lim }\frac{\mathrm{Tr}(E^n)}{{\lambda }_1^n}$$

$${\rho }_0=\sum _{\mathbf{j},\mathbf{k}}\mathrm{Tr}\left(R_0[1{]}_{k_1}^{j_1}{R}_0[2{]}_{k_2}^{j_2}\dots R_0[N{]}_{k_N}^{j_N}\right)|\mathbf{j}⟩⟨\mathbf{k}|,$$

$$X=\sum _{{\alpha }^{\prime },\alpha ,{\beta }^{\prime },\beta }{X}_{{\alpha }^{\prime },\alpha ,{\beta }^{\prime },\beta }|{\alpha }^{\prime }⟩⟨\alpha |\otimes |{\beta }^{\prime }⟩⟨\beta |.$$

$$x=\sum _{({\alpha }^{\prime },\alpha ),({\beta }^{\prime },\beta )}{x}_{({\alpha }^{\prime },\alpha ),({\beta }^{\prime },\beta )}|({\alpha }^{\prime },\alpha )⟩⟨({\beta }^{\prime },\beta )|.$$

$$x_{({\alpha }^{\prime },\alpha ),({\beta }^{\prime },\beta )}=\sum _{\gamma ,{\gamma }^{\prime }}{U}_{({\alpha }^{\prime },\alpha ),\gamma }{S}_{\gamma ,{\gamma }^{\prime }}{\left[V^{†}\right]}_{{\gamma }^{\prime },({\beta }^{\prime },\beta )},$$

$${\rho }_{\phi }^{\prime }=\sum _{n=1}^N\mathrm{i}\left[{\rho }_{\phi },h^{[n]}\right].$$

$${\rho }_{\phi }^{\prime }=\sum _{\mathbf{j},\mathbf{k}}\mathrm{Tr}\left(R^{\prime }[1{]}_{k_1}^{j_1}{R}^{\prime }[2{]}_{k_2}^{j_2}\dots R^{\prime }[N{]}_{k_N}^{j_N}\right)|\mathbf{j}⟩⟨\mathbf{k}|,$$

$$\{\begin{array}{ll}\left(\begin{array}{cc}\mathrm{i}({ϵ}_{k_1}-{ϵ}_{j_1})& 1\\ 0& 0\end{array}\right)\otimes R_{\phi }[1{]}_{k_1}^{j_1}& \text{for }n=1,\\ \left(\begin{array}{cc}1& 0\\ \mathrm{i}({ϵ}_{k_n}-{ϵ}_{j_n})& 1\end{array}\right)\otimes R_{\phi }[n{]}_{k_n}^{j_n}& \text{for }n\in [2,\dots ,N-1],\\ \left(\begin{array}{cc}1& 0\\ \mathrm{i}({ϵ}_{k_N}-{ϵ}_{j_N})& 0\end{array}\right)\otimes R_{\phi }[N{]}_{k_N}^{j_N}& \text{for }n=N.\end{array}$$

$$L=\sum _{\mathbf{j},\mathbf{k}}\mathrm{Tr}\left(S[1{]}_{k_1}^{j_1}S[2{]}_{k_2}^{j_2}\dots S[N{]}_{k_N}^{j_N}\right)|\mathbf{j}⟩⟨\mathbf{k}|,$$

$$S[n{]}_{k_n}^{j_n}=\overline{S[n{]}_{j_n}^{k_n}},$$

$$F({\rho }_0,L)=2\sum _{\alpha }{b}_{\alpha }S[2{]}_{\alpha }-\sum _{\alpha \beta }S[2{]}_{\alpha }{A}_{\alpha \beta }S[2{]}_{\beta }.$$

Full text

I Introduction

Quantum metrology Giovannetti et al. (2006); Paris (2009); Tóth and Apellaniz (2014); Demkowicz-Dobrzański et al. (2015); Schnabel (2017); Degen et al. (2017); Pezzè et al. (2018); Pirandola et al. (2018) 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 Feynman (1982). While small-scale problems are feasible via direct numerical study, even a slight increase in the number of elementary objects quickly makes such an approach intractable. This presents a serious conceptual roadblock for the development of quantum technologies—once inside the quantum enhanced regime one can no longer validate the performance of quantum devices using naive classical simulation Aolita et al. (2015). For example, just storing the tomographic reconstruction of a multiqubit quantum state quickly becomes unfeasible for more than 40 particles Cramer et al. (2010).

Despite the curse of dimensionality there has nevertheless been considerable progress in quantum metrology. This is because in the noiseless case Giovannetti et al. (2006), as well as in some noisy models including e.g. perfectly correlated dephasing Dorner (2012); Macieszczak et al. (2014), the description of metrological phase/frequency estimation protocols may, without loss of generality, be restricted to the fully symmetric subspace. Since the dimension of the symmetric subspace scales linearly with particle number it is possible to perform direct and efficient calculations in the large particle limit. Apart from that, a number of powerful methods have been designed to obtain fundamental precision bounds for uncorrelated noise models, circumventing the issue that the states involved in the protocols are no longer restricted to the fully symmetric subspace Escher et al. (2011); Demkowicz-Dobrzański et al. (2012); Kołodyński and Demkowicz-Dobrzański (2013); Knysh et al. (2014); Demkowicz-Dobrzański and Maccone (2014); Sekatski et al. (2017); Demkowicz-Dobrzański et al. (2017); Zhou et al. (2018). However, in cases when one deals with partially correlated noise metrological models, or whenever one wishes to study the performance of states outside the fully symmetric subspace, there are no efficient methods that can be applied and one is forced to restrict considerations to small-scale problems.

There are several ways in which correlated noise manifests itself in metrological problems. The first is when environmental memory effects are significant, leading to probe dynamics with a time-correlated noise. The simplest example here is that of the atomic clock stabilization problem, where the effective dephasing process of atoms is temporally correlated as a result of correlations of the frequency fluctuations of the local oscillator Ludlow et al. (2015). Because of this feature, identification of the optimally quantum clock stabilization strategies taking into account all the possible states of the atoms as well as quantum measurement and feedback strategies is a highly non-trivial task André et al. (2004); Borregaard and Sørensen (2013); Kessler et al. (2014). One of the approaches, addressing the optimality of the protocols to stabilise atomic clocks, was based on the concept of the quantum Allan variance (QAVAR) Chabuda et al. (2016). Unfortunately the computation of the QAVAR for even small systems becomes quickly unfeasible, as its computational complexity grows exponentially with the number of atomic clock interrogations. An even more challenging case involves models where time-correlated noise cannot be effectively described via some classical stochastic process Clerk et al. (2010); Szańkowski et al. (2017) and as such manifests non-Markovian features of quantum dynamics Chin et al. (2012), as e.g. in NV-center sensing models Rondin et al. (2014); Paz-Silva et al. (2017); Kwiatkowski and Cywiński (2018). A second natural setting exhibiting nontrivial noise correlations is that of many-body systems such as, e.g., spin chains. Here, typically, spatially correlated noise emerges. This is of crucial relevance for any models where the effective signal is obtained from spatially distributed probes in the mean field estimation Jeske et al. (2014) or e.g. in field-gradient metrology Altenburg et al. (2017); Apellaniz et al. (2018).

Although realistic quantum noise is unlikely to be strictly uncorrelated, it is not going to be arbitrarily complicated. On the contrary, 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. For typical short-range correlated noise processes the optimal performance of a metrological protocol is expected to be attainable with input probe states exhibiting entanglement within groups of a finite number of particles similarly as in the uncorrelated noise models Jarzyna and Demkowicz-Dobrzański (2013). This key physical insight is a strong hint for how to go beyond the extant methods for uncorrelated and Markovian noise: one needs a way to represent short-range correlations in many body systems.

The most successful approach for classically simulating short-range correlated many body systems is via the variational tensor-network state (TNS) ansatz (see, e.g., Bridgeman and Chubb (2017) for a recent review). Expressive ansatz classes such as the matrix product states (MPS) Fannes et al. (1992); Östlund and Rommer (1995), projected entangled-pair states (PEPS) Verstraete and Cirac (2004), or the multiscale entanglement renormalisation ansatz (MERA) Vidal (2007), have led to unparalleled insights into the physics of quantum many body systems ranging from the description of quantum phase transitions to new phases of matter such as topological order. Most relevant for the present work is the application of tensor networks to the study of dissipative physics via the matrix product operator (MPO) ansatz for density operators Zwolak and Vidal (2004); Verstraete et al. (2004a) and also their infinite particle limit known as infinite MPO (iMPO) Cincio and Vidal (2013).

With a few notable examples, the application of tensor-network methods to quantum metrology is still in its infancy. At the present time calculations of the QFI for mixed tensor-network states has been via approximate upper bounds, effectively amounting to a weighted sum of the QFI calculated on pure states Jarzyna and Demkowicz-Dobrzański (2013)— a task equivalent to calculating the variance of a local observable. A further problem is that quantities such as the Quantum Fisher Information (QFI) or the Bayesian cost are not familiar in the tensor-network literature and, until now, it was far from clear whether they can be easily calculated at all in the mixed-state setting.

In this Paper we develop tensor-network based methods to address quantum metrological problems, in particular that involving short-range quantum noise. We demonstrate how to directly calculate relevant metrological quantities, such as the QFI and Bayesian-type cost, using the MPO formalism. Moreover, we also show how to optimize input probe states for metrological problems using the same formalism. As a result, we provide an efficient iterative procedure to design optimal metrological protocols that remains within the MPO formalism and does not suffer from the curse of dimensionality.

The Paper is extensive as its aim is also to bridge the gap between the tensor network and quantum metrology communities and the outline is as follows. In Sec. II we review the main concepts of quantum metrology, including, the Cramér-Rao bound and the quantum Fisher information, and specify the optimization algorithm that leads to the identification of the optimal metrological protocols. We also define the class of metrological models with correlated (but short range) noise for which the developed MPO methods are expected to be efficient. In Sec. III we review the tensor networks formalism with a particular focus on the MPO and MPS construction and describe the whole metrological optimization algorithm in the language of MPO operations. We also present the iMPO approach where asymptotic performance of metrological protocols in the limit of large number of particles can be obtained in a direct way. This is followed, in Sec. IV, with a series of applications including: magnetic field sensing in presence of spatially correlated fluctuating field, calculation of the QAVAR for the atomic clock stabilisation problem taking into account temporal correlations of the local oscillator fluctuations, and finally demonstrate how the developed tools may be applied to problems outside the field of quantum metrology, namely, the calculation of the fidelity susceptibility of a finite-temperature thermal state of a many-body spin system. Finally, in Sec. V we conclude and provide future directions.

II Optimal quantum metrological protocols

II.1 General models

A typical problem in quantum metrology may be formulated as follows:

\rho_{0}$$\Lambda_{\varphi}$$\Pi_{x}$$x$$\widetilde{\varphi}(x)

.

Here a probe state $\rho_{0}$ is subject to a parameter-dependent quantum evolution, mathematically represented by a parameter-dependent quantum channel $\Lambda_{\varphi}$. A POVM type measurement $\{{\Pi}_{x}\}_{x}$ Nielsen and Chuang (2000) is then carried yielding a conditional probability distribution

[TABLE]

where

[TABLE]

Given the conditional probability distribution $p(x|\varphi)$ the objective is to estimate the value of the unknown parameter $\varphi$. To this end one employs an estimator function $\widetilde{\varphi}(x)$ to produce a given estimate for $\varphi$ given the measurement outcome $x$. The performance of an estimator function is limited by the average uncertainty

[TABLE]

where the expectation is over all measurement results $x$. The central goal of quantum metrology is to find the best input probe state $\rho_{0}$, the best measurement and estimator so as to minimize $\Delta^{2}\widetilde{\varphi}$. This is a challenging problem in general as both the optimal probe state and the measurement depend in a deeply nontrivial way on the channel $\Lambda_{\varphi}$.

Progress can be made by exploiting a fundamental result in quantum metrology, namely, the Cramér-Rao inequality Helstrom (1976); Holevo (1982); Braunstein and Caves (1994). This crucial result lower-bounds the average uncertainty of the best possible estimator $\widetilde{\varphi}(x)$ in terms of a quantity known as the QFI:

[TABLE]

where

[TABLE]

and $L$ is the symmetric logarithmic derivative (SLD) defined implicitly via

[TABLE]

where $\rho^{\prime}_{\varphi}$ is derivative of $\rho_{\varphi}$ with respect to $\varphi$. Instead of directly minimising $\Delta^{2}\widetilde{\varphi}$ from Eq. (3) one can instead maximise the QFI over input states $\rho_{0}$, an easier task in general as the optimization of the measurement $\Pi_{x}$ and the estimator is no longer required.

A calculation of the QFI requires solving Eq. (6) for the SLD operator $L$. This is a linear equation; its solution generally requires finding the eigenvalues and eigenvectors of $\rho_{\varphi}=\sum_{j}\lambda_{j}|\lambda_{j}\rangle\langle\lambda_{j}|$ Helstrom (1976); Holevo (1982); Braunstein and Caves (1994):

[TABLE]

Performing a full eigendecomposition is generally unfeasible for a system comprised of more than a small number of particles, and cannot be easily implemented using the efficient MPO description advocated in this paper. For this reason we do not use the above formula. Instead, it is much more numerically efficient to solve the linear equation (6) directly using standard linear equation solving methods in order to find $L$. Alternatively, one may reformulate the calculation of the QFI as the following maximization problem Macieszczak (2013); Macieszczak et al. (2014):

[TABLE]

where we will refer to $F(\rho_{0},L)$ as the figure of merit (FoM) for the problem. To see the equivalence of the definition Eq. (5) with Eq. (8) note that the above formulation is the supremum of a quadratic function of a Hermitian operator $L$. This optimization may be solved by formally taking the derivative with respect to $L$ and setting it to zero. The resulting extremum condition yields the equation for the SLD and hence the formula Eq. (5) for the QFI.

The QFI formula (8) has an advantage over the original (5), when one wants to additionally perform optimization of the QFI over the input states $\rho_{0}$ in order to find the optimal quantum metrological protocol. Using (8) this problem can be written as a double maximization problem:

[TABLE]

where the FoM is linear in $\rho_{0}$ 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 Demkowicz-Dobrzański (2011); Macieszczak et al. (2014) in the Bayesian estimation context, and then applied to the QFI FoM in Macieszczak (2013) (recently the method has been rediscovered in a slightly modified incarnation in Tóth and Vértesi (2018) and proved useful in studying metrological properties of PPT states). As we will see, each of these iterative steps may be performed efficiently using tensor networks.

Even though, in the above formulation, we focused solely on the QFI based approach, the above considerations are valid whenever the quantity to be optimized is given in the form (8): $\rho^{\prime}_{\varphi}$ 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 such as the QAVAR as discussed in Sec. IV.2.

II.2 Many-particle models with local parameter encoding and locally correlated noise

In this subsection we describe the class of metrological models currently challenging for state-of-the-art methods. These are the main motivation for the development of MPO-based techniques. We focus on systems comprised of $N$ distinguishable $d$-dimensional particles, so that the total Hilbert space is $\mathcal{H}=\bigotimes_{n=1}^{N}\mathbb{C}^{d}$. We assume that the parameter $\varphi$ is unitarily encoded in the output state $\rho_{\varphi}$ according to a product of unitaries given by the exponential of local generators (or Hamiltonians):

[TABLE]

where $h^{[n]}$ is the generator acting on the $n$th particle.

Most importantly for this paper, the noise, represented above by the operator $\Lambda$, is not assumed to be local, which makes the problem particularly challenging. Powerful methods capable of yielding fundamental metrological bounds in the large particle number regime effectively work only in case of uncorrelated noise models Escher et al. (2011); Demkowicz-Dobrzański et al. (2012); Kołodyński and Demkowicz-Dobrzański (2013); Knysh et al. (2014); Demkowicz-Dobrzański and Maccone (2014); Demkowicz-Dobrzański et al. (2017); Zhou et al. (2018) and cannot be directly used to study the effects of noise correlations. In many physically realistic situations, however, correlated noise tends to be only locally correlated, which gives one hope that more efficient methods to deal with such problems than simple brute force numerical optimization may be found.

We assume that $\Lambda$ may be effectively approximated as a product of single ($\Lambda^{[n]}$), two- ($\Lambda^{[n,n+1]}$), three- ($\Lambda^{[n,n+1,n+2]}$), etc. particle maps up to some cut-off point after which to is assumed that noise correlations do not extend beyond $r$ neighbouring particles. From a more physical perspective, consider a time-independent quantum master equation Breuer and Petruccione (2002) describing the noisy part of the evolution of an $N$-body quantum system:

[TABLE]

where

[TABLE]

In the above, the $L^{[n,\dots,n+k-1]}_{j}$ are noise operators acting on $k$ neighbouring particles (not to be confused with the SLD), and the combined effect of the $k$-particle noise acting on the subset $[n,\dots,n+k-1]$ of particles is represented by the $\mathcal{L}^{(k,n)}$ operator, where we neglect terms with higher range than $r$—note that a three-particle term may in particular represent a two-body interaction between next-nearest neighbours. The channel $\Lambda$ can now be obtained by integrating the evolution over some fixed time $t$:

[TABLE]

If all the operators in the above exponent commute we can immediately write $\Lambda$ as a product of single, two-, three- etc. maps acting on different subsets of particles, which will lead us immediately to an efficient MPO description of the dynamics, see Sec. III. Otherwise, one may approximate the evolution for a time $t$ as a product of short time steps, where in each time step we perform the Suzuki-Trotter decomposition Trotter (1959); Suzuki (1966, 1976).

In what follows we often use a vectorized density matrix notation $\rho\rightarrow|\rho\rangle$, which is useful in the MPO approach, and where clear from context, we switch between the cases where $\Lambda$ is understood as acting on $\rho$ or $|\rho\rangle$. For example, $|\Lambda(\rho_{0})\rangle=\Lambda|\rho_{0}\rangle$.

For definiteness, we assumed above that the noise acts before the unitary encoding. This entails no loss of generality if the noise commutes with the encoding, as is the case in the most popular metrological models of phase/frequency estimation in presence of dephasing or loss Dorner et al. (2009); Escher et al. (2011); Demkowicz-Dobrzański et al. (2012). If needed, one can extend our formalism to the case where the parameter dependence no longer commutes. This comes at the expense of a slightly higher complexity of formulas and numerics, since we are no longer able to write the derivative of $\rho_{\varphi}$ over the parameter as a commutator with the Hamiltonian. Instead, the entire channel structure $\Lambda_{\varphi}$ determines the derivative.

Note that, thanks to the assumptions of the model, if the input probe state $\rho_{0}$ is only locally correlated, which is a sufficient condition for an efficient MPO description, then it remains so under the above evolution. Moreover, the derivative of the output state with respect to the parameter $\varphi$—required for calculations of the QFI—reads $\rho^{\prime}_{\varphi}=\mathrm{i}\left[\rho_{\varphi},H\right]$ and, since $H$ is the sum of local Hamiltonians, is also be efficiently describable using MPO.

Most importantly, in essentially all realistic metrological protocols, the maximal achievable QFI scales linearly with $N$ in the limit of large particle numbers Escher et al. (2011); Demkowicz-Dobrzański et al. (2012), and the quantum-enhancement advantage appears in the form of a constant factor. This further implies that it is enough to consider input states $\rho_{0}$ with finite-range correlations to achieve almost optimal metrological performance Jarzyna and Demkowicz-Dobrzański (2013). As a result the MPO formalism is ideally suited to tackle this class of metrological problems and guarantees that the optimal values of QFI as well as the optimal probe states will be found via this approach.

III Matrix product operator approach

In this section we present the formalism of MPO (and also closely connected MPS) type tensor networks, adapted to quantum metrological problems, and show effective ways to use it in order to solve the optimizations formulated in the previous section.

III.1 Review of tensor networks

First we give a short review of the tensor-network formalism (for a more comprehensive recent review see Bridgeman and Chubb (2017)). To succinctly describe tensor networks we represent tensors diagrammatically: suppose that $T$ is a tensor with $N$ indices, each ranging from [math] to $d-1$, elements of this tensor we denote as $T_{j_{1}j_{2}\dots j_{N}}$. We depict such a tensor as a circle with $N$ legs, each labelled with an index:

T_{j_{1}j_{2}\dots j_{N}}=$$j_{1}$$j_{2}$$j_{3}$$j_{4}$$j_{5}$$j_{N}$$T$$\dots

.

Tensors in this paper are nothing more than lists of numbers; we don’t assume any particular transformation properties for our tensors. Some simple special cases include kets $|v\rangle$, bras $\langle v|$ ($N=1$) and operators $\mathcal{M}$ ($N=2$):

$v$

,         

$\overline{v}$ ,         

$\mathcal{M}$ ,

where the overline denotes complex conjugation.

Tensor contraction, whereby the components of two tensors $T$ and $W$ are multiplied and summed over repeated indices, is depicted by connecting the legs, e.g.,

j_{2}$$j_{3}$$j_{4}$$j_{5}$$j_{N}$$T$$\dots$$k_{1}$$k_{2}$$k_{3}$$k_{5}$$k_{M}$$W$$\dots

.

Here this tensor network depicts the contraction

[TABLE]

which is a tensor with $N+M-2$ legs. By combining tensors with three or more legs via tensor contraction we can build networks of arbitrary complexity. (Tensor networks involving tensors with one or two legs are necessarily a combination of lines and cycles.)

Among various tensor network classes, the most important for this paper are the matrix product state (MPS) and matrix product operator (MPO) tensor networks. MPS/MPO are a natural compact representation for states/operators with finite correlations so we expect that they form natural language to study quantum metrology problems with locally correlated noise. An MPS with open boundary conditions (OBC) is a state of the form

|\psi\rangle=$$\dots$$A[1]$$A[2]$$A[3]$$A[4]$$A[5]$$A[N]

.

Matrix product states with OBCs are now known to provide an excellent model for the ground states of one-dimensional quantum spin chains with (constant) spectral gap Hastings (2007).

We can accommodate periodic boundary conditions (PBCs) by joining the last tensor to the first via an additional “horizontal leg”:

|\psi\rangle=$$\dots$$A[1]$$A[2]$$A[3]$$A[4]$$A[5]$$A[N]

.

Such MPS with PBCs offer some numerical advantages for quantum spin systems on rings Verstraete et al. (2004b).

A matrix product operator (MPO) is a tensor network which is, in the PBCs case, a linear operator from $\left(\mathbb{C}^{d}\right)^{\otimes N}$ to itself parametrised according to

\dots$$A[1]$$A[2]$$A[3]$$A[4]$$A[5]$$A[N]

,

where the last horizontal leg is contracted with the first. Here the tensor $A[n]$ has four legs:

A[n]$$\alpha$$\beta$$k$$j

.

The horizontal-leg indices (also called the virtual indices) $\alpha$ and $\beta$ range from $1,2,\dots,D$, where the parameter $D$ is called the bond dimension and the vertical leg indices (also called the physical indices) $j$ and $k$ range from $0,1,\dots,d-1$, where $d$ is the physical dimension. Thus an MPO $\mu$ with PBCs is the following operator

[TABLE]

Exploiting channel/state duality we can bend the vertical legs upward ($|j\rangle\langle k|\rightarrow|j,k\rangle$) to write an MPO $\mu$ as an MPS $|\mu\rangle$ of $2N$ particles:

\dots$$A[1]$$A[2]$$A[3]$$A[4]$$A[5]$$A[N]

.

Defining a new vertical line

$=$

to range over a doubled index $(j,k)$, with $j,k=0,1,\dots,d-1$, we arrive at the equivalent tensor network for $|\mu\rangle$:

\dots$$A[1]$$A[2]$$A[3]$$A[4]$$A[5]$$A[N]

.

The multiplication of two MPOs $\mu$ and $\nu$, with bond dimensions $D_{1}$ and $D_{2}$, respectively, is given by

\dots$$B[1]$$B[2]$$B[3]$$B[4]$$B[5]$$B[N]$$\dots$$A[1]$$A[2]$$A[3]$$A[4]$$A[5]$$A[N]

.

Contracting the two tensors $A$ and $B$ vertically, and combining the two horizontal lines into a new horizontal line ranging from $1,2,\dots,D_{1}D_{2}$:

C[n]$$\gamma$$\delta$$k$$j$$=$$\alpha_{1}$$\beta_{1}$$j$$B[n]$$\alpha_{2}$$\beta_{2}$$k$$A[n]

,

results in a new MPO $\mu\nu$ with bond dimension $D_{1}D_{2}$:

\dots$$C[1]$$C[2]$$C[3]$$C[4]$$C[5]$$C[N]

.

The bond dimension $D$ is a refinement parameter limiting the correlations occurring in the MPO/MPS representation of an operator/state: when the correlations have a finite range there is a finite $D$ capable of representing operator/state accurately. As the cost of tensor network computations is polynomial in $D$, the MPO/MPS ansatz is a basis for powerful numerical methods. The bond dimension of an MPS is directly connected to the entanglement between a bipartition of the chain: when one computes the entanglement entropy of a contiguous collection $[j,k]=\{j,j+1,\dots,k\}$ of spins one may derive the bound

[TABLE]

on the entanglement between the region $[j,k]$ and the rest of the chain (see, e.g., Fannes et al. (1992); Vidal (2003, 2004) for an elaboration of this result amongst many others). Accordingly, if a quantum state has a large bipartite entanglement then a larger $D$ is required to represent it as an MPS, and hence the harder it is to approximate it numerically.

A central tool for tensor network manipulations is the singular value decomposition (SVD), according to which, for all operators $T$ there exist unitaries $U$ and $V$ and a diagonal matrix $S$ with non-negative real numbers on the diagonal called singular values, such that

[TABLE]

This may be diagrammatically represented as follows:

T$$=$$U$$S$$V^{\dagger}

.

So far we have discussed tensor networks for finite collections of particles. A crucial advantage of the tensor-network formalism is that we can easily extend the MPO (as well as MPS) ansatz to apply to infinite-sized systems; we simply allow the network to extend to infinity from either side:

\dots$$\dots$$\dots$$A[1]$$A[2]$$A[3]$$A[4]$$A[5]$$A[N]

.

In order to work with such networks it is expedient to assume translation invariance (TI), which is imposed by assuming the tensor does not vary from site to site, so for each $n$: $A[n]=A$. With this simple assumption it becomes possible to contract and evaluate infinite MPO (iMPO).

A key primitive operation for manipulations involving iMPOs is the trace. Diagrammatically the trace $\mathrm{Tr}(\mu)$ of an iMPO may be obtained by connecting the vertical legs:

\dots$$\dots$$A$$A$$A$$A$$A$$A

.

Define

=$$E$${A}

.

In this way we obtain for the trace $\mathrm{Tr}(\mu)$ a tensor network involving the infinite product of so-called transfer matrices $E$:

\dots$$\dots$$E$$E$$E$$E$$E$$E

or, in equations,

[TABLE]

To calculate this expression we note that when $E$ is diagonalizable (does not have to be Hermitian) we can decompose it: $E=\sum_{j}\lambda_{j}|r_{j})(l_{j}|$, where $|r_{j})$, $(l_{j}|$ are respectively right and left eigenvectors of $E$ which can be normalized that $(l_{i}|r_{j})=\delta_{ij}$. This decomposition in diagrammes looks like

E$$\displaystyle=\sum_{j}\lambda_{j}$$r_{j}$$l_{j}

,

and give us dominant contribution of $E^{n}$, determined by the (here assumed unique) leading eigenvalue $\lambda_{1}$:

[TABLE]

and $\mathrm{Tr}(E^{n})\sim\lambda_{1}^{n}$. Thus, the limit in an expression such as

[TABLE]

exists and is equal to $1$.

III.2 MPO representation of $\rho_{\varphi}$ and $\rho^{\prime}_{\varphi}$

In this subsection we exploit the tensor-network representation to write $\rho_{\varphi}$ and its derivative $\rho^{\prime}_{\varphi}$ as MPOs. The initial assumption we make to obtain our representation is that $\rho_{0}$ admits an efficient representation, with some finite bond dimension $D_{\rho_{0}}$, as an MPO:

[TABLE]

where for brevity we have introduced $\mathbf{j}=\{j_{1},\dots,j_{N}\}$. Exploiting channel/state duality to vectorize $\rho_{0}$ we obtain the MPS representation:

|\rho_{0}\rangle=$$\dots$$R_{0}[1]$$R_{0}[2]$$R_{0}[3]$$R_{0}[4]$$R_{0}[5]$$R_{0}[N]

.

Our first goal is to find the tensor network representation of $\rho_{\varphi}=\Lambda_{\varphi}(\rho_{0})$. We achieve this in two steps. First we exploit the local structure of the noise channel, Eq. (13), to build a tensor network for $\Lambda|\rho_{0}\rangle$—the result of this construction is not, in general, in MPS form. Then we exploit the singular value decomposition to put the resulting tensor network back into MPS form. Finally, we apply the local unitary parameter imprinting and return to the original MPO form.

For definiteness, we focus on the situation when the $\Lambda$ operator can be described by a subsequent action of singe-particle $\Lambda^{[n]}$ and two-particle terms $\Lambda^{[n,n+1]}$—which in the following are denoting by $Y$ and $X$, respectively—while retaining translation invariance. 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 $\Lambda$ operator therefore takes the following form:

\Lambda|\rho_{0}\rangle=$$\dots$$R_{0}[1]$$R_{0}[2]$$R_{0}[3]$$R_{0}[4]$$R_{0}[5]$$R_{0}[N]$$X$$X$$\dots$$X$$X$$Y$$Y$$Y$$Y$$Y$$Y

,

where we place $X$ operators in a skewed orientation in order to maintain a manifestly translation invariant model.

This state is no longer in MPS form. The operator $X$ is defined with respect to a product basis $|\alpha\rangle=|j,k\rangle$ for the doubled legs and it acts on two neighbouring subsystems $|\alpha\rangle$ and $|\beta\rangle$, i.e., it is the tensor

[TABLE]

By vectorizing each of the subsystems we can express it as a matrix:

[TABLE]

Applying the SVD to this matrix gives us the representation

[TABLE]

where $S_{\gamma,\gamma^{\prime}}=s_{\gamma}\delta_{\gamma,\gamma^{\prime}}$ with singular values $s_{\gamma}$ and $U$ and $V$ are unitary operators. Graphically this becomes:

=$$\alpha$$\alpha^{\prime}$$\beta$$\beta^{\prime}$$\alpha$$\alpha^{\prime}$$\beta$$\beta^{\prime}$$U$$S$$V^{\dagger}$$X

.

(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 $\sqrt{S}$ into the $U$ tensor from the right (respectively, into the $V^{\dagger}$ tensor from the left). Let $D^{(2)}$ (the upper index indicates two-particle nature of the noise) be the number of non-zero (or more practically non-negligible) singular values $s_{\gamma}$. Introducing $T=U\sqrt{S}$, $W=\sqrt{S}V^{\dagger}$ results in the following tensor network,

\dots$$T$$W$$T$$W$$T$$W$$T$$W$$T$$W$$T$$W

,

in place of the layer of $X$s.

The final step to obtain an MPO representation for $\rho_{\varphi}$ is to combine the tensors $R_{0}[n]$, $T$, $W$, $Y$ and a tensor $Z=e^{-\mathrm{i}h\varphi}\otimes(e^{\mathrm{i}h\varphi})^{\mathrm{T}}$—which represents the unitary phase encoding process—into a single new MPS tensor $R_{\varphi}[n]$:

R_{\varphi}[n]$$=$$R_{0}[n]$$W$$T$$Y$$Z

.

The doubled horizontal legs can be then combined into thicker horizontal legs to yield the MPS representation:

|\rho_{\varphi}\rangle=$$\dots$$R_{\varphi}[1]$$R_{\varphi}[2]$$R_{\varphi}[3]$$R_{\varphi}[4]$$R_{\varphi}[5]$$R_{\varphi}[N]

.

This representation can be put into MPO form by splitting the vertical legs back into the original two legs and bending:

\rho_{\varphi}=$$\dots$$R_{\varphi}[1]$$R_{\varphi}[2]$$R_{\varphi}[3]$$R_{\varphi}[4]$$R_{\varphi}[5]$$R_{\varphi}[N]

.

As a result we end up with an MPO with the bond dimension $D_{\rho}=D_{\rho_{0}}D^{(2)}$. The generalization of this derivation beyond the case of nearest-neighbour correlations will lead to an MPO representation of the density matrix $\rho_{\varphi}$ with bond dimension $D_{\rho}=D_{\rho_{0}}D_{r}$. Here $D_{r}=\prod_{k=2}^{r}D^{(k)}$ represents the contribution to the effective bond dimension of the output state resulting from the action of the correlated noise, where $D^{(k)}$ is the number of non-zero singular values that will appear when considering the $k$-particle noise term (the upper bound on $D^{(k)}$ is $d^{2(k-1)}$).

We finally move on to the task of writing the $\rho^{\prime}_{\varphi}$ operator as an MPO. It is possible to do this efficiently thanks to the fact that this operator is given as a commutator of $\rho_{\varphi}$ with $H$, where $H$ is a sum of local Hamiltonians

[TABLE]

As a result we can regard $\rho^{\prime}_{\varphi}$ as a sum of $N$ MPOs where each of them represents the original $\rho_{\varphi}$ MPO modified by an action of the Hamiltonian $h$ on consecutive particles. In what follows we assume that the basis $|j\rangle$ associated with the physical indices $j$ is chosen to be the eigenbasis of the local Hamiltonian $h$, $h=\sum_{j}\epsilon_{j}|j\rangle\langle j|$, where $\epsilon_{j}$ are the corresponding eigenvalues. With this choice of local basis, the MPO representation of $\rho^{\prime}_{\varphi}$ can be easily written at the cost of doubling the bond dimension ($D_{\rho^{\prime}}=2D_{\rho_{0}}D_{r}$):

[TABLE]

where $R^{\prime}[n]^{j_{n}}_{k_{n}}$ is equal to

[TABLE]

The $2\times 2$ matrices that appear in the above construction are responsible for the increase of the bond dimension but guarantee that the effect of trace in (26) is equivalent to that resulting from the sum of $\rho_{\varphi}$ MPO acted upon consecutively by the commutator of the local Hamiltonians corresponding to different particles.

III.3 Optimization of the QFI

As indicated in Sec. II.1, maximization of our FoM $F(\rho_{0},L)$ leading to the maximal possible QFI for a given metrological model is a two-step iterative process. In the first part of this subsection we show how to maximize the FoM over a Hermitian operator $L$ with fixed $\rho_{0}$, and in the second part we focus on the maximization over the input state $\rho_{0}$ with fixed $L$.

We search for the optimal $L$ in the form of an MPO

[TABLE]

with a finite bond dimension $D_{L}$. Without loss of generality, each $S[n]^{j_{n}}_{k_{n}}$ is assumed Hermitian in its physical indices,

[TABLE]

to ensure that $L$ is Hermitian. We call this condition a Hermitian gauge.

The bond dimension $D_{L}$ for $L$ is expected to be small for weakly correlated noise models. Indeed, in the limit of an uncorrelated product state $\rho_{\varphi}=\varrho_{\varphi}^{\otimes N}$, $L$ is a sum of local operators, $L=\sum_{n=1}^{N}l^{[n]}$. Here $l^{[n]}$ is the SLD for the single-particle problem applied to particle $n$. In a similar way as for $\rho_{\varphi}^{\prime}$ presented above, the sum can be represented by an MPO with a bond dimension $2$. Therefore, $D_{L}=2$ is the limiting value for uncorrelated noise models.

The FoM to be maximized is

2$$-$$\dots$$\dots$$R^{\prime}[1]$$S[1]$$R^{\prime}[2]$$S[2]$$R^{\prime}[3]$$S[3]$$R^{\prime}[N]$$S[N]$$\dots$$\dots$$\dots$$S[1]$$R_{\varphi}[1]$$S[1]$$S[2]$$R_{\varphi}[2]$$S[2]$$S[3]$$R_{\varphi}[3]$$S[3]$$S[N]$$R_{\varphi}[N]$$S[N]

.

Maximization of the FoM over $L$ is equivalent to a joint maximization over each tensor $S[n]$. We relax this optimization problem by iterating an optimization loop. In the loop we first find the optimal $S[1]$, then $S[2]$, and so on up to $S[N]$, after which we go back to $S[1]$. The optimization over each $S[n]$ 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 $S[n]$ and the optimal $S[n]$ is found as a solution to a linear equation. For definiteness, to explain the procedure, we focus on the generic example of $S[2]$. After vectorizing $S[2]\rightarrow|S[2]\rangle$,

S[2]$$\gamma$$\delta$$k_{2}$$j_{2}$$S[2]$$S[2]$$\alpha

,

we can represent the $S[2]$-FoM as

|b\rangle$$2$$-$$\dots$$\dots$$R^{\prime}[1]$$R^{\prime}[2]$$R^{\prime}[3]$$R^{\prime}[N]$$S[1]$$S[3]$$S[N]$$S[2]$$\alpha$$A$$\dots$$\dots$$\dots$$S[1]$$R_{\varphi}[1]$$S[1]$$S[3]$$R_{\varphi}[3]$$S[3]$$S[N]$$R_{\varphi}[N]$$S[N]$$S[2]$$R_{\varphi}[2]$$S[2]$$\beta$$\alpha

,

which can be written in a compact way as

[TABLE]

Here $b_{\alpha}$ are the elements of the vector $|b\rangle$, and $A_{\alpha,\beta}$ are the elements of the matrix $A$. Both $|b\rangle$ and $A$ describe the entire tensor network complementing the distinguished vector $|S[2]\rangle$ in the two respective terms of the $S[2]$-FoM. After taking a derivative with respect to $S[2]_{\alpha}$, we obtain a linear equation for the extremum:

[TABLE]

The $d^{2}D_{L}^{2}\times d^{2}D_{L}^{2}$ matrix $\widetilde{A}=\frac{1}{2}\left(A+A^{\mathrm{T}}\right)$ typically has a non-zero kernel and the linear equation does not have a unique solution. We use the Moore-Penrose pseudo-inverse, $\widetilde{A}^{+}$, to obtain a solution $|S[2]\rangle=\widetilde{A}^{+}|b\rangle$ that does not contain any zero modes of $\widetilde{A}$.

If the linear equation was non-singular, then its exact solution would satisfy the Hermitian gauge (29). For the typical singular case, using an SVD of $\widetilde{A}$ to construct its pseudo-inverse, we have to truncate singular values falling below a small but finite cut-off, set by $\kappa$ multiplied by the highest singular value. As the cut-off solution $|S[2]\rangle$ need not satisfy the Hermitian gauge condition exactly, we filter out its small anti-Hermitian part with the substitution:

[TABLE]

From experience, this substitution can improve numerical stability but is not necessary when all initial $S[n]$ are in the Hermitian gauge (29) and $\kappa$ is large enough to suppress the anti-Hermitian part of the solution. However, with too large a cut-off the final optimized $L$ does not reach the maximal possible value of the QFI. Therefore, we adjust $\kappa$ to obtain the highest QFI achievable without compromising the stability.

Now we move on to the maximization of the FoM over the input state $\rho_{0}$ for a fixed $L$. We start by rewriting $F(\rho_{0},L)$ as

[TABLE]

where by $\Lambda^{*}_{\varphi}(\cdot)$ we denote the channel which is dual to $\Lambda_{\varphi}(\cdot)$ (the evolution written in the Heisenberg picture). We can rewrite this as

[TABLE]

where we introduce $L^{*}_{\varphi}=\Lambda^{*}_{\varphi}(L)$, $L^{\prime*}_{\varphi}=\frac{\mathrm{d}L^{*}_{\varphi}}{\mathrm{d}\varphi}=\mathrm{i}\left[H,L^{*}_{\varphi}\right]$ and $L^{*}_{2,\varphi}=\Lambda^{*}_{\varphi}(L^{2})$. By analogy with the construction of the MPO representation for $\rho_{\varphi}=\Lambda_{\varphi}(\rho_{0})$ and $\rho^{\prime}_{\varphi}=\mathrm{i}\left[\rho_{\varphi},H\right]$ in Sec. III.2, we can easily construct the MPO representation of $L^{*}_{2,\varphi}$ and $L^{\prime*}_{\varphi}$ from the known MPO form of $L$. The tensors determining the MPO form of $L^{*}_{2,\varphi}$ and $L^{\prime*}_{\varphi}$ are denoted by $S_{2}[n]$ and $S^{\prime}[n]$, respectively, and their respective bond dimensions are $D_{L_{2}}=D_{L}^{2}D_{r}$ and $D_{L^{\prime}}=2D_{L}D_{r}$.

The quantity $F(\rho_{0},L)$ in (34) is maximal when $\rho_{0}$ is a projection on the eigenvector associated with the maximal eigenvalue of the Hermitian operator $2L^{\prime*}_{\varphi}-L^{*}_{2,\varphi}$. Hence, without loss of generality, we can assume a pure input state $\rho_{0}=|\psi\rangle\!\langle\psi|$ with $|\psi\rangle$ being an MPS with bond dimension $D_{\psi}$:

[TABLE]

The input state $\rho_{0}$ has bond dimension $D_{\rho_{0}}=D_{\psi}^{2}$ and its MPO tensors are $R_{0}[n]^{j_{n}}_{k_{n}}=P[n]^{j_{n}}\otimes\overline{P[n]^{k_{n}}}$.

The maximization of $F(\rho_{0},L)$ over the input state $|\psi\rangle$ is equivalent to the variational optimization for the ground state of a many-body “Hamiltonian” $L^{*}_{2,\varphi}-2L^{\prime*}_{\varphi}$, a problem widely discussed in the many-body physics MPS literature White (1992); Schollwöck (2011); Bridgeman and Chubb (2017). After reinterpreting our problem as a variational minimization of the “energy”

[TABLE]

we proceed iteratively in a similar way as in the case of the maximization of $F(\rho_{0},L)$ over $L$.

For example, in order to find the minimum over $P[2]$, we begin by vectorizing the tensor $P[2]\rightarrow|P[2]\rangle$ and expressing the “energy” (36) as a Rayleigh quotient

[TABLE]

Here the $dD_{\psi}^{2}\times dD_{\psi}^{2}$ matrices are

\alpha$$\dots$$\beta$$\dots$$\dots$$P[1]$$S_{2}[1]$$\overline{P[1]}$$P[3]$$S_{2}[3]$$\overline{P[3]}$$P[N]$$S_{2}[N]$$\overline{P[N]}$$S_{2}[2]$$\mathcal{F}_{\alpha\beta}=$$-$$\alpha$$\dots$$\beta$$\dots$$\dots$$P[1]$$S^{\prime}[1]$$\overline{P[1]}$$P[3]$$S^{\prime}[3]$$\overline{P[3]}$$P[N]$$S^{\prime}[N]$$\overline{P[N]}$$S^{\prime}[2]$$2

and

\alpha$$\dots$$\beta$$\dots$$P[1]$$\overline{P[1]}$$P[3]$$\overline{P[3]}$$P[N]$$\overline{P[N]}$$\mathcal{N}_{\alpha\beta}=

.

After taking the derivative of (37) we obtain the condition for the extremum:

[TABLE]

which is a generalized eigenvalue problem with eigenvalue $-F(\rho_{0},L)$. By multiplying it with a pseudo-inverse of matrix $\mathcal{N}$ 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 $|\psi\rangle$ can be brought to a canonical form, where the matrix $\mathcal{N}$ becomes an identity and Eq. (38) reduces to a standard eigenvalue problem. There is no need to pseudo-invert $\mathcal{N}$.

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

While running the algorithm, one has to choose the bond-dimensions for the input state, $D_{\psi}$, as well as for the SLD, $D_{L}$, 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, $D_{\psi}=1$, and optimized for $L$ with a minimal non-trivial $D_{L}=2$. Then we increased one of the bond dimensions, either $D_{\psi}$ or $D_{L}$, 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., $1\%$ and hence assumed that relative error of our method is around $1\%$.

III.4 Asymptotic limit

The previous subsection describes 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. Cincio and Vidal (2013). 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 $E^{N}$ 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, $E^{k}$ 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 $|\psi\rangle$ is TI then the final state $\rho_{\varphi}$ is TI as well. There is a problem, however, with the operator $\rho^{\prime}_{\varphi}$. Its construction as an MPO in Eq. (27) is not TI. Because of this, instead of calculating the derivative exactly, we approximate it by a difference of two TI iMPOs:

[TABLE]

with infinitesimal parameter $\varepsilon$. Motivated by the defining equation for the SLD (6), we can consider a similar expansion of the operator $L$:

[TABLE]

Here $\widetilde{L}$ and $\openone$ are solutions of Eq. (6) when $\rho^{\prime}_{\varphi}$ is replaced by, respectively, $\rho_{\varphi+\varepsilon}$ and $\rho_{\varphi}$. We search for the optimal $\widetilde{L}$ which is better suited for the TI formalism than the original operator $L$. Let us denote by

[TABLE]

the QFI per particle that we want to maximize:

[TABLE]

A TI iMPO is defined by only one tensor which we assign respectively as: $|\psi\rangle\rightarrow P$, $\rho_{0}\rightarrow R_{0}$, $\rho_{\varphi}\rightarrow R_{\varphi}$, $\widetilde{L}\rightarrow\widetilde{S}$, $\widetilde{L}^{*}_{\varphi}\rightarrow\widetilde{S}_{\varphi}$, $\widetilde{L}^{*}_{2,\varphi}\rightarrow\widetilde{S}_{2}$.

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 Corboz (2016). The main idea is to find the optimal tensor $A_{\mathrm{new}}$ at one site, treating all other tensors $A$ as fixed, and then rather then replacing all tensors by $A_{\mathrm{new}}$ perform a flexible substitution,

[TABLE]

with a mixing angle $\lambda$. 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 $\widetilde{L}$ which is equivalent to the determination of the optimal local tensor $\widetilde{S}$.

As explained in Sec. III.1, the trace of an operator represented as an iMPO is defined by its transfer matrix, so we start by introducing transfer matrices $E_{1}$ and $E_{2}$ associated with, respectively, $\rho_{\varphi+\varepsilon}\widetilde{L}$ and $\rho_{\varphi}\widetilde{L}^{2}$:

\widetilde{S}$$R_{\varphi+\varepsilon}$$E_{1}$$\displaystyle=\sum_{j}\lambda_{1,j}$$r_{1,j}$$l_{1,j}

\displaystyle\sum_{j}\lambda_{1,j}$$r_{1,j}$$l_{1,j}

,

\widetilde{S}$$R_{\varphi}$$\widetilde{S}$$E_{2}$$\displaystyle=\sum_{j}\lambda_{2,j}$$r_{2,j}$$l_{2,j}

\displaystyle\sum_{j}\lambda_{2,j}$$r_{2,j}$$l_{2,j}

,

which allow us to write $\mathrm{Tr}\left(\rho_{\varphi+\varepsilon}\widetilde{L}\right)=\mathrm{Tr}E_{1}^{N}$ and $\mathrm{Tr}\left(\rho_{\varphi}\widetilde{L}^{2}\right)=\mathrm{Tr}E_{2}^{N}$. Now using the fact that $E_{i}^{N}$ is determined by its leading eigenvalue (see Eq. (19)) we can write $\mathrm{Tr}E_{i}^{N}=\lambda_{i,1}^{N-1}(l_{i,1}|E_{i}|r_{i,1})$ and express $f(\rho_{0},\widetilde{L})$ as:

\widetilde{S}$$R_{\varphi+\varepsilon}$$r_{1,1}$$l_{1,1}$$\displaystyle\frac{2\lambda_{1,1}^{N-1}}{N\varepsilon^{2}}$$\widetilde{S}$$R_{\varphi}$$\widetilde{S}$$r_{2,1}$$l_{2,1}$$-\displaystyle\frac{\lambda_{2,1}^{N-1}}{N\varepsilon^{2}}$$-\displaystyle\frac{1}{N\varepsilon^{2}}

,

which after vectorizing $\widetilde{S}$ is equal to

\frac{2\lambda_{1,1}^{N-1}}{N\varepsilon^{2}}$$\widetilde{S}$$R_{\varphi+\varepsilon}$$l_{1,1}$$r_{1,1}$$|b\rangle$$\alpha$$-\frac{\lambda_{2,1}^{N-1}}{N\varepsilon^{2}}$$-\frac{1}{N\varepsilon^{2}}$$\widetilde{S}$$R_{\varphi}$$\widetilde{S}$$l_{2,1}$$r_{2,1}$$\beta$$\alpha$$A

,

or when written in as equation:

[TABLE]

The condition for an extremum reads:

[TABLE]

For $N\to\infty$ the powers of the eigenvalues may seem to pose a problem. Fortunately, however, this problem can be circumvented. For a given $\widetilde{L}$ we can calculate the associated value of our FoM per particle:

[TABLE]

but going back to the roots (Eq. (8)) we see that FoM per particle should have form $f(\rho_{0},\widetilde{L})=2f_{1}-f_{2}$ where $f_{1}$ and $f_{2}$ 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 $N\rightarrow\infty$, $\varepsilon\rightarrow 0$, and $N\varepsilon^{2}\rightarrow 0$, and remembering binomial expansion

[TABLE]

it is to be expected that the highest eigenvalues of the transfer matrices have the form:

[TABLE]

which after inserting into Eq. (46) and using binomial expansion to the first order give us exactly $f(\rho_{0},\widetilde{L})=2f_{1}-f_{2}$. Note that, it means that we can calculate value of FoM per particle in a simple way:

[TABLE]

It is clear now that for the purpose of solving Eq. (45) we can approximate $\lambda_{1,1}^{N-1}$ and $\lambda_{2,1}^{N-1}$ by ones and, hence, bring the condition for the optimal $\widetilde{S}$ to a simpler form:

[TABLE]

This equation is solved with a pseudo-inverse and its anti-Hermitian part is filtered out at every iteration step.

The SLD is always traceless in any unitary parameter estimation problem—see Eq. (7). Note also that $\langle\lambda|\rho^{\prime}_{\varphi}|\lambda\rangle=\mathrm{i}\langle\lambda|[H,\rho_{\varphi}]|\lambda\rangle=0$ for any $|\lambda\rangle$ which is an eigenstate of $\rho_{\varphi}$. We can ensure that solution $\widetilde{S}$ has proper normalization using the condition $\mathrm{Tr}L=0$ or, equivalently $\mathrm{Tr}\widetilde{L}=\mathrm{Tr}\mathds{1}=d^{N}$ which in the language of transfer matrices means that the highest eigenvalue of the transfer matrix,

$\widetilde{S}$

,

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 $\rho^{\prime}_{\varphi}$, we approximate the exact derivative of $L^{\prime*}_{\varphi}$ by its discrete version:

[TABLE]

After the expansion $L=(\widetilde{L}-\mathds{1})/\varepsilon$, our task becomes equivalent to minimization of the “energy density”:

[TABLE]

over $|\psi\rangle$. Using transfer matrices $E_{3}$, $E_{4}$ and $E_{5}$ associated with, respectively, $\langle\psi|\widetilde{L}^{*}_{2,\varphi}|\psi\rangle$, $\langle\psi|\widetilde{L}^{*}_{\varphi+\varepsilon}|\psi\rangle$ and $\langle\psi|\psi\rangle$,

\overline{P}$$\widetilde{S}_{2}$${P}$$=\displaystyle\sum_{j}\lambda_{3,j}$$r_{3,j}$$l_{3,j}

\overline{P}$$\widetilde{S}_{\varphi+\varepsilon}$${P}$$=\displaystyle\sum_{j}\lambda_{4,j}$$r_{4,j}$$l_{4,j}

\overline{P}$${P}$$=\displaystyle\sum_{j}\lambda_{5,j}$$r_{5,j}$$l_{5,j}

we can rewrite Eq. (52) in a diagrammatic form:

-f(\rho_{0},\widetilde{L})=$$\overline{P}$$P$$r_{5,1}$$l_{5,1}$$\displaystyle N\varepsilon^{2}$$\overline{P}$$P$$r_{5,1}$$l_{5,1}$$+\lambda_{5,1}^{N-1}$$\overline{P}$$\widetilde{S}_{\varphi+\varepsilon}$$P$$r_{4,1}$$l_{4,1}$$-2\lambda_{4,1}^{N-1}$$\overline{P}$$\widetilde{S}_{2}$$P$$r_{3,1}$$l_{3,1}$$\lambda_{3,1}^{N-1}

.

As previously, we expect that $\lambda_{i,1}=1+\varepsilon^{2}f_{i}$ and for the purpose of finding the optimal tensor $P$, we can approximate $\lambda_{3,1}^{N-1}$ and $\lambda_{4,1}^{N-1}$ by ones. After taking the derivative we obtain the condition for the extremum:

[TABLE]

where $g=-f(\rho_{0},\widetilde{L})N\varepsilon^{2}-\lambda_{5,1}^{N-1}$ is a generalised eigenvalue, whereas the matrices $\mathcal{F}$ and $\mathcal{N}$ are defined as:

\widetilde{S}_{2}$$r_{3,1}$$l_{3,1}$$\alpha$$\beta$$\mathcal{F}_{\alpha\beta}=$$\widetilde{S}_{\varphi+\varepsilon}$$r_{4,1}$$l_{4,1}$$\alpha$$\beta$$-2

,

r_{5,1}$$l_{5,1}$$\alpha$$\beta$$\mathcal{N}_{\alpha\beta}=$$r_{5,1}$$l_{5,1}$$\alpha$$\beta$$=$$\otimes$$\otimes$$\mathds{1}

.

The matrix $\mathcal{N}$ is a tensor product of three matrices, $r_{5,1}$, the identity, and $l_{5,1}$, hence its pseudo-inverse $\mathcal{N}^{+}$ can be obtained as a tensor product of (pseudo-)inverses of the smaller matrices. Applying $\mathcal{N}^{+}$ to Eq. (53) 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 $|\psi\rangle$ and require that $|\psi\rangle$ is normalized so that $\lambda_{5,1}=1$. Then we calculate the asymptotic value of the QFI per particle:

[TABLE]

Just as for the finite $N$, iterating the $\widetilde{L}$ and $|\psi\rangle$ optimization steps leads to the optimal solution with the maximal QFI per particle.

While performing the numerics one should choose $\varepsilon$ to be small but not too small as too small values may lead to numerical instabilities. Our general strategy in obtaining numerical results reported in the next section, was to lower the value of $\varepsilon$ 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 $\varepsilon\approx 10^{-3}-10^{-4}$ (instabilities started to appear for $\varepsilon<10^{-6}$). Notice that in the asymptotic iMPO approach described above we required $N\varepsilon^{2}$ to be small—on the order of the precision we expect from the numerical results. In other words setting the precision requirements to $10^{-2}$ this implies that $N\varepsilon^{2}\approx 10^{-2}$ and hence $N\approx 10^{4}-10^{6}$. What this physically means is that in our setup the QFI per particle does not change in any noticeable way for larger $N$ and hence the asymptotic behaviour in the actual $N\rightarrow\infty$ limit may be inferred from this results.

IV Applications

In this section we present three applications of our framework. The examples were 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. The first example – of magnetic field sensing with locally correlated magnetic field fluctuations from Sec IV.1 – falls into the general metrological model structure as outlined in Sec. II, and should be regarded as a typical representative of the models that can be dealt with efficiently using the MPO framework. In this case it is the standard QFI FoM that is being employed in the optimization process. However, in order to demonstrate that the utility of MPO based approach to metrology goes beyond the QFI optimization tasks, in Sec. IV.2 we show how the framework can be adapted to calculate the fundamental bound on the achievable Allan variance in the atomic clock stabilization problem. In this case not only the physical setup is different (temporal LO noise correlations affecting the atoms, rather than spatial noise correlations between the probes), but also the employed FoM. Instead of the QFI the QAVAR is optimized—a quantity in fact more closely related with the Bayesian variance rather than the QFI. Finally, in Sec. IV.3, the third example demonstrates that the utility of the techniques of calculating QFI for mixed state via the MPO formalism developed in this paper, goes beyond the realm of metrological applications and can be directly applied to exact calculation of fidelity susceptibility in many-body thermal states—a task deemed too hard for all the state-of-the art methods.

IV.1 Magnetic field sensing with locally correlated noise

Consider $N$ particles each with spin $s$ (in principle this might be an effective spin of some number of “sub-particles”) 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. We will go beyond the standard model of independent field fluctuations resulting in independent dephasing of atoms Huelga et al. (1997); Escher et al. (2011), by taking into account correlations between field fluctuations at the nearest neighbour particle sites. This will lead us to a model where we will be able to study the impact of correlations (and anti-correlations) in the effective dephasing process on the metrological potential of the system.

For a moment let us assume that the field $B$ is fixed. The Hamiltonian corresponding to the dynamics of an $n$-th spin in a static magnetic field is $-gS^{[n]}_{z}B$, where $S^{[n]}_{z}$ is the $z$ spin component of the $n$-th particle and $g$ is the gyromagnetic ratio of the particle. In order to stay consistent with the abstract notation introduced in Sec. II.2, we will identify $\varphi=gBt$, $h^{[n]}=S_{z}^{[n]}/\hbar$. The uncertainty of estimation of $B$ will be related with the uncertainty of a standard phase estimation problem via a simple proportionality relation $\Delta\widetilde{B}=\Delta\widetilde{\varphi}/(gt)$.

Now, taking into account the presence of fluctuations let us write the magnetic field at site $n$ as $B^{[n]}(t)=B+\delta B^{[n]}(t)$. Here we assume that fluctuations are Gaussian and have no relevant temporal correlations (white noise). The corresponding variance as well as the nearest-neighbour correlation functions of the fluctuating field read: $\langle\delta B^{[n]}(t)\delta B^{[n]}(t^{\prime})\rangle=\sigma^{2}\delta(t-t^{\prime})$, $\langle\delta B^{[n]}(t)\delta B^{[n+1]}(t^{\prime})\rangle=\chi\delta(t-t^{\prime})$, where $\chi$ 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.

Let $|\mathbf{j}\rangle=|j_{1},j_{2},\dots,j_{N}\rangle$, $j_{n}\in\{-s,\dots,s\}$ be the eigenbasis with well defined eigenvalues of local Hamiltonians $h^{[n]}=S^{[n]}_{z}/\hbar$ operators equal to $j_{n}$: $h^{[n]}|\mathbf{j}\rangle=j_{n}|\mathbf{j}\rangle$. Using this basis we can easily write the evolution of the density matrix under the action of noise by applying the standard technique of the cumulant expansion van Kampen (1981):

[TABLE]

where $\mathbf{j},\mathbf{k}$ are column vectors $\mathbf{j}=(j_{1},j_{2},\dots,j_{N})^{\mathrm{T}}$, and $C$ is the correlation matrix

[TABLE]

where $c_{1}=\sigma^{2}g^{2}t$, $c_{2}=\chi g^{2}t$. It is straightforward to extend the discussion here to deal with the more general longer range correlations covering up to $r$ neighbouring particles, in which case the matrix $C$ will be $2r+1$ diagonal.

For a better physical insight, let us provide the corresponding master equation in the form of Eq. 11. One may obtain it by simply differentiating (55) over $t$. As a result we get the master equation with single particle $L^{[n]}=\sqrt{\gamma_{1}}h^{[n]}$ and two particle $L^{[n,n+1]}=\sqrt{|\gamma_{2}|}(h^{[n]}+\mathrm{sgn}(\gamma_{2})h^{[n+1]})$ dephasing operators. The dephasing rates $\gamma_{1}$, $\gamma_{2}$ are related with field fluctuation properties as follows: $\gamma_{1}=(\sigma^{2}-2|\chi|)g^{2}$, $\gamma_{2}=\chi g^{2}$—note that $\sigma^{2}\geq 2|\chi|$ by virtue of positivity of the correlation matrix $C$ so the rate $\gamma_{1}$ is always positive.

In order to write the evolution manifestly in the MPO formalism, we will replace $|\mathbf{j}\rangle\langle\mathbf{k}|\rightarrow|\mathbf{j}\rangle|\mathbf{k}\rangle$ which forms a basis for the vectorized input density matrix $|\rho_{0}\rangle$. The action of $\Lambda$ on $\rho_{0}$ is identical to the action of the operator $e^{\Gamma}$ on $|\rho_{0}\rangle$, i.e., $|\Lambda(\rho_{0})\rangle=e^{\Gamma}|\rho_{0}\rangle$, where

[TABLE]

with

[TABLE]

and

[TABLE]

Note that the $\Upsilon^{[n]}$ and $\Xi^{[n,n+1]}$ mutually commute with each other, so that

[TABLE]

Denoting $Y^{[n]}=e^{-\frac{c_{1}}{2}\Upsilon^{[n]}}$ and $X^{[n,n+1]}=e^{-c_{2}\Xi^{[n,n+1]}}$ we finally obtain

[TABLE]

which is the form of evolution the same as discussed in Sec. III.2 guaranteeing efficient MPO description.

After evolution through quantum channel $\Lambda$, the phase is imprinted in our state through unitary evolution according to Eq. (10) with local Hamiltonians $h^{[n]}$—in the notation of Sec. III this is represented by the action of $\prod_{n=1}^{N}Z^{[n]}$, where $Z^{[n]}=e^{-\mathrm{i}h^{[n]}\varphi}\otimes(e^{\mathrm{i}h^{[n]}\varphi})^{\mathrm{T}}$. Written in the basis $|\mathbf{j}\rangle$, $\rho^{\prime}_{\varphi}=\mathrm{i}\left[\rho_{\varphi},H\right]$ reads:

[TABLE]

Now we are ready to make use of the methods described in Sec. III in order to find the optimal probe states and the corresponding maximal QFI. All the numerical results that follow correspond to the $N$ spin $1/2$ particle case, $h^{[n]}=\sigma_{z}^{[n]}/2$. We perform numerical calculations using MATLAB software with the help of ncon() function Pfeifer et al. (2014) for tensor contractions. We should stress that all the optimization algorithms are written from scratch with the metrological context in mind, since standard optimization procedures utilized by many-body physics community are not easily adapted to the optimization tasks that we face.

First, in Fig. 1(i) we present a comparison of results of the QFI optimization procedure for exemplary dephasing parameters $c_{1}=1$, $c_{2}=0.1$ (correlated noise) obtained using the finite number of particles $N$ MPO approach and the asymptotic value of the QFI per particle obtained using the iMPO approach—note that we plot $F/N$ as we expect asymptotic linear scaling of the QFI and hence convergence of $F/N$ to a fixed value. The results obtained via the two approaches are in very good agreement. When we plot QFI versus $N$ and fit to the data for large $N$ a straight line we obtain a slope of $0.500$. It compares well with the asymptotic limit obtained directly from the iMPO method as $0.504$. This is a numerical confirmation that indeed the iMPO approach which, as described in Sec. III.4 is much more conceptually involved, yields correct results. This is a highly relevant observation as it is numerically much more efficient to obtain asymptotic properties of the QFI using directly the iMPO approach rather than performing finite $N$ computations and extrapolating them to $N\rightarrow\infty$. While performing the optimization we have set the error tolerance on the level of $1\%$, which resulted in the maximal bond dimension required in the numerical procedure to be only $D_{\psi}=4$ for the input state and $D_{L}=2$ for the SLD. This demonstrates how efficient the MPO description is in this case. Crosses overlaid on the plot indicate the regime where direct calculations using the standard full Hilbert space description were possible on the same high performance PC on which the MPO algorithms were run. This regime is clearly very far from the one were we observe the convergence to the asymptotic linear scaling of the QFI with $N$, and is only accessible numerically using the MPO based methods. In Fig. 1(ii) we present the contour plot depicting the asymptotic value of the QFI per particle as a function of noise parameters and contrast with the achievable QFI when using product states, see Fig. 1(iii). In order to appreciate the amount of entanglement that is present in the optimal states, in Fig. 1(iv) we provide a contour plot of the results of calculation of the von Neuman entropy for the iMPO Cirac et al. (2011) corresponding to the reduced density matrix of the system when half of the particles is traced out (optimal QFI can be achieved by many different states so some fluctuations of entropy are to be expected). We see a clear relation, between the amount of entanglement and the increase in sensing precision.

Let us now ask the question, whether some insight into the problem could have been gained by ingeniously adapting the state-of-the-art methods of deriving fundamental bounds in quantum metrology developed with uncorrelated noise models in mind Escher et al. (2011); Demkowicz-Dobrzański et al. (2012); Demkowicz-Dobrzański and Maccone (2014); Demkowicz-Dobrzański et al. (2017); Zhou et al. (2018). These methods provide easily calculable asymptotic bounds, based on just the knowledge of the Kraus operators of elementary probe dynamics or alternatively noise jump operators appearing in the quantum master equation. Even though the noise in our problem is correlated we can formally divide the evolution as a collection of independent channels acting on two, three or more particles—a similar trick has been employed in a recent study of the impact of many-body effects in atomic interferometry Czajkowski et al. (2018). In the scheme below we indicate the possible constructions. First we formally decompose local dephasing gates $Y$ as products $Y=Y^{\prime}Y^{\prime\prime}$ of local dephasing gates with corresponding $c_{1}^{\prime}$ and $c_{1}^{\prime\prime}$ such that $c_{1}=c_{1}^{\prime}+c_{1}^{\prime\prime}$, and the unitary encoding gate $Z$ as a product $Z=Z^{\prime}Z^{\prime\prime}$, where $Z^{\prime}$, $Z^{\prime\prime}$ are phase gates with corresponding phases $\varphi^{\prime}$, $\varphi^{\prime\prime}$, such that $\varphi=\varphi^{\prime}+\varphi^{\prime\prime}$. We can now unravel the total dynamics as effectively composed of $N$ independent channels $\Lambda_{1}$, or group the gates into $N/2$ larger channels $\Lambda_{2}$, or $N/3$ channel $\Lambda_{3}$, etc.—see Fig. 2.

We may now apply the bounds derived for uncorrelated noise models Demkowicz-Dobrzański et al. (2012); Demkowicz-Dobrzański and Maccone (2014) using either decomposition of the dynamics into $\Lambda_{1}$, $\Lambda_{2}$ or $\Lambda_{3}$ channels—we can group the operations in any way we please as in this case all elementary evolutions commute. In order to obtain the tightest bound we numerically optimize the split of phases as well as noise contributions between gates $Y^{\prime}$, $Z^{\prime}$ and $Y^{\prime\prime}$, $Z^{\prime\prime}$, while making sure that the resulting $\Lambda_{i}$ is a legitimate quantum channel—note that the bare two qubit gate $X$ is not a proper quantum channel as it is not completely positive. The results obtained are depicted in the left inset of Fig. 1(ii). While the bounds are tight for the decorrelated noise model they are far from the actual achievable QFI in the correlated (or anti-correlated) noise regimes and as expected they improve when we increase the elementary channel size. Still, because this method scales badly with the elementary channel size we were forced to stop with $\Lambda_{3}$ channel. This 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.

As mentioned above, the use of the iMPO approach can greatly speed up process of calculating asymptotic value of QFI. In order to study the impact of noise correlations on the achievable QFI and the optimal states, for the rest of this section, we will therefore restrict ourselves to numerical results obtained using this approach. Using iMPO we have studied the optimal QFI for all noise parameters in the range $c_{1}\in[0.2,2]$, and $c_{2}\in(-c_{1}/2,c_{1}/2)$. We were able to impose the 0.3% relative error on the obtained results, by going with state bond dimensions $D_{\psi}$ up to $15$ while keeping $D_{\widetilde{L}}=1$. In interesting to note, that in the studied here cases higher $D_{\widetilde{L}}$ gave negligible gain in comparison to gain from increasing $D_{\psi}$—this is most probably related with the fact that local uncorrelated measurement are close to optimal, a fact often encountered in quantum metrology studies when the QFI is the only figure of merit to be optimised.

The main qualitative feature that clearly emerges from the Fig. 1(ii) is the decrease of the optimal QFI with the increase of correlated noise part parameter $c_{2}$. At the same time going into the anti-correlation regime (negative $c_{2}$) allows for a significant increase in the achievable QFI. This is to be expected as in the noise anti-correlation regime the noise operators responsible for correlations take the form of $L^{[n,n+1]}=\sqrt{|\gamma_{2}|}(h^{[n]}-h^{[n+1]})$ and are therefore linearly independent from the Hamiltonian operator which is a sum of $h^{[n]}$. This is related with the fact that in the extreme case of perfect anti-correlations (no independent local dephasing noise component) it is possible to employ quantum error-correction inspired protocols in order to preserve the Heisenberg scaling Layden and Cappellaro (2018). This is also reflected by the divergent behaviour of the state-of-the-art based bounds, which take into account all possible adaptive estimation strategies—note that there is no analogous divergence in the numerical results we obtain as we consider a parallel sensing strategy involving the most general entangled input states and most general measurement but no adaptive strategies. This again indicates, that if one aims at determining the potential of parallel entangled based strategies in presence of noise correlations, the state-of-the-art methods provide bounds which are very far from the actually achievable performance. The $c_{2}=0$ line on the Fig. 1(ii) corresponds to the strictly local dephasing model for which exact asymptotically saturable bound is known and reads $F/N=\eta^{2}/(1-\eta^{2})$ Escher et al. (2011); Demkowicz-Dobrzański et al. (2012), where $\eta=e^{-c_{1}/2}$. Our numerical results obtained using iMPO agree perfectly with this formula (see right inset on Fig. 1(ii)).

In case of purely local dephasing it is known that in the limit of large number of particles the fundamental bound, $F/N=\eta^{2}/(1-\eta^{2})$, can be saturated by protocols involving weakly spin squeezed states Ulam-Orgikh and Kitagawa (2001); Escher et al. (2011), e.g. one-axis twisted states Ma et al. (2011). Having obtained numerical optimal values of the QFI per particle using our MPO based methods, we want now to check whether the weakly spin squeezed strategy saturates the QFI in case of locally correlated noise model similarly as in the decorrelated noise scenario. For concreteness, consider the following one-axis squeezed state of $N$ particles:

[TABLE]

where $\theta$ is the squeezing strength. We follow the standard protocol Ma et al. (2011), where the above state is rotated to the equator of the Bloch sphere so that the $\langle\vec{S}\rangle$ points in the $x$ direction, in a way that the direction in which the angular momentum has minimal variance is $y$. The state is then subject to locally correlated dephasing evolution and is rotated by an unknown angle $\varphi$. Assuming we operate around $\varphi\approx 0$, we measure the $S_{y}$ observable as this is the optimal choice in this case, from which value we infer the value of $\varphi$. Using the standard linear error propagation formula the resulting uncertainty of estimating the phase reads: $\Delta\widetilde{\varphi}=\sqrt{\Delta^{2}S_{y}}/\left|\frac{\mathrm{d}\langle S_{y}\rangle}{\mathrm{d}\varphi}\right|$. In order to calculate the above quantity we move to the Heisenberg picture. Since we operate around $\varphi=0$, we can replace $\frac{\mathrm{d}\langle S_{y}\rangle}{\mathrm{d}\varphi}=\langle S_{x}\rangle$, and plug $\varphi=0$ everywhere. Under the locally correlated dephasing noise the relevant expectation values should be replaced according to the following rules:

[TABLE]

Taking the limit $N\rightarrow\infty$, $\theta\rightarrow 0$ in a way that $N\theta^{2}\ll 1$ we obtain that the relevant expectation values on the squeezed state read Ma et al. (2011) (assume $\hbar=1$): $\langle s_{x}^{[n]}\rangle\rightarrow\frac{1}{2}e^{-c_{1}}$, $\langle s_{y}^{[n]}\rangle=0$, $\langle s_{y}^{[n]}s_{y}^{[m]}\rangle\rightarrow-\frac{1}{4(N-1)}$, $\langle s_{x}^{[n]}s_{x}^{[m]}\rangle\rightarrow\frac{1}{4}$ ($n\neq m$). As a result $\langle S_{x}\rangle=\frac{N}{2}e^{-c_{1}}$ while

[TABLE]

This leads to the final formula for the asymptotic precision $\Delta\widetilde{\varphi}=\sqrt{(1-e^{-c_{1}}+2e^{-c_{1}}\sinh c_{2})/(Ne^{-c_{1}})}$ which can be related with the corresponding Fisher information per particle equal to:

[TABLE]

We have checked that this formula agrees with our numerical results up to the desired accuracy (${<1\%}$), and the representative comparison of the numerical data and this formula is provided in the left inset of Fig. 1(ii). This implies that similarly as in the uncorrelated dephasing models, weakly spin-squeezed states are asymptotically optimal. Note that this does not imply that the optimal MPO states we have obtained in our numerical procedure are these kind of states. Quite contrary, our MPO approach favours states with low bond dimension and local correlations, while the above states due to their fully symmetric nature correlate all the particles with each other irrespectively of their distance. Had we considered a product state strategy, the only modification in the above reasoning would be a substitution $\langle s_{y}^{[n]}s_{y}^{[m]}\rangle=0\ (n\neq m)$, which would lead to the corresponding QFI per particle:

[TABLE]

which also agrees perfectly with the numerical results we have presented in Fig. 1(iii).

Using the above formula, we may also go back to the original problem of magnetic field sensing. Utilizing the relation $\varphi=gBt$ we get the corresponding magnetic field sensing precision:

[TABLE]

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:

[TABLE]

which when optimized over $t$ reaches the minimal value when $t\rightarrow 0$ and yields:

[TABLE]

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 $N\rightarrow\infty$. In this case, following analogous calculations, we would arrive at the optimal magnetic field sensing precision of the form

[TABLE]

where $\chi_{k}$ represent magnetic field correlations for particles at distance $k-1$: $\langle\delta B^{[n]}(t)\delta B^{[n+k-1]}(t^{\prime})\rangle=\chi_{k}\delta(t-t^{\prime})$. Comparing this result with the performance of the GHZ states for the same model, see Layden et al. (2019), we notice that there is the $\sqrt{e}$ factor improvement in performance of the optimally spin-squeezed states over the GHZ states familiar from uncorrelated dephasing considerations Huelga et al. (1997); Ulam-Orgikh and Kitagawa (2001).

IV.2 Atomic clock stabilization

When one follows the Bayesian rather than the frequentist line of reasoning and focuses on minimization of the Bayesian variance, an apparently similar computational problem arises as in the cases of optimization of the QFI Macieszczak et al. (2014); Demkowicz-Dobrzański et al. (2015). The goal is then to minimize:

[TABLE]

where $p(x|\varphi)=\mathrm{Tr}({\rho}_{\varphi}{\Pi}_{x})$ and $p(\varphi)$ is a prior distribution for the parameter to be estimated (for simplicity we assume that the prior is centered at [math]: $\int\mathrm{d}\varphi\,p(\varphi)\varphi=0$). The minimal achievable quadratic Bayesian cost for the problem, optimized over all measurements and estimators, reads Helstrom (1976); Macieszczak et al. (2014):

[TABLE]

where $\Delta^{2}_{0}$ is the variance of the prior distribution, $\underline{\rho}=\int\mathrm{d}\varphi\,p(\varphi)\rho_{\varphi}$ is the output state averaged with the prior, while $\underline{\rho}^{\prime}=\int\mathrm{d}\varphi\,p(\varphi)\varphi\rho_{\varphi}$. It is clear from the above formula that the problem is computationally very similar to calculation of the QFI, as given in Eq. 8 (up to a replacement of the derivative of $\rho_{\varphi}$ with $\underline{\rho}^{\prime}$)

An important problem where the Bayesian line of reasoning is relevant, and which at the same time is very well suited for our tensor network framework is the atomic clock stabilization problem. A typical atomic clock operates in a feedback loop where the local oscillator (LO, e.g. laser) is stabilised 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 Ludlow et al. (2015). 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) Allan (1966); Riehle (2004):

[TABLE]

where $\left\langle\cdot\right\rangle$ represents averaging over frequency fluctuations of the LO described by some stochastic process (which from the Bayesian estimation perspective plays the role of the prior distribution), $\tau$ denotes averaging time, $\omega_{0}$ atomic reference angular frequency and $\omega(t)$ time-dependent angular frequency of the LO.

Fixing the physical properties of the atoms the goal is to optimize their initial states states, interrogations times, measurements and feedback corrections in order to minimize the AVAR. Performing such a comprehensive optimization is not feasible. In Chabuda et al. (2016) a lower bound on the achievable AVAR was introduced the quantum Allan Variance (QAVAR):

[TABLE]

where $\sigma_{\mathrm{LO}}^{2}(\tau)$ is the AVAR of free running LO, $\frac{1}{\omega_{0}^{2}}F_{\mathrm{A}}(\tau;\rho_{0},L,T)$ 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 $\underline{\rho}$ and $\underline{\rho}^{\prime}$ and refer the interested reader to Chabuda et al. (2016), but just note that they are analogs of $\underline{\rho}$, $\underline{\rho}^{\prime}$ as defined in the simple Bayesian estimation problem in Eq. (73). The important information is that, if the atoms with which the atomic clock interacts are described via states on some $d$ dimensional Hilbert space $\mathcal{H}$, then the $\underline{\rho}$ and $\underline{\rho}^{\prime}$ objects act on a tensor space $\mathcal{H}^{\otimes N}$, where $N=2(\tau/T)-1$ is the number of atomic cycles that need to be considered in order to calculate QAVAR. We assume here an idealized situation that one cycle of atomic clock operation lasts $T$—there is no dead time, and all the time is dedicated to interrogation of atoms by the LO. Typically we will think of the $d$ dimensional spaces of atomic states as being a fully symmetric subspace of $d-1$ two-level systems representing the relevant clock transition levels of the atoms—therefore we will use the notation for the states $|\psi\rangle=\sum_{k=0}^{d-1}a_{k}|k\rangle$, which correspond to the symmetric state, where $k$ atoms are in an excited state and $d-1-k$ atoms in the ground state. Interaction between the LO and the atomic sample has form of Ramsey interferometry and after a single interrogation step will effectively encode the phase in the above written state as $|\psi\rangle_{T}=\sum_{k=0}^{d-1}a_{k}e^{-ik\int_{0}^{T}\delta\omega(t)\mathrm{d}t}|k\rangle$, where $\delta\omega(t)$ is the detuning of the LO frequency from the atomic reference frequency $\omega_{0}$. The effect of LO fluctuations is equivalent to collective dephasing of atoms and hence the state will remain within the symmetric subspace—note that we are talking about the single Hilbert space $\mathcal{H}$ here which will be represented by single node in the MPO framework. The key feature from our perspective is the fact that the LO frequency fluctuations are temporally correlated, and hence the collective dephasing acting on atoms at different interrogation steps (in our representation different steps are formally represented as different product subsystems in the $\mathcal{H}^{\otimes N}$ space) will be correlated.

Assuming that LO fluctuations have finite correlations in time we can expect that QAVAR can be efficiently calculated using a tensor network in the form of a chain of length $N$ with $(d-1)$-dimensional physical indices on each site. Apart from a clear numerical efficiency advantages, the use of tensor networks approach also allows us to constrain the class of input states to be product (bond dimension $D_{\psi}=1$) which corresponds to the typical situation in which atomic samples in different time steps are independent of each other, as they prepared anew at the beginning of each interrogation step—note that we will still consider entanglement between physical atoms with which the LO interacts at a given interrogation time step.

LO fluctuations can be characterized by the autocorrelation function $R(t)$ which for the purpose of our example we choose to be a combination of Ornstein–Uhlenbeck (OU) process and white Gaussian frequency noise:

[TABLE]

where we can interpret parameters $\alpha$, $\beta$ as strength of respectively OU process and white noise and $1/\gamma$ as OU correlation range. We choose $\alpha=1\,(\mathrm{rad}/\mathrm{s})^{2}$, $\beta=0.1\,(\mathrm{rad}/\mathrm{s})^{2}\mathrm{s}$, $\gamma=2\,\mathrm{s}^{-1}$ for which noise correlations on the time scale that will correspond to the optimal interrogation time step $T$ (which in this case will happen to be around $1.5\,\mathrm{s}$) will be weak enough so that they will appreciable affect only the nearest-neighbour “time-step subsystems”—studying noise with further correlations is possible but would require more computational resources because of the larger bond dimensions. For such a noise the AVAR of the free running LO reads:

[TABLE]

and in the most interesting regime of large averaging times takes the form $\sigma_{\mathrm{LO}}^{2}(\tau)\simeq(2\alpha\gamma^{-1}+\beta)/(\tau\omega_{0}^{2})$. We expect that in this limit QAVAR also takes the form $\sigma_{\mathrm{Q}}^{2}(\tau)\simeq c/(\tau\omega_{0}^{2})$ with some constant $c$ which we will refer to as asymptotic coefficient.

These calculations have been attempted in Chabuda et al. (2016) 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. The tensor network framework, proposed in this paper, allows us to calculate QAVAR in previously inaccessible regime of large $\tau$. In Fig. 3(i) we present the exemplary results for an atomic clock operating on one two-level atom which also shows that completely neglecting noise correlations in analysing the clock performance is unjustified. We see that the QAVAR curve flatten for $\tau\gtrsim 50\,\mathrm{s}$, which when taking into account the optimal interrogation times which in this case approaches $\sim 1.3\,\mathrm{s}$, implies that for calculations of the QAVAR we would need to consider $\sim 80$ interrogation steps and hence if full Hilbert space description was used for this purpose would require $2^{80}$ dimensional space—clearly an impossible task. Similarly as in the previous example, we may directly access the asymptotic behaviour ($\tau\rightarrow\infty$) of the QAVAR function with the help of the iMPO approach. 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. 3(ii-iii). From this figures we see that the differences in QAVAR between cases with strictly local noise and when nearest neighbour 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. Note that, the definition (74) of the AVAR leads to an intrinsically not translationally invariant MPO for $\underline{\rho}^{\prime}$ in the expression for QAVAR. Because of this when implementing the iMPO approach we approximate AVAR by an asymptotically equivalent expression:

[TABLE]

which coincide with previous definition for $\tau$ much larger then the noise correlation length (which is exactly the regime in which iMPO approach operates). In order to keep relative errors below $1\%$ the numerical calculations required bond dimension $D_{L}=4$ for finite $N$ case and $D_{\widetilde{L}}=1$ in the case of the iMPO approach.

We confront the results (which are optimized over the input state) with values obtained using a NOON/GHZ states as an input, $|\psi\rangle=(|0\rangle+|d-1\rangle)/\sqrt{2}$. 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 visible in the Fig. 3(ii-iii), where we see that even though the optimal $T$ for the NOON state scales down (in fact as $(d-1)^{-2}$), the noise quickly destroys any gain from the feedback loop to AVAR of a free running LO. This is a manifestation of a generic poor performance of the NOON/GHZ states in realistic (noisy) scenarios with increasing particle number $N$ Huelga et al. (1997); Escher et al. (2011); Demkowicz-Dobrzański et al. (2012).

Our framework allow us also to easily study the optimal input states. In Fig. 3(iv) we plot the absolute values of the probability amplitudes $a_{k}$ for the optimal states (for LO noise which is strictly local or also includes the nearest neighbours correlations) alongside $a_{k}$ for coherent spin state (CSS) and the sine state Berry and Wiseman (2000) (which is the optimal input state for Bayesian estimation of phase with a flat prior distribution). We see that our optimal states are nonclassical, which we can quantify by calculating their spin-squeezing parameter $\xi^{2}=2\Delta^{2}J_{y}/\langle J_{x}\rangle$ Ma et al. (2011) which is $1$ for CSS, $0.522$ for sine state, $0.456$ and $0.378$ for the optimal states for LO noise which is respectively strictly local or also includes the nearest neighbours correlations (squeezing in this last case can be also observed from Husimi Q distribution on Bloch sphere in the centre of the Fig. 3(iv)).

IV.3 Fidelity susceptibility calculations for many-body thermal states

In condensed matter context, fidelity $\mathcal{F}(\varphi,\varphi+\varepsilon)=\mathrm{Tr}\sqrt{\sqrt{\rho_{\varphi}}\rho_{\varphi+\varepsilon}\sqrt{\rho_{\varphi}}}$ between many-body states $\rho_{\varphi}$ and $\rho_{\varphi+\varepsilon}$, that differ by a small variation of a parameter $\varphi$ in a Hamiltonian, is a mean to identify the location $\varphi_{c}$ of a phase transition Zanardi and Paunković (2006); Rams and Damski (2011). This is where the fidelity susceptibility $\chi_{\varphi}$, defined by $\mathcal{F}(\varphi,\varphi+\varepsilon)\approx 1-\frac{1}{2}\chi_{\varphi}\varepsilon^{2}$, has a maximum indicating a fundamental change in the state of the system. This concept was employed in Ref. Rams et al. (2018) to evaluate the usefulness of a quantum phase transition, that happens at zero temperature, for precise sensing of the parameter $\varphi$ in a realistic system at a finite temperature. QFI defines a metric in the space of quantum states (the Bures metric)Braunstein and Caves (1994) and is directly related with the fidelity susceptibility, namely $F=4\chi_{\varphi}$.

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 $\widetilde{\mathcal{F}}(\varphi,\varphi+\varepsilon)=\sqrt{\mathrm{Tr}\sqrt{\rho_{\varphi}}\sqrt{\rho_{\varphi+\varepsilon}}}$ defining a quasi-susceptibility, $\widetilde{\mathcal{F}}(\varphi,\varphi+\varepsilon)\approx 1-\frac{1}{2}\widetilde{\chi}_{\varphi}\varepsilon^{2}$, that provides bounds for the exact fidelity susceptibility Albuquerque et al. (2010); Sirker (2010):

[TABLE]

The Hamiltonian considered in Ref. Rams et al. (2018) was the spin-$\frac{1}{2}$ XX model

[TABLE]

with a quantum critical point at $\varphi_{c}=0$. Taking the MPOs studied in Ref. Rams et al. (2018), we bypass the tractability problem employing the part of our scheme with $\rho^{\prime}_{\varphi}=\left(\rho_{\varphi+\varepsilon}-\rho_{\varphi}\right)/\varepsilon$ to calculate the QFI $F=4\chi_{\varphi}$. This exact susceptibility for the chain with $64$ spins is shown in Fig. 4 together with the upper and lower bounds (80). The accuracy of the fidelity susceptibility is limited by the finite bond dimension $D_{L}$ as well the finite parameter difference $\varepsilon$. Nevertheless, we obtain satisfying results with relative error around $1\%$ for $\varepsilon=10^{-4}$ and $D_{L}=4$.

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 Zanardi and Paunković (2006); Rams and Damski (2011)—by now standard in condensed matter physics—to phase transitions in quantum many-body systems at finite temperature. Their thermal states can be represented either by MPO, when DMRG on a cylinder is employed Chen et al. (2018), or its two-dimensional generalization on an infinite lattice (iPEPO) Czarnik et al. (2016); Czarnik and Corboz (2019).

V Conclusions

We have provided a comprehensive framework for optimization of quantum metrological protocols using the MPO/MPS formalism. The potential to deal effectively with correlated noise models as well as directly access the asymptotic $N\rightarrow\infty$ is what makes this framework unique. We also expect that this framework may also be adapted to deal with even more challenging metrological problems including noisy multiparameter estimation Ragy et al. (2016); Szczykulska et al. (2016); Baumgratz and Datta (2016), waveform estimation Tsang et al. (2011); Berry et al. (2013) or the study of the effectiveness of adaptive metrological protocols including quantum error correction based schemes Arrad et al. (2014); Dür et al. (2014); Zhou et al. (2018); Layden and Cappellaro (2018); Gorecki et al. (2019). We also expect that this numerical framework may be crucial for understanding better metrological models with temporally correlated noise especially of non-Markovian nature Chin et al. (2012), where effective tools to find the optimal metrological protocols in such cases are yet to be developed.

Acknowledgements.

We would like to thank Marek M. Rams, David Layden, Maciej Lewenstein, Shi-Ju Ran, Piet O. Schmidt and Ian D. Leroux for fruitful discussions. We are also indebt to Marek M. Rams for sharing with us the data from Ref. Rams et al. (2018). KCh and RDD acknowledge support from the National Science Center (Poland) grant No. 2016/22/E/ST2/00559. Work of JD was funded by NCN together with European Union through QuantERA ERA NET program 2017/25/Z/ST2/03028. TJO was supported, in part, by the DFG through SFB 1227 (DQmat), the RTG 1991, and the cluster of excellence EXC 2123 QuantumFrontiers.