Quantum metrology for the Ising Hamiltonian with transverse magnetic field

Quantum metrology is one of the archetypical applications where quantum mechanics demonstrably exhibits a vast improvement over the best known classical strategies. The use of entangled input states of N qubits, such as Greenberger–Horne–Zeilinger (GHZ) states [1], is known to allow for a precise determination of an unknown parameter such as the relative phase in a Mach–Zender interferometer [2–4], or the frequency of an atomic transition [5, 6] with a precision that scales inversely proportional to N, the Heisenberg limit. In comparison the best classical estimation strategy, which employs separable states, gives a precision scaling inversely proportional to $\sqrt{N}$, the standard quantum limit [7, 8]. This observation has proven very useful in the development of ultra-precise atomic clocks [9–11], high resolution imaging [12], as well as the detection of gravitational waves [13, 14].

In the absence of any noise or decoherence effects the parameter of interest, λ, in a metrological scenario is imprinted onto the state of N probes via the unitary operator $U={{\rm{e}}}^{{\rm{i}}H(\lambda )t}$, where $H(\lambda )$ is the Hamiltonian describing the dynamics of the N probes. For the case of phase and frequency estimation—where the parameter to be estimated is $\lambda t$ and λ respectively—the parameter is a multicplicative factor of the Hamiltonian, $H(\lambda )=\lambda H$, and the Hamiltonian is local, i.e., $H={\displaystyle \sum }_{i=1}^{N}{h}^{(i)}$ where ${h}^{(i)}$ is the Hamiltonian describing the evolution of the $i\mathrm{th}$ probe. Indeed, quantum metrology using local Hamiltonians where the parameter of interest enters only as a multiplicative factor have been studied extensively both in the absence as well as in the presence of noise [15–23].

In stark contrast quantum metrology with more general Hamiltonians is only now beginning to attract attention. Some instances of quantum metrology with parameter dependent Hamiltonians concern the estimation of time-varying signals [24–26], the estimation of magnetic-field gradients along a spin chain [27–29], or the estimation of the anisotropy and/or decoherence along a spin-chain using non-equilibrium states [30]. The general problem of performing quantum metrology using parameter dependent Hamiltonians was treated in [31, 32] where it was shown that for parameter dependent local Hamiltonians Heisenberg scaling precision with respect to the number of probing systems is possible.

In this work we focus on noiseless quantum metrology of N interacting probe systems described by the Ising Hamiltonian

$$\begin{eqnarray}&&H(J,B)=J\underset{i=1}{\overset{N-1}{\displaystyle \sum }}{\sigma }_{x}^{(i)}{\sigma }_{x}^{(i+1)}+B\underset{i=1}{\overset{N}{\displaystyle \sum }}{\sigma }_{z}^{(i)},\end{eqnarray} \tag{ 1 }$$

where we are interested in determining the precision with which one can estimate either the strength of the transverse magnetic field, B, or the coupling interaction, J, provided the remaining quantity is known. The Hamiltonian of equation (1) is known to exhibit a phase transition [33], which have been discussed previously in the context of quantum metrology and were shown to be resourceful [34–37]. Moreover, as the Ising Hamiltonian is entanglement generating [38], it has found applications in ion-trap quantum computing architectures [39, 40], where either J or B can be controlled at will by modifying either the separation of the ions, or the global magnetic field strength.

Here we report the following results regarding precise estimation of parameters of the Ising Hamiltonian:

• 1.  
  The ultimate precision with which one can estimate either J or B, having complete knowledge of B or J, using N probe systems scales at the Heisenberg limit, and we provide an analytic expression for this achievable precision as well as the optimal states that achieve it. Our result answers the conjecture in [31] that Heisenberg scaling is the ultimate limit also for the case of N interacting probes for the case of nearest-neighbour interactions.
• 2.  
  We provide numerical evidence that the ultimate achievable precision strictly outperforms the optimal classical strategy, which deploys the N initial probes in a pure product state. For up to N = 11 our numerical study shows that the optimal product input state yields a precision that scales at the standard quantum limit.
• 3.  
  We analytically derive the precision achieved by the GHZ-type states, known to achieve the optimal precision when either J or B are equal to zero, and show that these states retain their Heisenberg scaling in precision (up to a constant factor) over the entire regime of parameters.

This paper is organized as follows. In section 2 we review the basics of quantum metrology as well as some important mathematical results regarding unitary operators generated by parameter dependent Hamiltonians. In section 3 we use the Jordan–Wigner transformation to determine the maximal possible precision with which one can estimate either J and B using N systems as probes. We compare this to the best possible precision that can be achieved by a separable state of N probes which we determine numerically for up to N = 11 qubits. In section 4 we analytically determine the performance of states that are known to be optimal at the extreme cases where J = 0 and B = 0. Finally, section 5 contains the conclusions of our investigation as well as some open questions for future work.

In this section we review the main results in quantum metrology and review some important facts pertaining to parameter-dependent Hamiltonians in general and to the Ising Hamiltonian (equation (1)) in particular. Specifically, we concentrate on noiseless quantum metrology and the quantum Fisher information (QFI); the central quantity of interest in quantum metrology. After introducing the QFI we provide a formula for calculating it for the case of general parameter-dependent Hamiltonians and then to the specific case of the Ising Hamiltonian. We note that the behaviour of the QFI, and in particular its scaling with time and number of probe systems, was investigated for general parameter dependent Hamiltonians in [31, 32].

A standard protocol in noiseless quantum metrology can be formulated as follows: N probes are prepared in a suitable state and undergo an evolution for some time, t, described by the unitary operator $U({\boldsymbol{\lambda }},t)={{\rm{e}}}^{-{\rm{i}}{tH}({\boldsymbol{\lambda }})}$, where $H({\boldsymbol{\lambda }})$ is the Hamiltonian describing the dynamics of the N probes and explicitly depends on the vector of parameters ${\boldsymbol{\lambda }}\equiv ({\lambda }_{1},\ldots ,{\lambda }_{M})$. Finally the N probes are measured and an estimate of $\hat{{\boldsymbol{\lambda }}}$ is obtained from the measurement statistics of ν repetitions of the above procedure. In what follows we shall assume that all other parameters except ${\lambda }_{i}$ are known, and shall be concerned with estimating ${\lambda }_{i}$ as precisely as possible. A lower bound on the error, $\delta {\lambda }_{i}={\lambda }_{i}-{\hat{\lambda }}_{i}$, for any unbiased estimator ${\hat{\lambda }}_{i}$ is given by the quantum Cramér-Rao bound [41]

