Towards an agnostic algorithm for sampling empirical structure models: The case of Uranus and Neptune

1. Problem Statement

Determining the internal density profiles of planets like Uranus and Neptune is a severely degenerate inverse problem.

• Data Scarcity vs. Model Complexity: Observational data (total mass, radius, and a few gravitational moments $J_{2n}$) are limited (fewer than 10 parameters), whereas a density profile requires orders of magnitude more data points (e.g., 1024 layers).
• Bias in Current Methods: Traditional approaches rely on Markov Chain Monte Carlo (MCMC) methods. To make these computationally feasible, researchers often parametrize density profiles using specific functional forms (e.g., polytropes, polynomials, or sigmoid functions). These parametric assumptions introduce bias, potentially excluding valid physical solutions or failing to explore the full solution space.
• Goal: The authors aim to develop an agnostic algorithm that samples the entire space of possible interior density profiles with minimal assumptions, avoiding the biases introduced by parametric simplifications or the slow convergence of MCMC.

2. Methodology

The authors propose a novel optimization-based gradient-descent approach that replaces MCMC.

Core Algorithm

• Free Parameters: Instead of parametric functions, the density profile is defined by $N=1024$ discrete density values at evenly spaced level surfaces. These $N$ values are treated as free parameters, constrained only by:
  1. No density plateaus or inversions (monotonically increasing density).
  2. Maximum density $\rho_{max}$ (set to 20,000 kg m$^{-3}$).
  3. Minimum density increase per step ($\Delta_{min}$).
• Optimization Scheme:
  • The algorithm minimizes a cost function based on the difference between calculated and observed gravitational moments ($J_2, J_4$) and the planetary mass.
  • It utilizes the Theory of Figures (ToF) to calculate gravitational moments from density profiles.
  • Gradient Descent: The authors use the Adam optimizer (Kingma & Ba, 2017) to navigate the parameter space.
  • Projection: After each step, parameters are projected to enforce the minimum density increase constraint.
• Unbiased Initialization: To ensure the starting profiles do not bias the search, the authors employ a binary subdivision algorithm. This recursively subdivides the density range to generate random decreasing functions, ensuring that the probability of generating a profile is uniform across the solution space volume (correcting for the non-uniform density of the simplex).
• Unbiased Search: The authors argue that their gradient descent method does not systematically favor global solutions over local ones, verified empirically by analyzing "discoverability" (Appendix C).

Simulation Setup

• Targets: Uranus and Neptune.
• Hardware: AMD EPYC 7702 CPU (96 cores).
• Scale: Hundreds of thousands of runs (approx. 700k for Uranus, 730k for Neptune).
• Success Criteria: Runs are considered successful if they respect $\rho_{max}$ and match observed gravitational moments within a multivariate 3-sigma threshold.

3. Key Contributions

1. Agnostic Sampling Algorithm: A new, computationally efficient method to sample the full space of empirical density profiles without relying on restrictive parametric forms (polytropes, etc.).
2. Comprehensive Solution Space Mapping: The first study to map the entire solution space for Uranus and Neptune, revealing the true degeneracy and diversity of possible interior structures.
3. Statistical Guidance on Discontinuities: A quantitative framework for modellers constructing models with fixed layers. The paper provides a statistical relation between the steepness threshold ($c_D$) used to define a "discontinuity" and the expected number of such discontinuities.
4. Benchmarking: Extensive comparison with previous literature (e.g., Nettelmann et al., Helled & Fortney, Movshovitz & Fortney), highlighting where previous models are biased or incomplete.

4. Key Results

Density Profiles and Solution Space

• Deep Interior Degeneracy: The innermost regions of the planets are the least constrained by current data. The 95% confidence interval for central density is $\sim 11$ g cm$^{-3}$ for Uranus and $\sim 13$ g cm$^{-3}$ for Neptune.
• Median Profiles: The median density profiles show a gradual increase consistent with composition gradients, maintaining a consistent slope toward the center (unlike some previous models that flatten out).
  • Uranus median central density: $\sim 6$ g cm$^{-3}$.
  • Neptune median central density: $\sim 7.5$ g cm$^{-3}$.
• Comparison to Literature:
  • Models with sharp discontinuities (e.g., Nettelmann et al.) tend to cluster around high central densities.
  • Models without discontinuities (e.g., Helled & Fortney) often cluster around the peak of the distribution but miss the median and the full diversity of high-density solutions.
  • The authors' ensemble captures a much wider spread, particularly for higher central densities.

Discontinuities

• Location: If discontinuities exist, they are concentrated around a normalized radius of $r/R \approx 0.65$ for Uranus and $r/R \approx 0.7$ for Neptune. This likely marks the transition from H/He-dominated outer layers to water/rock-dominated interiors.
• Frequency vs. Steepness:
  • If a discontinuity is defined as a density gradient $> 0.02$ kg m$^{-4}$, most solutions have at most one such discontinuity.
  • If the threshold is lower (e.g., $0.002$ kg m$^{-4}$), the number of discontinuities increases significantly.
  • Implication: Previous studies assuming a fixed number of discontinuities (e.g., exactly 3) exclude a significant portion of the valid solution space. The paper provides a lookup table (Table D.1) to guide modellers on expected discontinuity counts based on their chosen steepness criterion.

Derived Quantities

• Moments of Inertia: The average normalized moment of inertia ($\lambda$) is $0.22595$ for Uranus and $0.23597$ for Neptune, consistent with or slightly lower than previous findings.
• Flattening: The inferred flattening ratios agree well with observational data, validating the ToF implementation.

5. Significance and Conclusion

• Paradigm Shift: The paper demonstrates that progress in planetary science requires not just better data, but broader, more flexible interpretation methods. By embracing uncertainty rather than constraining it with parametric assumptions, the authors reveal the true complexity of planetary interiors.
• Robustness: The algorithm is highly efficient (approx. 20 seconds per run) and faithful to the data, capable of handling the high dimensionality of the problem without MCMC's slow convergence.
• Future Applications: This methodology can be applied to any planetary object with gravitational field data, including exoplanets. It provides a rigorous baseline for future compositional and thermal modeling, ensuring that subsequent physical interpretations are built upon a complete and unbiased understanding of the density solution space.

Limitations & Caveats:

• The generation of starting profiles introduces minor artifacts (kinks) at specific radii (inverses of multiples of two), though the optimization process largely removes these.
• The method relies on linear approximations for gradients; more complex problems with higher-order moments might require stochastic gradients or differentiable ToF models (e.g., via machine learning).
• Current results are limited by the precision of existing gravitational moment data (especially for Neptune).