diff --git a/crates/geopotential/src/coefficients.rs b/crates/geopotential/src/coefficients.rs index 9d85245..9a522b0 100644 --- a/crates/geopotential/src/coefficients.rs +++ b/crates/geopotential/src/coefficients.rs @@ -9,7 +9,11 @@ pub enum TideSystem { TideFree, /// Includes the permanent deformation (indirect) part, not the direct tidal potential. ZeroTide, - /// Includes both the direct and the indirect permanent parts. + /// Includes both the direct and the indirect permanent parts. The direct part, the tide-generating + /// potential itself, then sits in the Earth's own field, falling as `r⁻³` where it grows as `r²`: + /// right on the sphere `r = R` only ("there is no mean-tide C̄20", Mäkinen 2021, J. Geod. 95, + /// 106). Its gravity is +45.7 µGal at the equator against the true −30.5; take mean-tide gravity + /// from `gravimetry::tides::permanent`. MeanTide, /// Not stated; cannot be converted. Unknown, @@ -17,8 +21,9 @@ pub enum TideSystem { impl TideSystem { /// `C̄20` in this system minus `C̄20` tide-free. Uses the EGM96/EGM2008 convention - /// `ΔC̄20 = −3.11080e-8·k/√5` with `k = 0.3` (Rapp 1989), which gives the IERS 2010 value - /// −4.1736e-9 for zero-tide minus tide-free. + /// `ΔC̄20 = −3.11080e-8·k/√5` with `k = 0.3` (Rapp 1989), which gives −4.1736e-9 for zero-tide + /// minus tide-free, the value the IERS Conventions 2010 quote for EGM2008 (their own + /// `A₀H₀k₂₀`, Eq. 6.14 with `k₂₀ = 0.30190`, gives −4.2007e-9). fn c20_offset(self) -> Option { let direct = -3.11080e-8 / 5f64.sqrt(); let k = 0.3; diff --git a/crates/gravimetry/README.md b/crates/gravimetry/README.md index 55ed73e..3866134 100644 --- a/crates/gravimetry/README.md +++ b/crates/gravimetry/README.md @@ -133,6 +133,7 @@ maximum likelihood over many series at once. | `atom` | Phase, shot noise, and accelerometer and gradiometer noise of a cold-atom sensor, in orbit (double diffraction) or on the ground (single), its phase sensitivity and the gravitational Aharonov–Bohm rate ([Lévèque et al. 2021](https://arxiv.org/abs/2011.03382); [Stray et al. 2022](https://doi.org/10.1038/s41586-021-04315-3); [CARIOQA](https://arxiv.org/abs/2211.01215); [Trimeche et al. 2019](https://doi.org/10.1088/1361-6382/ab4548)) | | `gradiometry` | The error of each coefficient a gradiometer mission like GOCE achieves | | `local` | Gravity and gravity gradient of a buried sphere and horizontal cylinder, the shapes a ground survey fits, and the averaging time an instrument needs to see them | +| `tides` | The lunar and solar tide at a gravimeter, the gravimetric and diminishing factors of the Love numbers, and the permanent tide that separates mean-tide, zero-tide and tide-free gravity and geoids ([Longman 1959](https://doi.org/10.1029/JZ064i012p02351), as [LongmanTide](https://github.com/jrleeman/LongmanTide) codes it; [Mäkinen 2021](https://doi.org/10.1007/s00190-021-01541-5)) | | `sim` | A closed-loop simulation of a GRACE-type pair: simulate, estimate the field back, and compare its errors with the predicted ones | | `noise` | Noise with any spectrum, its whitening, and the FFT | | `technotes` | NASA's TN-14 (C̄20 and C̄30 from laser ranging, [Loomis et al. 2020](https://doi.org/10.1029/2019GL085488)) and TN-13 (degree 1, [Sun et al. 2016](https://doi.org/10.1002/2016JB013073)), from [PO.DAAC](https://archive.podaac.earthdata.nasa.gov/podaac-ops-cumulus-docs/gracefo/open/docs/) | @@ -161,6 +162,8 @@ errors of a real month.* files; `prompt_lead` gives how far ahead of the P-wave a light-speed gravity signal arrives. - **Site surveys and monitoring.** `local` gives the gravity and gradient of a buried void or dense body, and the averaging time a quantum gravimeter needs to resolve it. +- **Absolute gravimetry.** `tides` corrects a reading for the lunar and solar tide and moves it + between the zero-tide convention of gravity and the mean-tide convention of the IHRS. - **Mission design.** `missions`, `atom` and `sim` size the sensors and simulate what field a satellite pair recovers, including cold-atom instruments. @@ -179,6 +182,9 @@ errors of a real month.* equivalent mass sums to zero); on the 2025 Kamchatka fault model, its uplift matches the USGS's own surface displacement with r = 0.89. - **Point masses:** their potential matches Newton's law. +- **Tides:** LongmanTide's test vector to 3×10⁻¹¹ mGal; averaged over the 18.61-year nodal cycle, + Longman's tide is δ₂ times the permanent tide to 0.2 µGal; the permanent tide matches Mäkinen's + (2021) Eqs. 19–21 to their printed digits. - **Kalman smoother:** it matches the exact Gaussian-process solution to 1e-10, its fit recovers the settings of simulated data, and on the real record withheld months are predicted with calibrated errors. diff --git a/crates/gravimetry/examples/buried_tunnel.rs b/crates/gravimetry/examples/buried_tunnel.rs index 1c4bf6b..da5571e 100644 --- a/crates/gravimetry/examples/buried_tunnel.rs +++ b/crates/gravimetry/examples/buried_tunnel.rs @@ -1,30 +1,51 @@ //! Finding a void from the surface: the gravity and gradient of a buried tunnel, and how long a -//! quantum gravimeter must average to see it. The geometry is the Birmingham survey's -//! ([Stray et al. 2022](https://doi.org/10.1038/s41586-021-04315-3)), which found a tunnel at -//! about 1.9 m depth. +//! quantum gravimeter or gradiometer must average to see it. The geometry and the gradiometer are +//! the Birmingham survey's ([Stray et al. 2022](https://doi.org/10.1038/s41586-021-04315-3)), +//! which found a tunnel at about 1.9 m depth. //! //! ```text //! cargo run -p gravimetry --example buried_tunnel //! ``` +use gravimetry::atom::AtomInterferometer; use gravimetry::local::{EOTVOS, MICROGAL, averaging_time, cylinder}; fn main() { let (radius, depth, contrast) = (1.0, 1.9, -2000.0); // an air-filled tunnel in soil - println!("{:>8} {:>12} {:>12}", "offset/m", "Δg/µGal", "Γuu/E"); + + // Stray et al.'s gradiometer: two single-diffraction interferometers 1 m apart vertically, + // T = 85 ms, 10⁵ atoms at 0.7 Hz. It reads the difference in gravity between its sensors. + let stray = AtomInterferometer { t: 0.085, contrast: 1.0, atoms: 1e5, cycle: 1.0 / 0.7, diffraction: 1 }; + let baseline = 1.0; + let difference = |offset: f64| (cylinder(radius, depth, contrast, offset).0 - cylinder(radius, depth + baseline, contrast, offset).0) / baseline; + + println!("{:>8} {:>10} {:>10} {:>12} {:>12}", "offset/m", "Δg/µGal", "Γuu/E", "Δg/1 m in E", "phase/mrad"); for offset in [0.0, 1.0, 2.0, 3.0, 5.0] { let (gz, gradient) = cylinder(radius, depth, contrast, offset); - println!("{offset:>8.1} {:>12.1} {:>12.1}", gz / MICROGAL, gradient / EOTVOS); + let d = difference(offset); + println!("{offset:>8.1} {:>10.1} {:>10.1} {:>12.1} {:>12.2}", gz / MICROGAL, gradient / EOTVOS, d / EOTVOS, stray.phase(d * baseline) * 1e3); } let (peak, _) = cylinder(radius, depth, contrast, 0.0); println!("peak anomaly {:.0} µGal", peak / MICROGAL); // An Exail-class gravimeter at 500 nm/s²/√Hz needs this long per station for a 10:1 signal. let target = peak.abs() / 10.0; - println!("to resolve it 10:1 ({:.0} nm/s²): {:.0} s per station", target * 1e9, averaging_time(500e-9, target)); + println!("gravimeter, to resolve it 10:1 ({:.0} nm/s²): {:.0} s per station", target * 1e9, averaging_time(500e-9, target)); + + // The gradiometer at its projection limit, and at the 466 E/√Hz it reached outdoors. + let target = difference(0.0).abs() / 10.0; + for (what, asd) in [("projection limit", stray.gradient_asd(baseline)), ("outdoors", 466.0 * EOTVOS)] { + println!("gradiometer, {what} {:.0} E/√Hz: {:.0} s per station", asd / EOTVOS, averaging_time(asd, target)); + } + assert!((peak / MICROGAL + 44.1).abs() < 0.1); // The gradient changes sign away from the tunnel, which is what localises it. assert!(cylinder(radius, depth, contrast, 0.0).1 < 0.0 && cylinder(radius, depth, contrast, 3.0).1 > 0.0); + // Over a line mass the difference across the baseline is the point gradient at the geometric + // mean of the two depths, 2Gλ/(z₁z₂): about Stray et al.'s 150 E and 17.5 mrad. + let mid = cylinder(radius, (depth * (depth + baseline)).sqrt(), contrast, 0.0).1; + assert!((difference(0.0) / mid - 1.0).abs() < 1e-12); + assert!((difference(0.0) / EOTVOS + 152.2).abs() < 0.1 && (stray.phase(difference(0.0)) * 1e3 + 17.7).abs() < 0.1); } #[test] diff --git a/crates/gravimetry/src/lib.rs b/crates/gravimetry/src/lib.rs index 321a55f..967ec18 100644 --- a/crates/gravimetry/src/lib.rs +++ b/crates/gravimetry/src/lib.rs @@ -12,6 +12,7 @@ pub mod missions; pub mod noise; pub mod sim; pub mod technotes; +pub mod tides; /// Area of the world ocean, m²: 3.618×10⁸ km², so 361.8 Gt of water raise sea level 1 mm. pub const OCEAN_AREA: f64 = 3.618e14; diff --git a/crates/gravimetry/src/tides.rs b/crates/gravimetry/src/tides.rs new file mode 100644 index 0000000..08e7110 --- /dev/null +++ b/crates/gravimetry/src/tides.rs @@ -0,0 +1,224 @@ +//! The Earth tide at a gravimeter: Longman's lunar and solar tide, the Love-number factors that +//! scale a rigid Earth's tide to the observed one, and the permanent tide that separates the +//! mean-tide, zero-tide and tide-free conventions. +//! +//! ``` +//! use geopotential::julian_epoch; +//! use gravimetry::tides::{delta, longman, permanent}; +//! +//! // State College, Pennsylvania, at 2015-04-23 00:00 UT: 3.5 µGal of lunar and solar tide +//! // (LongmanTide's test vector), with Longman's own Love numbers. +//! let (lat, lon) = (40.7914f64.to_radians(), 282.1414f64.to_radians()); +//! let (moon, sun, total) = longman(lat, lon, 370.0, julian_epoch(2015, 4, 23, 0, 0)); +//! assert!((total / 1e-8 - 3.5347).abs() < 1e-4 && (moon + sun - total).abs() < 1e-18); +//! assert!((delta(2, 0.612, 0.303) - 1.1575).abs() < 1e-15); +//! // The permanent tide lowers gravity 30.49 µGal at the equator (Mäkinen 2021, Eq. 20). +//! assert!((permanent(0.0, 0.0).1 / 1e-8 + 30.49).abs() < 0.01); +//! ``` + +use core::f64::consts::TAU; + +use geopotential::geocentric; +use metrology::WGS84_A; + +/// The time average of the tide-generating potential in Mäkinen's (2021, Eq. 18) form +/// `W_T = A(r/a)²(sin²ψ − 1/3)`, m²/s²: −2.9166, from the IERS Conventions 2010, which Mäkinen +/// proposes for the International Height Reference Frame. Spectral expansions from Doodson's to +/// Kudryavtsev's agree to 0.0002 (his Table 1). +pub const PERMANENT_A: f64 = -2.9166; + +/// The gravimetric factor `δₙ = 1 + (2/n)hₙ − ((n+1)/n)kₙ` of degree `n` for Love numbers `h` +/// and `k`: the observed tide in gravity over a rigid Earth's. It sums the direct pull `nW/r` of +/// a tide potential `W ∝ rⁿ`, the free-air change `2hW/r` of the surface lifted `hW/g`, and the +/// pull `−(n + 1)kW/r` of the Earth's own redistributed mass. With the IERS Conventions' nominal +/// `h₂ = 0.6078` and `k₂₀ = 0.30190` it is 1.155, the conventional 1.16. +pub fn delta(n: u32, h: f64, k: f64) -> f64 { + let n = f64::from(n); + 1.0 + 2.0 * h / n - (n + 1.0) * k / n +} + +/// The diminishing factor `γ = 1 + k − h` for Love numbers `h` and `k` of one degree: the +/// observed tide in tilt, or in a potential difference measured between points that rise with +/// the ground, over a rigid Earth's. About 0.69 at degree 2 (Mäkinen 2021, Eq. 52). +pub fn gamma(h: f64, k: f64) -> f64 { + 1.0 + k - h +} + +/// The lunar and solar tide in gravity at geodetic latitude `lat` and east longitude `lon` (rad), +/// height `h` (m) and Julian epoch `epoch` (UT, as from [`geopotential::julian_epoch`]), m/s²: +/// `(moon, sun, total)`, positive upward. A tide that pulls up lowers gravity, so it is also the +/// correction to add to a reading. +/// +/// [Longman's (1959)](https://doi.org/10.1029/JZ064i012p02351) closed form from the mean +/// orbital elements: the lunar tide to degree 3, the solar to degree 2, both times +/// [`delta`]`(2, 0.612, 0.303)` = 1.1575, with Longman's own cgs constants and Love numbers, as +/// [LongmanTide](https://github.com/jrleeman/LongmanTide) (MIT) codes it. Its masses make the +/// lunar tide 0.09% and the solar 0.2% larger than today's `GM`s, and its degree-2 factor also +/// scales the degree-3 term, whose own `δ₃` is 1.07. Averaged over the 18.61-year nodal cycle it +/// is `δ₂` times the permanent tide of [`permanent`], to 0.2 µGal. +pub fn longman(lat: f64, lon: f64, h: f64, epoch: f64) -> (f64, f64, f64) { + // Longman's constants, cgs: G, the lunar and solar masses, the lunar eccentricity, the ratio + // of the mean motions of the sun and moon, the mean distances of the moon and sun, the + // Earth's equatorial radius, the lunar orbit's inclination to the ecliptic and the obliquity. + let (mu, moon, sun) = (6.673e-8, 7.3537e25, 1.993e33); + let (e, m, c, c1, a) = (0.05490, 0.074804, 3.84402e10, 1.495e13, 6.378270e8); + let (i, omega): (f64, f64) = (0.08979719, 23.452f64.to_radians()); + // Julian centuries from 1899 December 31, 12h, J2000.0 less one century, and the hour angle + // of the mean sun west of the station, from the days since J2000.0, 12h. + let tc = (epoch - 2000.0) / 100.0 + 1.0; + let day = (epoch - 2000.0) * 365.25; + let t = TAU * (day - day.floor()) + lon; + let (tc2, tc3) = (tc * tc, tc * tc * tc); + + // The moon: mean longitude, perigee, the sun's mean longitude and the ascending node. + let s = 4.72000889397 + 8399.70927456 * tc + 3.45575191895e-05 * tc2 + 3.49065850399e-08 * tc3; + let p = 5.83515162814 + 71.0180412089 * tc + 0.000180108282532 * tc2 + 1.74532925199e-07 * tc3; + let hs = 4.88162798259 + 628.331950894 * tc + 5.23598775598e-06 * tc2; + let node = 4.52360161181 - 33.757146295 * tc + 3.6264063347e-05 * tc2 + 3.39369576777e-08 * tc3; + let incl = (omega.cos() * i.cos() - omega.sin() * i.sin() * node.cos()).acos(); + let nu = (i.sin() * node.sin() / incl.sin()).asin(); + let chi = t + hs - nu; + let (cos_alpha, sin_alpha) = (node.cos() * nu.cos() + node.sin() * nu.sin() * omega.cos(), omega.sin() * node.sin() / incl.sin()); + let sigma = s - (node - 2.0 * (sin_alpha / (1.0 + cos_alpha)).atan()); + let l = sigma + 2.0 * e * (s - p).sin() + 1.25 * e * e * (2.0 * (s - p)).sin() + + 3.75 * m * e * (s - 2.0 * hs + p).sin() + 11.0 / 8.0 * m * m * (2.0 * (s - hs)).sin(); + + // The sun: perigee, the Earth's eccentricity and the sun's longitude. + let p1 = 4.90822941839 + 0.0300025492114 * tc + 7.85398163397e-06 * tc2 + 5.3329504922e-08 * tc3; + let e1 = 0.01675104 - 0.00004180 * tc - 0.000000126 * tc2; + let chi1 = t + hs; + let l1 = hs + 2.0 * e1 * (hs - p1).sin(); + + // Zenith angles, the station's radius and the two distances. + let (sl, cl) = lat.sin_cos(); + let cos_theta = sl * incl.sin() * l.sin() + cl * ((0.5 * incl).cos().powi(2) * (l - chi).cos() + (0.5 * incl).sin().powi(2) * (l + chi).cos()); + let cos_phi = sl * omega.sin() * l1.sin() + cl * ((0.5 * omega).cos().powi(2) * (l1 - chi1).cos() + (0.5 * omega).sin().powi(2) * (l1 + chi1).cos()); + let r = (1.0 / (1.0 + 0.006738 * sl * sl)).sqrt() * a + h * 100.0; + let (ap, ap1) = (1.0 / (c * (1.0 - e * e)), 1.0 / (c1 * (1.0 - e1 * e1))); + let d = 1.0 / (1.0 / c + ap * e * (s - p).cos() + ap * e * e * (2.0 * (s - p)).cos() + + 15.0 / 8.0 * ap * m * e * (s - 2.0 * hs + p).cos() + ap * m * m * (2.0 * (s - hs)).cos()); + let big_d = 1.0 / (1.0 / c1 + ap1 * e1 * (hs - p1).cos()); + + let gm = mu * moon * r / (d * d * d) * (3.0 * cos_theta * cos_theta - 1.0) + + 1.5 * mu * moon * r * r / (d * d * d * d) * (5.0 * cos_theta.powi(3) - 3.0 * cos_theta); + let gs = mu * sun * r / (big_d * big_d * big_d) * (3.0 * cos_phi * cos_phi - 1.0); + // Gal to m/s², times the gravimetric factor. + let k = 1e-2 * delta(2, 0.612, 0.303); + (gm * k, gs * k, (gm + gs) * k) +} + +/// The permanent tide at geodetic latitude `lat` (rad) and ellipsoidal height `h` (m): +/// `(W_T, g_T)`, the time-averaged tide-generating potential (m²/s²) and its part of gravity along +/// the ellipsoidal normal (m/s²), from [`PERMANENT_A`] (Mäkinen 2021, Eq. 18). +/// +/// [Mäkinen (2021)](https://doi.org/10.1007/s00190-021-01541-5): `g_T` is what to add to +/// zero-tide gravity (the IAG convention) for mean-tide gravity, −30.49 µGal at the equator and +/// +60.77 at the poles (his Eq. 20); `W_T/γ`, with `γ` normal gravity, is the height of the +/// mean-tide geoid above the zero-tide one, +99.4 mm at the equator and −196.4 mm at the poles +/// (Eq. 21). In Cartesian coordinates `W_T = (A/3a²)(2z² − p²)`, so this is exact for any height, +/// and his polynomials in `sin²φ` (Eqs. 19–21) agree with it to their printed digits. A tide-free +/// value differs from the zero-tide one by the Earth's response: `k` times `W_T`, and `δ − 1` +/// times `g_T` for the [`delta`] its tide correction used. +/// +/// Mean-tide gravity cannot be taken from a field whose `C̄20` was moved to +/// [`TideSystem::MeanTide`](geopotential::TideSystem::MeanTide): there `W_T` is a degree-2 +/// term of the Earth's own field, falling as `r⁻³` instead of growing as `r²`, which matches it on +/// the sphere `r = R` but gives +45.7 µGal at the equator instead of −30.5. +pub fn permanent(lat: f64, h: f64) -> (f64, f64) { + let (r, psi) = geocentric(lat, h); + let (p, z) = (r * psi.cos(), r * psi.sin()); + let k = PERMANENT_A / (3.0 * WGS84_A * WGS84_A); + // g_T = −∂W_T/∂h, along the normal (cos φ, sin φ) in the meridian plane. + (k * (2.0 * z * z - p * p), 2.0 * k * (p * lat.cos() - 2.0 * z * lat.sin())) +} + +#[cfg(test)] +mod tests { + use super::*; + use geopotential::{Coefficients, TideSystem, julian_epoch}; + use metrology::{GM_EARTH, WGS84_F}; + + const MICROGAL: f64 = crate::local::MICROGAL; + + #[test] + fn longmans_tide() { + // LongmanTide's test vector (tests/test_basics.py), mGal. The Julian epoch holds the time + // to about 4 µs, which is 3×10⁻¹¹ mGal. + let (lat, lon) = (40.7914f64.to_radians(), 282.1414f64.to_radians()); + let (moon, sun, total) = longman(lat, lon, 370.0, julian_epoch(2015, 4, 23, 0, 0)); + for (got, want) in [(moon, 0.0324029651226), (sun, -0.0288682178454), (total, 0.00353474727722)] { + assert!((got / 1e-5 - want).abs() < 1e-10, "{} against {want}", got / 1e-5); + } + } + + #[test] + fn averaged_over_the_nodal_cycle_it_is_the_permanent_tide() { + // The mean of the routine tide is not zero (Honkasalo 1964): over 18.61 years Longman's + // tide is δ₂ times −g_T, Mäkinen's average from the spectral expansions, to 0.2 µGal (0.19 + // at 45°). + let step = 1.37 / 24.0 / 365.25; // years, incommensurate with the day and month + let n = (18.61 / step) as usize; + for deg in [0.0, 45.0, 90.0] { + let lat = f64::to_radians(deg); + let mean = (0..n).map(|k| longman(lat, 0.0, 0.0, 2000.0 + k as f64 * step).2).sum::() / n as f64; + let want = -delta(2, 0.612, 0.303) * permanent(lat, 0.0).1; + assert!((mean - want).abs() < 0.25 * MICROGAL, "{deg}°: {} against {} µGal", mean / MICROGAL, want / MICROGAL); + } + } + + #[test] + fn love_number_factors() { + // Longman's own Love numbers, and the IERS Conventions 2010 nominal values: h₂ = 0.6078 + // (ch. 7), k₂₀ = 0.30190 (Table 6.3), h₃ = 0.292 (ch. 7), k₃ = 0.093 (Table 6.3). + assert!((delta(2, 0.612, 0.303) - 1.1575).abs() < 1e-15); + assert!((delta(2, 0.6078, 0.30190) - 1.15495).abs() < 1e-14); + assert!((delta(2, 0.6078, 0.30190) - 1.16).abs() < 0.01); // the conventional factor + assert!((delta(3, 0.292, 0.093) - 1.070_666_666_666_667).abs() < 1e-15); + // Mäkinen's Eq. 52, 1 + k − h ≈ 0.69, with k ≈ 0.30 and the nominal h. + assert!((gamma(0.6078, 0.30190) - 0.6941).abs() < 1e-12 && (gamma(0.6078, 0.30190) - 0.69).abs() < 0.005); + // A rigid Earth (h = k = 0) shows the tide as it is; a fluid one (h = 1 + k) shows no tilt. + assert_eq!((delta(2, 0.0, 0.0), gamma(0.0, 0.0), gamma(1.25, 0.25)), (1.0, 1.0, 0.0)); + } + + #[test] + fn makinens_permanent_tide() { + // Eq. 18's A follows from the IERS H₀ = −0.31460 m (Eq. 16): A = (3/2)H₀√(5/4π)g_e(a/R_e)². + let a_iers = 1.5 * -0.31460 * (5.0 / (4.0 * core::f64::consts::PI)).sqrt() * 9.79828685 * (WGS84_A / 6378136.55f64).powi(2); + assert!((a_iers - PERMANENT_A).abs() < 5e-5, "{a_iers}"); + // Eqs. 19–21 in sin²φ, to their printed digits, with GRS80 normal gravity (Moritz 2000). + let e2 = WGS84_F * (2.0 - WGS84_F); + let gamma0 = |s2: f64| 9.780_326_771_5 * (1.0 + 0.001_931_851_353 * s2) / (1.0 - e2 * s2).sqrt(); + for deg in (0..=90).step_by(5) { + let lat = f64::to_radians(f64::from(deg)); + let s2 = lat.sin().powi(2); + let (w, g) = permanent(lat, 0.0); + assert!((w - (0.9722 - 2.8841 * s2 - 0.0195 * s2 * s2)).abs() < 5e-5, "{deg}°: W_T {w}"); + assert!((g / MICROGAL - (-30.49 + 90.95 * s2 + 0.31 * s2 * s2)).abs() < 0.005, "{deg}°: g_T {}", g / MICROGAL); + assert!((w / gamma0(s2) * 1e3 - (99.40 - 295.41 * s2 - 0.42 * s2 * s2)).abs() < 0.005, "{deg}°: H_T"); + } + // The spherical form −2A/a·(sin²φ − 1/3) at the equator: −30.485 µGal; the mean-tide geoid + // 9.9 cm above the zero-tide one there. + assert!((permanent(0.0, 0.0).1 / MICROGAL + 30.485).abs() < 1e-3); + assert!((permanent(0.0, 0.0).0 / 9.780_326_771_5 - 0.0994).abs() < 1e-4); + // Eq. 19's height factor 1 + 2h/a, to 4×10⁻⁶ at 1 km. + let lat = 45f64.to_radians(); + assert!((permanent(lat, 1e3).0 / permanent(lat, 0.0).0 - (1.0 + 2e3 / WGS84_A)).abs() < 5e-6); + } + + #[test] + fn the_mean_tide_c20_is_right_only_on_its_sphere() { + // Moving C̄20 to the mean-tide system adds W_T as part of the Earth's own field. At r = R on + // the equator its potential is Mäkinen's to 2×10⁻⁴ (EGM's −3.11080×10⁻⁸/√5 against + // IERS's A₀H₀ = −1.39141×10⁻⁸); its gravity is 3W/r where the tide's is −2W/r. + let mut zero = Coefficients::zeros(GM_EARTH, WGS84_A, 2); + zero.tide_system = TideSystem::ZeroTide; + let mean = zero.to_tide_system(TideSystem::MeanTide).unwrap(); + let v = |c: &Coefficients, r: f64| c.partials(r, 0.0, 0.0).v; + let added = v(&mean, WGS84_A) - v(&zero, WGS84_A); + let (w, g) = permanent(0.0, 0.0); + assert!((added / w - 1.0).abs() < 2e-4, "{}", added / w); + let dr = 1.0; + let pull = -((v(&mean, WGS84_A + dr) - v(&zero, WGS84_A + dr)) - added) / dr; + assert!((pull / g + 1.5).abs() < 1e-3 && (pull / MICROGAL - 45.7).abs() < 0.1, "{}", pull / MICROGAL); + } +} diff --git a/crates/gravitics/README.md b/crates/gravitics/README.md index f949f09..ecc9b79 100644 --- a/crates/gravitics/README.md +++ b/crates/gravitics/README.md @@ -99,6 +99,7 @@ assert!((near * 1e3 - 13.28).abs() < 0.01); | `gem` | Gravitoelectromagnetism in three conventions with the spin-2 factor, the gravitomagnetic dipole as a `Field`, the extremal charge-to-mass ratio, and Blackett's magnetism-from-spin | | `anomalies` | CODATA 2022's sixteen G measurements with their correlations, their generalised-least-squares mean and Birge ratio; Anderson's flyby formula; column densities and shielding against Majorana and lunar laser ranging; Pioneer's thermal decay | | `modified` | The radial acceleration relation ([McGaugh et al. 2016](https://arxiv.org/abs/1609.05917)), Tully–Fisher velocity and MOND radius; Yukawa forces with the current limits ([Lee et al. 2020](https://arxiv.org/abs/2002.11761)), the dark-energy length and boson mass; Einstein–Cartan density; Weber-like perihelion and mass shift | +| `equivalence` | Tests of the equivalence principle: the Eötvös parameter η, the MICROSCOPE, Eöt-Wash and atom-interferometer results, the coupling to B, L or B − L that a bound on η allows for a vector or scalar boson, and the Yukawa form factor of a layered sphere ([Touboul et al. 2022](https://arxiv.org/abs/2209.15487); [Fayet 2025](https://arxiv.org/abs/2507.02723)) | | `quantum` | Bose–Marletto–Vedral entanglement phases ([Bose et al. 2017](https://arxiv.org/abs/1707.06050)); the Diósi–Penrose X-ray rate ([Donadi et al. 2021](https://arxiv.org/abs/2111.13490)) | | `hfgw` | High-frequency gravitational waves: the Li–Baker detector (flux conventions, graviton rate, the interference signal, the shot-noise time) and the Pais patents (radiated flux, surface field, magnetic conversion as he uses it and in SI) | | `superconductors` | The London moment; frame dragging of a ring in general relativity against Tajmar's report; the Cooper-pair density and gravitational London depth; a perfect shield's shadow; Podkletnov's pendulum and Marx bank; the measured bounds | @@ -152,6 +153,8 @@ Every test reproduces a published number from the paper's own inputs, mostly to - **Alternatives:** the RAR's published shape and deep-MOND limit; Eöt-Wash's 5.1 meV boson; Ruggiero & Tartaglia's Einstein–Cartan densities; the Weber perihelion ratios against a direct orbit integration. +- **Equivalence principle:** MICROSCOPE's final η; Fayet's coupling limits on B, L and B − L for + both spins, his form factors (homogeneous and five-shell Earth) and their weakening at 710 km. - **Quantum gravity:** Bose et al.'s phases; Donadi et al.'s germanium rate. The formulas were independently re-derived from their sources and every test value recomputed; diff --git a/crates/gravitics/src/equivalence.rs b/crates/gravitics/src/equivalence.rs new file mode 100644 index 0000000..e1f7710 --- /dev/null +++ b/crates/gravitics/src/equivalence.rs @@ -0,0 +1,214 @@ +//! Tests of the weak equivalence principle: the Eötvös parameter η, what the best experiments +//! measured, and how strongly a new long-range force may couple to baryon or lepton number +//! before MICROSCOPE would have seen it. +//! +//! A new force from a boson coupled to a charge `Q` (baryon number `B`, lepton number `L`, or +//! `B − L`) pulls on each body in proportion to its charge per unit mass, which differs between +//! elements, so two bodies of different composition fall differently. MICROSCOPE compared +//! titanium and platinum in orbit for two years and found no difference at 10⁻¹⁵ +//! ([Touboul et al. 2022](https://arxiv.org/abs/2209.15487)); in +//! [Fayet's (2025)](https://arxiv.org/abs/2507.02723) reading this bounds a `B − L` coupling at +//! 1.1×10⁻²⁵, 10²⁵ times weaker than electromagnetism. +//! +//! The identifiers say `eta`, since `gravimetry::local::EOTVOS` is the unit of gravity gradient. +//! +//! ``` +//! use gravitics::equivalence::{MICROSCOPE_LIMITS, Mediator, TESTS, coupling_bound}; +//! +//! // MICROSCOPE's η(Ti, Pt) = (−1.5 ± 2.3 ± 1.5)×10⁻¹⁵: 2.7×10⁻¹⁵ in all. +//! assert!((TESTS[0].combine() - 2.746e-15).abs() < 1e-18); +//! // The Earth's (B − L) per atomic mass unit, and Ti's less Pt's (Fayet 2025, Eqs. 6–7). +//! let g = coupling_bound(Mediator::Vector, 0.5138, -0.05625, MICROSCOPE_LIMITS); +//! assert!((g / 1.067e-25 - 1.0).abs() < 1e-3); // "|g_B−L| < 1.1×10⁻²⁵" +//! ``` + +use metrology::{ALPHA, AMU, C, G, HBAR}; + +/// The Eötvös parameter `η = 2(a_A − a_B)/(a_A + a_B)` of two bodies falling with accelerations +/// `a_a` and `a_b` (Touboul et al. 2022, Eq. 1). +pub fn eta(a_a: f64, a_b: f64) -> f64 { + 2.0 * (a_a - a_b) / (a_a + a_b) +} + +/// One measurement of η with its statistical and systematic errors (1σ). Where a paper quotes +/// one combined error it is in `stat` and `syst` is zero. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Test { + /// The experiment and the two materials. + pub name: &'static str, + /// η of the first material against the second. + pub eta: f64, + /// The statistical error. + pub stat: f64, + /// The systematic error. + pub syst: f64, + /// Where it comes from. + pub source: &'static str, +} + +impl Test { + /// The two errors in quadrature. + pub fn combine(&self) -> f64 { + self.stat.hypot(self.syst) + } +} + +/// Published tests of the equivalence principle, the tightest first: in orbit (MICROSCOPE's +/// Ti–Pt result and its Pt–Pt control), a torsion balance, and atoms in free fall. +pub const TESTS: [Test; 4] = [ + Test { name: "MICROSCOPE Ti–Pt", eta: -1.5e-15, stat: 2.3e-15, syst: 1.5e-15, source: "https://arxiv.org/abs/2209.15487" }, + Test { name: "MICROSCOPE Pt–Pt (control)", eta: 0.0, stat: 1.1e-15, syst: 2.3e-15, source: "https://arxiv.org/abs/2209.15487" }, + Test { name: "Eöt-Wash Be–Ti, towards the Earth", eta: 0.3e-13, stat: 1.8e-13, syst: 0.0, source: "https://arxiv.org/abs/0712.0607" }, + Test { name: "Stanford ⁸⁵Rb–⁸⁷Rb", eta: 1.6e-12, stat: 1.8e-12, syst: 3.4e-12, source: "https://arxiv.org/abs/2005.11624" }, +]; + +/// MICROSCOPE's one-sided 95% limits on η(Ti, Pt): below 4.5×10⁻¹⁵ if positive, below 6.5×10⁻¹⁵ +/// in size if negative. The result leans negative, so a positive η is the more constrained +/// (Fayet 2025, Eq. 16). +pub const MICROSCOPE_LIMITS: (f64, f64) = (4.5e-15, 6.5e-15); + +/// `α/(G u²) = e²/(4πε₀ G u²)`, 1.2536×10³⁶: the electric repulsion of two unit charges over the +/// gravitational attraction of two atomic mass units, the scale of Fayet's Eq. 7. +pub const ALPHA_OVER_G_U2: f64 = ALPHA * HBAR * C / (G * AMU * AMU); + +/// The spin of the boson that carries a new force. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Mediator { + /// Spin 1, like the photon: like charges repel. + Vector, + /// Spin 0: like charges attract. + Scalar, +} + +/// The η a new long-range force puts between two bodies falling towards the Earth (Fayet 2025, +/// Eq. 7): `η = ∓(α/(G u²))ε²(Q/A_r)_⊕Δ(Q/A_r)`, the upper sign for a vector mediator. `epsilon` +/// is the coupling in units of `e`, `earth` the Earth's charge per atomic mass unit and +/// `difference` the first body's less the second's. For a range much longer than the Earth's +/// radius; a finite range `λ` multiplies η by `(1 + r/λ)e^(−r/λ)` at distance `r` +/// ([`yukawa_factor`](crate::modified::yukawa_factor) less one) and by the Earth's form factor +/// (see [`sphere_form_factor`]) (Fayet's Eq. 94). +pub fn eta_from_coupling(mediator: Mediator, epsilon: f64, earth: f64, difference: f64) -> f64 { + let sign = match mediator { + Mediator::Vector => -1.0, + Mediator::Scalar => 1.0, + }; + sign * ALPHA_OVER_G_U2 * epsilon * epsilon * earth * difference +} + +/// The largest coupling `|g| = |ε|√(4πα)` that one-sided limits `(positive, negative)` on η +/// allow, for the same inputs as [`eta_from_coupling`]: the mediator and the sign of the charge +/// difference fix the sign of η, and so which limit applies. With [`MICROSCOPE_LIMITS`] and a +/// `B − L` charge it is 1.07×10⁻²⁵ for a vector and 1.28×10⁻²⁵ for a scalar, which Fayet rounds +/// to 1.1 and 1.3 (his Eqs. 24–25). A finite range raises it by `F^(−1/2)`, where `F` is the +/// factor that range puts on η. +pub fn coupling_bound(mediator: Mediator, earth: f64, difference: f64, (positive, negative): (f64, f64)) -> f64 { + let per_epsilon2 = eta_from_coupling(mediator, 1.0, earth, difference); + let limit = if per_epsilon2 > 0.0 { positive } else { negative }; + (limit / per_epsilon2.abs()).sqrt() * (4.0 * core::f64::consts::PI * ALPHA).sqrt() +} + +/// How much more strongly a homogeneous sphere of radius `R` sources a Yukawa force of range `λ` +/// than a point of the same charge, at `x = R/λ`: `φ(x) = 3(x cosh x − sinh x)/x³` (Fayet 2025, +/// Eq. 37). It is 1 for a long range and grows as `(3/2)eˣ/x²`, since the charge near the surface +/// is nearer. A layered body is the mass-weighted sum of its shells' terms (his Eq. 80). +pub fn sphere_form_factor(x: f64) -> f64 { + if x.abs() < 1.0 { + // Σ 3x²ᵏ/((2k + 1)!(2k + 3)): the closed form loses every digit to cancellation as x → 0. + let (mut term, mut sum) = (1.0, 0.0); + for k in 0..12 { + sum += 3.0 * term / f64::from(2 * k + 3); + term *= x * x / f64::from((2 * k + 2) * (2 * k + 3)); + } + sum + } else { + 3.0 * (x * x.cosh() - x.sinh()) / x.powi(3) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::modified::{boson_mass_ev, yukawa_factor}; + + #[test] + fn eotvos_parameters() { + // η is antisymmetric, and Δa/g to first order in Δa. + let g = 7.95; + assert!(eta(3.0, 1.0) == 1.0 && eta(1.0, 2.0) == -eta(2.0, 1.0) && eta(g, g) == 0.0); + assert!((eta(g + 1e-6, g) / (1e-6 / g) - 1.0).abs() < 1e-7); + // Touboul et al. (2022), Eq. 5 and the control: η(Ti, Pt) = [−1.5 ± 2.3 ± 1.5]×10⁻¹⁵, + // 2.746×10⁻¹⁵ in quadrature, a differential pull of 2.2×10⁻¹⁴ m/s² at 710 km. + let [ti_pt, pt_pt, eot_wash, stanford] = TESTS; + assert!((ti_pt.combine() / 2.745_906e-15 - 1.0).abs() < 1e-6); + assert!((ti_pt.combine() * metrology::GM_EARTH / (6371e3 + 710e3f64).powi(2) / 2.18e-14 - 1.0).abs() < 0.01); + assert!(ti_pt.eta.abs() < ti_pt.combine() && pt_pt.eta.abs() < pt_pt.combine()); + // Schlamminger et al. (2008), Eq. 2, and Asenbaum et al. (2020): 3.85×10⁻¹² in all. + assert!(eot_wash.eta.abs() < eot_wash.combine() && (stanford.combine() / 3.847e-12 - 1.0).abs() < 1e-3); + } + + #[test] + fn fayets_coupling_bounds() { + assert!((ALPHA_OVER_G_U2 / 1.2536e36 - 1.0).abs() < 5e-5); + // Eq. 6: Δ(Q/A_r) of Ti less Pt. Eq. 7 prints the combined coefficients ∓1.00×10³³ (B), + // ±3.623×10³⁴ (B − L) and ∓3.482×10³⁴ (L); they give the Earth's charges per mass. Its L + // and B − L add up to 1.0008 nucleons per u, between ¹⁶O's 1.0003 and ⁵⁶Fe's 1.0012; B's + // own coefficient and Δ(B/A_r) carry two or three digits, so its 1.010 is good to 1%. + let (db, dl): (f64, f64) = (0.00079, 0.05704); + let dbl = db - dl; + assert!((dbl + 0.05625).abs() < 1e-15); + let earth = |coefficient: f64, d: f64| coefficient / (ALPHA_OVER_G_U2 * d.abs()); + let (eb, ebl, el) = (earth(1.00e33, db), earth(3.623e34, dbl), earth(3.482e34, dl)); + assert!((ebl + el - 1.0008).abs() < 1e-4 && (eb - (ebl + el)).abs() < 0.01); + assert!((eta_from_coupling(Mediator::Vector, 1.0, ebl, dbl) / 3.623e34 - 1.0).abs() < 1e-12); + // Eqs. 24–25 (95% CL), and note 40's unrounded 1.07 and 1.28 for B − L. + for (mediator, q, d, want) in [ + (Mediator::Vector, ebl, dbl, 1.1e-25), (Mediator::Vector, eb, db, 7.7e-25), (Mediator::Vector, el, dl, 1.3e-25), + (Mediator::Scalar, ebl, dbl, 1.3e-25), (Mediator::Scalar, eb, db, 6.4e-25), (Mediator::Scalar, el, dl, 1.1e-25), + ] { + let got = coupling_bound(mediator, q, d, MICROSCOPE_LIMITS); + assert!((got / want - 1.0).abs() < 0.05, "{mediator:?}: {got:e} against {want:e}"); + } + let vector = coupling_bound(Mediator::Vector, ebl, dbl, MICROSCOPE_LIMITS); + let scalar = coupling_bound(Mediator::Scalar, ebl, dbl, MICROSCOPE_LIMITS); + assert!((vector / 1.067e-25 - 1.0).abs() < 1e-3 && (scalar / 1.283e-25 - 1.0).abs() < 1e-3); + // The two differ by √(6.5/4.5), the dissymmetry of the result. + assert!((scalar / vector - (6.5f64 / 4.5).sqrt()).abs() < 1e-12); + } + + #[test] + fn the_form_factor_of_a_sphere() { + // Fayet's Table IV, homogeneous sphere, to its printed digits. + for (x, want) in [(0.5, 1.0252), (1.0, 1.1036), (2.0, 1.4616), (3.0, 2.2428), (5.0, 7.1243), (10.0, 297.36), (20.0, 1.728e6), (30.0, 1.722e10), (50.0, 3.049e18), (100.0, 3.992e39), (200.0, 2.696e82)] { + assert!((sphere_form_factor(x) / want - 1.0).abs() < 2.5e-4, "x = {x}: {}", sphere_form_factor(x)); + } + // 1 + x²/10 + x⁴/280 + x⁶/15120 near zero, where the closed form has no digits left; the + // series and the closed form meet at x = 1. + for x in [0.0, 1e-8, 1e-3, 0.1] { + let x2: f64 = x * x; + let want = 1.0 + x2 / 10.0 + x2 * x2 / 280.0 + x2 * x2 * x2 / 15120.0; + assert!((sphere_form_factor(x) - want).abs() < 1e-8 * x2 * x2 * x2 + 2e-16, "{x}"); + } + let (below, above) = (sphere_form_factor(1.0 - 1e-15), sphere_form_factor(1.0)); + assert!((below - above).abs() < 4e-15); + assert!((sphere_form_factor(-2.0) - sphere_form_factor(2.0)).abs() < 1e-15); + } + + #[test] + fn a_layered_earth_at_satellite_height() { + // Fayet's five-shell Earth (Eq. 80): Φ(x) = Σ yᵢ φ(Rᵢx/R), Table IV's Φ(3) = 2.0021. + let shells = [(0.5, 1.0), (0.1333, 0.9953), (0.1715, 0.8948), (0.1929, 0.5462), (0.0023, 0.1916)]; + let earth = |x: f64| shells.iter().map(|&(y, s)| y * sphere_form_factor(s * x)).sum::(); + assert!((earth(3.0) / 2.0021 - 1.0).abs() < 1e-5 && (earth(0.0) - 1.0).abs() < 1e-15); + // Table V: at MICROSCOPE's 710 km, a boson of mass m weakens η by F = (1 + r/λ)e^(−r/λ)Φ, + // so the coupling limits grow by G = F^(−1/2). + let (radius, r) = (6371e3, 6371e3 + 710e3); + let hbar_c = HBAR * C / metrology::EV; // eV m + for (m, f, g) in [(1e-14, 0.957, 1.022), (2e-14, 0.867, 1.074), (5e-14, 0.574, 1.320), (1e-13, 0.281, 1.886), (2e-13, 0.0927, 3.28)] { + let lambda = hbar_c / m; + assert!((boson_mass_ev(lambda) / m - 1.0).abs() < 1e-12); + let big_f = (yukawa_factor(r, 1.0, lambda) - 1.0) * earth(radius / lambda); + // To the table's three digits. + assert!((big_f / f - 1.0).abs() < 2e-3 && (big_f.powf(-0.5) / g - 1.0).abs() < 2e-3, "m = {m:e}: F = {big_f}"); + } + } +} diff --git a/crates/gravitics/src/lib.rs b/crates/gravitics/src/lib.rs index 070b296..fef4ce7 100644 --- a/crates/gravitics/src/lib.rs +++ b/crates/gravitics/src/lib.rs @@ -3,6 +3,7 @@ pub mod anomalies; pub mod electrodes; pub mod electrogravity; +pub mod equivalence; pub mod gem; pub mod hfgw; pub mod modified;