$$\begin{eqnarray}&&\delta {\lambda }_{i}\geqslant \displaystyle \frac{1}{\sqrt{\nu {\mathcal{F}}({\rho }_{{\boldsymbol{\lambda }},t})}},\end{eqnarray} \tag{ 2 }$$

where ${\mathcal{F}}({\rho }_{{\boldsymbol{\lambda }},t})$ is the QFI of the state ${\rho }_{{\boldsymbol{\lambda }},t}$ describing the N probes after the unitary dynamics have acted. In the most general case the QFI can be computed as [41]

$$\begin{eqnarray}&&{\mathcal{F}}({\rho }_{{\boldsymbol{\lambda }},t})=\mathrm{Tr}\left\{{L}_{{\lambda }_{i}}{\rho }_{{\boldsymbol{\lambda }},t}{L}_{{\lambda }_{i}}\right\},\end{eqnarray} \tag{ 3 }$$

with

$$\begin{eqnarray}&&{L}_{{\lambda }_{i}}=2\underset{\alpha ,\beta }{\displaystyle \sum }\displaystyle \frac{\langle {\psi }_{\alpha }| {\partial }_{{\lambda }_{i}}{\rho }_{{\boldsymbol{\lambda }},t}| {\psi }_{\beta }\rangle }{\alpha +\beta }| {\psi }_{\alpha }\rangle \langle {\psi }_{\beta }| \end{eqnarray} \tag{ 4 }$$

the symmetric logarithmic derivative, where $\alpha \;(\beta )$ are the eigenvalues of${\rho }_{{\boldsymbol{\lambda }},t}$, $| {\psi }_{\alpha }\rangle ,(| {\psi }_{\beta }\rangle )$ the corresponding eigenvectors, and the sum in equation (4) is over all $\alpha ,\beta$ satisfying $\alpha +\beta \ne 0$. Here and in what follows ${\partial }_{x}\equiv \frac{\partial }{\partial x}$. We now review the case where the parameter of interest enters as a multiplicative factor of the Hamiltonian.

2.1. Parameter independent Hamiltonians

In the case of noiseless metrology, an easier expression for computing the QFI exists if one initializes the N probes in a pure state $\rho =| \psi \rangle \langle \psi |$. In this case the QFI can be shown to be

$$\begin{eqnarray}&&{\mathcal{F}}\left(| {\psi }_{{\boldsymbol{\lambda }},t}\rangle \right)=4\left(\left\langle{\partial }_{{\lambda }_{i}}{\psi }_{{\boldsymbol{\lambda }},t}| {\partial }_{{\lambda }_{i}}{\psi }_{{\boldsymbol{\lambda }},t}\right\rangle-{\left|\left\langle{\partial }_{{\lambda }_{i}}{\psi }_{{\boldsymbol{\lambda }},t}| {\psi }_{{\boldsymbol{\lambda }},t}\right\rangle\right|}^{2}\right),\end{eqnarray} \tag{ 5 }$$

where $| {\psi }_{{\boldsymbol{\lambda }},t}\rangle =U({\boldsymbol{\lambda }},t)| \psi \rangle$. If $H({\boldsymbol{\lambda }})={\lambda }_{i}H$ then a bit of algebra yields ${\mathcal{F}}\left(| {\psi }_{{\boldsymbol{\lambda }},t}\rangle \right)=4{t}^{2}{\Delta }^{2}\;H$, where ${\Delta }^{2}\;H\equiv \langle \psi | {H}^{2}| \psi \rangle -{| \langle \psi | H| \psi \rangle | }^{2}$ is the variance of H with respect to the state $| \psi \rangle$. If, in addition, $H={\displaystyle \sum }_{i=1}^{N}{h}^{(i)}$ where ${h}^{(i)}={h}^{(j)}=h,\forall i,j$ is the Hamiltonian acting on the ith probe system, then it can be shown that when $| \psi \rangle ={(\frac{| {\alpha }_{\mathrm{min}}\rangle +| {\alpha }_{\mathrm{max}}\rangle }{\sqrt{2}})}^{\otimes N}$, where $| {\alpha }_{\mathrm{min}(\mathrm{max})}\rangle$ are the eigenstates corresponding to the minimum and maximum eigenvalues of h, then ${\mathcal{F}}\left(| {\psi }_{{\boldsymbol{\lambda }},t}\rangle \right)={t}^{2}\;N{({\alpha }_{\mathrm{max}}-{\alpha }_{\mathrm{min}})}^{2}$ and give the standard quantum limit in estimation precision[8]. On the other hand if the probes are prepared in the GHZ state

$$\begin{eqnarray}&&| \psi \rangle =\displaystyle \frac{1}{\sqrt{2}}(\underbrace{{| {\alpha }_{\mathrm{min}},\ldots ,{\alpha }_{\mathrm{min}}\rangle }}_{N\;\mathrm{times}}+\underbrace{{| {\alpha }_{\mathrm{max}},\ldots ,{\alpha }_{\mathrm{max}}\rangle }}_{N\;\mathrm{times}}),\end{eqnarray} \tag{ 6 }$$

then ${\mathcal{F}}\left(| {\psi }_{{\boldsymbol{\lambda }},t}\rangle \right)={t}^{2}{N}^{2}{({\alpha }_{\mathrm{max}}-{\alpha }_{\mathrm{min}})}^{2}$, the Heisenberg scaling in estimation precision [8]. In general, Heisenberg scaling in precision implies an improvement in scaling with regards to the number of probe systems, over the optimally achievable precision which uses the same number of probe systems initialized in a separable state. Hence, in the case where the parameter to be estimated is a multiplicative factor of a local Hamiltonian, the use of highly entangled states leads to a quadratic improvement in scaling precision. Note that besides the GHZ states there exists a large class of pretty good states that scale at the Heisenberg limit up to a multiplicative factor [23].

