Method for Scanning an Object Using a Gravimeter

Abstract

Methods of scanning an object using a pendulum gravimeter are disclosed. The pendulum gravimeter may include an interferometer, such as a Sagnac interferometer, to determine a displacement on the pendulum by way of a mirror attached to the pendulum. Scanning of the object may be performed in 1D, 2D, or 3D, and may result in an image of the object. In another aspect, a mass may be tracked while in motion using a pendulum gravimeter to detect the gravitational attraction of the object.

Description

FIELD OF THE DISCLOSURE

The present disclosure relates to measuring gravitational attraction.

BACKGROUND OF THE DISCLOSURE

There is great interest in measuring gravitational fields nearly perpendicular to the Earth's gravitational field. A sensitive gravimetric camera could “see” a person walking around through hundreds of meters of concrete or monitor shipping crates for hazardous materials. Currently, precision gravimetry tends to have relatively low resolution and is used primarily for large scale structures. Gravimetry is used extensively in mapping the earth's local gravity, oil and gas exploration, mining, mapping temporal geological shifts, the determination of Newton's gravitational constant and gravitationally imaging opaque systems. Precision of the order of 1 mGal (1 Gal=1 cm/s²) to 1 μGal are often used for mapping geological variations. Both relative and absolute measurements are employed. A standard in the industry for absolute gravimetry is measuring interference fringes due to the free-fall of a corner cube in one arm of a Mach-Zehnder interferometer with a sensitivity of 100 μGal/Hz^1/2. Another competing gravimetric technology employs atomic interferometry achieving a resolution of 100 nGal after two days of integration. At the forefront, superconducting spheres suspended in the field of a superconducting coil have achieved 3 nGal resolution after one month integration or 1 nGal after one year. The current standard has been set by Kasevich's group (Stanford), which has achieved 1 nGal per shot (2 seconds) in their recently built 10 m atomic fountain and is working on a squeezed atom implementation with the potential for 20 dB improvement per shot.

BRIEF SUMMARY OF THE DISCLOSURE

Embodiments of precision gravimeters are disclosed herein. In a first embodiment, a gravimeter is based on a simple pendulum with a few 10's of ms response time. A shot noise limited acceleration sensitivity of 100 fm/s²/Hz^1/2 or 10 pGal/Hz^1/2 can be achieved with a few mW of low-coherence laser power. The gravimeter is a physical pendulum which is allowed to oscillate in only a single plane. The pendulum and interferometric readout measure relative gravitational fields perpendicular to the Earth's gravitational field using the fringes in the dark port of a displaced Sagnac interferometer operated in the inverse weak value regime. The system is small, lightweight, simple to build, easily portable, inexpensive, and requires no cryostat, or cold atom technology. This gravimeter can be scaled to several detectors to build a gravity camera. In this way, a one kg mass moving around at a few meters can be remotely tracked using only the gravitational attraction of the mass.

In a second embodiment, an gravimeter uses a long-period torsion pendulum. The shot noise limited predicted sensitivities are a function of the restoring torque or alternatively considered as the restoring period. A device according to this embodiment is able to remotely image precision gravimetric density profiles in nontrivial mass distributions (e.g., tumors in-situ). With relatively reasonable device parameters, the shot noise limited sensitivity of the torsion balance embodiment is sufficient to measure the center of mass of a 1 gram mass 10 cm away with picometer resolution. Used in a near-field configuration, the device may be used to, for example, measure the curvature of space due to the presence of 1 Watt of energy

$\left(i.e.,m=\frac{E}{c^2}10\mathrm{fg}\right)$

at a distance of a few mm.

DESCRIPTION OF THE DRAWINGS

For a fuller understanding of the nature and objects of the disclosure, reference should be made to the following detailed description taken in conjunction with the accompanying drawings, in which:

FIG. 1: In a first embodiment of the present disclosure, light from a fiber-coupled laser enters the input port of a 50/50 beam splitter of a displaced Sagnac interferometer. Clockwise and counterclockwise beams strike different heights of a mirror attached to a rigid pendulum. A tilt in the pendulum induces a relative phase between the two paths. The phase causes the shift of an interference fringe at the output of the interferometer. The interferometer is tuned to be close to the dark port with a single dark fringe of a bimodal interference pattern, which is referred to herein as the “inverse weak value regime.” This interference configuration leads to dark fringe that moves rapidly for small phase shifts.

FIG. 2: A Sagnac interferometer with a torsion pendulum integrated as a gravity sensor.

FIG. 3: A torsion pendulum according to another embodiment of the present disclosure is formed by attaching two mass m with a rigid, massless rod, of length 2 and fixing the center of mass in place, able to rotate only. A nearby mass M creates a torque on this torsion pendulum, causing it to rotate to an equilibrium angle δθ₀, which is detected optically. Lengths, angles, and mass labels are shown in the figure.

FIGS. 4A and 4B: We plot the static deflection angle θ (Eq. (9)) for the torque in the one-armed torsion pendulum (Eq. (18)) (4A), and the two-armed torsion pendulum (Eq. (19)) (4B), as shown on the colorbars, as a function of target distance r₀ (x-axis) and target mass M (y-axis). The plot, given as a log-log-log density plot, shows three decades of distance, and six of mass. The test mass M is placed at the point of optimal sensitivity for each device. The parameters in Table I are used for these plots.

FIGS. 5A and 5B: Contour plots of the angular displacement θ due to a static mass M=100 kg placed in the plane of a 1-armed (5A) and 2-armed pendulum (5B), as a function of r₀ and θ. The Figures show the entire angular dependence θ, and show values of r₀ ranging from 100 m to 5 km in the radial direction. The + shading legend denotes an angular displacement in the +θ (CCW) direction, while the − shading legend denotes dection in the −θ (CW) direction. Angular blindspots are at the juncture of the two shade gradients, where the deflection is zero, no matter the value of r₀ or M. Numerical values for the pendulum correspond to those shown in Table I. The − shading legend corresponds to the lower half (5A) and the upper left and lower right (5B), and the + shading legend corresponds to the upper half (5A) and the upper right and lower left (5B).

FIG. 6: Sample of a one dimensional image, given here as a piece-wise constant function.

FIG. 7: Predicted net displacement signal of the pendulum for 1-D mass as a function of position with 10% Gaussian noise added. The signal is given proportional to the net displacement angle. This is modeled as the convolution of the point spread function with the data given in FIG. 6, plus additive Gaussian noise to model environmental noise in the data collection process. The standard deviation is about 10% of the signal. Angular displacement is plotted versus relative position of the object to the imager.

FIG. 8: Result of the deconvolution algorithm to reconstruct the original mass distribution.

FIG. 9: The net torque of a point mass is proportional to g(,x)−g(−,x), which is plotted here versus x. This is also called the “point spread function” in Fourier imaging techniques. Note it is an odd function, and most sensitive at x=0.

FIG. 10: Sample of a two dimensional image, given here as a piece-wise constant function giving the letters “U R.”

FIG. 11: The signal is given proportional to the net displacement angle as the detector is scanned over the imaged distribution. This is modeled as the convolution of the 2D point spread function with the data given in FIG. 10, plus additive Gaussian noise to model environmental noise in the data collection process. The standard deviation is about 10% of the signal.

FIG. 12: Result of the deconvolution algorithm to reconstruct the original 2D mass distribution. Note the algorithm has reflected the image, which can be reversed by an appropriate transformation.

DETAILED DESCRIPTION OF THE DISCLOSURE

Precision gravimeters are disclosed herein, suitable for studying gravitational fields perpendicular to the earth's gravitational field. The gravimeters may be used for fundamental studies of gravity, high resolution bulk imaging, as a part of a gravity camera, etc. Shot noise limited measurements with 1 mW of laser power predict precision of the order of 10 pGal/Hz^1/2.

With reference to FIG. 1, the present disclosure may be embodied as a gravimeter 10 having a pendulum 20 and an interferometer 30 configured to measure a displacement of the pendulum 20. As is known in the art, the pendulum 20 may be considered to form a part of the interferometer 30, but for convenience the pendulum 20 will be generally referred to herein separate from the interferometer 30. The pendulum 20 may be a simple pendulum or a torsion pendulum—each of which is further described below. The pendulum 20 is configured to move in only a single plane. For example, where the pendulum 20 is a simple pendulum, the pendulum is configured to swing in a single vertical plane (see, e.g., FIG. 1). The pendulum 20 is damped, and in some embodiments the pendulum 20 is critically damped such that it moves to a new position in the shortest time possible without overshooting the new position. For example, in some embodiments, the pendulum comprises a swinging mass made from copper, and the pendulum is damped using a magnet (e.g., a rare-earth magnet, etc.)

