Experimental overview

A light-pulse atom interferometer is conceptually similar to an optical interferometer, with the roles of light and matter interchanged. Atoms, acting as matter waves, are subjected to a sequence of light pulses that impart momentum to them, acting analogously to mirrors and beamsplitters. Applying a light pulse for an appropriate length of time will cause a transition between the ground and excited states of an atom, accompanied by the absorption and stimulated emission of a photon. Such a pulse, commonly referred to as a π pulse, acts as an atom optic mirror owing to the momentum that is transferred. Similarly, tuning the light pulse such that it has only a 50% transition probability, commonly referred to as a π/2 pulse, acts as a beamsplitter through providing a momentum kick to only half of the atomic probability distribution. A matter-wave equivalent of the optical Mach–Zehnder interferometer can then be created through applying a π/2–π–π/2 pulse sequence with an evolution time, T, between the pulses. The resulting space-time area enclosed by the atomic trajectories (Extended Data Fig. 1b) is proportional to the local gravitational acceleration, which can then be measured from the relative population of the two atomic states after the final pulse.

A gravity gradiometer utilizes two such interferometers offset vertically and probed simultaneously with the same pulse sequence. This suppresses common-mode effects, such as noise from vibration or phase changes due to variations in tilt with respect to the gravity of the Earth, which are indistinguishable from the gravity anomalies of interest according to Einstein’s equivalence principle. Our device consists of two subunits (Extended Data Fig. 1a), a sensor head and a control system, with light and electrical signals transferred through a 5-m umbilical. The gradiometer is shown in Extended Data Fig. 2, with an overview of its size, weight and power characteristics.

The sensor head features a vacuum system with dual MOT preparation and interrogation regions in an hourglass configuration, with all light delivered to the atoms through on-axis counter-oriented telescopes. The light is delivered in each direction, with portions of the beam being redirected towards the atom-trapping region using in-vacuum mirrors, to form the radial cooling beams in each MOT. The central portion passes through, such that each input provides the vertical laser cooling beam in a given direction for both MOTs. This makes all fluctuations in intensity common for the radial cooling beams (preventing lateral offsets), and, through use of a Gaussian beam shape, provides a higher intensity for the vertical beams to better saturate the radiation pressure force in this direction. This results in a greatly improved stability and robustness of the laser cooling process, reducing fluctuations in temperature or atom cloud position (Fig. 1b) without the need for excessive laser powers that would inhibit field operation. In a comparable test system, this provided a reduction in average cloud centre-of-mass motion to (0.14 ± 0.09) mm as compared to (1.19 ± 0.86) mm over an hour in similar conditions with a six-beam MOT. Both MOT regions have two coils, each formed of 92 turns of 1-mm-Kapton-coated copper wire wound around an aluminium former (fixed using epoxy), with a slit to prevent eddy currents. The coils have a radius of 43 mm and separation of 56 mm, to produce a linear field gradient of 12.5 G cm⁻¹ at a driving current of 2.5 A. These are located around the vacuum system, such that the strong magnetic field axis of their quadrupole field is along the direction of travel of the cooling beam axis. In addition, two sets of rectangular coil pairs, each having 20 turns, are located around the MOT regions. These have a separation of 100 mm, and dimensions of 320 mm in the vertical and 90 mm in the horizontal, and can be used to compensate residual magnetic fields, or apply offsets. In practice, no compensation fields are used for the molasses phase. In the lower chamber, one coil pair is used to apply a 0.63 G field to adjust the atom cloud horizontal position by approximately 0.5 mm in the MOT phase, improving the interferometer contrast. A bias coil⁴² is positioned around the system to define a quantization axis and remove degeneracy between magnetic sublevels, with other coils being switched off after the magneto-optical trapping phase. This has a variable pitch shape to account for edge effects and improve field uniformity over the atom interferometry region. The system is enclosed in a magnetic shield that provides 25 dB attenuation of the external field. The in situ magnetic field profile is measured (through spectroscopy of the Raman transition) as being homogeneous to below 5% across the atom interferometry region, limited by internal magnetic field sources from vacuum pumps.

