Precision bounds for gradient magnetometry with atomic ensembles

We study gradient magnetometry with an ensemble of atoms with arbitrary spin. We calculate precision bounds for estimating the gradient of the magnetic field based on the quantum Fisher information. For quantum states that are invariant under homogeneous magnetic fields, we need to measure a single observable to estimate the gradient. On the other hand, for states that are sensitive to homogeneous fields, a simultaneous measurement is needed, as the homogeneous field must also be estimated. We prove that for the cases studied in this paper, such a measurement is feasible. We present a method to calculate precision bounds for gradient estimation with a chain of atoms or with two spatially separated atomic ensembles. We also consider a single atomic ensemble with an arbitrary density profile, where the atoms cannot be addressed individually, and which is a very relevant case for experiments. Our model can take into account even correlations between particle positions. While in most of the discussion we consider an ensemble of localized particles that are classical with respect to their spatial degree of freedom, we also discuss the case of gradient metrology with a single Bose-Einstein condensate.


I. INTRODUCTION
Metrology plays an important role in many areas of physics and engineering [1]. With the development of experimental techniques, it is now possible to realize metrological tasks in physical systems that cannot be described well by classical physics and instead quantum mechanics must be used for their modeling. Quantum metrology [2][3][4][5] is the novel field that is concerned with metrology using such quantum mechanical systems.
One of the basic tasks of quantum metrology is magnetometry with an ensemble of spin-j particles. Magnetometry with a completely polarized state works as follows. The total spin of the ensemble is rotated by a homogeneous magnetic field perpendicular to it. We would like to estimate the rotation angle or phase θ based on some measurement; this phase parameter can then be used to obtain the field strength. To determine the rotation angle, one needs, for instance, to measure a spin component perpendicular to the mean spin.
Up to now, it looks as if the total spin behaves like a clock arm and its position tells us the value of θ exactly. At this point one has to remember that we have an ensemble of N particles governed by quantum mechanics, and the uncertainty of the spin component perpendicular to the mean spin can never be zero. Hence, simple calculation shows that the scaling of the precision of the phase estimation is (∆θ) −2 ∼ N , which is called shot-noise scaling [2][3][4][5]. However, spin squeezing [6][7][8][9][10] can decrease the uncertainty of one of the components perpendicular to the mean spin and this can be used to increase the precision * iagoba.apellaniz@gmail.com † toth@alumni.nd.edu; http://www.gtoth.eu of the measurements [10]. While it is possible to surpass the shot-noise limit, for the case of a linear Hamiltonian [2][3][4][5], no quantum state can have a better scaling in the precision than (∆θ) −2 ∼ N 2 , called Heisenberg scaling.
In this paper, we compute precision bounds for the estimation of the magnetic field gradient (see Fig. 1). In general, in order to achieve these bounds, an estimate of the constant (homogeneous) part of the field is required. Hence, we have to use the formalism of multiparameter estimation. Magnetometry of this type can be realized with differential interferometry with two particle ensembles, which has raised a lot of attention in quantum metrology [15,[39][40][41][42][43][44]. Another possibility is considering spin chains, which can be relevant in trapped cold ions or optical lattices of cold atoms, where we have individual access to the particles [45][46][47].
Finally, gradient magnetometry can be carried out using a single atomic cloud, which is very relevant from the point of view of cold gas experiments. One can consider both atomic clouds of localized particles, as well as Bose-Einstein condensates. While most works in magnetometry with a single ensemble focus only on the determination of the strength and direction of the magnetic field, certain measurement schemes for the gradient have already been proposed and tested experimentally. Some schemes use an imaging of the ensemble with a high spatial resolution. They do not count as FIG. 1. Schematic representation of an atomic ensemble (blue cloud) placed in a magnetic field (green lines) in a Stern-Gerlach apparatus. From the final state the gradient of the field can be estimated. Note that the field intensity changes along the cloud of atoms, while the direction is always the same from north (N) to south (S) within the cloud. For an easier presentation, a setup is shown where the direction of the magnetic field is parallel to the cloud, however, this does not have always to be the case.
single-ensemble methods in the sense we use this expression in our paper since in this case not only collective observables are measured [18][19][20]. There is a method based on collective measurements of the spin length of a fully polarized ensemble given in Ref. [48]. There is also a scheme based on many-body singlet states described in Ref. [45].
We can also connect our results to entanglement theory [54][55][56]. We find that the shot-noise scaling cannot be surpassed with separable states, while the Heisenberg scaling can be reached with entangled states. However, the shot-noise scaling can be surpassed only if the particle positions are correlated, which is the case, for instance, if the particles attract each other.
Next, we present the main characteristics of our setup. For simplicity, as well as following recent experiments (e.g., Ref. [18]), we consider an ensemble of spin-j particles placed in a one-dimensional arrangement. The atoms are then situated along the x axis with y = z = 0. We assume that we have particles that behave classically with respect to their spatial state. That is, they cannot be in a superposition of being at two different places. On the other hand, they have internal degrees of freedom, their spin, which is quantum. This is a very good description to many of the cold gas experiments.
Based on these considerations, we assume that the state is factorizable into a spatial part and a spin part as where the internal state is decomposed in its eigenbasis as For the spatial part defined in the continuous Hilbert space, we assume that it can be modeled by an incoherent mixture of pointlike particles as where x = (x 1 , x 2 , . . . , x N ) is a vector which collects all the particle positions, P (x) is the spatial probability distribution function of the atoms, and dx denotes dx 1 dx 2 · · · dx N . Note that the spatial part (3) is diagonal in the position eigenbasis, which simplifies considerably our calculations (see Appendix B for more details). During the evolution of the state, correlations might arise between the internal and the spatial parts and the product form (1) might not be valid to describe the evolution of the system. First, we consider spin chains and two particle ensembles at different places. The gradient measurement with two ensembles is essentially based on the idea that the gradient is just the difference between two measurements at different locations. With these systems, it is possible to reach the Heisenberg scaling.
We also examine in detail the case of a single atomic ensemble. Since in such systems the atoms cannot be individually addressed, we assume that the quantum state is permutationally invariant (PI). We show that for states insensitive to the homogeneous magnetic field, one can reduce the problem to a one-parameter estimation scenario. Single-ensemble measurements have certain advantages because the spatial resolution can be higher and the experimental requirements are smaller since only a single ensemble must be prepared.
For completeness, we mention the case of Bose-Einstein condensates (BEC). The spatial state in this case is pure where |Ψ is the spatial state of a single particle. Hence, the spatial state is delocalized and it is not an incoherent mixture of various eigenstates of x. While we do not consider such systems in detail, our formalism could be used to model them.
We now outline the model we use to describe the interaction of the particles with the magnetic field. The field at the atoms is given as where we neglect the terms of order two or higher, and where O(ξ) is the usual Landau notation to describe the asymptotic behavior of a quantity, in this case for small ξ. We consider the magnetic field pointing in the z direction, hence, B 0 = B 0 (0, 0, 1) and B 1 = B 1 (0, 0, 1). For this configuration, due to the Maxwell equations, with no currents or changing electric fields, we have div B = 0, curl B = (0, 0, 0).
This implies l=x,y,z ∂B l /∂l = 0 and ∂B l /∂m−∂B m /∂l = 0 for l = m. Thus, the spatial derivatives of the field components are not independent of each other. In this paper, however, we consider an elongated trap. In the case of such a quasi-onedimensional atomic ensemble, only the derivative along the axis of the trap has an influence on the quantum dynamics of the atoms or a double-well experiment. We determine the precision bounds for the estimation of the magnetic field gradient B 1 . We calculate how the precision scales with the number of particles. We compare systems with an increasing particle number, but of the same size. As discussed later, if we follow a different route, we can obtain results that can incorrectly be interpreted as reaching the Heisenberg limit, or even a super-Heisenberg scaling.
The angular momentum of an individual atom is coupled to the magnetic field, yielding the following interaction term: where the operator B (n) z = B 0 + B 1x (n) acts on the spatial part of the Hilbert space andx (n) is the position operator of a single particle. Moreover, j (n) z is a single-particle spin operator, acting on the spin part of the Hilbert space. Finally, γ = gµ B where g is the gyromagnetic factor and µ B corresponds to the Born magneton, and we set = 1 for simplicity. We use the "ˆ" notation to distinguish the operatorx from the coordinate x. Later, we will omit it for simplicity. The Hamiltonian of the entire system is just the sum of all two-particle interactions of the type Eq. (7) and can be written as Equation (8) generates the time evolution of the atomic ensemble.
One could include also the kinetic energy in the Hamiltonian. Such an extra term causes that the gradient field pushes atoms in state |0 into one direction, while atoms in state |1 into the other direction. In our work, we do not take into account this effect. Moreover, we do not include in the model the initial thermal dynamics of the particles. Both of these effects are negligible in a usual setup, as shown in Appendix A.
We calculate lower bounds on the precision of estimating B 1 based on a measurement on the state after it passed through the unitary dynamics U = exp(−iHt), where t is the time spent by the system under the influence of the magnetic field. The unitary operator can be rewritten as where the b i = γB i t. The generator describing the effect of the homogeneous field is given as while the generator describing the effect of the gradient is We omit ⊗ and the superscripts (x) and (s) for simplicity, and use them only if it is necessary to avoid confusions. The operators H 0 and H 1 commute with each other. However, it is not necessarily true that the operators we have to measure to estimate b 0 or b 1 can be simultaneously measured. The reason for that is that both operators to be measured act on the same atomic ensemble. If the measurement operators do not commute with each other, then the precision bound obtained from the theory of QFI cannot necessarily be reached. For the particular cases studied in this paper, we prove that a simultaneous measurement to estimate both the homogeneous and the gradient parameter can be carried out (see Appendix E). On the other hand, in schemes in which the gradient is calculated based on measurements on two separate atomic ensembles or different atoms in a chain, the measuring operators can always commute with each other [14,15,46].
The paper is organized as follows. In Sec. II, general precision bounds for the estimation of the gradient of the magnetic field are presented. In Sec. III, we compute precision bounds for relevant spatial configurations appearing in cold atom physics such as spin chains and two ensembles spatially separated from each other. In Sec. IV, we consider a single atomic ensemble in a PI state and we calculate the precision bounds for various quantum states, such as the singlet spin state or the totally polarized state. In Sec. V, we consider Bose-Einstein condensates.

II. PRECISION BOUNDS FOR ESTIMATING THE GRADIENT
In this section, we show how the QFI helps us to obtain the bound on the precision of the gradient estimation. First, we discuss gradient magnetometry using quantum states that are insensitive to homogeneous fields. In this case, we need to estimate only the gradient and do not have to know the homogeneous field. Hence, this case corresponds to a single-parameter estimation problem.
Then, we discuss the case of quantum states sensitive to homogeneous fields. Even in this case, we are interested only in the gradient, and we do not aim at estimating the homogeneous field. In spite of this, gradient estimation with such states is a two-parameter estimation task. We introduce the basics of multi-parameter quantum metrology, and we adapt that formalism to our problem. We also show that the precision bound obtained does not change under spatial translation, which will be used later to simplify our calculations. In Appendix E, we show that even the precision bounds for states sensitive to the homogeneous field, appearing in this paper, are saturable.
Next, we summarize important properties of the QFI used throughout this paper (for reviews, see Refs. [2,4,53,[57][58][59][60][61][62]). Let us consider a quantum state with the eigendecomposition For two arbitrary operators A and B, and a state [Eq. (12)], the QFI is defined as [4,[49][50][51]53] where A k,k = k|A|k and B k,k = k|B|k . If the two operators are the same then, from Eq. (13), the usual form of the QFI is obtained: We list some useful properties of the QFI: (i) Based on Eq. (13), F Q [ , A, B] is linear in the second and third arguments This will make it possible to calculate the QFI for collective quantities based on the QFI for single-particle observables.
(ii) The QFI remains invariant if we exchange the second and the third arguments Equation (16) will help to simplify our calculations.
(iii) The following alternative form, is also useful since the correlation appears explicitly.
(v) The QFI is convex on the space of the density matrices, i.e., Hence, when maximizing the QFI, we need to carry out an optimization over pure states only.
In the following, we show the general form of the expressions giving the precision bounds for states insensitive to the homogeneous field, as well as for states sensitive to it. We also show that both bounds are invariant under the spatial translation of the system which makes the computing for particular cases much easier.

A. Precision bound for states insensitive to homogeneous fields:
Single-parameter dependence We will now consider quantum states that are insensitive to the homogeneous field. For such states, holds. Hence, the unitary time evolution given in Eq. (9) is simplified to and the evolved state is a function of a single unknown parameter b 1 .
When estimating a single parameter, the Cramér-Rao bound gives the best achievable precision as [4] It is always possible to find a measurement that saturates the precision bound, (22), which is indicated using the notation "| max =".
Observation 1. For states insensitive to the homogeneous fields, the maximal precision of the estimation of the gradient parameter b 1 is given as where the integral represents the correlation between the particle positions x n and x m . Moreover, Eq. (23) is translationally invariant, i.e., it remains the same after an arbitrary displacement d of the form of where d is the distance displaced and P x is the sum of all single-body momentum operators p (n) x in the x direction. Proof. We have to evaluate the right-hand side of Eq. (22). The state is a tensor product of the spatial and internal parts, and the spatial part is an incoherent mixture of position eigenstates, as in Eqs. (1) and (3). Hence, the eigenstates are |x, λ , where |x and |λ are defined in the spatial and internal Hilbert spaces, respectively. Then, the matrix elements of H 1 , which is diagonal in the spatial subspace, are obtained as Calculating Eq. (14) for A = H 1 , Eq. (22) leads to Eq. (23) (see Appendix C for details).
In the last part of the proof, we show that the precision (23) remains the same for any displacement of the system. We use the Heisenberg picture in which the operators must be transformed instead of the states. After the displacement, the operator H 1 describing the effect of the gradient is obtained as Hence, the unitary evolution operator of the displaced system is obtained as Using the commutation relation (20), we can see that Eq. (27) is equal to the time evolution given in Eq. (21).

B. Precision bound for states sensitive to homogeneous fields: Two-parameter dependence
We now show how to obtain the precision bounds for states sensitive to the homogeneous field. The homogeneous field rotates all the spins in the same way, while the field gradient rotates the spins differently depending on the position of the particles. Hence, in order to estimate b 1 , we have to consider the effect of a second unknown parameter b 0 . Note, however, that we are not interested to estimate b 0 precisely, we just need it to estimate b 1 .
In this case, the metrological performance of the quantum state is given by the 2 × 2 Cramér-Rao matrix inequality [4] where the covariance matrix is defined as The matrix elements of the quantum Fisher information matrix F Q are Unlike in the case of single-parameter estimation, Eq. (28) can be saturated only if the measurements for estimating the two parameters are compatible with each other [4,51,63]. Hence, we use " " instead of "| max " for the bounds for quantum states sensitive to the homogeneous fields.
Using the well-known formula for the inverse of 2 × 2 matrices, Eq. (28) yields for the precision of b 1 .

Observation 2.
For states sensitive to the homogeneous field, the expression to compute the precision bound for the gradient parameter takes the following form: Moreover, the bound, (31), similarly to Eq. (23), is invariant under spatial translations of the system.
Proof. To obtain the bound, (31), we need to consider the matrix elements of QFI one by one. First of all, we compute F 11 which has the same form as Eq. (23) Next, we have that H 0 , similarly to Eq. (25), is diagonal in the spatial |x basis, and its matrix elements in the |x, λ basis of the state are written as With this we obtain Note that Eq. (34) is not a function of the whole state but only of the internal (s) state. Finally, we compute F 01 and F 10 .
Let us now determine the bound on the precision for estimating the gradient on the translated system. We have to compute first the QFI matrix elements. We use the linearity of the last two arguments of F Q [ , A, B] given in Eq. (15), the fact that H 0 remains unchanged in the Heisenberg picture. We also use the formula (26) for the shifted H 1 operator. The diagonal element of the QFI matrix corresponding to the measurement of the homogeneous field is hence, it does not change due to the translation. For the diagonal element corresponding to the gradient measurement we obtain Finally, for the off-diagonal element, we get After determining all the elements of the QFI matrix, the bound for a displaced system can be obtained as The bound in Eq. (30) can be obtained from the right-hand side of Eq. (39) with straightforward algebra.

III. SPIN CHAIN AND TWO SEPARATED ENSEMBLES FOR MAGNETOMETRY
After presenting our tools in Sec. II, we start with simple examples to show how our method works. We calculate precision bounds for gradient metrology for spin chain and for two-particle ensembles separated by a distance.
Before considering the setups mentioned above, we introduce various quantities describing the distribution of the particles based on the probability distribution function appearing in Eq. (3). The mean particle position is The standard deviation of the particle positions, describing the size of the system, is computed as Finally, the covariance averaged over all particle pairs is The covariance is a large positive value if the particles tend to be close to each other, while it is negative if they tend to avoid each other. After presenting the fundamental quantities above, let us study concrete metrological setups. The first spatial state we consider is given by N particles placed equidistantly from each other in a one-dimensional spin chain, as shown in Fig. 2. Such a system has been studied also in the context of a single parameter estimation in the presence of collective phase noise [44]. The probability density function describing such a system is where a is the distance between the particles in the chain. For this system, the average position of the nth particle is whereas the two-point average (42) is The standard deviation defined in Eq. (41) is obtained as Next, we will obtain precision bounds for particles placed in a spin chain.
Observation 3. Let us consider a chain of N spin-j particles placed along the x direction separated by a constant distance, and a magnetic field pointing in the z direction. Then, for the spin-state totally polarized in the y direction, the precision bound is given by Here, σ ch denotes the standard deviation of the average position of the particles for the chain (ch). Proof. We use the precisions bound for states sensitive to the homogeneous field given in Eq. (31). We obtain Note that the bound can be saturated (see Appendix E). Here, for the last equality we used the definitions of the average quantites given in Eqs. (44) and (45) and we also used Eq. (18) giving the QFI for pure states. We can see that the standard deviation given in Eq. (46) coincides with a factor we have in Eq. (49), with which we conclude the proof. Note that the bound (49) seems to scale with the third power of the particle number N , and hence seems to overcome the ultimate Heisenberg limit. The reason is that the length of the chain increases as we introduce more particles into the system. We should compare the metrological usefulness of systems with different particle numbers, but of the same size. In our case, we use throughout the paper the standard deviation of the averaged particle positions as a measure of the spatial size of the system, and normalize the results with it. One can miss this important point since when only the homogeneous field is measured such a normalization is not needed. 1 After the spin chain, we consider estimating the gradient with two ensembles of spin-j atoms spatially separated from each other. Such systems have been realized in cold gases (e.g., Ref. [64]), and can be used for differential interferometry [15,39,44]. We will determine the internal state with the maximal QFI.
Let us assume that half of the particles are at one position and the rest at another one, both places at a distance a from the origin. The probability density function of the spatial part is Such a distribution of particles could be realized in a double-well trap, where the width of the wells is negligible compared to the distance between the wells. To distinguish the two wells we use the labels "L" and "R" for the left-hand side and right-hand side wells, respectively. Based on these, we obtain the single-point averages as The two-point correlation functions are x n x m P (x) dx = +a 2 if (n, m) ∈ (L,L) or (R,R), −a 2 if (n, m) ∈ (L,R) or (L,R).
For the average particle position we obtain µ = 0, while the standard deviation for the spatial state in the double well (dw) is Next, we calculate the achievable precision of the gradient estimation.
Observation 4. For the case of two ensemble of N spin-j particles, the state that maximizes the QFI is The best achievable precision is given as Equation (55) agrees with the results obtained in Ref. [39].
Proof. The state given in Eq. (54) is insensitive to the homogeneous field, hence we have to use the formula (23) to bound the precision. We obtain For the state (54), the equation above, (56), yields where we have used the definition of the QFI for pure states given in Eq. (18). A factor in Eq. (57) can be identified with the standard deviation, (46), from which the proof follows.
It is interesting to simplify the QFI for product states states |ψ (L) ⊗ |ψ (R) , where |ψ (L) and |ψ (R) are pure states of N/2 particles each. This approach is also discussed in Ref. [39]. Such states can reach the Heisenberg limit, while they are easier to realize experimentally than states in which the particles in the wells are entangled with each other.
Before obtaining the precision for the case above, we present a method to simplify our calculations. The system is at the origin of the coordinate system such that for mean particle position given in Eq. (40), holds. Thus, the second term in the expression for the bound for states sensitive to the homogeneous field (31) is zero since all F Q [ (s) , j (n) z , J z ] are equal considering product states of two equal permutationally invariant states, |ψ (L) ⊗ |ψ (R) . Hence, the bounds for states insensitive and sensitive to the homogeneous field, Eqs. (23) and (31), respectively, are the same in this case.
We now compute F Q [ρ, H 1 ] for the case when the state is sensitive to the homogeneous field, hence, we use the bound on the precision given in Eq. (31). Using the the probability density distribution function given in Eq. (50), and following steps leading to Eq. (57), we obtain where we used that F Q [|ψ (L) , J z ]. Note that our results concerning using product states for magnetometry can be interpreted as follows. In this case, essentially the homogeneous field is estimated in each of the two wells, and then the gradient is computed from the measurement results. The bounds for these type of states are also saturable (see Appendix E).
We will now present precision bounds for various well known quantum states in the two wells. We consider the TABLE I. Precision for differential magnetometry for various product states of the type |ψ (L) ⊗ |ψ (R) in the two ensembles. Note that there are NL = NR = N/2 particles in each ensemble. In the second column we show the QFI for the estimation of the homogeneous field appearing in the literature, for states with NL particles. The third column shows the result for the bounds obtained with Eq. (59).

IV. MAGNETOMETRY WITH A SINGLE ATOMIC ENSEMBLE
In this section, we discuss magnetometry with a single atomic ensemble. We consider a one-dimensional ensemble of spin-j atoms placed in a trap which is elongated in the x direction. The setup is depicted in Fig. 3. In the second part of the section, we calculate precision bounds for the gradient estimation with some important multiparticle quantum states, for instance, Dicke states, singlet states, and GHZ states.

A. Precision bound for an atomic ensemble
In an atomic ensemble of many atoms, typically the atoms cannot be individually addressed. We will take this into account by considering states for which both the internal state (s) and holds, where the summation is over all possible permutations P k of the variables x n . Hence, we do not need to sum over all possible n's in Eqs. (40) and (41), and neither to sum over all n's and m's in Eq. (42). All the terms in each sum are equal to each other due to the permutationally invariance of the probability distribution function (62). An interesting property of the covariance (42) is that it can only take values bounded by the variance in the following way: where both the lower and the upper bounds are proportional to the variance σ 2 . See Fig. 4 for examples on how different correlations are obtained in an atomic 1D lattice. Next, we present precision bounds for PI states. Observation 5. The maximal precision achievable by a single atomic ensemble insensitive to homogeneous fields is The precision given in Eq. (64) can be reached by an optimal measurement. Nevertheless, it is worth to note that the precision cannot surpass the shot-noise scaling because F Q [ (s) , j z ] cannot be larger than j 2 . Moreover, η cannot be smaller than −σ 2 /(N − 1) due to Eq. (63), which makes its contribution negligible for large N .
Proof. From the definition of the QFI for states insensitive to the homogeneous field [Eq. (23)] we obtain the bound for a single ensemble as Then, we have to use the fact that for states insensitive to the homogeneous fields F Q [ , J z ] = 0 holds, which implies Based on this, for such states the sum of QFI terms involving two operators can be expressed with the sum of QFI terms involving a single operator as Substituting Eq. (67) into Eq. (65), Observation 5 follows.
Observation 6. For states sensitive to homogeneous fields, the precision of estimating the gradient is bounded from above as which may surpass the shot-noise scaling whenever η is a positive constant.
Proof. We start from Eq. (31) and take into account that in this case the bound is saturable (see Appendix E). As explained in Sec. II B, if we move the system, the precision bounds do not change. We then move our system to the origin of the coordinate system yielding µ = 0, and making the second term appearing in Eq. (31) zero. Thus, we only compute the first term in Eq. (31) and obtain z ] to the last term and subtract it from the first term to make the expression more similar to Eq. (64).
Note that the second term on the right-hand side of Eq. (68) is new in the sense that it did not appear in the bound for states insensitive to homogeneous fields given in Eq. (64). Even if the first term cannot overcome the shot-noise limit, in the second term the covariance is multiplied by the QFI for estimating the homogeneous field and therefore this concrete term, for extremely correlated particle positions, allows to achieve the Heisenberg scaling.

B. Precision limit for various spin states
In this section, we present the precision limits for various classes of important quantum states such as the totally polarized state, the state having the best precision among separable states, the singlet state, the Dicke state (61), or the GHZ state (60). We calculate the precision bounds presented before, (64) and (68), for these systems.

Singlet states
A pure singlet state is a simultaneous eigenstate of the collective J z and J 2 operators, with an eigenvalue zero for both operators. We will now consider PI singlet states. Surprisingly, the precision bound is the same for any such state. PI singlet states are very relevant for experiments, since they have been experimentally created in cold gases [79,80] while they also appear in condensed matter physics [81].
Let us now see the most important properties of singlet states of an N -particle system. There are several singlets pairwise orthogonal to each other. The number of such singlets, D 0 , depends on the particle spin j and the number of particles N . It is the most natural to write the singlet state in the angular momentum basis. The basis states are |J, M z , D , which are the eigenstates of J 2 x + J 2 y + J 2 z with an eigenvalue J, and of J z with an eigenvalue M z . The label D is used to distinguish different eigenstates corresponding to the same eigenvalue of J and J z . Then, a singlet state can be written as where D p D =1.
Let us see some relevant single-particle expectation values for the singlet. Due to the rotational invariance of the singlet (s) singlet , we obtain that holds. We also know that for the sum of the second moments of the single particle angular momentum components holds. Hence, the expectation value of the second moment of the single-particle angular momentum component is obtained as After discussing the main properties of the singlet states, we can now obtain a precision bound for gradient metrology with such states.
Observation 7. For PI spin states living in the singlet subspace, i.e., states composed of vectors that have zero eigenvalues for J z and J 2 and all their possible statistical mixtures, the precision of the magnetic gradient parameter is bounded from above as Proof. First compute the QFI for the one-particle operator For that we need that when j (n) z acts on a singlet state, produces a state outside of the singlet subspace. Hence, for any pair of pure singlet states. Then, we use the formula (17) to compute the QFI. The second term of Eq. (17) is obtained as due to Eq. (75). It follows that the single-particle QFI for any singlet equals four times the second moment of the angular momentum component Note that Eq. (77) is true even though z ] for any n. Then, we have all the ingredients to evaluate the maximal precision given in Eq. (64), and with that we prove the Observation.
As mentioned earlier, singlet states are insensitive to homogeneous magnetic fields, hence determining the gradient leads to a single-parameter estimation problem. This implies that there is an optimal operator that saturates the precision bound given by Eq. (74). However, it is usually very hard to find this optimal measurement, although a formal procedure for this exists [4]. In Ref. [45], a particular setup for determining the magnetic gradient with PI singlet states was suggested by the measurement of the J 2 x collective operator. For this scenario the precision is given by In Appendix D, we show that this measurement provides an optimal precision for gradient metrology for all PI singlets.