The inferometer 30 includes a laser 32 for providing a laser beam A. A beam splitter 34 of the inferometer 30 is configured to split the laser beam A into a first beam B and a second beam C. Where a Sagnac inferometer is used, the first beam B and second beam C propagate in opposite directions and are each directed around a looped path to impinge on a mirror 22 affixed to the pendulum 20. The first beam B impinges on the mirror 22 at a first location to provide a first reflected beam D, and the second beam C impinges on the mirror 22 at a second location to provide a second reflected beam E. For example, where the pendulum is configured to move in a vertical plane, such as a simple pendulum, the first location and second location are arranged vertically with respect to each other. In embodiments where the pendulum is configured to move in a horizontal plane, such as a torsion pendulum, the first location and second location are arranged horizontally with respect to each other.

A photodetector 36 is configured to receive the first reflected beam D and the second reflected beam E. A controller 40 is in electrical communication with the photodetector 36. The controller 40 is configured to determine a displacement of the pendulum 20 based on interference between the first reflected beam D and the second reflected beam E. The controller 40 may be, for example, a programmable controller, programmed to determined the displacement of the pendulum.

In some embodiments, the pendulum is a torsion pendulum comprising a mass assembly. In some embodiments, the mass assembly is a rigid rod. In some embodiments, the mass assembly is made up of at least two masses rigidly coupled to each other, such as, for example, coupled using a rigid rod. The mass assembly is balanced on a fulcrum and configured to rotate about the fulcrum. The mass is biased to a first rotational position and configured such that a gravitational attraction (e.g., from an object) will rotate the mass assembly away from the first rotational position. For example, the mass assembly may be suspended by a wire, which acts as the fulcrum at the location where the wire is affixed to the mass assembly. In this way, the mass assembly is balanced while suspended by the wire, and able to rotate by twisting the suspension wire. Such twisting (rotational) motion also creates a restorative torque in the wire.

I. Weak Values

Precision gravimetry using weak value amplification techniques is discussed herein. Weak values were born through asking fundamental questions about quantum measurement limits. Unlike expectation values, weak values consider a normalized expectation of an operator (e.g., the Pauli operator Â=|++|−|−−|) using pre- and post-selected quantum states ψ_i,f

$\begin{array}{cc}{A}_w=\frac{〈{\psi }_fA{\psi }_i〉}{〈{\psi }_f|{\psi }_i〉}.& \left(1\right)\end{array}$

Because weak values can be much larger than their respective expectation values when ψ_f|ψ_i→0, they have been used to amplify small effects. Weak value amplification has been shown to be exceptionally valuable in suppressing technical noise in precision measurements. While these techniques do not beat the shot noise limit (with some exceptions), they can come close to reaching it because of the dramatically suppressed technical noise. Of particular interest is the recent inverse weak value work where an angular tilt measurement noise floor of 200 frad Hz^−1/2 was achieved. Remarkably, this sensitivity was for signals down to 1 Hz, where noise suppression can be incredibly difficult. This tilt corresponds to a displacement of less than a hair's breadth at the distance of the moon in one second of measurement time using only a few milliwatts of laser power. Using these techniques for precision gravimetry, even at the classical optical fundamental limits, they would push gravimetric sensitivity by several orders of magnitude beyond the state-of-the-art.

II. Inverse Weak Value Interferometry

A specialized weak values interferometer employs a laser beam of transverse width σ in a Sagnac interferometer. The laser beam enters a beamsplitter and propagates in opposite directions around a Sagnac interferometer. When the beams recombine, a small relative phase ϕ between the two returning beams causes a small transverse tilt, k, of the mirror attached to the pendulum to be amplified in the (nearly) dark port of the beamsplitter. The weak value limit occurs when ϕ»kσ. For this particular setup, the amplification of the small transverse tilt k shows up in the dark port beam as a spatial shift by kσ²/ϕ. In terms of weak values, the pre-selected state ψ_i is the field after passing through beamsplitter the first time. The post-selected state ψ_f is set by a combination of a phase shift and the second pass through the beamsplitter. The tilt of the mirror yields a small amount of which-path information of a photon in the interferometer, which is a weak measurement of σ_z.

Conversely, for an inverse weak values experiment, the parameters satisfy the inequality ϕ«kσ. In this case, we fix transverse tilt k and use the known k to amplify a small unknown phase. In this latter experimental regime, the interference pattern of the two beams in the dark port is now a bimodal distribution with a dark fringe at the center of the interference pattern for ϕ=0. The dark fringe moves rapidly with small changes in relative phase. These phase shifts are determined by measuring the relative intensity of the left versus right side of the interference pattern via a split detector. The amplification of the phase in this inverse weak value regime is given by the mean shift of the beam in the dark port, ϕ/k, which now is proportional to the inverse weak value A_w⁻¹.

A displaced Sagnac interferometer may be used to measure the relative phase shift ϕ for this inverse weak value regime. In a displaced Sagnac interferometer, two beams propagate with a transverse displacement (albeit parallel) in opposite directions. A small tilt of a mirror inside the interferometer causes a relative phase between the two paths since the path length increases for one path and decreases for the other. At the output port of the interferometer, the beams are brought back together to interfere with each other.

The shot noise limited angular resolution can be understood from a geometric argument. The relative phase between the two paths goes as δϕ=2√{square root over (2)}πLθ/λ, where √{square root over (2)} comes from impinging at 45 degrees, L is the distance between the centers of the beams propagating in opposite directions (see FIG. 2), λ is the wavelength of the laser light, and θ is the tilt angle of the mirror that we are interested in determining. Assuming the phase can be determined with shot noise limited sensitivity

$\mathrm{\Delta \phi }=\frac{1}{2\sqrt{N}},$

we find.

$\begin{array}{cc}{\mathrm{\Delta \theta }}_{\mathrm{SN}}=\frac{1}{4\sqrt{2\pi }}\frac{\lambda }{L\sqrt{N}},& \left(2\right)\end{array}$

where N is the number of detected photons. Using L=1 cm, approximately 3 mW of laser power, and a wavelength of 500 nm, we achieve a shot noise limited angular sensitivity of 30 frad/Hz^1/2. The inverse weak value method of readout for the optical phase ϕ can achieve this shot-noise limited sensitivity, up to a factor of √{square root over (π)}/2 associated with the resolution loss on the split detector.

This same inverse weak value setup with a displaced Sagnac is further described below. Further consideration is given to the physical limitations and sensing capabilities when the tilt mirror in the simple pendulum is replaced with a mirror rigidly connected to a torsion pendulum as shown in FIG. 2.

Simple Pendulum Gravimeter

An example of a pendulum gravimeter according to an embodiment of the present disclosure is shown in FIG. 1. A pendulum is in equilibrium with the static surrounding mass distribution. The introduction of another mass changes the equilibrium position of the pendulum (this is referred to as a relative gravimeter rather than an absolute gravimeter). If the pendulum is critically damped, it moves to its new equilibrium in the shortest time possible. The difficulty arises in determining the new relative angle of the pendulum due to the introduced mass.

A displaced Sagnac interferometer operated in the inverse weak value regime may be used to read out the tilt of a critically damped physical pendulum. A laser beam traverses two trajectories in the interferometer. One trajectory strikes the mirror on the pendulum at one height and the other beam strikes the mirror at a different height. The pendulum has a mass m and the height difference between the beams on the mirror is L. The pendulum design is such that it can only swing in a single plane. For convenience in discussion the operation of this embodiment, the pendulum was considered to be acted upon only by Earth's gravity and by the force of a distant remote object of mass M , which induces a force F perpendicular to Earth's gravity on the test mass of the pendulum. In the presence of the external object, the pendulum mass has an angle θ with respect to the vertical. The equilibrium condition is then F cos θ=mg sin θ. The small angle approximation was used to arrive at

$F=\frac{\mathrm{mg}\phi }{2\mathrm{kL}}.$

Rather than measure the angular displacement θ of the mass, the relative distance the mass shifts Lθ was measured. For normal reflection from the pendulum, the relative phase is given by ϕ=2kLθ, where k is the wavenumber. A force equation

$F=\frac{\mathrm{mg}\phi }{2\mathrm{kL}}$

was obtained. For shot noise limited measurements, an approximation

$\Delta \phi =\frac{1}{2\sqrt{N}}$