The laser system consists of telecom lasers that are frequency doubled to 780 nm, to be near the D2 line of rubidium-87 (refs. ^43,44). The light for laser cooling is generated by passing the laser output through an electro-optic modulator (EOM) and generating a sideband at a frequency of approximately 1.2 GHz output from the carrier. This is used to provide a locking signal using the ({|F}=3 > to |{F}^{{prime} }=4 > ) transition in rubidium-85, placing the carrier frequency such that it is tuneable around resonance with the ({|F}=2 > to |{F}^{{prime} }=3 > ) transition of rubidium-87 to provide the cooling light. A separate EOM is used to provide repumping light resonant with the ({|F}=1 > to |{F}^{{prime} }=2 > ) transition. Atom interferometry is realized through two-photon stimulated Raman transitions. The Raman laser used to drive these has a linewidth of 73 kHz and is locked with an offset of 1.9 GHz to the ({|F}=2 > to |{F}^{{prime} }=3 > ) transition. The second Raman frequency is generated using a pair of EOMs operating at 6.835 GHz. Performing the differential measurement suppresses phase noise that may arise owing to optical path-length changes between the two Raman beams (such as those due to vibration and thermally induced changes in the refractive index of fibres). This allows the two beams to be delivered independently without the need for a phase lock between them, facilitating an implementation in which the modulated spectrum is applied to only one of the input beams. This avoids parasitic Raman transitions that give rise to systematic offsets and dephasing when using conventional modulation-based schemes, such as those including a retro-reflected beam³¹. To realize a practical implementation of space-time area reversal³⁰, also known as wavevector reversal, the system has an EOM in each input direction of the Raman beams, and the modulation signal is applied to one arm in each measurement. This allows the direction of the momentum kick imparted to the atoms to be changed between measurements, by changing which arm the modulation signal is applied to using a radiofrequency switch (see Extended Data Fig. 1). The contributions to the interferometer phases due to acceleration under gravity are sensitive to the direction of the recoil imparted by the light, whereas those arising from many other effects, such as those due to magnetic fields, are not. This allows these effects to be removed when interleaved measurements are performed in the two recoil directions.

The light is delivered to the sensor head using polarization-maintaining optical fibres, with separate fibres for the cooling and Raman beams. These fibres deliver the light to optical telescopes that collimate the light at the desired beam size. The cooling beams have a waist of 24 mm, and contain a typical maximum power of 130 mW. These impinge on the in-vacuum mirrors, which are 15-mm right-angle prisms (Thorlabs, MRA15-E03), to deliver the horizontal cooling beams. The mirrors are mounted to a titanium structure (attached using Epo-Tek H21D adhesive) in a cross configuration such that there is a 15-mm aperture in their centre. The central portion of the cooling beams passes through these apertures to provide the sixth beam required for the opposite MOT. The Raman beams are overlapped with the cooling beams using a polarizing beamsplitter cube, such that they are then delivered along the same beam axis as the cooling light. The Raman beams, each containing a typical maximum power of 300 mW, have their waist set to 6.2 mm to limit aperturing and diffraction on the central aperture of the in-vacuum mirrors, allowing the Raman beams to pass through the system without being redirected by the prisms. Although aperturing is limited on the mirrors in the current instrument, it may be desirable to use a larger Raman beam than the aperture in more compact systems or those aiming to further reduce dephasing induced by laser beam inhomogeneity. Diffraction from the aperture would need to be given due consideration if pursuing this, as would the potential for further light shifts due to, in this case, one interferometer seeing extra light fields from mirror reflections. The polarization of the light is set to the appropriate configuration for cooling or driving Raman transitions through use of voltage-controlled variable retarder plates in the upper and lower telescopes used to deliver the light. The intensity of the Raman beams is actively stabilized using feedback from a photodiode to control acousto-optic modulators, which are also used to produce the laser pulses.

The experimental sequence starts by collecting approximately 10⁸ rubidium-87 atoms in each MOT from a background vapour over 1–1.5 s. Molasses cooling is then used to reduce the upper- and lower-cloud temperatures to (2.86 ± 0.09) μK and (3.70 ± 0.20) μK, respectively (see Fig. 1b). The differences in temperature arise from differences in local residual magnetic fields, arising primarily from the magnetic shield geometry, and small differences in optical alignment. Optical state and velocity selection is performed to select only atoms in the ({|F}=1,{m}_{{rm{F}}}=0 > ) magnetic sublevel and desired velocity class. This is achieved through application of π pulses and a series of blow-away pulses to remove atoms in undesired states and velocity classes. Atom interferometry is then performed with a pulse separation of T = 85 ms and π-pulse length of 4 µs. The interferometers are read out using bistate fluorescence detection to determine the atomic state population ratios of the |F = 2> and |F = 1> ground states, for which (2.7 ± 0.1) × 10⁵ and (1.7 ± 0.1) × 10⁵ atoms participate in the upper and lower interferometers, respectively, with a typical measurement rate of 0.7 Hz. The differential phase, from which the gravity gradient is derived, is extracted by plotting the upper interferometer outputs against the lower interferometer outputs, to form a Lissajous plot as shown in the inset of Fig. 2. In addition to random noise arising from vibration, we add a deliberate random phase value, from between 0 and 2π, to the final pulse of the interferometer. At ellipse phases that do not correspond to a circle, a clustering of points around the extremal points of the ellipse is visible even for uniform noise.

The quantum projection noise of the system based on the participating atom number is approximately 44 E/√Hz. The total noise budget includes contributions from further terms, and is shown in Extended Data Table 1, alongside relevant systematics observed during the survey. The noise budget was investigated through computer simulation of noise processes, compared to experimental data, and ellipse fitting.