Totally polarized state
The totally polarized state can easily be prepared experimentally.
It has already been used for gradient magnetometry with a single atomic ensemble [18,19]. For the gradient measurement as for the measurement of the homogeneous field, the polarization must be perpendicular to the field we want to measure.
We chose as before the totally polarized state along y axis, given in Eq. (47). The relevant variances for the state, (47), are We can see clearly that the precision scales as O(N ) for large N .
Let us now see which measurement could be used to estimate the field gradient with a totally polarized state. The homogeneous field rotates all spins by the same angle, while the gradient rotates the spin at different positions by different angles. Due to that, the homogeneous field rotates the collective spin, but does not change its absolute value. On the other hand, the field gradient decreases the absolute value of the spin since it has been prepared to be maximal, which has been used in Ref. [48] for gradient magnetometry (see Fig. 2). Hence, we can measure the spin length to estimate the field gradient.

Best separable state
We now turn our attention to the precision bound for all separable spin states. It is useful to obtain this value so we have a direct comparison on what the best classically achievable precision is. It turns out that for j > 1 2 , it is possible to achieve a precision higher than with the fully polarized state (47).
Let us consider general separable states, which are not necessarily PI. We do not know if the optimal separable state is sensitive or insensitive to the homogeneous field. The corresponding precision bounds for the gradient estimation are given in Eqs. (23) and (31), respectively. Since the probability density function (62) is PI, we have for all n. As explained in Sec. II B, by moving the ensemble, the precision bounds do not change. If we move the system to the origin of the coordinate system achieving µ = 0, we can make our calculations simpler since the second term appearing in Eq. (31) is zero. Thus, we only compute the first term in Eq. (31). Hence, the two bounds (23) and (31) are the same in this case and we arrive at where we already assume that the bound can be saturated (see Appendix E).
We now look for the separable state that maximizes the right-hand side of Eq. (82), which has to be a pure product state due to the convexity of the quantum Fisher information. Hence, we look for the pure product state maximizing ]. Based on Eq. (18), for product states we find that For all n, a state that maximizes Eq. (83) is for which the single-particle variances are maximal, i.e, (∆j (n) z ) 2 = j 2 . While we carried out an optimization over general, non-necessarily PI separable states, the optimal state is PI. Plugging the state (84) into the bound given in Eq. (82) leads to the precision bound for separable states as where we have used the definition of the variance of the particle positions (41) for a permutational invariant state. Note that the bound for the best separable state given in Eq. (85) is above the bound obtained for the singlet state (74), whereas the bound for the totally polarized state in Eq. (80) is below. Nevertheless, when the singlet state is used, the homogeneous magnetic field has no effect on the state. In contrast, the state (84) is sensitive to the homogeneous field.