is made, wnere N is me number of photons measured For relatively modest parameters, this system can measure very sensitive relative changes in gravity. For example, using visible light photons (k≈107 m⁻¹), approximately 1016 photons/second flux (approximately 1 mW of cw laser power) in the interferometer, a relative pendulum height L of 2 cm, an interferometric sensitivity on the order of 10⁻¹³ m/s² or 10 pGal after one second of integration is made. Note that the presently-disclosed pendulum gravimeter differs from several previous pendulum gravimeters, which were based on a Fabry-Perot etalons.

A system according to the present embodiment can be used to measure modest masses at relatively large distances. The force F is simply the gravitational force due to a mass M at a distance d. Solving for d, yields:

$\begin{array}{cc}d=\sqrt{\frac{4\mathrm{GkLM}\sqrt{N}}{g}}.& \left(3\right)\end{array}$

Using the same numbers as above with a mass of M=1 kg, yields d≈25 m. In other words, the device is sensitive enough to measure one kilogram mass up to 25 meters away with a shot-noise limited signal-to-noise ratio of 1 after one second of integration.

Achieving the shot noise limited sensitivity for gravimeters is exceedingly difficult owing to significant environmental noise that can couple into the gravimeter. The present disclosure draws upon ideas in weak values where technical noise can be suppressed, in some cases, by several orders of magnitude and phase stabilization can be made to a high degree of precision. A displaced Sagnac interferometer is operated in the inverse weak value domain. For the inverse weak value, k_tσ<1, where k_t is the transverse momentum shift and σ is the beam width. The operation of the interferometer in this configuration results in an average beam displacement given by ϕ/k_t, which can make the displacement a very sensitive function of phase for small k_t. However, this large shift does not change the ultimate sensitivity, because it reduces the signal strength at the detector by (k_tσ cos(ϕ/2))².

The displaced Sagnac interferometer allows the beams to sample two different heights of the reflecting mirror while still achieving Sagnac stability and using low coherence sources. The parity of the interferometer is determined by the number of mirrors, in addition to the beam splitter, inside the interferometer. The Sagnac configuration cancels translational vibration noise since both the upper and lower beams are translated in the same way, which preserves the phase stability. Further, noise sources outside the interferometer such as beam jitter or detector jitter are not amplified. The Sagnac technique was found to be advantageous over the Michelson technique, having improved stability by several orders of magnitude.

Unfortunately, weak values cannot suppress all types of noise. If the noise is in the device that is generating the signal, there is no attenuation of the noise. This means that environmental noise (e.g., acoustic, microseismic, laser heating, laser shot-noise torques) may couple into the pendulum and limit the resolution of the detector. Measures may be used to lower such noise and improve the performance of the device. For example, noise may be reduced by using active isolation, having the interferometer in high vacuum, increasing the effective mass of the device against acoustic noise, balancing the readout laser with respect to the oscillation axis, etc.

In the experimental embodiment shown in FIG. 1, light from a fiber-coupled, free-running, diode laser was sent into a 50/50 beamsplitter. The light entering the interferometer was tilted upward, so that the exit beam from the interferometer was significantly higher than the input beam. The interferometer mirrors were all 2.5 cm in diameter. The displaced Sagnac was set such that the height differential of the clockwise and counter-clockwise beams was approximately one cm on the swinging pendulum. The swinging mass of the pendulum was a rectangular copper block of mass 500 g. A long Neodymium magnet (longer than the length of the block) was placed below the copper block to dampen the vibrations via Eddy currents. The height of the block relative to the magnet was tuned until critical damping was achieved. The 1/e damping time was approximately 50 ms. The copper mass was suspended by a light aluminum rod fastened to a razor blade with approximately 2.5 cm pendulum length. The razor blade edges touch two flat metal platforms. The mirror was an aluminum coating deposited on an aluminum substrate to ensure that all parts of the pendulum were electrically connected and grounded. The other mirror(s) in the interferometer were used for alignment of the interferometer. One of the mirrors in the interferometer had a multi-axis piezo-actuator to ensure precise alignment and adjustment of the dark fringe and bi-modal distribution. The light leaving the interferometer was then measured on a quad cell which differenced the intensity on one side of the detector from the other. The table was grounded and the copper block was electrically connected through the all-metal pendulum. The pendulum, but not the remainder of the interferometer, was put into a low vacuum acrylic container so that the light could freely propagate around the interferometer through the vacuum walls. The apparatus was on top of a breadboard on top of sorbethane pads, which rested on a large concrete slab embedded in a sand base and decoupled from the building. The system had a broad (10 Hz) microseismic spectrum centered at 8 Hz and a narrow 30 Hz signal, but a low background below 2 Hz. An active, vibration isolation table may be used.

To readout the pendulum shifts, the interferometer was aligned very close to the perfect dark port until a bimodal interference pattern with a large dark fringe results. One axis of the piezo actuated mount was then used to place the dark fringe in the middle of the intensity pattern. The lobes of the bimodal intensity pattern were aligned to two halves of a quad detector so that their difference signal could be measured. The rapidly moving fringe then resulted in large changes in the amount of light striking each side of the quad detector. The sensitivity of the device was calibrated using two different techniques. The first technique used readouts of the sum and difference voltages on the quad cell along with the measurement of the beam waist to determine k_t and σ, respectively. The second technique was based on the careful characterization of the other piezo-actuated mirror in the interferometer. By applying a calibrated voltage on the other mirror and measuring the balanced signal, the mass-equivalent shift could be calibrated. Then the sensitivity of the device was calibrated with a known mass at a given distance while monitoring the signal voltage. This allowed for calibration of the signal response for the given test mass for that alignment. While the ultimate sensitivity of the inverse weak value and the standard homodyne technique resulted in the same shot noise limit, the inverse weak value was useful, because the small phase shifts that result from very small pendulum tilts can lead to large dark fringe shifts, which are easily determined from the detector.

Preliminary measurements of the exemplary device were performed using a lock-in amplifier. A 1.5 kg mass was placed on a 3 rpm turntable. The center of the turntable was 46 cm from the center of mass of the pendulum. The 1.5 kg test mass was placed 17 cm from the center of the turntable. The experiments were performed early in the morning when the acoustic and seismic noise was at a minimum. The lock-in amplifier was able to find the signal and hold the voltage to 1 part in 10 of the central voltage. The lock-in voltage was close to the simulated value for the given experimental parameters. This would imply that gravitational attraction was measured to a few nanoGal after a few 10's of minutes integration. Further tests are expected to demonstrate linearity in mass and the quadratic dependence on distance. The use of further damping, vacuum, and vibration isolation is expected approach the shot-noise limited sensitivity without the need for the lock-in amplifier.

Torsion Pendulum

In a second embodiment, shown schematically in FIGS. 2 and 3 a gravimeter is configured as a torsion pendulum. Deflection of the pendulum is read using a displaced Sagnac interferometer as described above, except that the displacement L of the readout beams is in the horizontal plane (rather than from beams at different heights as in the above-described simple pendulum). An advantageous feature of this type of balance is that the sensitivity of the device can be tuned by modifying a restoring torque of the system, such as, for example, a restoring torque provided by a wire suspending the two masses m joined by a connecting rod.

The torsion pendulum will allow ultra-precise near-field gravimetric measurements and far-field gravity gradiometry for high resolution imaging.

III. Geometry and Torque Calculation

Suppose there is a system comprising a point mass M that is detected via a torsion pendulum comprising two masses m, connected via a rigid massless rod of length 2. The mass M is located according to FIG. 3 in relation to the oriented torsion pendulum. Assume M is a point mass for the time being. Suppose that the motion of all those objects will be characterized by their moment of inertia I about the axis defined by the pivot. By way of example, that axis could be comprised of a wire attached to the rigid body made up of the masses m at both ends. An external torque exerted from the gravitational force of a massive object will disturb the equilibrium position of the oscillator, which will oscillate until its damps to the new equilibrium position, as described in the next subsection.

A. Converting Torque Into Angle

In order to detect this small torque, recall that the pendula have a linear restoring torque quantified by the torsion spring constant κ, such that:

τ_ext=−κδθ,   (4)

where δθ denotes the angular distance from its equilibrium position. Note that κ can be empirically found by finding the period of the oscillations. When the pendulum swings freely,

τ_ext=Iα=2m²{umlaut over (θ)}=−κδθ,   (5)

where I=2 m² is the moment of inertia. Putting this equation in the form {umlaut over (θ)}+ω₀²θ=0, the natural frequency of the pendulum is found to be ω₀=(1/)√{square root over (κ/2m))}, which may be inverted to find κ from a measurement of the frequency, or period of the pendulum ,

