Neural Symplectic Integrator with Hamiltonian
Inductive Bias for the Gravitational N -body Problem

Maxwell X. Cai
SURF Corporative / Leiden University
Science Park 140, 1098 XG Amsterdam / Niels Bohrweg 2, 2300 RA Leiden
The Netherlands
maxwell.cai@surf.nl
Simon Portegies Zwart
Leiden University
Niels Bohrweg 2, 2300 RA Leiden
spz@strw.leidenuniv.nl

Damian Podareanu
SURF Corporative
Science Park 140, 1098 XG Amsterdam
damian.podareanu@surf.nl

Abstract
The gravitational N -body problem, which is fundamentally important in astrophysics to predict the motion of N celestial bodies under the mutual gravity of eachother, is usually solved numerically because there is no known general analyticalsolution for N > 2. Can an N -body problem be solved accurately by a neuralnetwork (NN)? Can a NN observe long-term conservation of energy and orbital
angular momentum? Inspired by Wistom & Holman’s symplectic map, we present
a neural N -body integrator for splitting the Hamiltonian into a two-body part,solvable analytically, and an interaction part that we approximate with a NN. Our
neural symplectic N -body code integrates a general three-body system for 105
steps without diverting from the ground truth dynamics obtained from a traditional N -body integrator. Moreover, it exhibits good inductive bias by successfully
predicting the evolution of N -body systems that are no part of the training set.

1

Introduction

Deep neural networks (DNNs) have been widely used as universal function approximators to learn atarget distribution. Such a property is promising for physics-related applications, whose underlying
distributions are typically high-dimensional, computationally expensive, and/or difficult to derive
analytically. In astrophysics and particle physics, generative models are often used for MonteCarlo event generation [e.g., Erdmann et al., 2018, Zamudio-Fernandez et al., 2019], dimensionality
reduction [e.g., Portillo et al., 2020], classification [e.g., Fremling et al., 2021], and anomaly detection[e.g., Reis et al., 2018].
Despite the success of the aforementioned DNN applications on astrophysical data, the trainedmodels learn to encode representations from the training data, but they do not necessarily understandthe underlying physical laws. Making DNNs understand physical laws, therefore, requires special
formulations. Physics informed neural networks [PINNs, see Raissi et al., 2017a,b, Raissi et al.,2019], for example, take the form of differential equations and replace the differentiable mathematical
function f (which is sometimes unknown or expensive to calculate) by a differentiable neural network
fNN . In this way, PINNs are trained to approximate the function f , and thus differentiating the
PINNs will yield the evolution of the underlying system. Moreover, the outputs of fNN becomes
more interpretable. This idea is echoed in the recent work by Greydanus et al. [2019], in which theauthors replace the Hamiltonian of a dynamics system with a DNN and successfully apply it to solved
Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021).

problems such as pendulum, mass-spring, and the two-body systems. The resulting architecture,
namely Hamiltonian Neural Networks (HNNs), manages to learn and respect the conservation of totalenergy. Similarly, the Hamiltonian Generative Networks (HGNs) are capable of consistently learning
Hamiltonian dynamics from high-dimensional observations (such as images) without restrictive
domain assumptions, and the Neural Hamiltonian Flow (NHF), which uses Hamiltonian dynamics tomodel expressive densities [Toth et al., 2019]. Various improvements and adoptations have since been
made, see e.g., Chen et al. [2019], Sanchez-Gonzalez et al. [2019], Cranmer et al. [2020], DesmondZhong et al. [2019], Finzi et al. [2020].
All aforementioned works approximate the full Hamiltonian (or a general function f ) using a neuralnetwork. In this paper, however, we propose a formulation to improve the approximation. We will
see that, if the Hamiltonian can be rewritten as an analytically solvable part and a residual part, thenapproximating only the small residual part with a DNN will improve the prediction accuracy. As anapplication, we use this idea to construct a neural symplectic integrator for N -body simulations 1 ,based on the formulation by Wisdom and Holman [1991].

2

Methods

2.1

Theoretical Framework

Hamiltonian mechanics is a reformulation of Newtonian mechanics, defined as the sum of potential energy and kinetic energy, i.e., H = T + V , where T is the generalized kinetic energy andV is the generalized potential energy. The function H is called the Hamiltonian of the system,
defined in an abstract phase space s = (q, p), where q = (q1 , q2 , ..., qn ) is the position vector, andp = (p1 , p2 , ..., pn ) is the momentum vector. The subscript n indicates the dimensionality of thecoordinates. The dynamics can be obtained through the following set of differential equations:
q̇ ≡

