## Abstract

We have recently completed a measurement of the Newtonian constant of gravitation *G* using atomic interferometry. Our result is *G*=6.67191(77)(62)×10^{−11} m^{3} kg^{−1} s^{−2} where the numbers in parenthesis are the type A and type B standard uncertainties, respectively. An evaluation of the measurement uncertainty is presented and the perspectives for improvement are discussed. Our result is approaching the precision of experiments based on macroscopic sensing masses showing that the next generation of atomic gradiometers could reach a total relative uncertainty in the 10 parts per million range.

## 1. Introduction

This monographic issue clearly shows that measuring the Newtonian constant of gravitation *G* with a total uncertainty below 100 ppm is a formidable task. In the past 15 years, at least five independent groups reported measurements with a total uncertainly below 30 ppm [1–5]. Most measurements however are mutually incompatible according to the standard statistical tests and are scattered over a range of about 460 ppm.

Clearly, sources of systematic error must be present and this is acknowledged in the recommended value for *G*. In the last three versions of the CODATA database, the standard relative uncertainties moved from 150 ppm in 2002, to 100 ppm in 2006 and finally to 120 ppm in 2010 [6].

Measurements based on different techniques can then be useful in an attempt to detect systematic errors, even if their uncertainties are quite higher than the best reported ones. In this spirit, more than 10 years ago an experiment aiming at a measurement of *G* based on an atomic gradiometer with an uncertainty in the 100 ppm range was started in Florence [7]. After proof-of-principle measurements at 1% in Florence [8], 0.5% in Stanford [9] and finally 0.2% again in Florence [10], we have reported an uncertainty of 150 ppm [11] which is for the first time comparable with that of the current CODATA value.

Since atomic interferometry is a newcomer in the field of precision *G* measurements, special care must be taken in describing the process that leads from raw data to a value for *G*. In this article, after a schematic description of the operating principle of an atomic gradiometer we would like to discuss the evaluation of type A and type B uncertainties with emphasis on data analysis. A detailed description of the sources of instability and uncertainty that should limit our experimental apparatus below the 100 ppm level has already been published [12]. Finally, we will present the short-term perspectives for future high precision *G* measurements based on cold atoms.

## 2. Atomic gravity gradiometer

A comprehensive review on atomic interferometry is presented in [13] while the most recent results are summarized in [14]. Here, we will only describe the basic operating principle of an atomic gravity gradiometer and report the details that are necessary for the following analysis of the uncertainties. Our experimental apparatus has been described in [15], and its latest characterization can be found in [12]. More details are available in [16,17].

### (a) Principle of operation

Following the seminal work presented in [18], in our experiment a cloud of cold ^{87}Rb atoms can be launched along the vertical *z*-direction and prepared in the non-magnetic |*F*=1,*m*_{F}=0〉=|*a*〉 hyperfine state of the fundamental electronic level.

State preparation uses a combination of stimulated Raman transitions and blow-away optical pulses. The Raman pulses are induced by two counter-propagating laser beams (Raman beams) directed along the vertical *z*-axis that transfer atoms between the states |*a*〉 and |*F*=2,*m*_{F}=0〉=|*b*〉. The state |*b*〉 is also in the fundamental electronic level. The blow-away pulses use radiation pressure to remove from the cloud atoms that remain in the initial state after a Raman pulse.

Note that between the Raman beams and the atoms there is a small energy exchange of the order , where *ω*_{ab}≃2*π*×6.8 GHz is the frequency separation between |*a*〉 and |*b*〉, but a large momentum exchange , where *k*≃16×10^{6} m^{−1}, corresponding to twice the momentum of an optical photon at 780 nm. The associated recoil velocity for an atom is *v*_{r}≃12 mm s^{−1}. We can drive, at every launch, the transition choosing if the atomic recoil will be in the +*z* or −*z* direction. We will distinguish between these two cases, when necessary, by using the ↑ and ↓ subscripts respectively.

After preparation there are about 10^{5} atoms in the cloud which has a spatial Gaussian distribution with standard deviation *σ*_{s}≃3 mm.

During the parabolic flight the cloud is illuminated again by the Raman beams to create coherent superpositions of |*a*〉 and |*b*〉. Starting from a pure state |*a*〉 we can apply a *π*/2-pulse that creates a superposition of |*a*〉 and |*b*〉 with equal amplitudes, or a *π*-pulse that creates a pure state |*b*〉.

An ‘interferometric sequence’ here is a series of three pulses, *π*/2−*π*−*π*/2, separated by a time interval *T*. In our experiment *T* is 0.16 s. If a constant acceleration *a* in the *z*-direction is present, the probability of finding an atom in the state |*a*〉 at the end of the sequence can be written as where *ϕ*=*kaT*^{2}. Since *P*_{a} can be measured accurately by state selective fluorescence detection, a very sensitive acceleration sensor is available. To get a taste of all the subtleties actually involved in a high accuracy measurement of the local gravity acceleration, see [19]. Here we just want to stress that *ϕ* and *a* are linked by a conveniently large factor *kT*^{2} that is very easy to measure with high accuracy.