$\begin{array}{cc}\kappa =\frac{8{\pi }^2m{}^2}{2^{}}.& \left(6\right)\end{array}$

Adding in damping of the pendulum brings the dynamics into the form of Eq. (7).

For a vertical simple pendulum subject to Earth's gravitational field, the restoring force is simply the gravitational acceleration from the earth, causing a restoring torque of τ=−gmδθ for small angles. This also gives rise to dynamics of the form Eq. (7) but with a natural frequency given by ω₀=√{square root over (g/)}.

The source of the gravitational signal is a mass M near the pendulum, thereby applying an external torque τ. We wish to measure this signal. We assume that τ can be time-dependent in general. Damping of the oscillations will be critical for a quickly responding detector, so we also add a velocity-dependent damping term, to find the equation of motion of a damped/driven oscillator,

{umlaut over (θ)}+2ζω₀{dot over (θ)}+ω₀²θ=τ(t)/I   (7)

Here, ω₀=√{square root over (κ/I)}, and ζ is a dimensionless damping coefficient. The underdamped case corresponds to 0≤ζ≤1, whereas the overdamped case corresponds to ζ>1. The general solution for τ=0 (the homogeneous solution) can be expressed as

$\begin{array}{cc}{\theta }_{\mathrm{hom}\left(t\right)}=e^{-\zeta {\omega }_0t}\{{\theta }_0\left[e^{{\omega }_0t\sqrt{{\zeta }^2-1}}-\sinh \left({\omega }_0t\sqrt{{\zeta }^2-1}\right)\right]+\hspace{1em}\left[\frac{\zeta {\theta }_0+{\stackrel{.}{\theta }}_0\text{/}{\omega }_0}{\sqrt{{\zeta }^2-1}}\right]\sinh \left({\omega }_0t\sqrt{{\zeta }^2-1}\right)\},& \left(8\right)\end{array}$

where θ₀ and {dot over (θ)}₀ are initial conditions for the pendulum, and the terms inside the braces { . . . } describe decay for the overdamped case, and oscillations in the underdamped case.

If τ is fixed in time, for large damping ζ, then the oscillator will converge to its new equilibrium position exponentially in time with a rate ζω₀ according to the solution (8). After this time, the angular displacement can be approximated by the fixed point of (7), given by

$\begin{array}{cc}\stackrel{\_}{\theta }=\frac{\tau }{I{\omega }_0^2}=\frac{\tau }{\kappa }.& \left(9\right)\end{array}$

We wish to design the pendulum to respond sensitively to stimuli from the target objects, but do not want it to oscillate for a long time before returning to a new equilibrium position. There is a trade-off between sensitivity of the measurement and the speed of the response as will be explored in the following sections.

B. Pendulum Model

The pendulum is fixed with respect to its center of mass motion, and is allowed to only rotate about its center of mass in the plane of the figure. This geometry can be analyzed by computing the torque about the middle of the torsion pendulum.

Plane trigonometry dictates that the distances defined in the figure are given by:

r ₁ ²=(x−)²+d², r₂²=(x+)²+d², r₀²x²d².   (10)

The gravitational force between mass M and mass m_j, according to Newton, is given by:

$\begin{array}{cc}{F}_j=-\frac{{\mathrm{GM}}_jM}{r_j^2}{\hat{r}}_j.& \left(11\right)\end{array}$

The torque τ generated on mass j is given by:

τ_j=F_j cos θ_j,   (12)

where F_j is the magnitude of gravitational force of mass M on mass m_j.

For a Simple Pendulum (SP), we only have one of the masses in the pendulum of FIG. 3. The net torque is:

$\begin{array}{cc}{\tau }_{\mathrm{SP}}=\mathrm{GmM}\frac{d}{r_1^3}.& \left(13\right)\end{array}$

The above results are relevant for a single test mass other than the pendulum mass. In the following sections we will use either the earth's gravity as the restoring force (as usual pendulums do), or by orienting the pendulum perpendicular to the earth's field, can also use the restoring force of a rod to obtain longer periods.

For a balanced Torsion Pendulum (TP), we include torques that nearly counterbalance each other (the torque on mass 1 is positive in sign, and the torque on mass 2 is negative in sign). The net torque is given by:

$\begin{array}{cc}{\tau }_{\Sigma }=\sum _j{\tau }_j={\tau }_{\mathrm{TP}}=d\mathrm{GmM}\left(\frac{1}{r_1^3}-\frac{1}{r_2^3}\right),& \left(14\right)\end{array}$

where we have replaced cos θ_j=d/r_j.

The simple pendulum responds to the bare force on the sensing mass, and thus decays as 1/r² with respect to the test mass distance. The torsion pendulum balances the average force, and thus responds to the gradient of the field across the size of the torsion pendulum. This effect leads to a less sensitive response to objects far away; it may be beneficial since it efficiently screens out far away objects and allows the sensor to focus on nearby objects.

C. Limiting Case

In some experiments, this expression (Eq. (14)) can be further simplified because it is expected that «d, x, r₀. The expression for τ_Σ is proportional to the difference of the functions g()−g(−↑), where:

$\begin{array}{cc}g\left(\right)=\frac{1}{{\left(d^2+{\left(x-\right)}^2\right)}^{3/2}}.& \left(15\right)\end{array}$

Since is a small parameter, one can approximate g()−g(−)≈g′(0)(2). It is found that:

$\begin{array}{cc}{g}^{\prime }\left(0\right)=\frac{3x}{r_0^5},& \left(16\right)\end{array}$

so that, to a good approximation,

$\begin{array}{cc}{\tau }_{\Sigma }=\frac{6\mathrm{GmMd}{}^2x}{r_0^5}=\frac{6\mathrm{GmM}{}^2\cos \theta \sin \theta }{r_0^3},& \left(17\right)\end{array}$

where we approximate θ₁≈θ₂=θ, and write x=r₀ sin θ and d=r₀ cos θ. In this limit, the sensor does not respond to the net force, but rather to its gradient, as indicated by the r⁻³ law.

In this limit, the one-armed device (SP) equilibrates to an angle

$\begin{array}{cc}{\stackrel{\_}{\theta }}_{\mathrm{SP}}\approx \frac{\mathrm{GM}\cos \theta }{{\omega }_0^2r_0^2},& \left(18\right)\end{array}$

while its two-armed counterpart equilibrates to

$\begin{array}{cc}{\stackrel{\_}{\theta }}_{\Sigma }\approx \frac{3\mathrm{GM}\sin \theta \cos \theta }{{\omega }_0^2r_0^3},& \left(19\right)\end{array}$

In both cases we have used Eq. (9). We stress that in both cases, the sensing mass m only appears in the natural frequency, and in the case of the torsion pendulum, the length also drops out, indicating that small sensors work as well as large ones so long as their periods are the same.

These expressions can be applied to make an approximate survey of the sensitivity of the device to different objects. Specifically, we show the best-case angular response to a target M at distance r₀ in FIGS. 4A and 4B, and we plot the angular dependence of the sensing for each device in FIGS. 5A and 5B. Some example values for a small torsion pendulum are given in Table I, in reference to the geometry of FIG. 3.

TABLE I
Example parameter values
Pendulum mass	m	100	g
Wavelength of light	λ	500	nm
Pendulum length		5	cm
Period of oscillator		500	s
Torsion spring constant	κ	7.9 × 10−8 kg m2/s2
Length between beams on the mirror	L	1	cm

The above analysis may be extended to a continuous mass distribution by replacing the mass M by a differential element dM=ρ(x)dx, where a body with mass per unit distance ρ(x) distributed along the x direction is imagined. In that case, the next torque for such a mass distribution is given by:

τ_Σ=∫f(x)ρ(x)dx.   (20)

In the general case, ƒ(x)=dGm(r₁(x)⁻³−r₂(x)⁻³), whereas in the limiting case, it is given by ƒ=6²dGmx/(d²+x²)^5/2.

IV. Noise Considerations

It is noted that the pendulum will experience several kinds of noise that may be mitigated in order to reach the fundamental limits of angle detection that the system is capable of. This section focuses on three types of noise: thermal noise, measurement heating noise, and quantum noise of the oscillator.

A. Thermal Noise

Contributions of the thermal noise from the surrounding environment can be computed via the equipartition theorem assuming large temperatures. Both the mean kinetic energy and potential energy are given by the thermal energy for one degree of freedom each. In general,

$\begin{array}{cc}\frac{1}{2}\kappa 〈\delta {\theta }^2〉=\frac{\hslash \omega }{4}\coth \left(\hslash \omega \text{/}2k_BT\right).& \left(21\right)\end{array}$

