Nondestructive Probing of Means, Variances, and Correlations of Ultracold-Atomic-System Densities via Qubit Impurities

We show how impurity atoms can measure moments of ultracold atomic gas densities, using the example of bosons in a one-dimensional lattice. This builds on a body of work regarding the probing of systems by measuring the dephasing of an immersed qubit. We show this dephasing is captured by a function resembling characteristic functions of probability theory, of which the derivatives at short times reveal moments of the system operator to which the qubit couples. For a qubit formed by an impurity atom, in a system of ultracold atoms, this operator can be the density of the system at the location of the impurity, and thus, means, variances, and correlations of the atomic densities are accessible.


I. INTRODUCTION
The impressive quantum control and versatility achieved in cold quantum gas experiments [1] exemplify their use for quantum technologies, such as quantum simulation [2,3] and quantum information processing [4]. Typical methods for measuring these systems, such as time-of-flight imaging of the momentum distribution [1,[5][6][7][8] and in situ imaging of density [9][10][11][12], require destroying the system. Each measurement is of a single shot of the system state, and hence finding averages requires repeated trapping and cooling of the system.
In this work, we propose an alternative method where trapped atomic impurities forming qubit probes are immersed in the system, entangling system and probes. Each qubit is dephased, and measuring this dephasing reveals information about the system [13][14][15][16][17][18]. The dephasing relates to the system operator to which the probes couple, and for typical lowenergy interactions between cold atoms this is the density of the system at the location of the impurity. We show that this enables measurement of not only the mean density but also density variances and correlations between different locations. This exploits the similarity of the dephasing function to characteristic functions of probability theory, allowing moments of the density to be determined. A major advantage of this scheme is that it is potentially nondestructive, leaving the system intact and allowing for repeated measurements.
As our example, we consider a gas of bosons in a onedimensional lattice, described by the Bose-Hubbard model. We demonstrate, by simulating the gas and impurity atoms for typically accessible experimental parameters, that our protocol can faithfully capture, even in the presence of errors, the behavior of density-related properties, over a range of phases of the system, thus providing an additional tool for measuring *  ultracold atoms. Moreover, since we develop the protocol generically, independent of the particular system, our work contributes to and furthers the ongoing development of methods for characterizing systems using coupled nonequilibrium dynamics, in cold atom and other setups, by observing the dephasing of a qubit in an environment [19][20][21][22][23][24][25][26].

A. Extracting the dephasing function
We consider a qubit with Hamiltonian H q = ω q |1 1| immersed in a system with Hamiltonian H sys . Here ω q is the energy difference between the qubit states {|0 ,|1 }. State |1 couples with strength κ to system operator H int , which we call the interaction Hamiltonian, such that the combined system is described by the Hamiltonian where I is the identity operator. We define the Hamiltonians H 0 = H sys and H 1 = H sys + κH int , describing the system evolution for each qubit state [27]. In a Ramsey-like scheme ( Fig. 1), we initialize the qubit at times t < 0 in the noninteracting state |0 and the system in state ρ sys . At t = 0, the qubit is suddenly switched to the state (|0 + |1 )/ √ 2 [28]. For times t 0, the system and qubit then evolve according to the Hamiltonian (1), such that at time t the qubit state is 1 .
Here we define the dephasing function (sometimes called the overlap function) L(t) = e iH 1 t e −iH 0 t , with expectations taken for the initial system state ρ sys . Values of the dephasing function L(t) are related to the expectations of the Pauli operators: σ x = Re[e iω q t L(t)] and σ y = Im[e iω q t L(t)].
The dephasing function has been investigated previously to study properties such as the orthogonality catastrophe in ultracold fermions [19,20], the Luttinger parameter [13], and superfluid excitations [18], as a method of thermometry [14][15][16][17], and to extract work distributions [23][24][25]. In this work, we proceed by investigating the derivatives of the dephasing function, motivated by the similar structure it shares with the characteristic function of a probability distribution [29].

B. Derivatives of the dephasing function
In a strong-coupling limit, where energies associated with the interaction Hamiltonian κH int are much larger than those of the system Hamiltonian, at sufficiently short times that e iH sys t ≈ 1 the dephasing function tends towards L strong (t) = e iκH int t . This is the characteristic function of the interaction Hamiltonian. Thus, like such functions in probability theory, from this all moments of H int can be obtained directly from derivatives: The drawback of this limit is that it may not be easily accessible in experiment for some systems, as it requires a large coupling strength to be engineered between system and probes. Increasing the coupling strength also decreases the time scales over which L(t) evolves, and hence requires sufficiently quick switching of the qubit state or ramping of κ, and sufficiently fine time resolution of the qubit measurements that the derivatives at t = 0 are resolvable. Positively, it is possible to extract the first two moments of H int for arbitrary κ. Considering the derivatives of L(t) directly, and where we have used that the commutator [H 0 ,H 1 ] is anti-Hermitian and thus has a purely imaginary expectation value. We now focus on these first two moments of the interaction Hamiltonian H int , as they can be extracted from derivatives of L(t) for any κ, allowing more experimental flexibility.

C. Estimation protocol and errors
At short times, the first and second derivatives Eqs. (2) and (3) dominate the imaginary and real parts of the dephasing function derivatives, respectively, and thus it is possible to extract the first two moments H int and H 2 int by fitting linear and quadratic functions to the initial behavior of Im[L(t)] and Re[L(t)], respectively, using that L(0) = 1. More precisely, we first obtain estimatesL(n t) of values L(n t) at discrete times n t, for integer n up to a maximum N corresponding to time t = N t. We then perform two least-squares estimations, minimizing with respect to α and β, and take the minimizing valuesᾱ and β as our estimates of κ H int and κ 2 H 2 int , respectively [30]. The choice of t should be such that it optimizes the number of points used in the fit without being so large that the imaginary and real parts of L(t) stop behaving approximately linearly and quadratically, respectively [31].
Each estimateL(n t) of L(n t) is obtained by estimating σ μ (μ = {x,y}), by averaging outcomes of N repeated measurements of σ μ . The resulting estimate will be unbiased, and for enough measurements will be Gaussian with a variance that can be calculated from L(n t), given by (1 + L)(1 − L)/ 4N , and hence attenuates to a few percent for a reasonable number of measurements. In our examples, we will treat this finite number of measurements as the main source of error in estimating L(n t), though other noise-based errors (such as imprecision in the time at which measurements are made and stochastic imperfections in the gate implementation) will behave in the same manner, and so can be considered equivalently. We will also implicitly account for the error arising from the discrete and finite nature of the time steps at which measurements are made, and in the backaction of the system-probe interaction on the system state. Other errors that we will not directly account for in our simulations are the finite time required to implement gates (which may be neglected when this occurs on time scales much shorter than that at which the measurements are made) and systematic imperfections in the gate implementation, which result in a deterministic multiplicative factor to the dephasing function [17], and may be accounted for by the corrective factor needed to ensure L(0) = 1.

III. IMPLEMENTATION IN ULTRACOLD ATOMIC SYSTEMS
We now consider a possible implementation of the protocol for probing ultracold (single species) atomic gases. A probe qubit is formed by two internal states of an impurity atom of a different species [20], trapped deeply such that its spatial wave function ψ q (x) is fixed [32]. The necessary gates may be applied to the impurity qubit using a Rabi laser pulse, and measurements performed using gates and state-dependent fluorescing. With interparticle interactions suppressed within the impurity gas, multiple probes can be active simultaneously, as illustrated in Fig. 2; such mixtures of quantum gases have already been achieved experimentally [33][34][35][36]. Repeat measurements do not require the gas and impurities to be retrapped, and can hence be nondestructive, though the measurement will perturb the system and is hence not nondemolition.
The impurity qubit and atoms comprising the system interact through s-wave scattering, potentially controlled through Feshbach resonances [8,37]; hence we have an interaction of the form where n(x) is the density of atoms in the system, andñ(x) is the system density course-grained over the impurity density. Our protocol thus allows probing of the moments of this coursegrained density. When the impurity is highly localized around x i ,ñ(x) is then given by the gas density at this point, n(x i ). In this case, our protocol will estimate the moments of the density at this point, n(x i ) and n(x i ) 2 [Eqs. (2) and (3)]. Additionally, since impurity atoms can be localized to regions smaller than the wavelength of light, the spatial resolution is potentially higher than for in situ imaging.
If the impurity is in a superposition of being localized to two distinct locations x i and x j (for example, by placing it in a lattice potential, and using tunneling as a beamsplitter operation [38]), thenñ(x) = [n(x i ) + n(x j )]/2, and our protocol will then estimate n(x i ) + n(x j ) and (n(x i ) + n(x j )) 2 . Using these and the previous results, an estimate for the correlation function n(x i )n(x j ) is obtained. Alternatively, this can also be achieved using two qubits localized at x i and x j with entangled internal states (|00 + |11 )/ √ 2, thus behaving as a single effective qubit with a density equal to the sum of that of the individual probes.

IV. SIMULATION FOR ULTRACOLD ATOMS IN A LATTICE
We now simulate an example application of the protocol to a Bose gas trapped in an optical lattice, obeying the onedimensional Bose-Hubbard Hamiltonian [2,39] Here b i (b † i ) annihilates (creates) a boson localized at site i (with number operator n i = b † i b i ), and i,j denotes a sum over nearest-neighbor sites. The energies J and U parametrize hopping between neighboring sites and on-site interactions, respectively. We simulate a one-dimensional system with 101 lattice sites and unit filling factor, and we examine the system for U/J in the interval 0.1 : 6, spanning the whole phase diagram of the model at this filling factor [2]. We are concerned entirely with investigating the ground state ρ sys of the system.
We use parameters that assume a realistic time-resolution in the measurements of the probe qubit dephasing. We choose κ = J , and t = 0.05J , allowing us to take N 20 measurement points in our fit. For typical J ∼ 10 3 Hz [39], this corresponds to a measurement interval t ∼ 50 μs [40]. With these parameters, we simulate application of our protocol to the Bose-Hubbard model, by calculating L(n t), adding stochastic noise to these values to simulate how estimatesL(n t) would be obtained in a real experiment, then estimate one-and two-site correlations of site occupation numbers from these simulated values. These estimated values are then compared to exactly calculated equivalents. For both parts, calculating dephasing function L(t) and ground-state expectation values, we use the time-evolving block decimation (TEBD) algorithm [41][42][43] (see the Appendix for details).
We begin our analysis by considering an impurity localized at some point x i near the central site i = 51 such that H int = n i . We plot the real and imaginary parts of the dephasing function L(t) for U/J = 3 in Fig. 3(a), together with their simulated noisy estimatesL(n t) obtained from N = 10 4 measurements of the qubit. We apply our fitting procedure to estimate the derivatives of the dephasing function and thus the moments of the interaction Hamiltonian, giving estimates of n i and n 2 i . In Fig. 3(b) we analyze the choice of t = N t to use in estimating the expected occupation n i , by averaging over many noisy trajectories. As expected, the accuracy initially increases with N , benefiting from an increase in the number of points used in the fit due to a corresponding decrease in random error. However, for larger t the accuracy decreases with N as times are included for which higher-order terms in L(t) beyond the linear fit begin to play a significant role and thus introduce a systematic error. The same is shown for n 2 i in Fig. 3(c), with the slight decrease in accuracy highlighting the increasing difficulty in estimating higher-order derivatives and thus moments. We find t J = 0.2 to be approximately optimal for estimating the moments.
With this t we assume the protocol is repeated with the impurities configured so as to probe H int = n j and n i + n j . We then estimate the mean n i and variance n 2 i − n i 2 of occupation at the central site, and the correlations C i,j = n i n j − n i n j between this central site and neighboring site j = i + 1. The results are plotted in Fig. 3(d). The markers show the value calculated from a single run of a noisy trajectory at each U/J , and the dashed lines show a cubic smoothing spline fit [44] over varying interaction strength for these data, while solid lines show the exact values. We observe that, though individual data points display a noticeable error, fitting over the whole parameter range the quantitative values can be obtained faithfully to a high degree of accuracy.

V. SUMMARY AND OUTLOOK
We have investigated a method of using a qubit to probe properties of a system to which it is coupled. In particular, we have introduced a way to calculate moments of a so-called interaction Hamiltonian through the derivatives of a dephasing function that can be obtained through measurements of solely the qubit state, even for a weakly interacting qubit. Further, we have discussed how this protocol could be implemented for ultracold atomic systems, to reveal properties of the gas density. We trialed the protocol by applying it successfully to simulations of atoms trapped in an optical lattice in one dimension. A direct application of this would be to probe phase transitions in the Bose-Hubbard model; for example, meanfield treatments, valid in higher dimensions, have found b i to be a suitable order parameter [45,46] for the superfluid-Mott insulator transition. The density fluctuations reveal whether this is nonzero, and hence also map out the phase diagram. Since our protocol is valid in any number of dimensions, this suggests that it could potentially be used as a way to nondestructively probe the transition. Our protocol could be extended beyond the current system, perhaps by manipulating the interaction between probe and system to allow other quantities to be probed, not just those related to density. For example, for bosons with spin, if κ is sensitive to the spin (i.e., κ ∝ S z ), then magnetization properties could also be probed. Alternatively, as the protocol has been derived generically, it could also be applied to other systems outside of cold atoms where similar probe-system interactions can be achieved, such as trapped ions [23] and nuclear magnetic resonance spins [25].

ACKNOWLEDGMENTS
The authors thank Stephen Clark and Dieter Jaksch for useful discussions. T.J.E. thanks the Engineering and Physical Sciences Research Council (Doctoral Training Account) for financial support. T.H.J. thanks the National Research Foundation and the Ministry of Education of Singapore for support.

APPENDIX: TEBD ALGORITHM
TEBD allows for the efficient classical simulation of the dynamical evolution of a pure quantum system on a one-dimensional lattice [41][42][43]. The algorithm operates by storing the state as a matrix product state (MPS) [47]. An arbitrary quantum state may require an exponential scaling in the bond dimension χ of this MPS. However, ground states and low-lying states (including short-time quenches from a ground state, as considered here) of local one-dimensional Hamiltonians may be very accurately represented by an MPS with a small bond dimension O(1 − 10 2 ) [48], resulting in a tractable representation of the state. The MPS is then evolved for each discrete time step δt according to a Hamiltonian H formed of single-site and two-site nearest-neighbor terms, after performing a Suzuki-Trotter decomposition [49] of the evolution operator e −iH δt . The algorithm can also be used to find the ground state of a Hamiltonian, by performing the evolution in imaginary time. Expectation values of single-site and two-site operators can be calculated for an MPS after decomposing them into a matrix product representation.
The one-dimensional Bose-Hubbard model is well suited to such simulation, as the Hamiltonian (4) consists only of oneand two-site nearest-neighbor operators. We have here used imaginary time evolution to obtain the ground-state MPS, and evolved this state according to each of H 0 and H 1 . The overlap of these evolved states then gives us an expression for the dephasing function (up to a known phase). We also found the expectation values for the observables n i n j (n i = b † i b i ),  for the purpose of providing exact values for comparison to the results obtained from the protocol. In our simulations we limited the maximum occupation per site to four bosons, with χ = 50 and time step δtJ = 10 −3 .