2.2. Parameter dependent Hamiltonians

We now use equation (5) to compute the QFI for Hamiltonians of the form

$$\begin{eqnarray}&&H({\boldsymbol{\lambda }})={\lambda }_{1}{H}_{1}+{\lambda }_{2}{H}_{2},\end{eqnarray} \tag{ 7 }$$

where $[{H}_{1},{H}_{2}]\ne 0$. We make no assumptions on the structure of ${H}_{1},{H}_{2}$; in particular we do not assume that they are local Hamiltonians. Notice that the Ising Hamiltonian of equation (1) is a special case of equation (7), whose QFI was studied in detail in [31] with ${\lambda }_{1}{H}_{1}=J{\displaystyle \sum }_{i=1}^{N}{\sigma }_{x}^{(i)}{\sigma }_{x}^{(i+1)}$ and ${\lambda }_{2}{H}_{2}=B{\displaystyle \sum }_{i=1}^{N}{\sigma }_{z}^{(i)}$. From equation (5) we need to compute $| {\partial }_{{\lambda }_{i}}{\psi }_{{\boldsymbol{\lambda }},t}\rangle ={\partial }_{{\lambda }_{i}}U({\boldsymbol{\lambda }},t)| \psi \rangle$. As $[{H}_{1},{H}_{2}]\ne 0$ and using

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{e}}}^{-{\rm{i}}\;t({\lambda }_{1}{H}_{1}+{\lambda }_{2}{H}_{2})}}{{\partial }_{{\lambda }_{i}}}=\underset{N\to \infty }{\mathrm{lim}}\displaystyle \frac{\partial {\left({\mathbb{1}}-\frac{i\;t}{N}({\lambda }_{1}{H}_{1}+{\lambda }_{2}{H}_{2})\right)}^{N}}{{\partial }_{{\lambda }_{i}}},\end{eqnarray} \tag{ 8 }$$

equation (5) reads

$$\begin{eqnarray}{\mathcal{F}}\left(| {\psi }_{{\boldsymbol{\lambda }},t}\rangle \right) & = & 4\left(\langle \psi | {U}^{\dagger }({\boldsymbol{\lambda }},t){{\mathcal{O}}}_{i}^{2}({\boldsymbol{\lambda }},t)U({\boldsymbol{\lambda }},t)| \psi \rangle \right.\\ & & \left.-{\left|\langle \psi | {U}^{\dagger }({\boldsymbol{\lambda }},t){{\mathcal{O}}}_{i}({\boldsymbol{\lambda }},t)U({\boldsymbol{\lambda }},t)| \psi \rangle \right|}^{2}\right)\\ & = & 4{\Delta }^{2}{{\mathcal{O}}}_{i}({\boldsymbol{\lambda }},t),\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{{\mathcal{O}}}_{i}({\boldsymbol{\lambda }},t)\equiv {\displaystyle \int }_{0}^{t}\;{\rm{d}}s\;U({\boldsymbol{\lambda }},s){H}_{i}{U}^{\dagger }({\boldsymbol{\lambda }},s)\end{eqnarray} \tag{ 10 }$$

and the variance of ${{\mathcal{O}}}_{i}$ is computed with respect to the evolved state $| {\psi }_{{\boldsymbol{\lambda }},t}\rangle$.

Thus, the QFI is maximized by the states that are linear superpositions of the eigenstates corresponding to the minimum and maximum eigenvalues of ${{\mathcal{O}}}_{i}({\boldsymbol{\lambda }},t)$. An interesting question is whether the presence of ${\lambda }_{1}{H}_{1}$ (${\lambda }_{2}{H}_{2}$) can help boost the precision of estimation of parameters for ${\lambda }_{2}$ (${\lambda }_{1}$) respectively. It was shown in [31] that this is not the case. Specifically, if one has control over either ${\lambda }_{1}{H}_{1}$ or ${\lambda }_{2}{H}_{2}$ and wishes to estimate ${\lambda }_{2}$ or ${\lambda }_{1}$ respectively, then the optimal strategy is to set the dynamics under our control to zero.

In the next section we use the expressions in equation (9) to determine the optimal QFI for either the magnetic field B or interaction strength J of the Ising Hamiltonian using separable and entangled states respectively.

We are interested in determining the optimal precision in estimating either the magnetic field strength, B, or interaction strength, J, of the Ising Hamiltonian (equation (1)). In particular, we will show that the optimal precision in estimating either B or J scales at the Heisenberg limit, up to a constant factor which depends only on the ratio of J and B, and is achievable by states that are linear superpositions of the vacuum and fully occupied states of free fermions of a suitable type. Furthermore, we numerically determine the best achievable precision using a separable state for up to 11 qubits and show that the optimal precision scales, to within best fit errors, linearly with N. Hence, our results provide strong evidence that the entanglement generated by the Ising Hamiltonian when acting on an initially pure separable state is not enough to boost the precision in estimation from the SQL to the Heisenberg limit.

We begin by first determining the optimal precision in estimation of either B, or J, and the corresponding optimal states. To do so we note that via the use of the Jordan–Wigner transformation [42] the Ising Hamiltonian of equation (1) can be expressed as a quadratic Hamiltonian in fermionic creation and annihilation operators which can be suitably diagonalized. Specifically the mapping

$$\begin{eqnarray}{a}_{j} & \equiv & \left(\underset{k=1}{\overset{j-1}{\displaystyle \otimes }}{\sigma }_{z}^{(k)}\right)\otimes {\sigma }_{-}^{j}\\ {a}_{j}^{\dagger } & \equiv & \left(\underset{k=1}{\overset{j-1}{\displaystyle \otimes }}{\sigma }_{z}^{(k)}\right)\otimes {\sigma }_{+}^{j}\end{eqnarray} \tag{ 11 }$$

and its inverse