In the limit of high temperatures, the equipartition of potential energy indicates that

1/2κδθ²=1/2k_BT,   (22)

which gives the typical rms noise of the torsion pendulum,

δθ_rms=√{square root over (k_BT/τ.)}  (23)

The value can be estimated using the values given in Table I, and room temperature k_BT=4.1×10⁻²¹ J to find, δθ_rms=2.3×10⁻⁷ rad. In order to access below the picoradian regime, it will therefore be necessary to either cool the oscillator, or to time-average the signal for some time. One could also increase the value of κ, but that would also decrease the angular precision.

B. Measurement Heating

As demonstrated in the first section, it is advantageous to apply as much optical power as possible to the interferometer to maximize the precision of the angle measurement. However, because the gravitational sensor is freely moving, it is possible that the sensing laser may drive excitations, effectively heating the torsion pendulum. We will now calculate the effect of this heat, which may put a bound on the sensing power.

The displaced Sagnac geometry (see FIG. 2), causes two laser beams to strike the sensing mirror at a lever arm

$\pm \frac{L}{2}$

from the axis of rotation. We can imagine that photons each carrying energy hc/λ=k₀c arrive at the mirror at discrete times, either at

$+\frac{L}{2}\mathrm{or}-\frac{L}{2}.$

On average, since the intensity of the light on the left and right side of the mirror is the same from the 50/50 beamsplitter, we have

$〈N_{+}〉=〈N_{-}〉=\frac{N}{2}.$

The time-averaged optical power being applied to the pendulum is P=Nk₀c/T, which could be understood as arising from

$\begin{array}{cc}\stackrel{\_}{P}=\frac{1}{T}{\int }_0^T\hslash k_0c\left(\sum _n^{N_{\pm }}\delta \left(t-t_n\right)\right)\mathrm{dt},& \left(24\right)\end{array}$

provided that the duration T is much longer than the typical time between individual photons, which arrive at the t_n, and where P=P₊+P₋. Each photon which strikes the mirror imparts an impulse Δp=k₀/√{square root over (2)}, where the factor √{square root over (2)} arises because the photons strike the mirror at an angle. The time-averaged torque may be expressed using the same logic as above, i.e.,

$\begin{array}{cc}{\stackrel{\_}{\tau }}_{\pm }=\frac{1}{T}{\int }_0^T\frac{L}{2}\Delta p\left(\sum _n^{N_{\pm }}\delta \left(t-t_n\right)\right)\mathrm{dt},& \left(25\right)\end{array}$

or for T»|t_n−t_n−1|,

$\begin{array}{cc}〈\stackrel{\_}{\tau }〉=〈{\stackrel{\_}{\tau }}_{+}-{\stackrel{\_}{\tau }}_{-}〉=\frac{{\mathrm{Lk}}_0\hslash }{T\sqrt{2}}〈N_{+}-N_{-}〉.& \left(26\right)\end{array}$

Since the number of photons striking on either side of the point of rotation is equal on average, we see that τ=0; however, there will be fluctuations from the coherent states of light, which depend on

(N₊−N₋)²=N₊²−N₋²−2N₊N₋=N.   (27)

The term N₊N₋=0 because the ± photon beams are uncorrelated. Furthermore, given the geometry and using statistical properties of coherent states,

$〈N_{+}^2〉=〈N_{+}〉=\frac{N}{2}\mathrm{and}〈N_{-}^2〉=〈N_{-}〉=\frac{N}{2}.$

Therefore the torque fluctuations averaged over a given time interval T are given by

$\begin{array}{cc}〈{\stackrel{\_}{\tau }}^2〉=\frac{L^2k_0^2{\hslash }^2}{2T^2}〈{\left(N_{+}-N_{-}\right)}^2〉=\frac{L^2k_0^2{\hslash }^2}{2T^2}〈N〉=\frac{L^2{\stackrel{\_}{P}}^2}{2c^2〈N〉},& \left(28\right)\end{array}$

where we have re-written the duration of the experiment in terms of the power P, and used the relationships above. This leads to an increased variance of the angle, defined by δθ²=τ²/κ², or taking the square root, fluctuations in the coherent states of light lead to

$\begin{array}{cc}\delta \stackrel{\_}{\theta }=\frac{L\stackrel{\_}{P}}{\kappa c\sqrt{2N}}.& \left(29\right)\end{array}$

Equating with the shot noise limited value (2), we find that the measurement heating makes a smaller contribution to the noise than shot noise, when the condition

$\begin{array}{cc}\frac{L^2\stackrel{\_}{P}k_0}{\kappa }\le \frac{c}{2}& \left(30\right)\end{array}$

is satisfied. This implies that for the values in Table I, the noise contribution due to measurement heating would overtake the shot noise when the measurement laser power exceeds P=9.3 mW.C. uantum Noise

An intriguing aspect of the oscillators is the fundamental limitation of sensitivity due to quantum noise. As will be shown, the shot noise limited resolution is approximately equal to the quantum ground state uncertainty of the oscillator when the integration time is approximately equal to the period of the pendulum. From a quantum mechanical perspective, ground state quantum noise limitations are quite interesting in light of the large masses used in these experiment. Such studies may be valuable in probing quantum gravity. On the other hand, this also places fundamental noise limits on the resolution.

To determine the ground state angular uncertainty, we set the mean potential energy of the oscillator to the ground state energy of the oscillator

(1/2)κ(δθ)²=(1/4)ω,   (31)

which follows from Eq. (21) when k_bT«ω. Solving for the angle we obtain

δθ_rms=√{square root over (ω/2κ.)}  (32)

Using values listed in Table I, we find δθ_rms=2.9×10⁻¹⁵ rad. This resolution can be achieved when the integration time is roughly equal to the period of the pendulum using a few milliwatts of laser power.

V. Limits of Resolution

The results of the previous sections can now be combined to give the sensitivity limits of the SP and TP to forces, which can be translated into either mass or range uncertainty. Using Eq. (9) and the angular uncertainty Eq. (2), we find at the optimally sensitive response point (x=0 so d≈r, θ=0), the resolution of (usual) simple pendulum acceleration a, relative to the gravitational acceleration g of the SP to be

$\begin{array}{cc}\frac{\delta a}{g}=\mathrm{\delta \theta }& \left(33\right)\end{array}$

Consequently, the acceleration uncertainty in units of the accelerations due to gravity near the surface of the earth is simply the same as the angular uncertainty.

If instead, we consider a one-armed torsion pendulum with a torsion constant of κ, oriented perpendicular to the gravitational field of the earth, then the period of the oscillation can be much longer. The angular resolution is given by Eq. (18) so the acceleration uncertainty is reduced to

$\begin{array}{cc}\frac{\delta a}{g}=\frac{\mathrm{\kappa \delta \theta }}{gm}.& \left(34\right)\end{array}$

For the parameters in Table I, this reduces the acceleration uncertainty by a factor of 1.6×10⁻⁶, leading to 60 zepto g Hz^−1/2. Remarkably, the speed is only a thousand times slower, because of the inverse square relationship of Eq. (6). In either geometry, the acceleration is given by a=GM/r₀², so the sensitivity of the acceleration to a change in test mass δM at fixed r₀, or a change in the distance r₀ for fixed test mass M is given by

δa=GδM/r₀²−2GMδr₀/r₀³,   (35)

from which the mass or distance uncertainty is easily found. The response of the one-armed torsion pendulum is plotted in FIG. 4A for different values of test mass M and range R=r₀.

For a balanced torsion pendulum, a test mass far from the pendulum will respond according to Eqs. (16, 18). Setting θ=π/4 for maximum sensitivity, the angular response to a gravitating body will be

$\begin{array}{cc}\mathrm{\delta \theta }=\frac{3\mathrm{GmM}{}^2}{{\kappa r}_0^3}.& \left(36\right)\end{array}$

The r₀⁻³ law gives a smaller sensitivity, but also screens off distant objects. This cannot be directly translated into acceleration of a single mass, but gives the response of the detector to the gradient of the gravitational field. The torsion pendulum response is plotted in FIG. 4B for different values of test mass M and range R=r₀.