Survey practice and processing of the measurement data

For each measurement on the survey, 600 runs of the atom interferometer were typically taken with the sensor head in one location (with the horizontal position being measured using a total station, Leica TS15, and the vertical position from the road surface being approximately 0.5 m for the lower sensor and 1.5 m for the upper sensor), giving twelve 25-point ellipses in each of the interferometer directions and therefore 12 separate estimates of the gravity gradient. Repeat measurements were taken on each measurement position, with typically three points on each position. A measurement was taken at a base station between each measurement point, with the final base-station measurement for one location used as the first for the next. The quality of fitting to each ellipse was identified using the error metric, (varepsilon ), defined as

in which N is the number of data points, L is the minimum distance between each data point and a point on the ellipse, and a and c correspond to an ellipse defined parametrically by equations (x=a{rm{sin }}theta +b) and (y=c{rm{sin }}left(theta +varphi right)+d), respectively. Errors in the ellipse fitting are sensitive to changes in the ellipse opening angle⁴⁷. On the basis of numerical simulations, we estimate this effect to be less than a few parts in one thousand; therefore, a 100 E change would be subject to an error of less than 0.5 E. Such errors are therefore small compared to other errors. Such a 100 E change in gradient would correspond to an 11.6 mrad change in the ellipse shape. This phase shift can be compared to a 2π measurement range, meaning that measurement range of the instrument in this configuration is relevant to the majority of practical features of interest (these being typically below 400 E).

Ellipse fits found to have (varepsilon > 0.05) were automatically discarded. This resulted in 98.4% of all data being usable in normal operation, representing a favourable data up time compared to that of similar conventional geophysical devices.

To process the data, a straight line was fitted to the base-station points, with this line then being subtracted from all data points. This is standard practice to remove drift in geophysical surveys. The leading source of drift is believed to be due to the a.c. Stark shift, with this also being relevant owing to the difference in the temperature of the two clouds. The gravity gradient value is then taken as the average of the measurement points, resulting in an estimate of the difference in gradient between the measurement location and the base station. Furthermore, the variations in the data points are used to make an estimate of the error in the difference value. When multiple measurements from the same location were combined, a weighted average was used, giving less weight to measurements with greater errors. The weighting factor is proportional to the reciprocal of the variance of each measurement⁴⁸. The data, as shown in Fig. 3a, are not corrected for terrain or effects such as tides. Tidal effects are not corrected, being negligible through the differential measurement of the gravity gradient.

The average of the gravity gradient error found across the measurement positions of the survey is 17.9 E. Comparing this to an approximate signal size of 150 E gives an approximate signal-to-noise ratio of 8.

Inference from gravity gradiometer data

Bayesian inference is a framework within which prior beliefs can be updated with information contained in data. For a model parameter vector ((theta )) and a data vector ((d))

in which (pleft(d|theta right)) is the likelihood, (pleft(theta right)) is the prior, (pleft(dright)) is a normalization constant and (pleft(theta |dright)) is the posterior distribution.

The likelihood function provides the misfit between the measured data, (d), and the modelled data values calculated from the model parameter vector, (theta ). The model used here is that of a three-dimensional cuboid³⁵; the free model parameters are shown in Extended Data Fig. 3, along with the functional form of the respective prior distributions. The rationale behind the chosen prior distributions is detailed in Extended Data Table 2. The total uncertainty for each measurement point is calculated using the Pythagorean sum of the standard error and the model uncertainty random variable multiplied by the average of the standard error across all of the measurement positions.

The probabilistic Python package pymc3 (ref. ⁴⁹) is used to implement the cuboid model, define the model parameter prior distributions and sample the posterior distribution, using a no U-turn sampler⁵⁰. Extended Data Fig. 4 shows the Bayesian posterior distribution for select model parameters.

The parameter posterior distributions represent the updated beliefs about the model parameters, given the measurement data. To aid interpretation of the posterior distribution, the POE³⁶ is calculated, which represents the spatial probability of the anomaly underground, given the model and prior distributions (as shown in Fig. 3c). The horizontal position of the tunnel centre is determined as (0.19 ± 0.19) m along the survey line, with the distribution being approximately Gaussian. The depth from the origin, defined in the vertical using the lowest point on the survey line, to the centre is (1.7 −0.59/+2.3) m. At the horizontal position of the tunnel, the distance to the surface from the origin is approximately 0.19 m, meaning that the total distance from the surface to the tunnel centre is (1.89 −0.59/+2.3) m. From the tunnel geometry, this places the top of the tunnel at approximately 0.89 m depth from the surface.

The signals arising from local features are used to create a distinct site model. This is used to provide an estimate of the expected shape of the gravity gradient signal over the site, for comparison with the inference output. These features include the tunnel of interest, basements from nearby buildings, walls and a drain. They are shown in the scale drawing of Fig. 3b.