$$\begin{eqnarray}{\sigma }_{-}^{(j)} & \equiv & \mathrm{exp}\left({\rm{i}}\pi \underset{k=0}{\overset{j-1}{\displaystyle \sum }}{a}_{k}^{\dagger }{a}_{k}\right){a}_{j}\\ {\sigma }_{+}^{(j)} & \equiv & \mathrm{exp}\left({\rm{i}}\pi \underset{k=0}{\overset{j-1}{\displaystyle \sum }}{a}_{k}^{\dagger }{a}_{k}\right){a}_{j}^{\dagger },\end{eqnarray} \tag{ 12 }$$

where $\{{\sigma }_{x},{\sigma }_{y},{\sigma }_{z}\}$ are the Pauli matrices² , with ${\sigma }_{\pm }^{(j)}\equiv {\sigma }_{x}^{(j)}\pm i{\sigma }_{y}^{(j)}$, and ${a}_{j},{a}_{j}^{\dagger }$ are the fermionic annihilation and creation operators for mode j respectively and satisfy the anti-commutation relations $\{{a}_{j}^{\dagger }\;{a}_{k}^{\dagger }\}=\{{a}_{j}\;{a}_{k}\}=0$, $\{{a}_{j},{a}_{k}^{\dagger }\}={\delta }_{{jk}}{\mathbb{1}}$. Substituting equation (12) into equation (1) yields

$$\begin{eqnarray}&&H(J,B)=J\underset{j=1}{\overset{N}{\displaystyle \sum }}({a}_{j}^{\dagger }-{a}_{j})({a}_{j+1}^{\dagger }+{a}_{j+1})+2B\underset{j=1}{\overset{N}{\displaystyle \sum }}{a}_{j}^{\dagger }{a}_{j}.\end{eqnarray} \tag{ 13 }$$

Any quadratic Hamiltonian in the fermionic operators can be brought to the diagonal form

$$\begin{eqnarray}&&\tilde{H}(J,B)=2\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\sqrt{{\alpha }_{k}^{2}+{\beta }_{k}^{2}}\quad {c}_{k}^{\dagger }{c}_{k},\end{eqnarray} \tag{ 14 }$$

where ${\alpha }_{k}=J\mathrm{cos}\left(\frac{2\pi k}{N}\right)+B$, ${\beta }_{k}=J\mathrm{sin}\left(\frac{2\pi k}{N}\right)$, and ${c}_{k}^{\dagger }$, c_k are the creation and annihilation operators of free fermions in mode k, with ${c}_{N}={c}_{0}$. The eigenstates of $\tilde{H}(J,B)$ are fermionic Fock states, $| {\boldsymbol{k}}\rangle$, where ${\boldsymbol{k}}$ is an N-bit binary string indicating which modes are occupied by fermions. Without loss of generality we may set the vacuum state, $| {\bf{0}}\rangle$ to have energy$\sqrt{{\alpha }_{0}^{2}+{\beta }_{0}^{2}}=0$. The maximally occupied fermionic Fock state, $| {\bf{1}}\rangle$ has energy equal to ${\displaystyle \sum }_{k=0}^{N-1}\sqrt{{\alpha }_{k}^{2}+{\beta }_{k}^{2}}$.

As equation (14) is of great importance in the remainder of this work, we now discuss the steps required for obtaining it. Starting from the quadratic Hamiltonian of equation (13) one first performs the Fourier transformation

$$\begin{eqnarray}&&{b}_{j}=\displaystyle \frac{1}{\sqrt{N}}\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\;{{\rm{e}}}^{-{\rm{i}}\frac{2\pi {jk}}{N}}{a}_{k},\end{eqnarray} \tag{ 15 }$$

of the mode operators a_k. After substituting equation (15) into equation (13) the Hamiltonian can be written in matrix form as

$$\begin{eqnarray}\tilde{H}(J,B)=\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\left({b}_{k}^{\dagger },{b}_{N-k}\right)\left(\begin{array}{cc}{\alpha }_{k} & {\rm{i}}{\beta }_{k}\\ -{\rm{i}}{\beta }_{k} & -{\alpha }_{k}\end{array}\right)\left(\begin{array}{c}{b}_{k}\\ {b}_{N-k}^{\dagger }\end{array}\right).\end{eqnarray} \tag{ 16 }$$

This block-diagonal structure of the Hamiltonian in terms of the mode operators b_k and ${b}_{N-k}$ makes it particularly useful for calculating the operator ${\mathcal{O}}$ of equation (10). Finally, from equation (16) one can go to the diagonalized Hamiltonian of equation (14) by performing a suitable Bogoliubov transformation on the two-dimensional block of mode operators ${b}_{k},{b}_{N-k}$

$$\begin{eqnarray}{c}_{k} & = & \mathrm{cos}{\theta }_{k}\;{b}_{k}-{{\rm{e}}}^{{\rm{i}}{\phi }_{k}}\mathrm{sin}{\theta }_{k}\;{b}_{N-k}^{\dagger }\\ {c}_{N-k} & = & {{\rm{e}}}^{{\rm{i}}{\phi }_{k}}\mathrm{sin}{\theta }_{k}\;{b}_{k}^{\dagger }+\mathrm{cos}{\theta }_{k}\;{b}_{N-k}\end{eqnarray} \tag{ 17 }$$

with ${\theta }_{k}=\frac{1}{2}{\mathrm{tan}}^{-1}\left(\frac{J\mathrm{sin}\left(\frac{2\pi k}{N}\right)}{J\mathrm{cos}\left(\frac{2\pi k}{N}\right)+B}\right)$ and ${\phi }_{j}=\frac{\pi }{2},\forall j$.

We now determine the optimal achievable precision in estimating j given that b is known.

3.1. Estimating the Interaction strength

We now proceed to estimate the interaction strength J for the Ising Hamiltonian. The QFI is given by equation (9) with

$$\begin{eqnarray}&&{{\mathcal{O}}}_{J}(J,B,t)={\displaystyle \int }_{0}^{t}U(J,B,s){H}_{1}{U}^{\dagger }(J,B,s){\rm{d}}s,\end{eqnarray} \tag{ 18 }$$

where ${H}_{1}={\displaystyle \sum }_{i=1}^{N}{\sigma }_{x}^{(i)}{\sigma }_{x}^{(i+1)}$. Writing the latter in terms of the fermionic operators b_k one obtains