When two atomic clouds are launched sequentially and move with the same velocity, separated by a constant distance *d* along *z* (in our experiment *d*≃0.3 m), it is possible to apply the interferometric sequence to both clouds simultaneously and measure each *P*_{a}. In the following, we will use the superscripts ‘u’ and ‘l’ to denote quantities related to the ‘upper’ and ‘lower’ cloud, respectively. Both *ϕ*^{u} and *ϕ*^{l} are subjected to noise added by parasitic accelerations if special care is not taken in insulating the experiment from seismic or acoustic vibrations. The noise is however mostly in common mode on *ϕ*^{u} and *ϕ*^{l} so we can measure the difference *φ*=*ϕ*^{u}−*ϕ*^{l} with high accuracy by plotting the experimental values of versus obtained in many measurements and by fitting an ellipse to the experimental points [20].

If a constant acceleration gradient *γ* is present we have that *φ*≃*kdγT*^{2} where we have taken only the first order in *γT*^{2}. When using the numbers for our experiment and for *γ* the Earth gravitational gradient at the surface i.e. *γ*≃−3×10^{−6} s^{−2}, we obtain |*φ*|≃0.37 rad, a value sufficient to produce in the experimental versus plot a clearly non-degenerate ellipse. The conversion factor accuracy between *γ* and *φ* is typically limited by *d* while the magnitude is eventually constrained by the vertical size of the apparatus that fixes *d* and *T*.

With an atomic gradiometer, it is possible to measure *G* using a version of the Cavendish experiment where two source masses are moved vertically ‘inside’ and ‘outside’ the positions that the two clouds will reach at the apogees, as shown schematically in figure 1, in order to apply known gravitational forces along *z*.

In the following, we will label *φ*_{C} and *φ*_{F} the angles measured by the gradiometer when the masses are, with reference to figure 1, in the (*a*) (Close) and (*b*) (Far) configurations, respectively. Clearly, the angle *φ*_{C}−*φ*_{F} is proportional to *G*. Since the contribution to *φ* due to the gravitational field is sensitive to the recoil direction while that due to some spurious effects, namely the first-order light shift and the effect of inhomogeneous magnetic fields, is not, we alternate between *k*_{↑} and *k*_{↓} Raman transitions. Adjusting also the velocities of the clouds we then take two consecutive measurements at the same spatial position but exchanging |*a*〉 and |*b*〉 along the two interferometer's paths. We obtain then a phase shift *Φ*, averaged over the two *k*-directions, defined as
2.1

The protocol that leads to a *Φ* measurement will be described in detail later.

### (b) From phase shifts to *G*

The simple relation between the phase shift for a single cloud *ϕ* and *a* holds only for the ideal case of a uniform acceleration.

In the general case, there are different theoretical frameworks at increasing levels of accuracy and complexity that can be applied to evaluate *ϕ* starting from the gravitational potential *V* .

In the simplest case, the external degrees of freedom of the wave function are treated as momentum eigenstates and *ϕ* is evaluated using the semiclassical approximation [21] as the sum of three terms labelled as propagation phase *ϕ*_{p}, laser phase *ϕ*_{l} and separation phase *ϕ*_{s} [22]. At the same time, the length *τ* of the *π*/2 pulses (and 2*τ* of the *π* pulses) is considered negligible when compared with *T*.

When this last approximation is removed the corrections are expressed as a power series in *τ*/*T*. Since in our experiment *τ*=12 μs the first-order correction is sufficient for 10 ppm accuracy.

A more sophisticated treatment, where Gaussian wave packets are considered both in coordinates and momenta and the dynamics during the interferometric pulses are evaluated, is summarized in [23] while a detailed account can be found in [24].

At an accuracy level in *ϕ* of about 10 ppm a semiclassical approach is adequate. Moreover, we adopt a perturbative treatment: to the linear potential *V* _{0}=*mgz*, where *m* is the atomic mass, we add as a perturbation *V* _{1} which includes the Earth gravity gradient and the source masses potential *V* _{s}. The effect of Earth rotation does not need to be taken into account since during the interferometric sequence it is compensated by counter-rotating the Raman beams as described in [25,26]. Details on our implementation can be found in [12,17]. Since compensation is not perfect it affects the error budget for *G*, as discussed later. Perturbation theory must be adopted since in both the theoretical frameworks the exact solution is derived under the hypothesis of a Hamiltonian at most quadratic both in the momenta and the coordinates while *V* _{s} contains higher order terms.

In our experiment, *V* _{s} is comparable to the potential due to the Earth gradient so *m*|*γ*|*δz*^{2}≃|*V* _{1}| over a distance *δz*≃*gT*^{2}. In the perturbative series, we have terms scaling as *V* _{1}⋅(*V* _{1}/*V* _{0})^{n} and since *V* _{1}/*V* _{0}≃*γT*^{2}≃10^{−7} we can stop at the *n*=0 term. We then can write
2.2
where I and II are the two classical paths of the interferometer evaluated considering only *V* _{0} as a potential.

Moreover, in a perturbative calculation at the first order *ϕ*_{l} and *ϕ*_{s} do not contribute to *Φ* since the endpoint of a classical path does not change, as explained in appendix F in [27].

To include the effect of finite *τ*, we note that the approximation *V* _{s}=const. holds well during the pulses since Δ*V* _{1}/*V* _{1}∼*τ*^{2}/*T*^{2}. We can then evaluate the phase shift during the interferometric pulses by integrating the optical Bloch equations during the pulses, noting that the Doppler effect is compensated by changing linearly with time the frequency difference of the Raman beams, as in [28]. The result is very simple: the *π* pulse does not contribute to *Φ* so we can simply not integrate over the time of the *π* pulse in equation (2.2). A correction should also be applied to the contribution to *ϕ*_{p} due to the *π*/2 pulses but this is negligible at our level of approximation.