We now briefly discuss the angular response of both types of pendula as the test mass is placed at different angles relative to the axis of rotation. The one-armed torsion pendulum has blind spots at θ≈0 and θ≈π, where a target mass applies no torque, and its sensitivity is maximized at θ≈π/2 and θ≈3π/2. The two-armed torsion pendulum has four blind-spots, as illustrated in FIG. 5B. Notice that the scaling of the deflection angle in terms of the one-armed pendulum's construction parameters really depends only on , and that the the smaller we make , the larger the deflection angle will get (the moment of inertia in the denominator wins out over the greater torque with greater arm length). The dependence cancels out entirely from the two-armed device, except for its appearance in the natural frequency.

Near Field Measurements

To illustrate the sensitivity of the device, assume that the mass of interest is in the near-field of the pendulum (i.e., a distance of the mass from one pendulum mass is much closer than the separation between the two pendulum masses), the minimum detectable acceleration of a single mass close to the test mass is:

$\begin{array}{cc}a=\frac{4{\pi }^2l}{{\mathrm{kLT}}^2\sqrt{N}},& \left(37\right)\end{array}$

where l is the radius of the pendulum masses, k is the wavenumber, L is the distance between the beams reading out the deflection of the pendulum, N is the number of photons measured during an integrated time interval and T is the period of the pendulum. For a critically damped pendulum, the damping time is of the order of T/2. For r=0.2 m, L=0.1 m, T=1000 s and a 3 mW beam, the shot noise limit for acceleration is 10⁻²² m/s². To demonstrate the force sensitivity of the device, such a gravimeter is sufficient to measure the effective mass (m=E/c²) of 1 Watt (m=10¹⁷ kg) of energy at a few mm away from the pendulum, assuming all other effects due to the energy could be ameliorated (for example, one can imagine a resistor dissipating 1 W close to the pendulum with an ultra-broadband mirror reflecting the radiation away from the gravimeter in between). Using a pendulum mass of 100 g, the pendulum could measure forces with a resolution of 10 yocto Newtons, which could be used in almost any force measurement from protein folding to the Watt Balance. Interestingly, a ns, visible, single-photon normally reflecting from a mirror, imparts 10⁻¹⁹ N. There may be many limitations on the resolution of the device including the fact that shot noise in the read-out laser will cause torques on the mirrors.

GravCam, Solving the Inverse Problem

More than one of the devices can be arrayed. For example, with four such devices, a three-dimensional determination of the position and mass of an object can be determined assuming a large enough signal-to-noise ratio using the correlation of the independent detectors. Movement of the object may also be determined. Because the damping times are still relatively long, the objects must move relatively slowly, but with a knowledge of the velocity of an object (velocity can be determined from temporal responses of the pendula), additional structure information can be gathered. Further, with a large array of sensors allowed to pivot about their tilt-axis or translated, even static environments can be gravitationally imaged.

The sensitivity of embodiments of the disclosed device is expected to be sufficient to gravitationally map density variations within the body. As a simple proof-of-principle simulation, a uniform smaller sphere of known density, but unknown size and position inside a larger sphere of known radius and uniform density was considered. For simplicity, the position along a single axis was estimated. With two unknowns, size and single-axis position, only two gravity measurements were needed. Two gravity detectors separated by a distance of one meter with the same z=0 cm, y=0 cm coordinate were used. The center of each sphere was assumed to have a z=30 cm and y=0 m, but unknown x for the smaller sphere, albeit within the larger sphere. For the simulation, the density of the larger sphere was assumed to be ρ_L=1000 kg/m³ (the density of water) and only a slightly higher density ρ_S=1100 kg/m³ for the smaller sphere. The radius of the larger sphere was 10 cm and the center of the large sphere was half way between the two detectors. A precision of 10 pGal was assumed and using a smaller sphere radius of 3 cm. Solving the inverse problem, the simulation was found to predict the position of the smaller sphere to within 0.5 cm along the x axis and the simulation determined the mass to within 1%. The density problem is easier in Airport scanners where hidden objects typically have much higher density contrasts. In some embodiments, the detectors are moved such that even static environments can be imaged. Several pendula may be used in a gravity camera along with the torsion pendulum discussed later.

In another example, a torsion pendulum may be used to obtain high-resolution images, for example, when using the gravitational effects of both masses of the pendulum (see FIGS. 2 and 3). For simplicity, assume that the object of interest is close to centered between the masses (a small distance x from the symmetry line), but is a distance d=l from the rotation axis. In other words, along the z-axis it is a distance r and along the x-axis it is a distance x. In this small distance limit, the net torque on the pendulum is given by

$\sum \tau =\frac{2\mathrm{GMmx}}{r^2}.$

Setting this equal to the restoring torque and taking the differential, an advantageous spatial resolution is found to be:

$\begin{array}{cc}\Delta x=\frac{2{\pi }^2l^4}{{\mathrm{kLGMT}}^2\sqrt{N}}.& \left(38\right)\end{array}$

Using the numbers from above, along with a mass of 1 g, a 2 picometer spatial resolution is achievable at 20 cm away. This would indicate that the gravimetric gradiometer may be used for performing high resolution imaging well beyond the diffraction limit deep inside a bulk object.

Suppose the torsion pendulum is used as a sensor to detect an unknown mass density profile ρ(x). Suppose the mass distribution is moving at a constant speed v along the x direction. This will give a time-dependent torque τ_Σ(t) (proportional to the equilibrium angle δθ₀(t)). Then, given this time dependent signal, the mass distribution ρ(x) can be found.

The time-dependent torque is given by:

τ_Σ(t)=∫_−∞^∞f(x)ρ(x−vt)dx.   (39)

Defining a coordinate y=vt, the central insight is to interpret this integral as a convolution (ƒ*ρ)(y) of two functions ƒ and ρ (where * denotes convolution). By the convolution theorem, the unknown function ρ(x) is given by:

$\begin{array}{cc}\rho \left(x\right)={\int }_{-\infty }^{\infty }\frac{\mathrm{dk}}{2\pi }e^{\mathrm{ikx}}\left({\tilde{\tau }}_{\Sigma }\left(k\right)/\tilde{f}\left(k\right)\right).& \left(40\right)\end{array}$

where {tilde over (ƒ)}(k) and {tilde over (τ)}_Σ(k) are defined as the Fourier transforms of the functions ƒ(x), which is a known function, and τ_Σ(y), which is the measured data. It is noted that both the Fourier transform of the measurement data and the inverse Fourier transform to find the solution can be made with numerical fast Fourier transforms. Consequently, this procedure can be done very quickly, once the functional form of Fourier transform of ƒ(x) is known. The analytic form of the function ƒ is calculated below, giving the final step in the solution of the problem.

Fourier Transform of ƒ(x)

The Fourier transform of ƒ (in the general case) is proportional to the difference of the Fourier transform of g(x,)−g(x,−), where:

$\begin{array}{cc}\tilde{g}\left(k,\right)={\int }_{-\infty }^{\infty }\mathrm{dx}\frac{e^{-\mathrm{ikx}}}{{\left[{\left(x-\right)}^2+d^2\right]}^{3/2}}.& \left(41\right)\end{array}$

The integral is nontrivial because of the branch cut. The result is found as:

$\begin{array}{cc}\tilde{g}\left(k,\right)=e^{-\mathrm{ikl}}\frac{2k}{d}K_1\left(dk\right),& \left(42\right)\end{array}$

where K₁(z) is the first modified Bessel function. This function vanishes at k=0 since K₁(0) is finite, and decays exponentially at large k. Taking the difference of {tilde over (g)}(k,) at and −, as well as restoring the prefactors, the following result is found,

{tilde over (ƒ)}(k)=(−4i)Gm sin(k)|k|K₁(|k|d).   (43)

Together with the main result Eq. (40), the problem is solved. It is noted that the velocity v is unimportant, so long as it is slow enough to allow the torsion pendulum to come to equilibrium, and the optics to readout the angle to suffciently high precision.

An example of a 1-D mass distribution is shown in FIG. 6, which shows two step square functions of differing spatial extent and amplitude. FIG. 7 shows the predicted static equilibrium of the pendulum as a function of the position of the mass relative to the pendulum. The net displacement signal is shown with 10% Gaussian noise. Using a Wiener filter and a total variation reconstruction, FIG. 8 shows the reconstructed image.

Generalization to Two and Three Dimensions

Suppose that a 2D or 3D object is to be gravitationally image. There will be a mass distribution ρ(x,y) or ρ(x,y,z). One may proceed by assuming the torsion pendulum can rotate only about one axis (the z-axis), and consequently only the z-component of the torque is computed. This is found as:

τ_z={circumflex over (z)}·({right arrow over (r)}_i×{right arrow over (F)}_i),   (44)

where {right arrow over (r)}_i is a vector from the origin to mass m_i, and {right arrow over (F)}_i is the force vector acting on mass i from the differential mass region at point (x,y,z). The net z-component of the torque is found as,