dq
∂H
=
,
dt
∂p

ṗ ≡

dp
∂H
=−
.
dt
∂q

(1)

Since H is the total energy, it is a conserved quantity in systems where external perturbations do notexist or can be safely ignored. For a gravitational N -body systems, each celestial body experiences
the gravity exerted by the other N − 1 bodies, and the corresponding Hamiltonian is
H=

N
−1
X

N
−2
N
−1
X
X
p2i
mj
mi
−G
,
2mi
r
i=0
i=0
j=i+1 ij

(2)

where pi ≡ |pi | is the magnitude of the momentum vector for particle i, and rij ≡ |qi − qj | isthe distance between particles i and j. While there is no known analytic solution for the N -bodyproblems with N > 2, Wisdom and Holman [1991] point out that the Hamiltonian of an N -bodyproblem can be naturally split into two parts:
H = Hkepler + Hinter .

(3)Here Hkepler represents the Kepler motion of the bodies with respect to the center of mass, andHinter represents the perturbation of the minor bodies on each other. This splitting is rooted in theperturbation theory, and this formulation allows the integrator to reduce numerical errors by taking
the advantage that Hkepler can be solved analytically. Working in Jacobi coordinates, Hkepler can be
written as
!
N
−1
X
p̃2i
M̃i m̃i
Hkepler =
−G
.
(4)
2m̃i
r̃i
i=1
P

Pi
i−1
Here m̃i =j=0 mj /j=0 mj mi is the normalized mass for body i, p̃i = m̃i ṽi is the momentum of body i in Jacobi coordinates, and r̃i is the position vector of body i in Jacobi coordinates. Theinteracting Hamiltonian Hinter takes the form
Hinter = G
1

N
−1
X

N
−2
N
−1
X
X
m0 mi
mj
−G
mi
.
r̃i
r
i=1
i=0
j=i+1 ij

(5)

Source code available at https://github.com/maxwelltsai/neural-symplectic-integrator.

2

The Wisdom-Holman (WH) integrator is a symplectic mapping that employs a kick-drift-kick (ordrift-kick-drift) strategy. During the drift step, the integrator assumes that the secondary body orbits
the central body (or the center of mass) without experiencing any perturbations. This assumption
allows one to propagate the state analytically by solving Kepler’s equation. The error for leaving the
perturbations from other bodies unaccounted for is subsequently corrected by the kick step, in which
the mutual interactions between each body (e.g., planets) are regarded as velocity kicks. A velocity
kick ∆v = ap ∆t accounts for the cumulative effect of all planets due to the perturbing acceleration
ap over a time interval ∆t. Since a kick step only changes the velocity vector without touching theposition vector, it ensures that the evolution of position vectors is smooth and differentiable.
2.2

Neural Interacting Hamiltonian (NIH)

Hkepler is an O(N ) operation and Hinter is a more expensive O(N 2 ) operation. The former canbe solved analytically, but the latter has to be solved numerically. We solve Hkepler analyticallybut Hinter using a neural network. Following Greydanus et al. [2019], we use a simple multi-layer
perceptron (MLP) backbone network Hinter,θ to serve as a function approximator of Hinter . The
network, which we refer to as the neural interacting Hamiltonian (NIH), outputs a single scalar value
analog to the total energy of the interaction part. This network is differentiable, and by taking its
in-graph gradients with respect to q and p, we obtain the velocities and the accelerations, respectively.
The accelerations are used to calculate the velocity kicks in the WH mapping. Therefore, the networks
takes the state vector (q, p) of the underlying N -body system as an input, and generates the time
derivative (q̇θ , ṗθ ) (every quantity related to the neural network is hereafter denoted with the subscript
θ). The loss is calculated as
L == kq̇ − q̇θ k2 + kṗ − ṗθ k2 + kq̇θ × p − q × ṗθ k2 .
| {z } | {z } |{z}
L2 (v)

L2 (a)(6)

L̇residual