Given the initial conditions of the classical motion for the atoms, equation (2.2) is applied four times to obtain a value for *Φ* according to equation (2.1). An average value *Φ*_{s} for *Φ* then is obtained by a Monte Carlo simulation that picks the initial conditions according to the measured initial positions and velocities distributions of the atoms in the clouds.

Note that *Φ*_{s} is linear in *V* _{1} so we can just let *V* _{1}=*V* _{s} in equation (2.2) since the Earth gradient contribution is cancelled anyway when evaluating *Φ*_{s}. The linearity also implies *Φ*_{s}∝*G*.

We usually simulate 10^{5} atoms per cloud, i.e. a number comparable with the actual one. Increasing the number of simulated atoms by a factor 10 does not change significantly the result. The uncertainties in the initial distributions and in other parameters entering in the calculation are the main source of type B uncertainties.

## 3. Type A uncertainties

Before discussing statistical uncertainties, we need to describe briefly how data are acquired.

In our experiment, we can take a measurement every 1.9 s. We shall call this a ‘point’. Starting with the masses in the close position (figure 1), 720 points are acquired alternating between *k*_{↑} and *k*_{↓} Raman transitions, then the masses are moved to the far position and 720 more points are recorded. The whole process takes about 1 h to complete. By plotting the points, we obtain four ellipses that can be fitted to obtain the angles *φ*_{C↑},*φ*_{C↓},*φ*_{F↑} and *φ*_{F↓} that are combined, using equation (2.1), in a single *Φ* measurement with a typical relative error in the 1000 ppm range. An example of these ellipses relative to the *φ*_{C↑} and *φ*_{F↑} angles is shown in figure 2*a*.

Our measurement of *G* is based on a total of 107*Φ* measurements, shown in figure 2*b*, taken between 5 and 12 July 2013. This corresponds to more than 1.5×10^{5} points. The error bar on each *Φ* measurement is obtained from the fit to the four ellipses. The *χ*^{2} of our *Φ* measurements is 145.4 with 106 degrees of freedom. We attribute this to an underestimate of the error bar by the fit to the ellipses so we rescale the error bars by a factor 1.17 in order to bring the reduced *χ*^{2} to 1. After rescaling, the average measured value for *Φ* is *Φ*_{m}=0.547870(63) *rad*, with a relative uncertainty of 116 ppm.

According to the analysis presented in [12], we control many experimental parameters at a level which is sufficient for 100 ppm accuracy on *Φ*_{m} but probably not good enough for a significant improvement, so there is no guarantee that the statistical error will keep decreasing as expected, when taking more data, unless better control on some experimental parameters is achieved. For this reason, we did not attempt to further reduce the statistical error.

## 4. Type B uncertainties

Type B uncertainties stem mainly from the input parameters to the Monte Carlo simulation linking *Φ*_{s} to *G*. There are also contributions associated with the source masses not included in the simulation that have been evaluated separately. We will start discussing the effects on *Φ*_{s} due to the source masses then consider the trajectories and sizes of the clouds, and finally we will analyse the detection sensitivity and other minor effects before summarizing the results.

### (a) Source masses

The source masses and the positioning system have already been described in detail in [29]. Here, we discuss only what is relevant for evaluating *Φ*_{s}.

The masses can be divided into three groups: the 24 cylinders that, with a total mass of about 516 kg, generate about 0.97*Φ*_{s}; the supports for the cylinders, with a mass of about 70 kg, responsible for 0.028*Φ*_{s}; other minute parts, dominated by the mass of the four stainless steel sledges inside the vertical translators moving the supports, that amount to a total of about 7 kg and contribute to 0.002*Φ*_{s}. This hierarchy establishes the accuracy required to characterize each group of masses.