Unpolarized Dicke states |DN and |DN x
Next, we compute precision bounds for entangled states. In this section, we consider unpolarized Dicke states, which play an important role in quantum optics and quantum information science. The Dicke state |D N l [Eq. (61)], with a maximal J 2 x + J 2 y + J 2 z and J l = 0 for any l ∈ x, y, z is particularly interesting due to its entanglement properties and its metrological usefulness [71,72]. This state has been created in photonic experiments [73][74][75] and in cold atoms [76,77], while a Dicke state with J z > 0 has been created with cold trapped ions [78]. The Dicke state |D N is an eigenstate of J z so it is insensitive to a homogeneous magnetic field pointing into the z direction. Thus, the precision bound can be saturated by some measurement. The Dicke state |D N x is sensitive to the homogeneous field. Moreover, it is very useful for estimating the homogeneous field as it has been shown in Ref. [76]. Here, we consider large particle numbers, to make the results simpler.
Let us now see the most important properties of Dicke states. For the expectation values of the single-particle angular momentum components j (n) l = 0 (86) hold for l = x, y, z for all n. The second moments of the collective angular momentum components are given as Let us now see two-body correlations. Since the Dicke state is PI, we have for all m = n and l = x, y, z. Hence, the collective second moments are connected to the single particle and two-particle operator expectation values as for l = x, y, z. Considering the symmetry under rotations around z axis, we also have (j (1) . Based on these and using Eq. (72) for j = 1/2, we arrive at [45] (j (n) for l = x, y, z.
After discussing the main properties of the Dicke states, we can now obtain a precision bound for gradient metrology with such states.
Observation 8. For large N , the precision bound for the Dicke state |D N is For the Dicke state |D N x , the precision is bounded as which allows in principle a Heisenberg limited behavior due to the second term on the right-hand side.
Proof. Let us prove first Eq. (91). Since |D N is a pure state, the QFIs appearing in Eq. (64) are simply four times the corresponding variances of j l , respectively, we obtain From Eq. (93) and the bound for states insensitive to the homogeneous field (64), the precision bound for the Dicke state |D N follows.
We prove now the bound for the |D N x states given in Eq. (92). The second moments (j (n) z ) 2 for |D N x can be obtained from the second moments computed above for |D N by relabeling the coordinate axes. Since |D N x is a pure state, the QFI again equals four times the corresponding variance. Hence, we obtain and using the bound for states sensitive to homogeneous fields given in Eq. (68) we have all we need to prove Observation 8.

GHZ state
The GHZ states are defined for qubits in Eq. (60). Such states are very sensitive to the homogeneous field. GHZ states are highly entangled and play an important role in quantum information theory [65]. They have been created experimentally in photonic systems [66][67][68] and trapped ions [69,70].
Let us see first the relevant expectation values for GHZ states. Direct calculation shows that Moreover, for the second moments hold.
Let us now calculate the precision bound. We recall that for pure states the QFI is given as Eq. (18). Using the bound for states sensitive to homogeneous fields given in Eq. (68), we obtain From (97) follows that we can reach the Heisenberg limit with such states, but only in cases where η is positive, i.e., when the particles are spatially correlated.

Summary of results
Finally, we summarize the precision bounds obtained for various quantum states in Table II. In Fig. 5, we show the mean values and variances of the collective angular momentum components for these states. Note that for these PI states the optimal estimators for the homogeneous field and the gradient field are compatible (see Appendix E). It means that the two parameters can be estimated at once even for the states sensitive to the homogeneous fields.

V. GRADIENT MAGNETOMETRY WITH A BOSE-EINSTEIN CONDENSATE
In this section we study the case when our external state is a Bose-Einstein condensate instead of an incoherent mixture of pointlike particles. We can write the spatial state of a BEC [Eq. (4)] as where we define the state |0 as the pure product state representing the BEC. Since all particles are in the same spatial state, several important quantities describing the ensemble can easily be computed. For such a quantum state, for the average particle position defined in Eq. (40) we obtain for all n. For the variance given in Eq. (41), holds for all n. Finally, there is no correlation between particle positions, i.e., x (n) x (m) = x (n) x (m) if n = m. Hence, the covariance, (42), is zero Finally, as explained in Sec. II B, the precision bounds do not change if we translate the system. We move the atomic ensemble to the origin of the coordinate system such that This will make our calculations much simpler. Based on Eq. (22), for states insensitive to the homogeneous field we obtain Based on Eq. (30), for states sensitive to the homogeneous field we obtain Here we used that is zero due to Eq. (102). The bounds needed in Eqs. (103) and (104) are equal to each other and can be obtained as follows. We will compute a bound on F 11 on pure states. Straightforward algebra leads to One can see that the optimal spin state for gradient estimation is the state totally polarized in the z direction which is separable. Hence, the precision is bounded for spin-j particles as This is quite surprising, since under the dynamics coupling to the z component of the spin and hence it rotates around the z axis. One would naively expect that the optimal state is the state totally polarized in the y direction (47) studied in Sec. IV B 2 for the case of cold atomic ensembles. Due to the convexity of the quantum Fisher information, the bounds are also valid for the case of a mixed spin state. Based on Eq. (108), we see that the Heisenberg scaling cannot be reached in this case. Interestingly, this is true for any spatial wave function. For instance, if a single BEC is in a double-well potential, it still cannot have a scaling better than the shot-noise scaling in gradient estimation. In contrast, in Sec. III we have seen that a Heisenberg scaling is possible in a double well, if two independent BECs are in the two wells.

VI. CONCLUSIONS
In this work, we investigated the precision limits of measuring the gradient of a magnetic field with atomic ensembles arranged in different geometries and initialized in different states. We were particularly interested as to how the best achievable precision scales with the number of particles. For spin chains and the two-ensemble case, the precision of the estimation of the gradient can reach the Heisenberg limit. For a single ensemble with localized particles, the shot-noise limit can be surpassed and even the Heisenberg limit can be achieved if there is a strong correlation between the particle positions. We also studied the case of a single Bose-Einstein condensate, and found that the shot-noise limit can not be surpassed in this case. However, even if the Heisenberg limit is not reached, single-ensemble methods can have a huge practical advantage compared to methods based on two or more atomic ensembles since using a single ensemble makes the experiment simpler and can also result in a better spatial resolution. Independently from our work, Ref. [82] studied gradient metrology for different configurations of N particles distributed on a line. split in two. Moreover, the more distance between the two final subensembles, the larger the gradient of the field. Hence, surprisingly, taking into account the movement of the particles induced by the gradient reduces the error in the estimation, so neglecting it, our bounds on the precision are still valid.
Nevertheless, the gradient induces a force which depends on the spin state of the atoms. The force is constant, thus the position will change quadratically in time. On the other hand, the spin state changes linearly. Hence, for small enough evolution times the displacement of the particles can be neglected.
Moreover, in a typical experiment for sensing the gradient of the magnetic field, sensitivities of the order of 1 pT/µm can be reached, for a gradient of the magnetic field of 100 nT/µm [48,80,84]. Hence, the classical acceleration due to the gradient of the magnetic field is a ≈ g F µ B B 1 /m, where m and g F are the mass and the gyromagnetic g-factor of a 87 Rb atom, respectively, m ≈ 87 u and g F ≈ 0.5. This results in an acceleration of the order of 3 × 10 −2 m/s 2 . After 0.5 ms of evolution [48], the atom travels a distance of the order of 10 nm, which is irrelevant compared with the size of these systems.
Next, let us consider the thermalization of the state which introduces random displacements of the particles potentially blurring the signal. A typical cigar-shape ensemble of 87 Rb atoms used for gradientometry is a couple of millimeters long and temperatures around 20 µK [80,84]. We use the formula that connects the mean-root-average velocity of the particles and the temperature,v = 3k B T /m. Note that not all the particles move towards the same direction but randomly in any direction. Hence, we compute the average of the modulus of the projection of the velocity parallel to the direction of the cloud as |v | =v/2. We conclude that the atoms are displaced by around 19 µm along the axis to the cloud, which again is irrelevant for clouds of the size of millimeters [48,84].
Moreover, the displacement due to the gradient and thermal dynamics can be clearly neglected in the cases of the spin chain, the two ensembles, and the BEC, which are discussed in Secs. III and V. Hence, the precision bounds computed in this paper can be used as a tool to characterize different states.
Concerning the sensitivity of our magnetometer, we can say the following. Assuming N = 8.5 × 10 6 atoms, trap length σ = 3 mm, and for the completely polarized state discussed in Sec. IV B 2, we obtain ∆B 1 ≈ 3 pT/mm, which is similar to the state of the art of other cold gas magnetometers [48]. The precision can be considerably improved if we use entangled states and we have correlation between the particle positions. There are other setups that work at much lower length scales, however, it is difficult to compare them to our system since they would not work at mm length scales [48].
Appendix B: Spatial state of thermally distributed pointlike particles We discuss the spatial state represented by Eq. (3). For that, let us introduce the position operator aŝ where x is a vector of the particle positions, and |x denotes a spatial state in which the pointlike particles are at given positions with the usual normalization as expected. Based on Eq. (B1), we see that Thus, |x is an eigenstate of the operatorx. In order to obtain a quantum state that represents N pointlike particles placed in the locations determined by the x vector, we have to normalize it as From Eq. (B4) and using that there is a probability distribution function, P (x), and defining P (x) as the probability to find particles at a given position x, we arrive at Eq. (3).
Appendix C: Calculation of the QFI matrix elements for pointlike particles In this appendix, we show how to compute the QFI F Q [ , H i , H j ] if the spatial part of the state is written as Eq. (3).
Let us write first the density matrix in its eigenbasis as where P (x)p λ / x|x are the eigenvalues. Based on Eq. (13), the QFI matrix elements are written as [P (x)p λ − P (y)p ν ] 2 P (x)p λ + P (y)p ν × (H i ) x,λ;y,ν (H j ) y,ν;x,λ dx dy. (C2) Note that x|x ≡ y|y and that the integral is over 2N variables, x and y.
We now use the fact that the generators H 0 and H 1 are diagonal in the spatial basis [see Eqs. (33) and (25)]. Hence, the matrix elements can be rewritten as for i = 0, 1, where H i is a shorthand for  In this appendix, we prove that the precision limits for gradient metrology can be saturated for singlet states if we measure J 2 x . Observation 9. Let the initial spin state of an atomic ensemble be an arbitrary PI singlet state (s) singlet . Consider the experimental setup when b 1 is obtained by measuring J 2 x . The precision of estimating b 1 , which is given by the error propagation formula, is optimal in the short-time limit, i.e., where J k x (t) = U † (t)J k x U (t), the time-evolution unitary operator is of the form U (t) = e −ib1H1 , and H 1 is defined in Eq. (11).
Proof. Since for any pure singlet We see that both the numerator and denominator of the right-hand side of Eq. (D1) go to zero as t → 0, thus the l'Hospital rule can be used applying the derivative ∂ b1 in both the denominator and the numerator, which yields . (D4) However, here the numerator and the denominator are again zero at t = 0, so we employ the l'Hospital rule once again and obtain lim t→0 (∆b 1 ) −2 = = lim where we simplified the expectation values that are 0 and we used the Heisenberg equation of motion twice for the second derivatives and simplified the result, using Eq. (D2) and the definition of the commutator, to rewrite the equation. Next, we will compute the numerator and the denominator in Eq. (D5). First of all using the angular momentum commutation relation [j One can see by inspection that the operators given in Eqs. (E5a) and (E5b) commute with each other. Hence, all our bounds for PI states in Sec. IV can be saturated, and all PI states discussed in other sections as well.
Finally, we have to discuss the states appearing in Table I. They are product states of two PI states of N/2 particles each. Thus, in terms of "L" and "R," we have the following expression z is the z projection of the total angular momentum of the "L" subsystem, as J (R) z is for "R." Clearly, the operator L( , H 1 ) commutes with L( , H 0 ), which is given in Eq. (E5a).
With this, we conclude this appendix which let us demonstrate that all bounds in this paper can be saturated.