Here q̇ and ṗ are ground truth velocities and accelerations obtained from a traditional WH integrator.
The loss function consists of three parts: the first part L2 (v) aims to minimize the difference between
the predicted velocities and the ground truth velocities; the second part L2 (a) aims to minimize
the difference between the predicted accelerations and the ground truth acceleration. The third partL̇residual is essentially the time derivative of the orbital angular momentum L:
dL
d (q × p)
dq
dp
∂Hθ
∂Hθ
≡
=
×p+q×
=
×p−q×
dt
dt
dt
dt
∂p
∂q

(7)

In the absence of external perturbations, there is no external torque and the total angular momentumof an N -body system is conserved, hence dL/dt = 0. Therefore, minimizing the residual of dL/dt
is equivalent to forcing Hinter,θ to obey the law of angular momentum conservation.
In an attempt to improve the backbone network, we also tested a sparsely-connected MLP model[see e.g., Constantin Mocanu et al., 2017, Pieterse and Constantin Mocanu, 2019] where 90% of its
neurons are randomly zeroed out as a form of regularization. This approach forces a more efficientuse of neurons. The sparse MLP model performs slightly better than the vanilla MLP backbone. For
simplicity, hereafter we use only the vanilla MLP backbone to carry out all the tests.
2.3

Simulations, Activation Function, and Training:

The SymmetricLog Activation Function Since the neural N -body integrator is expected to handle
real astrophysical data like any traditional integrators, we need to make sure that Hinter,θ is capableof handling raw data without the need to perform data pre-processing. Given the vast dynamic ranges
in astrophysical data, however, commonly used activation functions, such as ReLU, sigmoid, andtanh have insufficient dynamic range. For example, the ReLU activation function does not activate
for negative values, and the tanh activation function saturates to 1 or -1 for absolute value > 3. As
such, we design a custom activation function SymmetricLog:
f (x) = tanh(x) log [x tanh(x) + 1] .

(8)

As shown in Fig. 1, the SymmetricLog activation function behaves similarly as tanh for values
close to zero, which delivers strong gradients for backward propagation; for larger values, it behaves
like a logarithmic function, which does not saturate or explode the gradients; it is well defined for
3

6
4

10 1

SymmetricLog
log(x + 1)
tanh(x)

10 2
10 3
dE/E0

y

2
0
2

10 4
10 5
10 6

4

10 7

6

SymmetricLog
tanh

10 8
1000 750

500

250

0
x

250

500

750 1000

0

250

500

750 1000 1250 1500 1750 2000
t [time steps]

Figure 1: Left: The activation function used in the neural interaction Hamiltonian, as compared toy = log(x + 1) and y = tanh(x). Right: Relative integration energy error dE/E0 as a function oftime. The SymmetricLog activation performs ∼ 3 orders of magnitude better in terms of energy
conservation compared to tanh.

any real values, making it possible for the neural network to handle raw astrophysical data (unlessthe quantities are smaller than 10−3 or larger than 104 , in which case data normalization is neededto achieve optimal prediction accuracy). Its symmetry with respect to zero makes it particularly
suitable for using in N -body simulations because the coordinates of celestial bodies are usuallyquasi-symmetric with respect to the center of the reference frame. We perform a test of planetarysystems with 10 planets orbiting a solar-type star. The semi-major axes of these planets vary from 0.5
au to 10 au, which means that any planets with a semi-major axis greater than ∼ 3 au will saturate thetanh activation function. ReLU and SoftPlus cannot be used because they cannot generate negative
activation outputs. As shown in Fig. 1, SymmetricLog performs ∼ 3 orders of magnitude better interms of energy conservation.
Training We deliberately train the NIH with a small number of simulations (see Appendix A.1
for the initial conditions). In doing so, we wanted to test whether the sub-optimally trained neuralnetwork has good inductive bias in transferring the physics laws it learn from the training examples
to distinctively different N -body systems. We train the neural network with only N = 3 systems(there is no NIH for N = 2), but the same network is subsequently asked to deal with problems with
arbitrary N . For brevity, we refer to a traditional Wisdom-Holman integrator as WH, and the modified
Wisdom-Holman integrator with neural interacting Hamiltonian as WH-NIH. We implement WH-NIH
using PyTorch, and use the Adam [Kingma and Ba, 2014] optimizer to train the NIH. Given that thetraining dataset is deliberately made small (roughly 300 MB of float32 data), full training can bedone on a single GPU or a multi-core CPU within one hour.

3

Results

3.1

Two-body Problems