$$\begin{eqnarray}&&{H}_{1}=\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\left(\begin{array}{l}{b}_{k}^{\dagger }{b}_{N-k}\end{array}\right)\underbrace{{\left(\begin{array}{l}\mathrm{cos}\left(\frac{2\pi k}{N}\right)\mathrm{isin}\left(\frac{2\pi k}{N}\right)\\ -\mathrm{isin}\left(\frac{2\pi k}{N}\right)-\mathrm{cos}\left(\frac{2\pi k}{N}\right)\end{array}\right)}}_{{M}_{k}}\left(\begin{array}{l}{b}_{k}\\ {b}_{N-k}^{\dagger }\end{array}\right).\end{eqnarray} \tag{ 19 }$$

In order to calculate the operator ${{\mathcal{O}}}_{J}(J,B,t)$ of equation (18) we need to determine the action of $H(J,B)$ on the fermionic operators b_k, i.e., we need to determine ${b}_{k}(s)=U(J,B,s){b}_{k}(0){U}^{\dagger }(J,B,s)$. This is simply the Heisenberg equation of motion for the mode operators b_k. Due to the block-diagonal structure of $\tilde{H}(J,B)$, when written in terms of fermionic operators b_k, the solution to the Heisenberg equation of motion can be seen to be

$$\begin{eqnarray}&&\left(\begin{array}{l}{b}_{k}(s)\\ {b}_{N-k}^{\dagger }(s)\end{array}\right)=\underbrace{{\mathrm{exp}\left(-{\rm{i}}\;s\left(\begin{array}{cc}{\alpha }_{k} & {\rm{i}}{\beta }_{k}\\ -{\rm{i}}{\beta }_{k} & -{\alpha }_{k}\end{array}\right)\right)}}_{{R}_{k}(s)}\left(\begin{array}{l}{b}_{k}(0)\\ {b}_{N-k}^{\dagger }(0)\end{array}\right).\end{eqnarray} \tag{ 20 }$$

Henceforth, we drop the explicit time dependence of the mode operators b_k for convenience.

Substituting equation (20) into equation (18) yields

$$\begin{eqnarray}&&{{\mathcal{O}}}_{J}(J,B,t)=\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\left({b}_{k}^{\dagger }\;\;{b}_{N-k}\right)\left({\displaystyle \int }_{0}^{t}\;{\rm{d}}{{sR}}_{k}^{\dagger }(s){M}_{k}{R}_{k}(s)\right)\left(\begin{array}{l}{b}_{k}\\ {b}_{N-k}^{\dagger }\end{array}\right).\end{eqnarray} \tag{ 21 }$$

Performing the integration over s gives

$$\begin{eqnarray}{{\mathcal{O}}}_{J}(J,B,t)=\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\left(\begin{array}{l}{b}_{k}^{\dagger }\;\;{b}_{N-k}\end{array}\right)\left(\begin{array}{ll}{\Omega }_{k} & {\Delta }_{k}\\ {\Delta }_{k}^{*} & -{\Omega }_{k}\end{array}\right)\left(\begin{array}{l}{b}_{k}\\ {b}_{N-k}^{\dagger }\end{array}\right),\end{eqnarray} \tag{ 22 }$$

where

$$\begin{eqnarray}{\Omega }_{k} & = & \frac{{\alpha }_{k}\left({\alpha }_{k}\mathrm{cos}\left(\frac{2\pi k}{N}\right)+{\beta }_{k}\mathrm{sin}\left(\frac{2\pi k}{N}\right)\right)t+\displaystyle \frac{{\beta }_{k}\left({\beta }_{k}\mathrm{cos}\left(\frac{2\pi k}{N}\right)-{\alpha }_{k}\mathrm{sin}\left(\frac{2\pi k}{N}\right)\right)\mathrm{sin}\left(2{\omega }_{k}t\right)}{2{\omega }_{k}}}{{\omega }_{k}^{2}}\\ {\Delta }_{k} & = & \left({\rm{i}}{\beta }_{k}\left({\alpha }_{k}\mathrm{cos}\left(\frac{2\pi k}{N}\right)+{\beta }_{k}\mathrm{sin}\left(\frac{2\pi k}{N}\right)\right)t+\left({\beta }_{k}\mathrm{cos}\left(\frac{2\pi k}{N}\right)\right.\right.\\ & & -\left.\left.{\alpha }_{k}\mathrm{sin}\left(\frac{2\pi k}{N}\right)\right)\left({\mathrm{sin}}^{2}\left({\omega }_{k}t\right)-\frac{{\rm{i}}{\alpha }_{k}\mathrm{sin}\left(2{\omega }_{k}t\right)}{2{\omega }_{k}}\right)\right)/({\omega }_{k}^{2}),\end{eqnarray} \tag{ 23 }$$

where we have separated the linear and oscillatory behaviour of ${\Omega }_{k}$ and ${\Delta }_{k}$ (see also [32]) and ${\omega }_{k}=\sqrt{{\alpha }_{k}^{2}+{\beta }_{k}^{2}}$. As equation (22) is of the same form as equation (16) it can be brought to the diagonal form

$$\begin{eqnarray}&&{{\mathcal{O}}}_{J}(J,B,t)=2\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\sqrt{{\Omega }_{k}^{2}+{\left|{\Delta }_{k}\right|}^{2}}\quad {d}_{k}^{\dagger }{d}_{k}\end{eqnarray} \tag{ 24 }$$

by a suitable Bogoliubov transformation (see equation (17)). Note that the free fermions corresponding to ${d}_{k}^{\dagger },{d}_{k}$ are different from those of equation (14), and we may choose without loss of generality the positive square square root of ${\Omega }_{k}^{2}+{\left|{\Delta }_{k}\right|}^{2}$, which corresponds to choosing the vacuum state for free fermions to be zero.

The optimal precision in estimating the interaction strength J is now easy to determine. As the latter is inversely proportional to the square root of the variance of ${{\mathcal{O}}}_{J}(J,B,t)$, we simply need to determine the maximum achievable variance for the operator in equation (24). This is achieved by preparing the equally weighted superposition of the vacuum state and the state where all N modes are occupied by fermions. The variance with respect to this state is simply given by

