Within the currently emerging quantum era, quantum metrology [1, 2, 3, 4, 5, 6, 7, 8] has proven itself to be a powerful tool for the accurate estimation of even very small physical parameters, such as gravitational wave amplitudes [9], or to limit the attainable measurement accuracy of fundamental constants like the speed of light [10]. As a result, quantum metrology promises to revolutionize the existing technologies of measurement.

While quantum metrology has frequently been employed in the quantum optical context [11, 12, 13, 14], more recently the corresponding experiments and theory are also exploring coherent matter waves cf., e.g., Refs. [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. Photons freely propagating in the quantum vacuum are to a very good approximation noninteracting particles and are well described by plane waves of definite momentum. Matter waves forming Bose-Einstein condensates at very low temperatures are, however, interacting by the scattering of their elementary atomic or molecular constituents, and are spatially confined (trapped) by arbitrary scalar potentials. In what follows, we show that the full self-consistency of the quantum many-body evolution of such a system needs in general to be taken into account, to yield reliable parameter estimation. We demonstrate that the interplay of Fock space amplitudes and time-dependent field operator modes ($\coloneqq$ orbitals), representing self-consistent many-body evolution, is crucial. This interplay is not obtained when fixing the orbitals’ shape, thereby significantly restricting the associated Hilbert space.

We take as an archetypical model system and for concreteness a tilted double well, where the parameter to be estimated is the linear slope $p_{4}$ which could, e.g., represent exposing the gas to constant gravitational acceleration (see Fig. 1). To facilitate comparison with conventional interferometry, we operate within a (continuously monitored) two-mode approximation (TMA), corresponding to two interferometric arms. We consider a simple (Mach-Zehnder type) experiment which counts the final number of particles on the left and right. It is demonstrated that while a non-self-consistent evolution yields a null result for $p_{4}$ [zero classical Fisher information (CFI)], a self-consistent quantum many-body evolution gives finite CFI, enabling $p_{4}$ estimation. We therefore conclusively show that self-consistency is crucial for the correct interpretation of parameter estimation data in a trapped interacting many-body system. We note that Heisenberg-limit number scaling of estimation precision with $1/N$ [instead of the shot noise (standard quantum) limit $\propto 1/\sqrt{N}$], for cat state distributions in Fock space [28], or an enhanced $N$-scaling relying on $k$-body interactions [31, 32], are not our aim in the present work. The focus is on the fundamental imprint of the self-consistency of many-body evolution on the accuracy of parameter estimation from the metrological protocol.

To determine the many-body evolution self-consistently, we perform a multiconfigurational time-dependent Hartree (MCTDH) analysis. This is a well tested method to describe the dynamics of spatially confined quantum systems, such as the one under investigation here, for which self-consistency is crucial [33, 34, 35, 36, 37]. The general $N$-body state reads

where $\sum_{\vec{n}}|C_{\vec{n}}|^{2}=1$ for state normalization and $\vec{n}$ denotes the set of occupation numbers $\{n_{i}\,|\,i=1,2,\cdots,M\}$ in each mode (orbital), with $\sum_{i=1}^{M}n_{i}=N$. The time-dependent Fock basis state $|\vec{n}(t)\rangle$ indicates that the orbitals change in time as a result of finite $M$ (see, Fig. 1). Their dynamics follows the system of nonlinear coupled integrodifferential Eqs. (S1) in the supplement [38]. Numerical solution enables the determination of Fock space coefficients $C_{\vec{n}}$ and orbitals in a self-consistent manner at any time [39].

The quantum metrological approach to parameter estimation, see for example Refs. [40, 41, 42, 43, 44], proceeds essentially as follows. An initial state $|\psi\rangle$ experiences a dynamical evolution, e.g., $e^{-i\hat{H}_{X}t}$, during the time $t$ and the final state $|\psi_{X}\rangle$ contains the information of the parameter $X$. One chooses an appropriate measurement on $|\psi_{X}\rangle$ to estimate $X$. Previous studies on quantum metrology with ultracold atoms [15, 16, 42, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30] have focused on the coefficients $C_{\vec{n}}(X;t)$ in Eq. (1), and have calculated the quantum Fisher information (QFI) ${\mathfrak{F}}_{X}$ from $C_{\vec{n}}(X;t)$ only. However, since the orbitals also evolve by Eqs. (S1) and the time evolution relies on $X$, the $|\vec{n}(X;t)\rangle$ must be considered in the calculation of both the QFI and the CFI. To evaluate the sensitivity of a quantum mechanical state to a parameter change thus requires full exploitation of the information encoded in the state. Here, we make full use of the parameter dependence of the state, reflected in both coefficients and orbitals. As a result, we establish numerically exact many-body parameter estimation for trapped interacting quantum gases.

A scalar bosonic gas with contact interactions, trapped in a quasi-one-dimensional (quasi-1D) double-well potential, is described by the Hamiltonian

where $V$ is the trap potential. We render $\hat{H}$ in a dimensionless form, fixing a unit length $L$; the unit of time is then $L^{2}$ for $\hbar=m=1$ ¹¹1In ${}^{87}\text{Rb}$, $L=1\,\,\mu$m yields a time unit $\Delta t=1.366$  msec.. The quasi-1D interaction coupling $g$ is controllable by either Feshbach or geometric scattering resonances [46]. Exact solutions of the Schrödinger equation associated to Eq. (2) exist for homogeneous gases under periodic boundary conditions or in a box trap; their metrology was explored in [47].

We assume a trap potential of the form $V(x)=V_{\rm dw}(x)+p_{4}x=\frac{1}{2}p_{1}x^{2}+p_{2}\exp[-x^{2}/(2\,p_{3}^{2})]+p_{4}\,x$, where $p_{4}=0$ initially, cf. Fig. 1. The remaining parameters, $p_{1}=0.5$, $p_{2}=50$, and $p_{3}=1$, are fixed throughout the evolution. The two lowest-energy single-particle states are symmetric and antisymmetric with respect to the origin, respectively, and their addition and subtraction results in two well-localized orbitals: left $\phi_{\rm L}(x)$ and right $\phi_{\rm R}(x)$, see Fig. 1. These orbitals, approximate ground states of each well, furnish the two modes [48].

The estimation process of $p_{4}$ proceeds with an initial state in the form of Eq. (1). With $\phi_{\rm L}(x)$ and $\phi_{\rm R}(x)$, we employ two coefficient distributions: A NOON (cat) state, $|\psi_{0}\rangle=(|N,0\rangle+|0,N\rangle)/\sqrt{2}$, and a spin-coherent state $|\psi_{0}\rangle=\sum_{k=0}^{N}\sqrt{\frac{N!}{k!\,(N-k)!}}\cos^{N-k}(\frac{\pi}{4})\sin^{k}(\frac{\pi}{4})|N-k,k\rangle$. Then a small $p_{4}$ tilt is switched on. The state then evolves according to Eqs. (S1); finally, the number of particles in each well is measured and $p_{4}$ is estimated from the set of outcomes.

As the relative interaction strength $gN$ (typical ratio of interaction over single-particle energies), increases, more modes than two are required to correctly reproduce the many-body dynamics [37]. Keeping $gN$ fixed when taking the limit $N\rightarrow\infty$ has been demonstrated to reproduce the Gross-Pitaevskii ground state energy [49, 50]. We thus keep in the following $gN$ fixed when varying $N$, to remain close to the TMA. In order to adequately compare self-consistent (SC) evolution to conventional SU(2) two-mode interferometry (TMI), which operates with changing Fock space coefficients only, we maintain TMA validity throughout the time evolution. This can be assessed by evaluating $\rho_{\rm tm}\coloneqq(\rho_{1}+\rho_{2})/N$, cf. Fig. 2, where $\rho_{j}$ is the $j$th largest eigenvalue of the reduced one-body density matrix, $\rho^{(1)}(x,x^{\prime};t)$; see (S2) in the supplement [38]. After diagonalization, $\rho^{(1)}(x,x^{\prime})=\sum_{j}\rho_{j}(t)\,\phi_{j}^{({\rm no})\ast}(x^{\prime},t)\,\phi_{j}^{({\rm no})}(x,t)$, with the natural orbitals $\{\phi_{j}^{({\rm no})}(x,t)|\,j=1,2,\cdots\}$. When a single $\rho_{1}=O(N)$ (in the formal limit $N\rightarrow\infty$), a Bose-Einstein condensate is obtained. For several $\rho_{j}$ of $O(N)$, we have a fragmented condensate [51]. The validity of the TMA (two-fold fragmented condensate) depends on whether $\rho_{1}\simeq\rho_{2}\simeq N/2$ and $\rho_{\rm tm}\simeq O(1)$ hold. In the supplement [38], we provide further evidence for the appropriateness of using the TMA, by demonstrating the rapid convergence of our results with increasing $M$.

We define an initial state using four orbitals ($M=4$), which are including $\phi_{\rm L}(x)$ and $\phi_{\rm R}(x)$, and then monitor the natural occupations $\rho_{j}$ under self-consistent evolution at nonzero $p_{4}$, for both cat and spin-coherent states. Fig. 2 shows $\rho_{\rm tm}$, where we observe it is close to unity. Also, both $\rho_{1}$ and $\rho_{2}$ are macroscopically occupied during the evolution, with negligible occupations $\rho_{3}$ and $\rho_{4}$. This however also depends on the parameter regime used. When $gN=0.1$, two modes are sufficient, but when $gN=1$, $\rho_{\rm tm}$ discernibly dips below unity. Increasing $p_{4}$ further, the TMA fails. An appropriate regime of parameters where $\rho_{\rm tm}\simeq 1$ is obtained by fixing $gN=0.1$ and $p_{4}=0.1$. Also, even though natural orbitals are used in the discussion above, we can apply the two-mode criterion to the left/right orbitals or their time-evolved forms, i.e., $\phi_{1}(x,t)$ and $\phi_{2}(x,t)$; there always exists a unitary transformation such that $\phi_{j}(x,t)=\sum_{jk}U_{jk}\,\phi_{k}^{({\rm no})}(x,t)$.

Expanding the second-quantized form of Eq. (2) with $\hat{\Psi}(x)=\hat{b}_{\rm L}\,\phi_{\rm L}(x)+\hat{b}_{\rm R}\,\phi_{\rm R}(x)$, a two-site single-band Bose-Hubbard Hamiltonian is obtained. In terms of the usual SU(2) Pauli matrices $\hat{J}_{x}=\frac{1}{2}(\hat{b}_{\rm L}\hat{b}_{\rm R}+\hat{b}^{\dagger}_{\rm R}\hat{b}_{\rm L})$ and $\hat{J}_{z}=\frac{1}{2}(\hat{b}_{\rm L}\hat{b}_{\rm L}-\hat{b}^{\dagger}_{\rm R}\hat{b}_{\rm R})$,

where $\tau$ is tunneling amplitude, $\epsilon$ an energy offset between wells, and $U\propto g$ an interaction coupling, all depending on integrals involving the two orbitals. The implementation of quantum metrological protocols using the above Hamiltonian was carried out, e.g., in Refs. [19, 24, 44]: An initial state $|\psi_{0}\rangle$, as defined by the distribution of coefficients $C_{\vec{n}}$, evolves as $\exp(-i\hat{H}t)|\psi_{0}\rangle$, and the parameter of interest, e.g., $\epsilon$, is estimated from the population imbalance [19, 24]. In our setup, the strong barrier renders the initial $\tau$ exponentially small compared to $\epsilon$ and $U$, and $\epsilon\neq 0$, as the symmetry of $V(x)$ is broken by $p_{4}$. Hence the TMI time evolution operator is, to very good accuracy, $\exp(-i(\epsilon\hat{J}_{z}+U\hat{J}_{z}^{2})t)$, and the QFI can be analytically calculated by $\mathfrak{F}_{\epsilon}=4\,\langle\psi_{0}|(\Delta\hat{J}_{z})^{2}|\psi_{0}\rangle$. When the cat state is used, $\mathfrak{F}_{\epsilon}=N^{2}\,t^{2}$, which is denoted as the Heisenberg limit. For the spin-coherent state, $\mathfrak{F}_{\epsilon}=N\sin^{2}(\theta)\,t^{2}$, representing the shot-noise (standard quantum) limit. Any nonzero $\tau$ deteriorates the $N$-scaling of $\mathfrak{F}_{\epsilon}$, confirmed by numerically calculating the QFI [19]. Note that here only the change of Fock space coefficients has been considered, while the orbitals are fixed in TMI. Because of the latter fact, the still exponentially small $\tau$ and $U$ are kept constant during the evolution, and $\epsilon$ is abruptly switched on at $t=0$. In our setting, $\epsilon\in(-0.7,-0.6)$, and $U\in(0.002,0.03)$, with concrete values determined by $gN$ and $N$, on which in turn the initial orbitals $\phi_{\rm L}(x)$ and $\phi_{\rm R}(x)$ depend. Recall that our target parameter is $p_{4}$, not $\epsilon$, thus by using the chain rule, $\mathfrak{F}_{p_{4}}=\mathfrak{F}_{\epsilon}\times(h_{1}-h_{2})^{2}$, where the single-particle energies $h_{i}\coloneqq\int dx\,\phi_{i}^{\ast}(x)\left[-\frac{1}{2}\frac{\partial^{2}}{\partial x^{2}}+V(x)\right]\phi_{i}(x)$.

First, we compare the QFIs of the SC approach and the conventional TMI. The QFI with respect to a given parameter $X$ inscribed onto a pure state $|\psi_{X}\rangle$, is $\mathfrak{F}_{X}=4\big{(}\langle\partial_{X}\psi_{X}|\partial_{X}\psi_{X}\rangle-|\langle\psi_{X}|\partial_{X}\psi_{X}\rangle|^{2}\big{)}$, and insertion of Eq. (1) into $|\psi_{X}\rangle$ gives Eq. (S4), which facilitates calculation of the QFI using the ingredients of MCTDH from the expansion in Eq. (1), also cf. [38]. The first row in Fig. 3 shows QFI versus time $t$ and particle number $N$, respectively. For each initial state, the SC method reproduces very well the QFI predicted by TMI. The influence of self-consistency thus plays a subdominant role for $\mathfrak{F}_{X}$. The latter completely depends on the final state, and self-consistency (changing orbitals) essentially represents fitting that state more exactly. When the bosons weakly interact ($gN=0.1$) and the disturbance to the system is small ($p_{4}=0.1$), conventional TMI therefore approximates well the QFI [52].

We now turn to the CFI [53, 54], associated to a concrete measurement, for which, as we show, the impact of self-consistency becomes manifest. Counting the number of bosons in each well constitutes our measurement, so that the CFI is defined as

where $P(\vec{\mathbf{n}}|p_{4})$ is the probability distribution (likelihood) $\vec{\mathbf{n}}=(n_{\rm L},n_{\rm R})$, given $p_{4}$, and $n_{\rm L}$ and $n_{\rm R}$ are the numbers of particles that reside in the left and the right well, respectively. An appropriate $P(\vec{\mathbf{n}}|p_{4})$ has to be constructed to calculate the CFI and the conventional TMI approach considers that $P(\vec{\mathbf{n}}|p_{4})\coloneqq|\langle\vec{n}(t)|\Psi(t)\rangle|^{2}$. Then, the CFI always exactly vanishes, irrespective of the initial state, as the Hamiltonian contains only $\hat{J}_{z}$ and $\hat{J}_{z}^{2}$. TMI time evolution therefore just changes the phases of the coefficients: $C_{k}(t)=\exp(-i\epsilon\frac{N-2k}{2}t)\,\exp(-iU(\frac{N-2k}{2})^{2}t)\,C_{k}(0)$, where the state $|\psi(t)\rangle=\sum_{k=0}^{N}C_{k}(t)|N-k,k\rangle=\exp(-i(\epsilon\hat{J}_{z}+U\hat{J}_{z}^{2})t)|\psi(0)\rangle$ and $\hat{J}_{z}|N-k,k\rangle=\frac{N-2k}{2}|N-k,k\rangle$. Then $P(\vec{\mathbf{n}}=(N-k,k)|p_{4})=|C_{k}(t)|^{2}=|C_{k}(0)|^{2}$, independent of $p_{4}$, yielding vanishing CFI, cf. [38].

On the other hand, when both orbitals and Fock space coefficients evolve in the SC framework, the initial interpretation of the orbitals cannot be maintained throughout the time evolution. Initially well-localized orbitals, i.e., $\phi_{\rm L}(x)=\phi_{1}(x,0)$ and $\phi_{\rm R}(x)=\phi_{2}(x,0)$, correspond to particles being found in the left and right well, respectively. However, the orbitals change with time during the many-body evolution, so $\phi_{1}(x,t)$ and $\phi_{2}(x,t)$, do not necessarily imply left or right localization at the time of measurement. Therefore, simply computing $P(\vec{\mathbf{n}}|p_{4})=|\langle\vec{n}(t)|\Psi(t)\rangle|^{2}$, as in the TMI approach, is not applicable. In the bosonic field operator $\hat{\Psi}(x)=\sum_{j}\hat{b}_{j}(t)\,\phi_{j}(x,t)$, it is clear that the bosonic annihilation operator $\hat{b}_{j}(t)$ corresponds to the time-evolving orbital $\phi_{j}(x,t)$, which delocalizes with increasing $t$, see Fig. 1. Thus the Fock state in Eq. (1) $|\vec{n}(t)\rangle=\frac{\big{(}\hat{b}_{1}^{\dagger}(t)\big{)}^{n_{1}}\big{(}\hat{b}_{2}^{\dagger}(t)\big{)}^{n_{2}}\cdots\big{(}\hat{b}_{M}^{\dagger}(t)\big{)}^{n_{M}}}{\sqrt{n_{1}!n_{2}!\cdots n_{M}!}}|0\rangle$ cannot by itself project the quantum state into any of $|\vec{\mathbf{n}}\rangle$ and $\langle\vec{n}(t)|\Psi(t)\rangle$ cannot be interpreted as probability amplitude for each measurement outcome as in TMI. In other words, $|\vec{n}(t)\rangle$ and $|\vec{\mathbf{n}}\rangle$, while initially identical, become different due to SC evolution. Hence we need to re-establish a connection between time-evolving orbitals and $|\vec{\mathbf{n}}\rangle$ that corresponds to a measurement outcome, resulting in the proper distribution $P(\vec{\mathbf{n}}|p_{4})=|\langle\vec{\mathbf{n}}|\Psi(t)\rangle|^{2}\neq|\langle\vec{n}(t)|\Psi(t)\rangle|^{2}$.

Within SC time evolution after a trap tilt, $\phi_{1}(x,t)$ and $\phi_{2}(x,t)$ remain well-localized in the case of a cat state. That is, $\phi_{1}(x,t)$ ($\phi_{2}(x,t)$) begins from $\phi_{1}(x,0)=\phi_{\rm L}(x)$ [$\phi_{2}(x,0)=\phi_{\rm R}(x)$] and their absolute value remains nearly identical except slightly wider (narrower) width and shorter (taller) height, respectively. For the spin-coherent state, however, $\phi_{1}(x,t)$ and $\phi_{2}(x,t)$ spread into opposite wells while they evolve under nonzero $p_{4}$; see Fig. 1. Then even when a particle resides in $\phi_{1}(x,t)$ or $\phi_{2}(x,t)$, to assign it to the left or right well is ambiguous. One can however still define mathematically “left” or “right” by integrating the orbitals from $-\infty$ to the center point $x=0$ of $V(x)$ and from the center point of $V(x)$ to $\infty$, respectively: $P_{\text{Left}}=\int_{-\infty}^{0}\!\!|\phi_{j}(x,t)|^{2}\,dx$ and $P_{\text{Right}}=\int_{0}^{\infty}\!\!|\phi_{j}(x,t)|^{2}\,dx$. One can then construct $P(\vec{\mathbf{n}}|p_{4})$ by computing the permanent of a special matrix composed from $P_{\text{Left}}$ and $P_{\text{Right}}$. The computable $N$ range is limited, though, due to rapidly increasing algorithmic complexity when calculating a matrix permanent, see [38].

The second row in Fig. 3 shows the CFI versus $t$ and $N$. Fixed orbitals results in vanishing CFI, as expected. Also, even within SC evolution, the cat state shows almost vanishing CFI, which is attributed to the fact that the orbitals stay localized in each well during the whole evolution time and the probabilities, i.e., $P_{\text{Left}}$ and $P_{\text{Right}}$, remain nearly constant (for small $p_{4}$). Thus under the given measurement the change of $P(\vec{\mathbf{n}}|p_{4})$ with respect to $p_{4}$ is negligible. However, the spin-coherent state displays a significant change in orbitals and increasing CFI during the early stage. Bottom right in Fig. 3 shows the $N$-scaling of the CFI; the SC approach with spin-coherent state shows an almost linearly increasing CFI. The complex fluctuation pattern appears because of the short time $t=1.77$ after a nonzero $p_{4}$ is suddenly applied, and for increasing $t$ these fluctuations smoothen out. Thus a SC approach may yield drastically different metrological predictions from a TMI based method.

Fig. 4 shows our primary result: The implementation of parameter estimation at the final stage of the metrological protocol. The maximum likelihood estimator (MLE) [6, 55] is used as a concrete example, because it asymptotically saturates the Cramér-Rao bound for an infinite number of measurements ($\nu\rightarrow\infty$, see also [38])

Here, the parameter $X$ is our $p_{4}$, $X_{\text{est}}$ is an estimator of $X$, and $\delta X_{\text{est}}:=X_{\text{est}}/\big{|}\partial\langle X_{\text{est}}\rangle_{X}/\partial X\big{|}-X$, with average $\langle\cdots\rangle_{X}$ taken with respect to $P(\vec{\mathbf{n}}|p_{4})$.

In Fig. 4, the likelihood function $P(\vec{\mathbf{n}}|p_{4})$ is displayed, supposing, for concreteness, that the measurement outcome is $\vec{\mathbf{n}}=(7,3)$, where $n_{\rm L}$ and $n_{\rm R}$ denote number of particles in left and right well, respectively. Red solid line shows the maximum of $P(\vec{\mathbf{n}}|p_{4})$ at $p_{4}\simeq 0.109$, thus the estimate of $p_{4}$, given $\vec{\mathbf{n}}=(7,3)$, is $X_{\text{est}}\simeq 0.109$. Similarly, every single outcome, $11$ in total, is connected to an estimate of $p_{4}$. The orange dashed line obtained by conventional TMI with a spin-coherent state stays flat, which means that the information given by the measurement outcome $\vec{\mathbf{n}}=(7,3)$ is zero, and an estimate of $p_{4}$, in accordance with $F_{p_{4}}=0$, is not possible.

When $\langle X_{\text{est}}\rangle_{X}=X$ holds, an estimator is unbiased. The present MLE is not unbiased except infinitesimally close to $p_{4}=0$. The larger $p_{4}$, the more bias its estimation acquires. For instance, in Fig. 4, the true value of $p_{4}$ is $0.1$ and the bias of the MLE is on the level of 10 %, $\langle X_{\text{est}}\rangle_{X=0.1}\simeq 0.11$. To compensate for this bias, we calibrate the MLE by obtaining $\langle X_{\text{est}}\rangle(X)$ from our MCTDH simulations, see [38]. Also, $\partial\langle X_{\text{est}}\rangle_{X}/\partial X|_{X=0.1}\simeq 0.6$, which then provides the mean-square deviation $\langle(\delta X_{\text{est}})^{2}\rangle_{X=0.1}$. The inset in Fig. 4 verifies that, asymptotically, $\langle(\delta X_{\text{est}})^{2}\rangle_{X=0.1}$ approaches the Cramér-Rao lower bound for large $\nu$.

In conclusion, we have found, using a self-consistent many-body approach, that many-body metrology utilizing interacting trapped bosons needs to conform to the final self-consistently computed many-body state. The QFI completely depends on the parameter dependence of the final state itself, and is thus relatively unaffected by self-consistency (in the weakly interacting regime). However, even in the latter regime, the CFI for a parameter estimation experiment is strongly affected by self-consistency due to its sensitive dependence on the orbitals’ time evolution. As a particularly notable example, fitting the outcome of a number-statistics experiment in a double well to conventional TMI gives a null result for estimating the slope parameter $p_{4}$. The SC approach we employ, however, enables $p_{4}$ estimation. Metrology with trapped ultracold quantum gases thus in general requires self-consistency of dynamical evolution, to correctly predict the estimation precision that can be accomplished in a given metrological protocol.

This work has been supported by the National Research Foundation of Korea under Grants No. 2017R1A2A2A05001422 and No. 2020R1A2C2008103.