We test the neural integrator with a two-body problem. Since the integrator is trained with three-bodysimulations, we want to see whether the NIH understands that a two-body problem can be seem as athree-body problem without a perturbing third body. To compare with existing works, we use thesame initial condition as Greydanus et al. [2019]. We carry out the simulation for 1,000 time steps,
and WH-NIH yields a relative energy error of dE/E0 ∼ 10−9 by the end of the simulation. This result
is about 7 orders of magnitude more accurate comparing to Greydanus et al. [2019].
3.2

The TRAPPIST-1 Exoplanet System

T RAPPIST-1 [Gillon et al., 2016] is an extrasolar planetary system with seven Earth-size planets.
Because of its very special orbital architecture (compact super-Earth system with seven known planets
orbiting a spectral type M8V star), it has since attracted considerable interests concerning its origin
and potential for habitability. We test the WH-NIH integrator using the orbital parameters observedin T RAPPIST-1. The initial conditions of its orbital parameters are obtained from the exoplanets
4

0.06

0.01

0.04
0.02
0.00
0.02

b
c
d
e
f
g
h

0.04
0.06
0.050 0.0250.000 0.025 0.050
x

50

100

150

200
b
c
d
e
f
g
h

0.06
0.05
a [au]

0.06

0

0.04
0.03

0

0.015
0.010
0.000

200

200
b
c
d
e
f
g
h

0.020

0.005
100 150t [years]

150

0.025

0.01

50

100

0.030

0.02
0

50

0.000500.00048

b
c
d
e
f
g
h
0

50

dE/E0

89.90
89.85
89.80
89.75
89.70
89.65
89.60
89.55

100

150

89.8
89.7

100 150t [years]

200

50

0

50

100

150

200

100 150t [years]

200

0.000150.000100.00005

89.5
50

0

0.0002089.6

00.00042

200

b
c
d
e
f
g
h

89.9

0.000460.00044

dE/E0

0.02
0.050 0.0250.000 0.025 0.050

y

0.03

0.04

b
c
d
e
f
g
h

I [deg]

0.02

0.04

0.016
0.014
0.012
0.010
0.008
0.006
0.004
0.002

I [deg]

0.00

0.05

e

y

0.02

b
c
d
e
f
g
h

0.06

e

b
c
d
e
f
g
h

0.04

a [au]

0.06

0

50

100 150t [years]

200

0.00000

Figure 2: Evolution of the T RAPPIST-1 system for 200 years (equivalent to 2 × 105 time steps).
First row: results from a Wisdom-Holman N -body integrator with neural interacting Hamiltonian;second row: ground truth results from a traditional N -body integrator. First column: orbit plots ofJupiter and Saturn; second column: orbital semi-major axes (a) as a function of time t; third column:
orbital eccentricities e as a function of t; fourth column: orbital inclinations I (with respect to theheliocentric plane) as a function of t; last column: relative energy error (defined as the energy drift
normalized to the initial energy) as a function of t. The legends “b”, “c”, ..., “h” represent the first,
second, ..., and the seventh planet of the T RAPPIST-1 system, respectively.

database.2 The innermost planet T RAPPIST-1b has an orbital period of 1.51 days, and its orbital
semi-major axis of a = 0.0011 au is well below the minimum semi-major axes (0.4 au) in the training
dataset, so the WH-NIH integrator has never been trained to handle such an extreme system. For thisparticular system, we will use an integration time step of h = 0.365 days.
As shown in Fig. 2, we integrate the system for 200 years, which corresponds to 2 × 105 integration
time steps or 4.9 × 104 orbits of the innermost planet. Interestingly, even though the T RAPPIST-1
is well beyond the regime of the initial training data set for WH-NIH, the neural symplectic N -body
integrator handles it nicely. The evolution of the Cartesian coordinates and semi-major axis areindistinguishable from the results of the numerical integration. Concerning the t − e and t − I plots,
the WH-NIH integrator conserves the total angular momenta but fails to capture the rich dynamics
of secular resonance between multiple planets. Nevertheless, we demonstrate that it is possible tolearn the behavior of physical laws using a neural network and that the WH-NIH can apply the learned
physical laws to entirely different systems.

4

Conclusions