The cylinders have been individually weighed, taking into account air buoyancy, with an accuracy of 10 mg with a precision mass comparator (Mettler XP64003L, courtesy of INRIM), and measured with a contact three-dimensional scanner with an accuracy of 4 μm (DEA Scirocco 100907 with PH10 Renishaw probe). The diameters *d* and heights *h* are *d*=(99.895±0.003) mm and *h*=(150.107±0.007) mm, respectively, where the quoted uncertainties are the standard deviations of the sample.

Each group of 12 cylinders is placed around the axis of the Raman beams (figure 1) in a hexagonal pattern with the nominal horizontal positions of the centres lying on two circumferences with radii *d* and , respectively.

The actual positions of the cylinders relative to the centre of the Raman beams have been measured before and after acquiring the data used to determine *G* with a laser tracker (Leica AT 901 LR) with an accuracy of 10 μm on the horizontal plane [17]. In the vertical direction, the absolute calibration was also provided by the lasers tracker but motion between the close and far configuration was monitored with a pair of optical rulers (Heidenhain LS 603) with a reproducibility of 2 μm RMS.

A destructive test with a spare cylinder from the same batch has been performed in order to evaluate the density homogeneity. As expected for this material (Inermet180, a sintered alloy of 95% W, 3.5% Ni and 1.5% Cu), density has a minimum at the centre due to the production technique even if the parts where exposed to a hot isostatic pressing treatment that makes the average density of the batch quite uniform [29]. Indeed from the mass and the geometric measurements an average density *ρ*=(18263±2) g cm^{−3} is obtained. The uncertainty is the standard deviation of the sample.

The centre of gravity of another spare cylinder has been determined in the radial direction by measuring the oscillation period on an air bearing and in the axial direction with a precision scaler and a lever. Since the distance from the geometrical centre was comparable with the laser tracker resolution and in our experiment we have 24 cylinders randomly oriented, we assume the centres of gravity to be located, on the average, at the geometrical centre of the cylinders. The uncertainty on *Φ*_{s} due to the cylinders is 21 μrad.

For the supports, the nominal machining resolution of 0.1 mm was assumed and, where possible, verified with a Vernier caliper, while only one of the two supports has been taken apart after completing the *G* measurement and weighed. Its total mass is *M*_{s}=(34476±7) g. The other is assumed to have the same mass. The positions of the supports have also been measured with the laser tracker. The uncertainty on *Φ*_{s} due to the supports is of the order of 3 μrad.

The remaining moving parts have been weighed with a scale with an accuracy of 2 g and the distance of the centre of gravity from the Raman beam axis determined with a Vernier caliper with a resolution of 0.5 mm. The uncertainty on *Φ*_{s} is 3 μrad.

According to simulations, the density inhomogeneity of the cylinders is the dominant source of uncertainty for *Φ*_{s} and will be evaluated later in this section after discussing how *V* _{s} has been calculated. The uncertainty on the positions of the source masses amounts to 21 μrad dominated by the uncertainty on the horizontal positions of the cylinders. Another small correction arises from the non-zero value of the air density.

#### (i) Evaluation of the gravitational potential

The gravitational potential of the source masses *V* _{s} must be evaluated not only with sufficient accuracy but also with reasonable computational efficiency since obtaining *Φ*_{s} requires computing *V* _{s} about 10^{8} times due to the need to numerically evaluate equation (2.2) for each simulated particle.

An analytical expression for the gravitational potential *V* _{c} of a homogeneous cylinder of mass *M*, height 2*H* and radius *R* exists but it involves evaluating an integral of elliptic functions [30].