$$\begin{eqnarray}&&{\Delta }^{2}{{\mathcal{O}}}_{J}{(J,B,t)}_{\mathrm{max}}={\left(\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\sqrt{{\Omega }_{k}^{2}+{\left|{\Delta }_{k}\right|}^{2}}\right)}^{2}.\end{eqnarray} \tag{ 25 }$$

As the sum includes N summands, the variance of ${{\mathcal{O}}}_{J}(J,B,t)$ scales as N² up to some constant factor that depends solely on the ratio between J and B as we now explain.

Using equation (23) one can easily show that

$$\begin{eqnarray}&&{\Omega }_{k}^{2}+{\left|{\Delta }_{k}\right|}^{2}\\ &&\quad =\displaystyle \frac{{\left({\beta }_{k}\mathrm{cos}\left(\frac{2\pi k}{N}\right)-{\alpha }_{k}\mathrm{sin}\left(\frac{2\pi k}{N}\right)\right)}^{2}\left(1-\mathrm{cos}\left(2{\omega }_{k}t\right)\right)+2{t}^{2}{\omega }_{k}^{2}{\left({\alpha }_{k}\mathrm{cos}\left(\frac{2\pi k}{N}\right)+{\beta }_{k}\mathrm{sin}\left(\frac{2\pi k}{N}\right)\right)}^{2}}{2{\omega }_{k}^{4}}.\end{eqnarray} \tag{ 26 }$$

For long interaction times, i.e., $t\to \infty$ the term quadratic in t in equation (26) completely dominates. Moreover, for N large the sum in equation (25) can be replaced, to a good approximation, by an integral resulting in the following simple expression for the variance

$$\begin{eqnarray}{\Delta }^{2}{{\mathcal{O}}}_{J}{(J,B,t)}_{\mathrm{max}} & = & {N}^{2}{t}^{2}\;G\left(\displaystyle \frac{B}{J}\right)\\ G(g) & = & {\left(\displaystyle \frac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }\sqrt{\displaystyle \frac{{(1+g\mathrm{cos}(x))}^{2}}{1+{g}^{2}+2\;g\mathrm{cos}(x)}}{\rm{d}}x\right)}^{2}.\end{eqnarray} \tag{ 27 }$$

Thus, the variance scales as ${N}^{2}{t}^{2}$, i.e., at the Heisenberg limit, up to an overall constant factor that only depends on the ratio $B/J$. The function G is plotted in figure 1. Notice that $G(B/J)$ exhibits a phase transition at J = B.

Zoom In Zoom Out Reset image size

We now proceed to determine the optimal precision with which one can estimate J, knowing B, if we restrict the input state of N systems to be a product state. As it is not immediately evident what separable states look like in the basis that diagonalizes the operator ${{\mathcal{O}}}_{J}(J,B,t)$, we work directly with equation (10). However, due to the form of equation (10) it is highly non-trivial to perform an analytical optimization over all possible separable states of N qubits. To that end we perform a brute-force optimization of the variance of equation (25) over all possible input product states for up to N = 11 qubits. The results are shown in figure 2 where one can already see a difference in scaling between the optimal product state strategy and the optimal quantum strategy even for small system sizes. Moreover, within a high margin of certainty the scaling of the QFI using the optimal product input state is at most linear in N.

Zoom In Zoom Out Reset image size

3.2. Estimating the field strength

We now proceed to estimate the magnetic field B, given we know J exactly. The procedure is identical to that of section 3.1. Writing ${H}_{2}={\displaystyle \sum }_{i=1}^{N}{\sigma }_{z}^{(i)}$ in terms of the fermionic operators b_k we obtain

$$\begin{eqnarray}{H}_{2}=\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\left({{lb}}_{k}^{\dagger }\;\;{b}_{N-k}\right)\left(\begin{array}{cc}1 & 0\\ 0 & -1\end{array}\right)\left(\begin{array}{l}{b}_{k}\\ {b}_{N-k}^{\dagger }\end{array}\right).\end{eqnarray} \tag{ 28 }$$

Substituting equations (20) and (28) into equation (10) and performing the integration over s yields

$$\begin{eqnarray}{{\mathcal{O}}}_{B}(J,B,t)=\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\left(\begin{array}{l}{b}_{k}^{\dagger }{b}_{N-k}\end{array}\right)\left(\begin{array}{ll}{A}_{k} & {B}_{k}\\ {B}_{k}^{*} & -{A}_{k}\end{array}\right)\left(\begin{array}{l}{b}_{k}\\ {b}_{N-k}^{\dagger }\end{array}\right),\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}{A}_{k} & = & \displaystyle \frac{{\alpha }_{k}^{2}t}{{\omega }_{k}^{2}}+\displaystyle \frac{{\beta }_{k}^{2}\mathrm{sin}(2{\omega }_{k}t)}{2{\omega }_{k}^{2}}\\ {B}_{k} & = & \displaystyle \frac{{\rm{i}}{\alpha }_{k}{\beta }_{k}t}{{\omega }_{k}^{2}}+\displaystyle \frac{2{\beta }_{k}{\omega }_{k}{\mathrm{sin}}^{2}({\omega }_{k}t)-{\rm{i}}{\alpha }_{k}{\beta }_{k}\mathrm{sin}(2{\omega }_{k}t)}{2{\omega }_{k}^{2}},\end{eqnarray} \tag{ 30 }$$

where we have again separated the linear and oscillatory parts in t. As equation (29) is of the same form as equation (16) it can be brought to the diagonal form

$$\begin{eqnarray}&&{{\mathcal{O}}}_{B}(J,B,t)=2\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\sqrt{{A}_{k}^{2}+{\left|{B}_{k}\right|}^{2}}\quad {f}_{k}^{\dagger }{f}_{k}\end{eqnarray} \tag{ 31 }$$

by a suitable Bogoliubov transformation (see equation (17)). The optimal achievable precision for estimating the field strength B is given by the optimal variance of ${{\mathcal{O}}}_{B}(J,B,t)$ in equation (31) which can be easily computed to be