$\begin{array}{cc}{\tau }_{z,\Sigma }=\mathrm{GmMy}\left(\frac{1}{{\left(r_1^{\prime }\right)}^3}-\frac{1}{{\left(r_2^{\prime }\right)}^3}\right),& \left(45\right)\end{array}$

where r′_1,2=√{square root over ((x∓)²+y²+z²)} has been modified. Here, the basic effect of the z direction is seen as contained in the mass distribution itself, and the additional fall off in the denominator the further the mass is from the pendulum. In the y direction, the constant value d now can change as y, as one gets closer and further away from the pendulum. The two dimensional case corresponds to fixing the y variable.

The continuum limit can be reached by introducing a mass density ρ, with a net torque of:

$\begin{array}{cc}{\tau }_{z,\Sigma }=\int \mathrm{dxdydz\rho }\left(x,y,z\right)f\left(x,y,z\right),\mathrm{where}\text{:}& \left(46\right)\\ f\left(x,y,z\right)=\mathrm{Gmy}\left(\frac{1}{{\left(r_1^{\prime }\right)}^3}-\frac{1}{{\left(r_2^{\prime }\right)}^3}\right).& \left(47\right)\end{array}$

In 2D, one only integrates over x and z. The inversion in 2D and 3D works in principle as the 1D case. Focusing here on the 2D case of imaging in a plane, the torque is only mapped out as the mass density is being raster scanned. This gives the torque as a function of the relative coordinates τ_z,Σ(x′,z′). As in the 1D case, this is the same as Eq. (46), with ρ replaced by ρ(x−x′,z−z′). This is now a two-dimensional convolution. Consequently, ρ(x,z) can be represented with the inverse Fourier transform ,

ρ(x,z)=[{tilde over (τ)}_z,Σ({right arrow over (k)})/{tilde over (ƒ)}({right arrow over (k)})],   (48)

where the known function {tilde over (ƒ)}=Gm({tilde over (g)}()−{tilde over (g)}(−)) is given by the multi-dimensional Fourier transform,

$\begin{array}{cc}\tilde{g}\left(\overrightarrow{k}\right)=\int d^2{\mathrm{xe}}^{-i\overrightarrow{k}·\overrightarrow{x}}\frac{y}{{\left({\left(x-l\right)}^2+y^2+z^2\right)}^{3/2}}.& \left(49\right)\end{array}$

Finite Sampling in 3D

In the previous sections, a continuous data curve was assumed as the imaged object was displaced in time. It is likely more realistic to consider a set of data points taken at discrete locations around the object in question. One may consider finite data sets from such an array of measurements at discrete points. For simplicity, consider a uniform grid of object points in three dimensions, {x_i,y_j,z_i}, where i=1, . . . , p. These points approximate a three dimensional object by shrinking the mass within a lattice unit cell into a point of mass m_i. This mass array is sensed with a detector sampling the gravitational field at a series of discrete points {x_j^d,y_j^d,z_j^d}, with j=1, . . . , q. These points can be taken around the outside of the body that is being imaged, with several detectors, or a single detector that is displaced (or kept stationary, with the body being imaged moved around). Since the gravitational force is proportional to the two attractive masses, and is additive as a vector, a general relation can be derived between the measured quantities θ_j and the system masses m_i, given as a linear relationship,

θ=R·M,   (50)

where θ is a vector of dimension q with elements θ_j, and M is a vector of dimension p with elements m_i. The matrix R is a q×p matrix, which encodes both the geometric relationships between detector and system, as well as the type of detector in use (a free particle, a simple pendulum, a torsion pendulum, etc.). For each of those detectors, the measured quantities may be different (acceleration, displacement, an angle, etc.)

If the dimensions of the system and detector degrees of freedom match, then the matrix has a unique inverse (assuming is has a nonzero determinant). Otherwise, the solution is either over or under determined in general. The Moore-Penrose pseudoinverse, R⁺ can be used to give the best fit solution to this linear equation,

M _bƒ =R ⁺·θ,   (51)

where the pseduo-inverse solution minimizes the norm of all solutions M. The pseduo-inverse may be constructed in a simple way with the singular value decomposition, R=UΣV^†, where U,V are unitary matricies, and Σ is a diagonal matrix,

R ⁺ =VΣ ⁺ U ^†.   (52)

The pseduo-inverse of the diagonal matrix Σ, denoted as Σ⁺, may be found by inverting the nonzero elements of Σ and taking the transpose.

Matrix R is now computed for the special case of a torsion pendulum detector discussed above. Adapting the geometry shown in FIG. 3 to three dimensions, and aligning the torsion pendulum in the y direction, one considers sampling the gravitational field in a plane of constant z. In that case, the net torque about the pendulum axis is

$\begin{array}{cc}{\tau }_{\mathrm{ij}}={\mathrm{Gm}}_im_j\left(\frac{z_j^d-z_i}{{\left[{\left(z_j^d-z_i\right)}^2+{\left(y_j^d-y_i\right)}^2+{\left(x_j^d-x_i-l\right)}^2\right]}^{3/2}}-\frac{z_j^d-z_i}{{\left[{\left(z_j^d-z_i\right)}^2+{\left(y_j^d-y_i\right)}^2+{\left(x_j^d-x_i-l\right)}^2\right]}^{3/2}}\right).& \left(53\right)\end{array}$

Assuming the masses of the torsion pendulum are the same at every point, and therefore set the masses of the pendulum m_j=m for all j. The net displacement angle of the pendulum located at position i, for a restoring torsion constant κ, is given by

$\begin{array}{cc}{\theta }_j=\sum _i{\tau }_{\mathrm{ij}}/\kappa =\sum _iR_{\mathrm{ji}}m_i,& \left(54\right)\end{array}$

where the last equality defines the matrix R_ji.

In order to switch to another imaging plane, the torque given in Eq. (53) can be modified accordingly. For example, to align the pendulum in the z direction, and sample in the plane of constant y, the numerator would involve y_j^d−y_i, and the denominator's dependence would be in the y degree of freedom. The matrix R′_ji obtained from that plane would then be appended to the matrix R_ji from the first measurement plane, to obtain a generally rectangular matrix of all the data.

Once the response matrix R is defined, it is then straightforward to change the mass distribution within the region, and to use the same reconstruction matrix, which will still be valid so long as the geometric relationships between the measurement planes and grid points are the same.

Numerical Examples

The problem to solve is a deconvolution problem: Given a known “point spread function,” ƒ({right arrow over (x)}), and a measured torque function versus scan coordinate {right arrow over (x)}′, the mass density profile ρ({right arrow over (x)}) is to be reconstructed where there is added white Gaussian noise v from the measurement process,

τ=ƒ*ρ+v   (55)

Here * indicates the convolution operation. The simple inverse of the Fourier transform (in the absence of noise) is problematic, since the Fourier transform of the point spread function rapidly goes to zero, so its inverse diverges rapidly, leading to instability in the deconvolution. There exist a suite of deconvolution techniques, the classic one being the Wiener filter, defined as:

$\begin{array}{cc}W\left(k\right)=\frac{{\left(ℱ\left[f\right]\left(k\right)\right)}^{*}}{{ℱ\left[f\right]\left(k\right)}^2+1/}.& \left(56\right)\end{array}$

Here is the signal-to-noise ratio of the data collection, which cuts off the divergence in inverting the small values of the transform of ƒ. The reconstructed approximation is given by:

92 _rc =τ*W   (57)

Note that because [ƒ](k=0), the overall scale of ρ cannot be directly reconstructed, so the inverse is defined only up to an overall factor; that is, Cδ(k) can be added to the Fourier transform, or a constant function to ƒ(x).

One having skill in the art will understand that the Wiener filter is the oldest example in what is now a lengthy number of possible methods to help in the deconvolution problem. The present disclosure should not be limited to the disclosed methods—e.g., the non-limiting, exemplary use of the Wiener filter.

To obtain a 2D scan, the device may be raster scanned or multiple devices may be used. In embodiments where a 3D scan is obtained, multiple detectors can be incorporated to work together.

Resolution of the Imaging System

With respect to the sensitivity of such a detector to fine-grain details of an image, in an optical context, this resolution is limited by the Rayleigh Criterion. However, in the gravitational case, there is no known “wavelength,” so in principle there is no limit to the resolution within Newtonian gravity.