For computational purposes, it is then much more convenient expanding the potential in multipoles. Owing to the symmetry of the problem only the *m*=0,*l* even spherical harmonics *Y* _{0l} appear in the expansion and since where *P*_{l} are the Legendre polynomials, we can write
4.1
where *r* and *θ* are spherical coordinates in a reference frame with the origin at the centre of the cylinder and *θ* is measured from the cylinder axis, *α* is *H*/*R* and
4.2

Since *P*_{2n}(*x*) contains only even powers in *x* up to the degree 2*n* the integrand in equation (4.2) is actually a polynomial in *u*=*z*/*R* and *v*=*r*/*R* so *Q*_{2n}(*α*) can be easily evaluated for every *n* to a polynomial of degree 2*n* in *α*. To speed up convergence, it is clearly better to have *α*<1. In our cylinders, *α*≃1.5 so it is convenient, for computation, to split each cylinder in two and compute *V* _{c} for 48 cylinders with *α*=0.75. We truncate the series for *V* _{c} at *n*=20 obtaining a relative error always better than 10^{−11} along the interferometer paths.

The supports were modelled with a series of seven rings with density *ρ*, internal radius *R*_{i}, external radius *R*_{e} and height *H*. The potential *V* _{r} of a ring in cylindrical coordinates can be written as
4.3
where all the lengths in the integrand have been rescaled by dividing them by *R*_{i} while *α*=*H*/*R*_{i} and *β*=*R*_{e}/*R*_{i}. Again, for computational efficiency, instead of applying the results given in [30] we use an approximate expression. We are interested to *V* _{r} in the region close to the ring axis so we expand the integrand in a power series of *r*. After integrating in *θ* only even powers of *r* remain. We truncate the expansion at the *r*^{2} term, neglecting contributions to *V* _{r} of order *r*^{4}. In our experiment, the atomic clouds have a size of the order of 5 mm while *R*_{e} is about 5 cm for the smallest rings so *r*≃0.1 and the neglected terms introduce a change in *Φ*_{s} of less than 1 μrad. The integrations in *v* and *u* are straightforward but the resulting formula for *V* _{r} is quite long so we will not explicitly report it here.

Finally, the other masses are treated like point masses. They are at a distance *r* of about 30 cm from the axis of the Raman beams and can be enclosed in a sphere of radius *R*≃5 cm. The dominant contribution comes from the four sledges mentioned above which are reasonably symmetric and therefore have a negligible dipole term. The first neglected term in the multipole expansion is then the quadrupole so the error on *Φ*_{s} is again of the order of 1 μrad.

#### (ii) Density inhomogeneity

According to the measurements reported in [29], the density of our cylinders is not uniform but rather modelled by
4.4
where *u* and *v* are the rescaled coordinates as in equation (4.2), *k*_{u}=(1.5±0.1)×10^{−3} and *k*_{v}=(2.3±0.1)×10^{−3}. The errors have been estimated by assigning to fluctuations in *k*_{u} and *k*_{v} the measured standard deviation of the density of the 24 cylinders. The choice of a linear dependence is forced by the small number of samples in which the spare cylinder has been divided.

Instead of deriving *V* _{c} by expanding in multipoles a cylinder with a density given by equation (4.4), we have evaluated the effect of inhomogeneity by adding a negative point mass *m* at the centre of a homogeneous cylinder of mass *M*. By adjusting *M* and *m*, it is possible to match the first two coefficients of the correct multipole expansion. If *M*_{c} is the cylinder mass the monopole coefficients match if *M*_{c}=*M*+*m* while for the quadrupole, to the first order in *m*/*M*_{c}, we obtain
4.5
which gives *m*≃−8.8 g. The effect on *Φ*_{s} is (−50±10) μrad where the uncertainty has been estimated to 20% of the correction due to the uncertainty in the model for *ρ*(*u*,*v*).

#### (iii) Air density correction

It should be noted, with reference to figure 1, that when moving the masses from (*a*) to (*b*) an equivalent volume of air is also moved from (*b*) to (*a*). The effect on *Φ*_{s} of the air masses has been estimated in (−33±3) μrad. The uncertainty has been conservatively set to 10% of the correction since during the *G* measurement we did not record the ambient pressure.

### (b) Clouds dimensions

According to equation (2.2), we need to determine the classical paths of the atoms to evaluate *Φ*_{s}. For a single atom, we need to know the local value of *g* and the position and velocity at the time *t*_{0} of launch. For simulating a cloud, the initial conditions are randomly chosen from the phase space density distribution *p*(**r**,**v**,*t*_{0}).

The value of *g* has been measured with an absolute gravimeter with an uncertainty of about 10^{−8} *m* s^{−2}, more than adequate for our purposes. The result has been reported in [31]. Since we do not correct for tidal effects the uncertainty on *g* is actually of the order of 10^{−6} m s^{−2}. The effect on *Φ*_{s} is still negligible.

If we assume that the phase space density at *t*=*t*_{0} can be factored as *p*(**r**,**v**,*t*_{0})=*g*_{r}(**r**)*g*_{v}(**v**) where *g*_{r} and *g*_{v} are Gaussians, the temporal free evolution of *p*(**r**,**v**,*t*) for *t*>*t*_{0} is very simple: the averages **r**_{0}(*t*) and **v**_{0}(*t*) follow the classical equation of motion for a point particle, the variance does not change while . We can then measure the clouds at *t*>*t*_{0}, when it is more convenient, and propagate back in time to *t*_{0}. Actually, the effects of the state preparation and the Raman pulses must be taken into account making the computations slightly more involved. The clouds in the horizontal plane are reasonably symmetric so the variances can be described with just two values, in the vertical direction and in the horizontal plane.

An essential feature of our source masses is their high density. This greatly helps in reducing the required trajectory accuracy. If *V* _{s} were not large enough to be comparable with the term due to *γ*, a vertical position fluctuation *δz* in one of the two clouds centres would induce a fluctuation in the gradiometer output *δφ*=*kγδzT*^{2} which amounts to 1.2 rad m^{−1} in our experiment. If instead *V* _{s} can be made large enough to locally generate a gradient comparable with *γ* but with opposite sign, we can choose classical paths and mass positions that minimize the sensitivity of *Φ* to *δz*.