$$\begin{eqnarray}&&{\Delta }^{2}{{\mathcal{O}}}_{B}{(J,B,t)}_{\mathrm{max}}={\left(\underset{k=0}{\overset{N-1}{\displaystyle \sum }}\sqrt{{A}_{k}^{2}+{\left|{B}_{k}\right|}^{2}}\right)}^{2},\end{eqnarray} \tag{ 32 }$$

and is again achieved by the equally weighted superposition of the vacuum state and the state where all N fermionic modes are occupied. We note that because the Bogoliubov transformation diagonalizing equation (29) explicitly depends on the coefficients of equation (30) the fermions described by modes f_k here are different than those described by modes d_k and c_k in equations (14) and (24) respectively. As the sum in equation (32) includes N summands, the maximum variance of ${{\mathcal{O}}}_{B}(J,B,t)$ scales as N² up to some factor.

Using equation (30) one can easily show that

$$\begin{eqnarray}&&{A}_{k}^{2}+| {B}_{k}{| }^{2}=\displaystyle \frac{{\beta }_{k}^{2}\left(1-\mathrm{cos}\left(2{\omega }_{k}t\right)\right)+2{t}^{2}{\alpha }_{k}^{2}{\omega }_{k}^{2}}{2{\omega }_{k}^{2}}.\end{eqnarray} \tag{ 33 }$$

For $t\to \infty$ and N very large the quadratic term in t dominates and the maximal variance can be given explicitly as

$$\begin{eqnarray}&&{\Delta }^{2}{{\mathcal{O}}}_{B}{(J,B,t)}_{\mathrm{max}}={N}^{2}{t}^{2}G(J/B),\end{eqnarray} \tag{ 34 }$$

where G is the function given in equation (27) and figure 1.

That the optimal precision in estimating either B or J is asymptotically given by the same expression can be understood via the duality of the one-dimensional Ising chain [43]. The duality is associated with whether one adopts the spin degrees of freedom on the chain or the kink degrees of freedom—associated with the links of the chain—as qubits. If adjacent spins are parallel then there is no kink, else there is a kink. One can then recast the Ising Hamiltonian of equation (1) in terms of the kink degrees of freedom as

$$\begin{eqnarray}&&H(J,B)=J\underset{i=1}{\overset{N}{\displaystyle \sum }}{\tau }_{x}^{(i)}+B\underset{i=1}{\overset{N-1}{\displaystyle \sum }}{\tau }_{z}^{(i)}{\tau }_{z}^{(i+1)},\end{eqnarray} \tag{ 35 }$$

where $\{{\tau }_{\alpha }:\alpha \in (x,y,z)\}$ have the same commutation relations as the Pauli matrices for the spin degree of freedom. Up to a basis change, the Hamiltonian in equation (35) is identical to that of equation (1), except that the roles of J and B are reversed. This is the reason why for the estimation of J the constant factor is given by $G(B/J)$ whereas for the estimation of the field it is given by $G(J/B)$. We note that the duality works well for spins and kinks in the bulk of the one-dimensional chain but is problematic with spins and kinks close to the edges of the chain. However, for large enough N the effects at the boundaries of the chain can be safely neglected. Also note that the quaternionic representation of the operators in the kink degrees of freedom is different than that of spins in one crucial way. The degeneracies of eigenstates in the kink representation are not the same as the ones for spins. Indeed, one can easily show that the ordered spectrum of eigenvalues of ${\displaystyle \sum }_{i=1}^{N}{\sigma }_{z}^{(i)}$ is $\{{\lambda }_{m}=\frac{N}{2}-m| \;m\in (0,\ldots ,N)\}$ with eigenvalue m having a degeneracy of $\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{N}{m}\right)$, whereas the ordered spectrum of ${\displaystyle \sum }_{i=1}^{N}{\tau }_{x}^{(i)}$ is given by $\{{\lambda }_{m}=\frac{N}{2}-(m+\frac{1}{2})| \;m\in (0,\ldots ,N-1)\}$ with corresponding degeneracies given by $2\displaystyle \left(\displaystyle \genfrac{}{}{0em}{}{N-1}{m}\right)$. Notice however that the duality is not exact for finite N.

This difference in the number of degenerate states can be exploited, in the case of finite N, to improve the precision with which one can estimate parameters in the Bayesian estimation regime.

The optimal precision for estimating B, given we know J, using a separable strategy is again numerically calculated for up to N = 11 qubits with the results shown in figure 3. Just as in the case of estimating J, one can already observe a difference in scaling of the QFI between the optimal product state strategy and the corresponding optimal quantum strategy. Moreover, within a high margin of certainty, the scaling of the QFI using the optimal product state is at most linear in N.

Zoom In Zoom Out Reset image size

In the previous section we showed that the optimal precision in estimating either the interaction strength or magnetic field is achieved by states that are linear superpositions of the vacuum and all N sites occupied by fermions. Whereas such states can be prepared efficiently, i.e., with a quantum circuit that grown polynomially with the number of qubits [44], preparing such states in practice may still be quite challenging due to the optimal states dependence on both time, and J (B). An important question, then, is whether there exist states that are easy to prepare and yield Heisenberg scaling in precision for all time and all values of J (B). For example in ion-trap set-ups one can prepare the GHZ state using a single Sø rensen-Mø lmer gate [45]. In this section we analytically determine the performance of GHZ-type states for estimating either the interaction strength or magnetic field. We will show that the GHZ-type states allow for Heisenberg scaling in precision for all values of J (B), with only a constant factor difference from the optimal precision achievable.

Let us first determine the ultimate precision of estimating B using the GHZ state

$$\begin{eqnarray}&&| \mathrm{GHZ}\rangle =\displaystyle \frac{1}{\sqrt{2}}\left({| 0\rangle }^{\otimes N}+{| 1\rangle }^{\otimes N}\right).\end{eqnarray} \tag{ 36 }$$

The GHZ state yields the ultimate precision in estimating B for the case where the latter is imprinted via the unitary operator $U=\mathrm{exp}(-{\rm{i}}{{BtH}}_{2})$. In order to calculate its performance for the Ising Hamiltonian of equation (1) we will use the Jordan–Wigner transformation and the discrete Fourier transform to express the GHZ state in terms of the fermionic creation operators ${b}_{k}^{\dagger }$.