We propose an efficient neural N -body solver based on the splitting the symplectic Hamiltonian map,
as proposed by Wisdom and Holman [1991]: instead of approximating the full Hamiltonian H, weapproximate only a part of it and solve the other part analytically. In doing so, the error made by
the DNN is constrained. We also develop a new activation function SymmetricLog that allows theneural integrator to work directly on the raw data without the need of pre-processing. We develop acustom loss function that takes into account the conservation of energy and angular momentum. We
have demonstrated that our methods allows a DNN to solve the N -body problem for ∼ 105 steps
without diverging from the ground truth physics.
The proposed neural N -body solver is experimental, and there are some important limitations. If
an N -body system has different interacting Hamiltonian, then the neural integrator will have to beretained. In addition, the neural integrator has been unable to understand the subtlety of secularangular momentum exchanges, and therefore it is not yet capable of modeling scenarios such as theeccentric Lidov-Kozai mechanism [see, e.g., Kozai, 1962, Naoz, 2016].

2

http://exoplanet.eu

5Acknowledgments and Disclosure of Funding
It is a pleasure to thank Sam Greydanus, Javier Roa, Adrian Hamers, Vikram A. Saletore, Capsar van
Leeuwen, Valeriu Codareanu, and Veronica Saz Ulibarrena for insightful discussions. This researchis partially sponsored by SURF Open Innovation Lab (SOIL) and Intel Corporation.

References
Z. Chen, J. Zhang, M. Arjovsky, and L. Bottou. Symplectic Recurrent Neural Networks. arXiv e-prints, art.
arXiv:1909.13334, Sept. 2019.
D. Constantin Mocanu, E. Mocanu, P. Stone, P. H. Nguyen, M. Gibescu, and A. Liotta. Scalable Training of
Artificial Neural Networks with Adaptive Sparse Connectivity inspired by Network Science. arXiv e-prints,
art. arXiv:1707.04780, July 2017.
M. Cranmer, S. Greydanus, S. Hoyer, P. Battaglia, D. Spergel, and S. Ho. Lagrangian Neural Networks. arXiv
e-prints, art. arXiv:2003.04630, Mar. 2020.
Y. Desmond Zhong, B. Dey, and A. Chakraborty. Symplectic ODE-Net: Learning Hamiltonian Dynamics with
Control. arXiv e-prints, art. arXiv:1909.12077, Sept. 2019.
M. Erdmann, L. Geiger, J. Glombitza, and D. Schmidt. Generating and refining particle detector simulations
using the Wasserstein distance in adversarial networks. arXiv e-prints, art. arXiv:1802.03325, Feb. 2018.
M. Finzi, K. A. Wang, and A. G. Wilson. Simplifying Hamiltonian and Lagrangian Neural Networks via Explicit
Constraints. arXiv e-prints, art. arXiv:2010.13581, Oct. 2020.
C. Fremling, X. J. Hall, M. W. Coughlin, A. S. Dahiwale, D. A. Duev, M. J. Graham, M. M. Kasliwal, E. C.
Kool, A. A. Miller, J. D. Neill, D. A. Perley, M. Rigault, P. Rosnet, B. Rusholme, Y. Sharma, K. M. Shin,
D. L. Shupe, J. Sollerman, R. S. Walters, and S. R. Kulkarni. SNIascore: Deep Learning Classification of
Low-Resolution Supernova Spectra. arXiv e-prints, art. arXiv:2104.12980, Apr. 2021.
M. Gillon, E. Jehin, S. M. Lederer, L. Delrez, J. de Wit, A. Burdanov, V. Van Grootel, A. J. Burgasser, A. H. M. J.
Triaud, C. Opitom, B.-O. Demory, D. K. Sahu, D. Bardalez Gagliuffi, P. Magain, and D. Queloz. Temperate
Earth-sized planets transiting a nearby ultracool dwarf star. Nature, 533(7602):221–224, May 2016. doi:
10.1038/nature17448.
S. Greydanus, M. Dzamba, and J. Yosinski. Hamiltonian Neural Networks. arXiv e-prints, art. arXiv:1906.01563,
June 2019.
D. P. Kingma and J. Ba. Adam: A Method for Stochastic Optimization. arXiv e-prints, art. arXiv:1412.6980,
Dec. 2014.
Y. Kozai. Secular perturbations of asteroids with high inclination and eccentricity. AJ, 67:591–598, Nov. 1962.
doi: 10.1086/108790.
S. Naoz. The Eccentric Kozai-Lidov Effect and Its Applications. ARA&A, 54:441–489, Sept. 2016. doi:
10.1146/annurev-astro-081915-023315.
J. Pieterse and D. Constantin Mocanu. Evolving and Understanding Sparse Deep Neural Networks using Cosine
Similarity. arXiv e-prints, art. arXiv:1903.07138, Mar. 2019.
S. K. N. Portillo, J. K. Parejko, J. R. Vergara, and A. J. Connolly. Dimensionality Reduction of SDSS Spectrawith Variational Autoencoders. AJ, 160(1):45, July 2020. doi: 10.3847/1538-3881/ab9644.
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics Informed Deep Learning (Part I): Data-driven Solutionsof Nonlinear Partial Differential Equations. arXiv e-prints, art. arXiv:1711.10561, Nov. 2017a.
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics Informed Deep Learning (Part II): Data-driven Discovery
of Nonlinear Partial Differential Equations. arXiv e-prints, art. arXiv:1711.10566, Nov. 2017b.
M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-informed neural networks: A deep learning frameworkfor solving forward and inverse problems involving nonlinear partial differential equations. Journal ofComputational Physics, 378:686–707, 2019.
I. Reis, D. Poznanski, D. Baron, G. Zasowski, and S. Shahaf. Detecting outliers and learning complex structureswith large spectroscopic surveys - a case study with APOGEE stars. MNRAS, 476(2):2117–2136, May 2018.
doi: 10.1093/mnras/sty348.