In a simulation, first we fix the mass positions in the close configuration and vary the initial conditions of the clouds' centres until *ϕ*_{p} is stationary for each cloud. We then fix the initial conditions, move the masses to the far configuration and adjust their vertical coordinates until *ϕ*_{p} is stationary again. It is not actually possible to have *ϕ*_{p} stationary for both *k*_{↑} and *k*_{↓} interferometers simultaneously and, for computational limits, we do not simulate a whole cloud in this minimization but just a single trajectory corresponding the cloud's centre. Nevertheless, eventually evaluating numerically *δΦ* as a function of *δz* the linear coefficient is reduced by at least a factor 5. The price to pay for this is that *δz* or a change in the vertical dimension of the cloud will add a bias to *Φ* due to the quadratic coefficient. In the horizontal plane, there is a similar situation since *V* _{s} has a parabolic minimum at *r*=0. The curvatures in the vertical direction and in the horizontal plane have opposite sign so, when varying the clouds sizes, there is a partial cancellation of the bias induced on *Φ* by the finite sizes of the clouds.

#### (i) Horizontal plane

In the horizontal plane, destructive fluorescence images of the clouds have been acquired with a CCD camera from two orthogonal directions when the clouds cross the detection beams. The images have been integrated in the *z* coordinate in order to obtain the horizontal profiles and these have been fitted with Gaussians. By measuring the clouds both on their way up and down in their parabolic flight, we obtain enough data to determine the horizontal components of **r**_{0} and **v**_{0}, i.e. **r**_{0r} and **v**_{0r}, respectively, and their standard deviations *σ*_{r}=(3.0±0.2) mm, *σ*^{l}_{vr}=(16±3) mm s^{−1} and *σ*^{u}_{vr}=(22±4) mm s^{−1} at *t*_{0}. The measured *r*_{0r} and *v*_{0r} have been taken into account in the simulation.

The values for *σ*_{vr}, corresponding to horizontal temperatures of about 3 μK and 5 μK for the lower and upper clouds, respectively, were also verified and confirmed by observing the Doppler width of counter-propagating Raman transitions in the horizontal plane.

#### (ii) Vertical direction

In the vertical direction, measurements for *σ*_{z} are obtained directly from the widths of the fluorescence detection signals recorded during the *G* measurement. After deconvolving for the effects of the finite size of the detection beams and the limited bandwidth of the photodiode, propagation at *t*_{0} gives four distinct values for *σ*_{z}, one for each of the combinations obtained from upper and lower clouds and *k*_{↑} and *k*_{↓} momenta, varying between (2.4±0.1) mm and (3.1±0.1) mm. For *σ*_{vz}, the value (2.84±0.10) mm s^{−1} has been assumed *a priori* since after the velocity selection and the interferometric pulses this result is weakly dependent from the initial values of *σ*_{vz}(*t*_{0}) for the two clouds, which have been measured and found to be consistent with the horizontal values. The uncertainty on *σ*_{vz} arises on the actual value of *σ*_{vz}(*t*_{0}). The assumptions on *σ*_{vz} are further discussed in the next section.

The initial position *r*_{0z} and velocity *v*_{0z} have been determined before the *G* measurement when adjusting the launch parameters to place the clouds on the optimum trajectories determined by simulation. The procedure has been repeated after the *G* measurement. The stability of *r*_{0z} and *v*_{0z} has also been monitored during the *G* measurement by checking the fluctuations in the times of arrivals of the clouds at detection.

The initial velocity has been determined from the time of detection and the time *t*_{a} at which the clouds reach the apogees of their parabolic flights since the vertical coordinates of the detection beams are known. To measure *t*_{a}, we apply only a single *π* pulse at a time *t*_{π} close to *t*_{a}. In these conditions, both *k*_{↑} and *k*_{↓} transitions are induced simultaneously so the detection signal for |*b*〉 splits in two peaks, separated at detection by a time *t*_{s}. It is easy to see that if *t*_{π}=*t*_{a} then *t*_{s}=2*v*_{r}/*g*. By scanning *t*_{π} and fitting to a straight line *t*_{s} versus *t*_{π} we can determine *t*_{a} with an uncertainty below 10^{−4} s. Once *v*_{0z} and *t*_{a} are known *r*_{0z} can be evaluated.

#### (iii) Simulation

In the simulation, the vertical velocity selection and the horizontal dependence on the intensity of the Raman beams must be taken into account. The effective Rabi frequency *Ω*_{R} of the Raman transition can be written as
4.6
where *Ω*_{0} is the unperturbed Rabi frequency, dependent on the Raman beams intensities *I*_{1} and *I*_{2} and Δ is the detuning from the Raman transition frequency. The probability *P*_{ab}(*t*) for a transition induced by a square pulse of length *t* of the Raman beams is [32]
4.7

For a *π* pulse, we have *t*=2*τ* and *Ω*_{0}2*τ*=*π* so
4.8

Equation (4.8) shows that only a narrow class of vertical velocities for which , i.e. |*v*|<*π*/(2*τk*)≃8.2 mm s^{−1} in our case, is effectively addressed by the *π* pulse. The complete sequence of Raman pulses, consisting of three velocity selection *π* pulses and the interferometer sequence, results into a velocity distribution well approximated by a Gaussian with *σ*_{vz} given as an *a priori* value in the previous section. Instead of starting with a the initial measured value of *σ*_{vz} and consider the change in the velocity distribution after each pulse, we model an equivalent expansion where the pulses are ignored and *σ*_{vz} is constant.