We begin by noting that the GHZ state can be written in terms of the fermionic operators $\{a,{a}^{\dagger }\}$ as

$$\begin{eqnarray}&&| \mathrm{GHZ}\rangle =\displaystyle \frac{1}{\sqrt{2}}\left({\mathbb{1}}+\underset{{\rm{k}}=1}{\overset{{\rm{N}}}{\displaystyle \prod }}{{\rm{a}}}_{k}^{\dagger }\right)| {\bf{0}}\rangle .\end{eqnarray} \tag{ 37 }$$

Using equation (15) to transform the fermionic operators, $\{{a}_{k}^{\dagger }\}$ to $\{{b}_{k}^{\dagger }\}$, and noting that the Fourier transform leaves the vacuum invariant, yields

$$\begin{eqnarray}&&| \mathrm{GHZ}\rangle =\displaystyle \frac{1}{\sqrt{2}}\left(| {\bf{0}}\rangle +{{\rm{e}}}^{{\rm{i}}\phi }| {\bf{1}}\rangle \right),\end{eqnarray} \tag{ 38 }$$

where ϕ is a relative phase that arises due to the Fourier transform.

Using equation (38) to calculate the variance of the operator in equation (29), one finds that

$$\begin{eqnarray}\langle \mathrm{GHZ}| {{\mathcal{O}}}_{B}({\rm{J}},{\rm{B}},{\rm{t}})| \mathrm{GHZ}\rangle & = & 0\\ \langle \mathrm{GHZ}| {{\mathcal{O}}}_{B}^{2}({\rm{J}},{\rm{B}},{\rm{t}})| \mathrm{GHZ}\rangle & = & {\left(\underset{k=0}{\overset{N-1}{\displaystyle \sum }}{A}_{k}\right)}^{2},\end{eqnarray} \tag{ 39 }$$

and consequently ${\Delta }^{2}{{\mathcal{O}}}_{B}(J,B,t)=\langle \mathrm{GHZ}| {{\mathcal{O}}}_{B}^{2}({\rm{J}},{\rm{B}},{\rm{t}})| \mathrm{GHZ}\rangle$. Using equation (30), denoting $\theta =\frac{2\pi k}{N}$, and letting $N\to \infty$ the variance of ${\Delta }^{2}{{\mathcal{O}}}_{B}(J,B,t)$ explicitly reads

$$\begin{eqnarray}&&{\Delta }^{2}{{\mathcal{O}}}_{B}(J,B,t)=\displaystyle \frac{{N}^{2}}{4{\pi }^{2}}{\left({\displaystyle \int }_{0}^{2\pi }\displaystyle \frac{{(1+g\mathrm{cos}\theta )}^{2}}{1+{g}^{2}+2\;g\mathrm{cos}\theta }{\rm{d}}\theta \right)}^{2},\end{eqnarray} \tag{ 40 }$$

where g = J/B. From equation (40) it follows that the variance of ${{\mathcal{O}}}_{B}(J,B,t)$ with respect to the GHZ state always scales quadratically with N, i.e., at the Heisenberg limit, with a prefactor F(g) slightly lower than that for the optimal states as shown in figure 1.

We now determine the precision with which one can estimate J using the state

$$\begin{eqnarray}&&| \psi \rangle =\left(\displaystyle \frac{{| 00\rangle }^{\otimes \frac{N}{2}}+{| 01\rangle }^{\otimes \frac{N}{2}}}{\sqrt{2}}\right),\end{eqnarray} \tag{ 41 }$$

which is, up to a basis change, the state that yields the ultimate precision in estimating J for the case where the latter is imprinted via the unitary dynamics $U=\mathrm{exp}(-{\rm{i}}{{JH}}_{1})$. We note that since H₁ has a doubly degenerate eigenspace for both its minimum and maximum eigenvalue, the optimal state for estimating J is not unique.

If we express this state in terms of the kink degrees of freedom, as opposed to the spin degrees of freedom then what we obtain is, up to an overall Hadamard transformation, the GHZ state, i.e, a linear superposition of zero kinks, and N kinks. As the Ising Hamiltonian expressed in the kink degrees of freedom is given by equation (35), up to an overall Hadamard transformation, it follows that the variance of ${{\mathcal{O}}}_{J}(J,B,t)$ with respect to the state $| \psi \rangle$ is given by equation (40), with g = B/J.

In this work we investigated precision limits for noiseless quantum metrology in the presence of parameter dependent Hamiltonians, and in particular the Ising Hamiltonian. We showed that the ultimate limit in estimating the interaction or magnetic field strength scales quadratically with the number of probe systems N, i.e., at the Heisenberg limit, and that the states that achieve this precision are linear superposition of the vacuum and fully occupied states of N free fermions. Moreover, whereas the Ising Hamiltonian generates entanglement this entanglement does not help in boosting the precision scaling with respect to N that can be achieved with product states. In addition, we showed that the achievable precision in estimating either J or B for the Ising Hamiltonian using GHZ-type states also scales at the Heisenberg limit.

Whereas we have shown that the entanglement generating properties of the Ising Hamiltonian do not boost the precision in estimation for product states, it may still be the case that we can exploit this property of the Ising Hamiltonian to reduce the amount of entanglement required in the initial input state of the N probes. This would be of high interest for practical realizations of quantum metrology where the creation of highly entangled states of N systems remains a challenge.

In addition, our analysis deals with optimal states and bounds in the absence of noise. It would be interesting to investigate the achievable precision bounds in the presence of several physical noise models, such as uncorrelated dephasing or depolarizing noise, as well as spatial and temporal correlated noise. Furthermore, it would be interesting to determine which of these types of noise can we readily combat via the use of error-correcting techniques, or by dynamical decoupling [46–48].

The authors would like to thank the anonymous referee for his valuable comments and suggestions during the reviewing process of this article. This work was supported by the Austrian Science Fund (FWF): P24273-N16 and the Swiss National Science Foundation grant P2GEP2_151964.