6

A. Sanchez-Gonzalez, V. Bapst, K. Cranmer, and P. Battaglia. Hamiltonian Graph Networks with ODE
Integrators. arXiv e-prints, art. arXiv:1909.12790, Sept. 2019.
P. Toth, D. Jimenez Rezende, A. Jaegle, S. Racanière, A. Botev, and I. Higgins. Hamiltonian Generative
Networks. arXiv e-prints, art. arXiv:1909.13789, Sept. 2019.
J. Wisdom and M. Holman. Symplectic maps for the N-body problem. AJ, 102:1528–1538, Oct. 1991. doi:
10.1086/115978.
J. Zamudio-Fernandez, A. Okan, F. Villaescusa-Navarro, S. Bilaloglu, A. Derin Cengiz, S. He, L. Perreault
Levasseur, and S. Ho. HIGAN: Cosmic Neutral Hydrogen with Generative Adversarial Networks. arXiv
e-prints, art. arXiv:1904.12846, Apr. 2019.

A

Appendix

A.1

Initial Conditions of the Simulations

To generate training data for the NIH, we carry out 50 three-body simulations with random initialconditions. For simplicity, we define the central body as a 1M 3 star. The masses of the two orbiters
vary from 1 M (two equal-mass stars orbit around each other, also known as a binary stellar system),10−3 M (a Jupiter-mass planet orbiting around the Sun), 10−6 M (an Earth-mass planet orbiting
around the Sun) to 10−8 M (a moon-mass object orbiting around the Sun). The orbital semi-majoraxes ranges from 0.4 au4 (comparable to Mercury’s orbit) to 100 au (comparable to 3 times the sizeof Neptune’s orbit). For convenience, we set G = 4π 2 as the value of the universal gravitational
constant,5 unless otherwise noted.

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [Yes] See the discussion in eachsubsection in Section 3.
(c) Did you discuss any potential negative societal impacts of your work? [N/A](d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section 2
for the assumptions of the underlying theoretical framework.
(b) Did you include complete proofs of all theoretical results? [Yes] See Section 2 for thederivation of the equations.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] The sourcecode and data analysis scripts will be hosted on a GitHub repository.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they
were chosen)? [Yes] See Section 2.
(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [No] No, the result presentation is example-based, but theneural network predictions are compared directly with the ground truth. See Fig 2.
(d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [Yes] See Section 2, where it states thatthe training can be done on a single GPU or CPU within 1 hour.
1 M ≈ 1.9985 × 1030 kg is the mass of the Sun in our Solar System.
1 au ≈ 1.496 × 108 km is the time-averaged distance from the Sun to the Earth.
5
When G = 4π 2 the corresponding mass unit is M , length unit is au, and time unit is years.

347

4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [N/A](b) Did you mention the license of the assets? [N/A](c) Did you include any new assets either in the supplemental material or as a URL? [N/A](d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]8