Since, however, *Ω*_{0} depends on *I*_{1} and *I*_{2}, there is a spatial dependence of *P*_{ab} on the horizontal coordinate *r*. We have measured *P*_{ab}(*r*) by scanning an iris across the Raman beams profile and by measuring the transfer efficiency of a pulse duration *t*=2*τ* as a function of the iris position. This function has been inserted in the simulation to obtain a relative contrast for each simulated trajectory. The contrasts have then been used as weights to evaluate *Φ*_{s}.

The uncertainties on the parameters describing *p*(**r**,**v**,*t*_{0}) give the largest contribution to the uncertainty on *Φ*_{s} by adding 39 μrad.

### (c) Detection

The detailed implementation of the state selective detection has already been reported in [17], and an analysis of the stability required on physical parameters such as lasers powers and frequencies can be found in [12]. Here, we consider only the effects of the detection geometry on *Φ*_{s}.

Detection of the atoms in the |*a*〉 and |*b*〉 states is implemented with two horizontal beams of resonant light 15×4 mm^{2} wide. The beams are vertically separated by about 20 mm.

Each beam has its own collection optics. The relative sensitivity in the horizontal plane of the collection optics has been modelled with a Gaussian *S*_{o}(*x*,*y*) with *σ*_{x}=(15±1) mm and *σ*_{y}=(6.3±0.1) mm centred with the intersection of the axis of the Raman beams within 1 mm both in the *x*- and *y*-directions.

The *x*- and *y*-axes are, respectively, perpendicular and parallel to the collection optics axis. The function *S*_{o}(*x*,*y*) has been determined by measuring the fluorescence of the thermal Rb vapour induced by the Raman beams, when tuned in resonance with the *D*2 transition, with the same scanning iris mentioned in the previous section.

The detection beams form an angle of 45° with the *x*-axis. Using a scanning blade the sensitivity, including saturation effects, as a function of the transverse coordinate *x*′, has been determined to be
4.9
with *c*=2.38±0.06 and *σ*_{x′}=(3.9±0.5) mm. The total detection sensitivity *S*_{d}=*S*_{o}(*x*,*y*)*S*_{b}(*x*′) has been used to compute *Φ*_{s} as a weighted average over individual trajectories as described above for the contrast dependence on the horizontal coordinates. The total uncertainty on *Φ*_{s} due to detection amounts to 8 μrad.

### (d) Elliptic fit bias

As explained in [12], the detection of the atoms in the |*a*〉 and |*b*〉 states does not have the same efficiency. Instead of performing a calibration which in our set-up would be difficult to obtain with the required accuracy, we have chosen to treat the efficiencies ratio as a fit parameter under the assumption that the correct ratio gives the smallest error on the fit to the ellipses. Moreover, our fitting routine implements a simple least-squares method and this is known to introduce a bias in *φ*. If the noise model is known, a Bayesian estimator can be used to provide more accurate values for *φ* [33].

We have checked our fitting procedure by generating four sets of 10^{5} ellipses, to which random noise was added, with parameters (number of points, RMS noise, contrast, *φ* and detection efficiencies ratio) as close as possible to the measured data. The average result for *Φ* as defined in equation (2.1) was at (7±2) μrad above the expected value so the same correction must be applied to *Φ*_{s} before comparing it to *Φ*_{m}.

### (e) Earth rotation compensation

In the laboratory reference frame, the Coriolis acceleration adds to the gravimeter phase shift *ϕ* a term
4.10
where *Ω* is the angular velocity of the Earth, and **v** is the cloud velocity. This affects also *φ* if there is a component *V* _{E} of **V**=**v**^{u}−**v**^{l} along the East direction. A bias can be added to *Φ* if *V* _{E} or the direction of **k** change when moving the masses between the close and far configurations due to a deformation in the apparatus.

As reported in [17], we have measured Δ*V* _{E}=(*V* _{EC}−*V* _{EF})=(17±11) μm s^{−1} while the direction of **k** has been actively stabilized, reducing its variation to negligible levels during the *G* measurement, improving on [12]. The effective component of *Ω* along the North direction seen by the atoms can be cancelled by counter-rotating **k** during the interferometric sequence. We estimate for our experiment . The combined effect of non-zero Δ*V* _{E} and leads to an uncertainty on *Φ* of 20 μrad.

### (f) Other effects

Here, we include the uncertainties in the linear acceleration gradient *γ*, the wavenumber of the Raman transition *k* and the main frequency reference of the experiment.