In practice, however, there will be limit set by how precisely the detector can measure the angle shift of the torsion pendulum. Consider Δθ as the standard deviation of the measured angle. As can be seen from FIG. 9, the response of the torsion pendulum is most sensitive at x=0. The position resolution of the detector can be estimated by finding the standard deviation Δx corresponding to the angular uncertainty. Physically, this corresponds to the ability to resolve when a point mass is shifted by a small amount Δx in space.

The resolution Δx may be found by error propagation,

$\begin{array}{cc}\Delta x={\frac{\partial \theta \left(x\right)}{\partial x}}^{-1}\mathrm{\Delta \theta },& \left(58\right)\end{array}$

where one considers when the slope of θ is steepest with respect to x—i.e., x=0.

The equilibrium angle for a point mass M is given by θ(x,,y)=C(g(x,y,)−g (x,y,31)), where C=yGmM/κ. The slope of the function at the origin S₀ is computed,

$\begin{array}{cc}{S_0=\frac{\partial \theta }{\partial x}}_{x=0}=\frac{6\mathrm{C}}{{\left({}^2+y^2\right)}^{5/2}}=\frac{6\mathrm{GmM}}{\kappa }\frac{{\mathrm{yl}}^2}{{\left({}^2+y^2\right)}^{5/2}}.& \left(59\right)\end{array}$

In order to made the uncertainty on Δx as small as possible, the slope S₀ is to be as large as possible. Therefore, S₀ versus is maximized, keeping y fixed. The maximal slope is found to be when =√{square root over (2/3y)}≈0.8615y, or y≈1.22, which corresponds to

$\begin{array}{cc}{S}_{0,\max }\approx \frac{1.115\mathrm{GmM}}{\kappa y^2}& \left(60\right)\end{array}$

where again, y is assumed to be fixed.

Estimating the position uncertainty δx, m and M are assumed to be about 1 kg, and the precision with which the object can be located is considered. Using the same torsion constant as before, κ of 10⁻³ Nm/rad, let y=1 m and G≈6×10⁻¹¹ Nm²/(kg)₂. For the angular uncertainty, a value found in previous work was used, Δθ=400 frad=4×10⁻¹³ rad. Putting this all together, for a 1 kg mass at 1 meter,

Δx≈6 μm,   (61)

which shows very good results. If the test mass is reduced to a mass of 1 mg, the resolution decreases by three orders of magnitude to about 6 mm, which is still good.

Torsion Pendulum Rod Embodiment

Considering another geometry where there is just a rod pendulum of mass m and length L. A point mass is considered, located at perpendicular position y from the rod, and lateral position x from the center of mass of the rod.

Integrating the net torque produced by a mass element dm over the length of the rod, the torque is found as,

$\begin{array}{cc}{\tau }_{\Sigma }=\frac{\mathrm{GMmy}}{L}\left(\frac{\mathrm{Lx}/\left(2y^2\right)-x^2/y^2-1}{\sqrt{{\left(L/2-x\right)}^2+y^2}}+\frac{\mathrm{Lx}/2+x^2/y^2+1}{\sqrt{{\left(L/2+x\right)}^2+y^2}}\right)& \left(62\right)\end{array}$

The moment of inertia is I=(1/12)mL². It is noted that this moment of inertia of the bar is down by a factor of ⅓ compared to the dumbell (two masses interconnected by a rigid rod) of the same length and total mass (m→2 m, L→2). This gives the period T since τ_Σ=iα is of the pendulum. Here, the torsion constant κ is expressed, in terms of T, as:

$\begin{array}{cc}\kappa =\frac{{\pi }^2mL^2}{3T^2}.& \left(63\right)\end{array}$

This gives the angle θ=τ_Σ/κ.

The function Eq. (62), divided by κ, is the new “point spread function” for imaging a mass distribution with the rod imaging system.

Exemplary Methods

In an embodiment, a method of scanning an object includes determining a first displacement of a pendulum caused by the object. The displacement may be determined using a Sagnac inferometer as described herein. The object is moved with respect to the pendulum. For example, either of the object or the pendulum can be moved, or both can be moved. A second displacement of the pendulum (i.e., displacement of the pendulum caused by the object) is measured. A profile of the object is determined. The profile may include, for example, a position (such as, for example, a 1D position), a mass, and/or a mass distribution (e.g., gravimetric density profile). As described above, additional measurements of the pendulum displacement may be made while the object is in various locations.

In another embodiment, a method of scanning an object comprises providing a pendulum gravimeter. The pendulum gravimeter may be of any variety described herein. For example, in an embodiment, the pendulum gravimeter has a pendulum with a mirror. The pendulum is configured to have a restoring force, for restoring the pendulum to a neutral position. The pendulum gravimeter may further include an interferometer, such as, for example, a Sagnac interferometer. The interferometer is configured to measure a displacement of the pendulum from the neutral position. The method includes measuring the gravitational field of the object at a first measurement location and measuring the gravitational field of the object at a second measurement location. Measuring the gravitational field may include measuring a displacement of the pendulum from the neutral position, and determining the gravitational field based on the measured displacement. The mass distribution of the object is calculated based on the measured gravitational fields. The mass distribution of the object may be calculated in any way described herein.

For example, the object may be measured at a plurality of additional measurement locations for a total of q measurement locations. The object may be scanned in three dimensions and the measured gravitational field of the object θ is represented by a matrix of the measurements θ_j at each measurement location j having coordinates {x_j^d,y_j^d,z_j^d}, where j=1, . . . , q and q is the total number of data points. The mass distribution of the object M may be determined as a matrix of point masses m_i over a set of grid points {x_i,y_j,z_i}, where i=1, . . . , p and p is the total number of grid points. The mass distribution M of the object may be determined according to M=R⁺·θ, where R⁺ is the pseudoinverse of a predetermined scanner matrix R.

The calculated mass distribution may be used to create, for example, an image of the object. In such a case, scanning the object may be considered as imaging the object, and the gravimeter may be considered to be a “gravimetric camera.”

In some embodiments, the object is measured at the second measurement location using a second pendulum gravimeter. In some embodiments, the object is measured at the second measurement location by moving the pendulum gravimeter (e.g., moving the gravimeter and subsequently obtaining the measurement). In some embodiments, the object is measured at the second measurement location by moving the object with respect to the pendulum gravimeter (e.g., moving the object and obtaining the reading when the gravimeter is at the second measurement location with respect to the object).

Although the present disclosure has been described with respect to one or more particular embodiments, it will be understood that other embodiments of the present disclosure may be made without departing from the spirit and scope of the present disclosure.

Claims

1. A method of scanning an object, comprising: providing a pendulum gravimeter comprising: a pendulum having a mirror, the pendulum configured to have a restoring force for restoring the pendulum to a neutral position; andan interferometer configured to measure a displacement of the pendulum from the neutral position;measuring the gravitational field of the object at a first measurement location;measuring the gravitational field of the object at a second measurement location; andcalculating the mass distribution of the object based on the measured gravitational fields.

2. The method of claim 1, wherein measuring the gravitational field comprises: measuring a displacement of the pendulum from the neutral position; anddetermining the gravitational field based on the measured displacement.

3. The method of claim 1, wherein the object is measured at the second measurement location using a second pendulum gravimeter.

4. The method of claim 1, wherein the object is measured at the second measurement location by moving the pendulum gravimeter to the second measurement location.

5. The method of claim 1, wherein the object is measured at the second measurement location by moving the object with respect to the pendulum gravimeter such that the pendulum gravimeter is at the second measurement location.

6. The method of claim 1, wherein the object is scanned in one dimension (1D).

7. The method of claim 1, wherein the object is scanned in two dimensions (2D).

8. The method of claim 1, wherein the object is scanned in three dimensions (3D).

9. The method of claim 1, wherein the object is measured at a plurality of additional measurement locations for a total of q measurement locations.

10. The method of claim 9, wherein the object is scanned in three dimensions and the measured gravitational field of the object θ is represented by a matrix of the measurements θj at each measurement location j having coordinates {xjd,yjd,zjd}, where j=1, . . . , q and q is the total number of data points; and wherein the mass distribution of the object M is determined as a matrix of point masses mi over a set of grid points {xi,yj,zi}, where i=1, . . . , p and p is the total number of grid points.

11. The method of claim 10, where the mass distribution M of the object is determined according to M=R+·θ, where R+ is the pseudoinverse of a predetermined scanner matrix R.

12. A method of scanning an object, comprising: determining, using a Sagnac interferometer, a first displacement of a pendulum caused by the object;moving the object with respect to the pendulum;measuring a second displacement of the pendulum caused by the object; anddetermining a profile of the object based on the displacement of the pendulum.