The gradient *γ* has been measured using the gradiometer itself after unloading the cylinders. The result is (−3.11±0.01)×10^{−6} s^{−2}, within 1% of the nominal value of −3.09×10^{−6} s^{−2}. This produces a negligible correction on *Φ*_{s} of (1±0.5) μrad.

It is very simple to measure *k* with a relative uncertainty smaller than 10^{−7}. Since the relative uncertainty induced on *G* is the same, it can be neglected. The change on *k* due to Earth rotation compensation is also negligible.

Even if our 10 MHz main frequency reference was not disciplined by a GPS receiver, we checked its value before and after the *G* measurement and found that its drift can be safely neglected.

### (g) Summary

A simulated value *Φ*_{s}=0.548105(1) rad, where the uncertainty is due only to the statistics on the number of simulated trajectories, is obtained when assuming for *G* the present CODATA value. After applying the corrections with their uncertainties discussed above and summarized in table 1, we obtain *Φ*_{s}=0.548029(51) rad.

By comparing *Φ*_{s} with the measured value *Φ*_{m} given above, we can finally extract *G*=6.67191(77)(62)×10^{−11} m^{3} Kg^{−1} s^{−2} where the first uncertainty, about 116 ppm, arises from the statistical error on *Φ*_{m} while the second uncertainty, about 92 ppm, is due to the simulation and the source masses.

## 5. Outlook

In this section, we would like to discuss, based on the experience gained with this experiment, what are the prospects for an improvement on our present uncertainty implementing only technologies that have already been demonstrated.

Clearly, an improvement on the statistical uncertainty will require very long acquisition periods so it is necessary to increase the reliability of the gradiometer. Our experiment is fully computer-controlled and presently can run unattended for up to 40 h at a time in the best days. A second generation apparatus should be capable of self-operating for weeks rather than for days. The improvements required might be tedious to implement but rely on well-consolidated techniques. Increasing the cycle rate of the experiment or the number of atoms in each cloud will also allow to reduce the integration time. All these improvements will still reduce the relative uncertainty on *Φ*_{m} as , where *N* is the number of detected atoms.

Trying to improve beyond the standard quantum limit towards the 1/*N* Heisenberg limit is an active field of research (see [34] for a recent result). Presently, however it seems a long-term experimental goal to beat the standard quantum limit when the number of atoms in each cloud exceeds 10^{5}.

Reducing type B uncertainties will obviously require a tighter control on some experimental parameters, especially laser powers and frequencies and compensation of the Earth rotation. They are expected to be time consuming but straightforward improvements. The uncertainty related to the cloud sizes and trajectories can be better controlled by using colder atoms and smaller clouds as implemented in [35].

The problem of density inhomogeneity for the source masses can be tackled adopting different materials.

Interferometers exchanging at every pulse a momentum in excess of have been demonstrated [36] even if at the expense of contrast. Recent results however show that pulses exchanging up are possible with excellent contrast [37,38]. Since *Φ*_{s} is proportional to the momentum transfer, implementing these beam splitters will allow to obtain the same signal and signal-to-noise ratio while reducing the source masses. Smaller source masses are indeed expected to reduce the inhomogeneity problem and allow a better positioning.

To summarize, present-day technology is already enabling the design of an experiment based on atomic interferometry that could significantly reduce the main sources of uncertainty presently limiting our set-up.

## 6. Conclusion

We have measured *G* with a gradiometer based on cold atoms and, for the first time, obtained a value with an uncertainty comparable with the one reported by CODATA. An improvement on the total uncertainty of at least a factor of five is still necessary to reduce the uncertainty to a level comparable with the one quoted by experiments based on macroscopic sensing masses.

Reducing the statistical error by this amount will require an increase of at least a factor 25 in the number of detected atoms since our experiment is essentially limited by quantum projection noise. Reducing type B uncertainties at the tens of ppm level will also demand a tighter control on many experimental parameters, especially atomic clouds sizes and trajectories and density homogeneity of the source masses.

We are confident that all the improvements necessary to measure *G* with an uncertainty in the 10 ppm range with atomic interferometry are feasible in the near future. This can help to solve the *G* puzzle.

## Funding statement

This experiment has been funded by the Istituto Nazionale di Fisica Nucleare (INFN).

## Acknowledgements

The authors gratefully acknowledge the contributions of A. Cecchetti and B. Dulach (INFN-LNF) to the design of the mechanical structure of the source masses supports and the source masses positioning system, and of A. Peuto, A. Malengo and S. Pettorruso (INRIM) to density tests and mass measurements on the cylinders. G.M.T. acknowledges seminal discussions with M. A. Kasevich and J. Faller and useful suggestions by A. Peters in the initial phase of the experiment. The authors would also like to thank all those that worked at this experiment in an early stage particularly, in chronological order, J. Stuhler, T. Petelski, M. Fattori, M. de Angelis, G. Lamporesi, A. Bertoldi and Y. H. Lien.

## Footnotes

One contribution of 13 to a Theo Murphy Meeting Issue ‘The Newtonian constant of gravitation, a constant too difficult to measure?’

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.