The JPL Planetary and Lunar Ephemerides DE440 and DE441

The planetary and lunar ephemerides called DE440 and DE441 have been generated by fitting numerically integrated orbits to ground-based and space-based observations. Compared to the previous general-purpose ephemerides DE430, seven years of new data have been added to compute DE440 and DE441, with improved dynamical models and data calibration. The orbit of Jupiter has improved substantially by fitting to the Juno radio range and Very Long Baseline Array (VLBA) data of the Juno spacecraft. The orbit of Saturn has been improved by radio range and VLBA data of the Cassini spacecraft, with improved estimation of the spacecraft orbit. The orbit of Pluto has been improved from use of stellar occultation data reduced against the Gaia star catalog. The ephemerides DE440 and DE441 are fit to the same data set, but DE441 assumes no damping between the lunar liquid core and the solid mantle, which avoids a divergence when integrated backward in time. Therefore, DE441 is less accurate than DE440 for the current century, but covers a much longer duration of years −13,200 to +17,191, compared to DE440 covering years 1550–2650.


Introduction
Modern-day planetary ephemerides are computed by fitting numerically integrated orbits to various types of ground-based and space-based observations (Folkner et al. 2014;Pitjeva & Pitjev 2018;Fienga et al. 2020). The Jet Propulsion Laboratoryʼs (JPL) planetary and lunar ephemerides Development Ephemeris (DE) series includes the positions of the Sun, the barycenters of eight planetary systems, the Moon, the Pluto system barycenter, and lunar libration angles, as well as their associated velocities. The high-precision orbits and lunar rotations around the three axes have a wide range of practical and fundamental applications (Thornton & Border 2003;Park et al. 2020b;Vallisneri et al. 2020;U.S. Nautical Almanac Office & Her Majestyʼs Nautical Almanac Office 2018). Without an update, however, the errors in orbits grow for several reasons. For Jupiter (Juno 2016-2020 < 12 yr period), Saturn (Cassini 2004-2018, <30 yr period), Uranus, Neptune, and Pluto the high-quality data are available for less than one orbit. For Mars, the range data is of high quality, but the mainbelt asteroid masses are a limitation partly due to the large number of asteroids and partly due to long-period perturbations (Folkner et al. 2014). In the future, these errors grow nonlinearly with time.
The planetary and lunar ephemerides DE440 replaces DE430 released in 2014 (Folkner et al. 2014) and its precursors. Since the DE430 release, several interim ephemerides have been released. Each interim DE file was for a specific flight project, which has been tuned for the flight projectʼs target body. For example, the last release was DE438 in 2018 for the Juno mission (Bolton et al. 2017). DE440 has updated all bodies using all available data, including the Moon, which has not been updated since DE430.
Since DE430, several updates have been made to the dynamical model used to integrate DE440. Perturbations from 30 individual Kuiper belt objects (KBOs) and a circular ring representing the rest of the Kuiper belt, modeled as 36 point masses with an equal mass located at 44 au, have been added to the model (Pitjeva & Pitjev 2018). The Lense-Thirring (LT) effect from the Sunʼs angular momentum has also been added (Park et al. 2017). For the orientation of Earth, the Vondrak precession model has been used (Vondrak et al. 2011), which, according to Vondrak et al. (2009), is more accurate for integrations beyond ±1000 years than the Lieske precession model (Lieske 1979) used for DE430. For the Moon, the effect of geodetic precession on lunar librations has been added as well as the solar radiation pressure force on the Earth-Moon system orbits.
Compared to DE430, new data spanning over about 7 years have been added to compute DE440. The shapes of Mercury, Venus, and Mars orbits are determined mainly by the radio range data of the MErcury Surface, Space ENvironment, GEochemistry, and Ranging (MESSENGER), Venus Express, and Mars-orbiting spacecraft, respectively. The orientations of inner planet orbits are tied to the International Celestial Reference Frame (ICRF) via Very Long Baseline Interferometry (VLBI) of Mars-orbiting spacecraft (Folkner & Border 2015;Park et al. 2015). The orbit of the Moon is determined from laser ranging to lunar retroreflectors. The orbit accuracy of Jupiter has improved substantially by fitting the orbit to the radio range and Very Long Baseline Array (VLBA) data of the Juno spacecraft. The orbit of Saturn is determined by the radio range and VLBA data of the Cassini spacecraft, with improved spacecraft orbits used for processing the radio range data. The orbits of Uranus and Neptune are determined by astrometry and radio range measurements to the Voyager flybys. The orbit of Pluto is now mainly determined by stellar occultations reduced against the Gaia star catalog (Gaia Collaboration et al. 2018;Desmars et al. 2019).
For the Moon, viscous damping between the liquid core and the solid mantle are observed in the lunar laser ranging (LLR) data. This implies an excitation of the relative motion of a lunar core and mantel in the past, possibly due to a spin/orbit resonance that occurred in geologically recent times (Rambaux & Williams 2011). Both DE440 and DE441 have been fit to the same data set, but DE441 assumed no damping between the lunar liquid core and the solid mantle. In this way, a divergence can be avoided when integrated backward for thousands of years. As a consequence, DE441 is less accurate than DE440 for the current century, but more accurate over longer past time spans. DE441 covers years −13,200 to +17,191, compared to DE440 covering years 1550-2650. The difference in the orbits of the planets between DE440 and DE441 are less than 1 m over the 1100 years of the DE440 time span. The difference in the orbit of the Moon between DE440 and DE441 is less than 2 m during the time span of the LLR data, i.e., 1970-2020, but increased over a longer time span, especially in the along-track direction (i.e., velocity direction). Figure 1 shows the difference in the lunar orbit relative to Earth between DE440 and DE441 over 200 years (i.e., DE441 minus DE440). In general, DE440 is recommended for analyzing modern data while DE441 is recommended for analyzing historical data earlier than the modern range data.

Inertial Reference Frame
The inertial coordinate frame of the planetary and lunar ephemerides is connected to the International Celestial Reference System (ICRS). The current ICRS realization is achieved by VLBI measurements of the positions of extragalactic radio sources (i.e., quasars) defined in the Third Realization of the International Celestial References Frame (ICRF3; Charlot et al. 2020), which is adopted by the International Astronomical Union (IAU). The orbits of the inner planets are tied to ICRF3 via VLBI measurements of Mars-orbiting spacecraft ) with respect to quasars with positions known in the ICRF. Overall, the orientations of inner planet orbits are aligned with ICRF3 with an average accuracy of about 0.2 mas (Folkner & Border 2015;Folkner et al. 2014;Park et al. 2015). The orbits of Jupiter and Saturn are tied to ICRF3 via VLBA measurements of Juno and Cassini spacecraft (Jones et al. 2020), respectively.

Solar System Barycenter
The solar system barycenter (SSBC) is defined as (Estabrook 1971) å å Here, the summation is over all bodies with finite mass, r i , is the position of body i, and m i * is defined as where GM i is the mass parameter of body i, c is the speed of light, v i is the barycentric speed of body i, and =r r r ij i j | | is the distance between bodies i and j.
For DE440, the bodies used for computing the SSBC were the Sun, barycenters of eight planetary systems, the Pluto system barycenter, 343 asteroids, 30 KBOs, and a KBO ring representing the main Kuiper belt. The 343 asteroids were the same set of asteroids used in DE430, which consist of ∼90% of the total asteroid-belt mass. The mass of the 30 largest known KBOs were from Pitjeva & Pitjev (2018). The circular KBO ring was modeled as 36 point masses with equal mass located in the ecliptic plane with a semimajor axis of 44 au, with the ring mass estimated. Figure 2 shows the motion of the SSBC relative to the Sun for 100 years (2000-2100), which is sometimes called the solar inertial motion (SIM). Compared to DE430, SSBC has shifted by ∼100 km, which is mainly due to the addition of KBOs. It is important to note that, to the first order, Earth orbits around the Sun, not around the SSBC. This point is reflected in Figure 3, which shows the time history of the closest (e.g., perihelion) and farthest (e.g., aphelion) points of the Earth-Moon barycenter (EMB) relative to the Sun and the SSBC. The near-constant distance of the perihelion and aphelion of the EMB relative to the Sun indicates that SIM does not affect the orbit of Earth relative to the Sun.

Ephemeris Coordinate Time
JPLʼs DE series are integrated using the barycentric dynamical time (TDB), which is defined relative to the barycentric coordinate time (TCB; Petit & Luzum 2009). All of the data used to compute DE440 and DE441 had the intrinsic time tag in the coordinated universal time (UTC), which differs from the international atomic time (TAI) by leap seconds (i.e., TAI = UTC + leap seconds). In order to process these UTC-tagged measurements, the conversion from UTC to TDB would be needed (Soffel et al. 2003;Petit & Luzum 2009). Once the TAI time is computed, 32.184 s are added to compute the terrestrial time (TT; i.e., TT = TAI + 32.184 s). The where L G = 6.969290134 × 10 −10 defines the rate of TT with respect to geocentric coordinate time (TCG), L B = 1.550519768 × 10 −8 defines the rate of TDB with respect to TCB, T 0 is 2443144.5003725 Julian days, = -´-TDB 65.5 0 10 86400 6 days, v E is the velocity of the Earth, v E is the velocity vector of the Earth, r E is the position vector of the Earth, and r S is the position vector of a measurement station. The potential term w 0E is defined as with the summation over all bodies other than Earth. The potential due to external oblate figures w LE is defined as where J 2e is the unnormalized second-degree gravitational zonal harmonic of the Sun, R e is the solar radius, and j E,e is the heliocentric ecliptic latitude of Earth. The term w iE is defined as  where the summation is over all bodies other than Earth. Lastly, Δ E is defined as where a i is the acceleration of body i and the summation is over all bodies.

Orientation of the Moon
LLR measures the round-trip light time of a laser pulse between an Earth LLR station and a retroreflector on the Moon. Thus, LLR data are not only sensitive to where the Moon is, but they are also sensitive to its orientation.
The orientation of lunar exterior (mantle and crust, hereafter referred to as the mantle) is defined by the principal axes (PAs) of the undistorted lunar mantle, and thus its moment of inertia matrix is diagonal. The directions of the PAs are taken from analyses of the Gravity Recovery and Interior Laboratory (GRAIL) data (Konopliv et al. 2013;Lemoine et al. 2013). The Euler angles that define the rotation from the PA frame to the inertial ICRF3 frame are: f m , the angle from the X-axis of the inertial frame along the XY plane to the intersection of the mantle equator; θ m , the inclination of the mantle equator from the inertial XY plane; and ψ m , the longitude from the intersection of the inertial XY plane with the mantle equator along the mantle equator to the prime meridian. The rotation from the lunar PA frame (i.e., lunar mantle frame) to ICRF3 is given as The Euler angles f m , θ m , and ψ m are also known as lunar libration angles, stored in the DE440 and DE441 files. The rotation matrices use the right-hand rule and are defined as where ω m,x , ω m,y , and ω m,z represent the components of the lunar mantle angular velocity ω m expressed in the mantle frame. The time derivatives of ω m are given in Section 4. The lunar orientation model includes a fluid core. The orientation of the core with respect to the ICRF is represented by the Euler angles f c , θ c , and ψ c , which are also numerically integrated. Since the shape of the core/mantle boundary is modeled as fixed to the frame of the mantle, it is more convenient to express the core angular velocities with respect to the mantle frame. The time derivatives of the core Euler angles are then given by where the core angular velocity ω c is related to the angular velocity w c † in a frame defined by the intersection of the core equator with the inertial XY plane by The time derivatives of ω c are given in Section 4. Most of lunar cartographic products are defined relative to the DE421 mean-Earth/mean-rotation (MER) frame. The MER frame is defined by the X-axis pointing toward the mean-Earth direction and the Z-axis pointing toward the mean-rotation axis direction. The rotation from the DE440 PA frame to the DE421 MER frame is estimated by comparing the coordinates of the lunar retroreflectors estimated in the DE440 PA frame and the retroreflector coordinates in the DE421 MER frame, which yield Table 1 shows the five lunar retroreflector positions in the DE440 PA frame and the corresponding lunar retroreflector positions in the DE421 MER frame.

Orientation of Earth
Only the long-term change of the Earth orientation is modeled in the ephemeris integration. The Earth orientation model used for DE440 and DE441 is based on a long-term precession model (Vondrak et al. 2011) and a modified The Earth pole unit vector in the inertial frame, p E , can be computed by the following steps.
First, the mean longitude of the ascending node of the lunar orbit measured on the ecliptic plane from the mean equinox of date is computed by where the mean obliquity ē is given by ( ) The pole unit vector in the inertial frame p E is computed by precessing the pole of date to inertial coordinates using the long-term precession model (Vondrak et al. 2011) plus an estimated frame offset in x and y rotations, where Φ x and Φ y are the estimated offsets of the EME2000 pole from the ICRF pole and the precession matrix  V is given by =´´ n k n k n n k n n k n n where k is the ecliptic pole vector and n is the mean equatorial pole vector, derived from polynomial fits to a numerically integrated long-term orientation of Earth (Vondrak et al. 2011). The |·| operator represents the norm of a vector.

Translational Equations of Motion
This section presents the dynamical models of the planetary and lunar ephemerides, including changes and updates made compared to DE430. Some materials from DE430 (Folkner et al. 2014) are repeated so that this paper can be self-contained and the results can be reproduced.

Point-mass Acceleration
The point-mass interaction between planetary bodies is governed by the parameterized post-Newtonian (PPN) formulation (Will & Nordtvedt 1972;Moyer 2003) where the summations are over all bodies, and β and γ are the Eddington-Robertson-Schiff parameters representing the measure of nonlinearity in the superposition law for gravity and the amount of space curvature produced by a unit rest mass, respectively, and are constrained to unity as predicted by the general theory of relativity (GTR). DE440 integrated the same set of bodies (i.e., Sun, barycenter of eight planets, the Moon, Pluto barycenter, and 343 asteroids) used in DE430, but also included perturbations from 30 KBOs and a KBO ring discussed in Section 2.2. The key mass parameters used in DE440 are shown in Table 2, and all other relevant parameters are given in the comment blocks of the DE440 and DE441 files.

Point-mass Interaction with Extended Bodies
Nonspherical gravitational interaction has been modeled using a spherical harmonic expansion. For Earth, the interaction of the zonal harmonics up to the fifth degree and the point masses of the Moon, Sun, Mercury, Venus, Mars, Jupiter, and Saturn have been modeled. For the Moon, the interaction of a  Park et al. 2014) degree and order 6 gravity field and the point masses of Earth, Sun, Mercury, Venus, Mars, Jupiter, and Saturn have been modeled. For the Sun, the interaction of the second-degree zonal harmonic with all other bodies has been modeled. The acceleration due to an extended body can be represented as where the ξηζ coordinate system is defined such that the ξ-axis is defined outward from the extended body to the point mass, the ξζ-plane contains the figure spin-pole of the extended body, and the η-axis completes the triad. Here, r is the distance between the two bodies; n 1 and n 2 represent the maximum degrees of the zonal and nonzonal spherical harmonic coefficients, respectively; P n and P nm represent the unnormalized degree-n Legendre polynomial and associated Legendre function with degree-n and order-m, respectively; ¢ P n and ¢ P nm represent the derivative of P n and P nm with respect to j sin , respectively; J n represents the degree-n zonal harmonic coefficient; C nm and S nm represent the nonzonal spherical harmonic coefficients for the extended body; R represents the reference radius of the extended body; and λ and j represent the longitude and latitude of the point mass in the extendedbody fixed coordinate system. Once the accelerations are computed in the body fixed frame, they are transformed into the inertial frame for integration.
There is also an interaction between the figure of an extended body and a point mass, often called the indirect acceleration. Given a i,figi−pmj , which denotes the acceleration of the extended body i interacting with the point-mass external body j expressed in an inertial frame, the corresponding indirect acceleration of the point mass, a j,figi−pmj , is For the Moon, the second-degree spherical harmonic coefficients vary with time due to distortion by tides and rotation. These coefficients are computed from the moment of inertia tensors as a function of time and the detailed description is given in Section 4. The second-degree spherical harmonic coefficients are given by where I ij,T represent the elements of the total lunar moment of inertia matrix (see Section 4), m M is the lunar mass, and R M is the lunar radius.

Acceleration of the Moon from Earth Tides
The lunar orbit is affected by the tides raised on Earth by the Sun and Moon. The tidal distortion of Earth can be modeled using the second-degree gravitational Love numbers, k 2j,E , where the order j is 0, 1, and 2 for long-period, diurnal, and semidiurnal responses, respectively. See Section 2.2 in Williams & Boggs (2016) for more information.
A time-delay tidal model has been applied to account for the tidal dissipation. The distorted response of Earth is delayed with respect to the tide-raising forces from the Moon or Sun. The appropriate time delay depends on the period of each tidal component. Consequently, different time delays have been employed for each order j. To allow for time delays shifting across the diurnal and semidiurnal frequency bands, separate time delays are associated with Earthʼs rotation and the lunar orbit.
The acceleration of the Moon due to the Earth tides is evaluated separately for the tides raised by the Sun and the tides raised by the Moon. The Earth tides depend on the position of the tide-raising body with respect to Earth, r T , where T can denote either the Sun or the Moon. The position of the tide-raising body is evaluated at an earlier time t -¢ t j , where t ¢ j denotes the orbital time lag, for long-period ( j = 0), diurnal ( j = 1), and semi-diurnal ( j = 2) responses. The rotational distortion of Earth is delayed by a rotational time lag τ j , so that the distortion leads the direction to the tide-raising body by an angle qt j  , where q  is the rotation rate of Earth. The long-period zonal tides ( j = 0) do not depend on the rotation of Earth, so τ j = 0. The acceleration of the Moon due to the distorted Earth depends on the position of the Moon with respect to Earth (r) and on the modified position vector for the tide-raising body (r j *) that is given for each order j by where  E rotates the time-delayed position of the tide-raising body with respect to Earth from the inertial frame to the Earth fixed frame and  z here means a right-handed rotation of the vector t -¢ r t T j ( )by the angle qt j  about Earthʼs rotation axis.
For evaluation of the acceleration of the Moon, the vectors r (Moon with respect to Earth) and r j * (time-delayed position of a tide-raising body with respect to Earth) are expressed in cylindrical coordinates with the Z-axis perpendicular to Earthʼs equator so that r = ρ + z and r = + r z j j j The acceleration of the Moon due to the tide raised on Earth by each tide-raising body (Sun or Moon), a M,tide is given by where m T is the mass of the tide-raising body. The tidal acceleration due to tidal dissipation is implicit in the above acceleration. Tides raised on Earth by the Moon do not influence the motion of the EMB. The effect of Sun-raised tides on the barycentric motion is not considered.
The tidal bulge leads the Moon and its gravitational attraction accelerates the Moon forward and retards Earthʼs spin. Energy and angular momentum are transferred from Earthʼs rotation to the lunar orbit. Consequently, the Moon moves away from Earth, the lunar orbit period lengthens, and Earthʼs day becomes longer. Some energy is dissipated in Earth rather than being transferred to the orbit.

LT Acceleration
In DE440, the acceleration of each body other than the Sun due to the gravito-magnetic effect of GTR, also known as the LT effect, has been implemented (Moyer 2003;Park et al. 2017): where the LT angular velocity vector Ω I is given by where G is the universal gravitational constant, γ is the Eddington-Robertson-Schiff parameter from Section 3.1, c is the speed of light, r i is position of the body with respect to the Sun, M e is the Sunʼs mass, C e is the Sunʼs polar moment of inertia divided by M R , R e is the Sunʼs equatorial radius (696,000 km), ω e the is Sunʼs rotation rate (14.1844 deg/day), and p  is the unit spin-pole direction of the Sun (R.A. of 286°. 13 and decl. of 63°.87; Archinal et al. 2018).
The LT effect is small, but it is important for fitting the MESSENGER range data (Park et al. 2017). Overall, the LT effect causes Mercuryʼs perihelion to precess at the rate of about −0 0020/Julian century, which is about 7% of the precession caused by the solar oblateness.

Solar Radiation Pressure
Photons carry energy and momentum so there is a very small force directed away from the Sun. Like Newtonian gravity, solar radiation pressure depends on the inverse-square distance of the Sun. A simple solar radiation pressure model has been implemented for Earth and the Moon, with accelerations given by where the acceleration due to solar radiation pressure is a very small fraction of gravitational acceleration (Vokrouhlický 1997), i.e., ε E,srp = 2 × 10 −14 for Earth and ε M,srp = 1.44 × 10 −13 for the Moon. Here, r SE is the inertial Sun-to-Earth position vector and r SM is the inertial Sun-to-Moon position vector.

Rotational Dynamics of the Moon
The Moon is modeled as an anelastic mantle with a liquid core. The orientations of the core and mantle are numerically integrated for the core and mantle angular velocities. The angular momentum vectors of the mantle and core are the product of the angular velocities and the moments of inertia. The angular momentum vectors change with time due to torques and due to the distortion of the mantle.

Rate of Change of Lunar Angular Velocities
In a rotating system, the change in angular velocity ω is related to torques N by where I is the moment of inertia tensor. The second term on the right-hand side puts the time derivative into the rotating system. The total lunar moment of inertia I T , which is the sum of the moment of inertia of the mantle I m and the moment of inertia of the core I c , is proportional to the mass m M times the square of the radius R M . Because the fractional uncertainty in the constant of gravitation G is much larger than that for the lunar mass parameter Gm M , Equation (42) is evaluated in the integration with both sides multiplied by G.
The components of vectors can be given in the inertial frame, mantle frame, or other frames. Since the moment of inertia matrices are nearly diagonal in the mantle frame, there is great convenience to inverting matrices and performing the matrix multiplications in the mantle frame. The resulting vector components can then be rotated to other frames if desired. The

 
where I m is the lunar mantle moment of inertia matrix, N M,figM−pmj is the torque on the lunar mantle due to the pointmass body j, N M,figM−figE is the torque on the lunar mantle due to the interaction between the extended figure of the Moon and the extended figure of Earth, ν gp is the angular rate due to geodetic precession, and N cmb is the torque due to the interaction between lunar mantle and core. The geodetic precession rate is defined as where L is the rotation matrix from the mean J2000 frame to the lunar body fixed (i.e., selenographic) frame, v EM is the inertial Earth-to-Moon velocity vector, r EM is the inertial Earthto-Moon position vector, and r SM is the inertial Sun-to-Moon position vector. The fluid core is assumed to be rotating like a solid and constrained by the shape of the core-mantle boundary at the interior of the mantle, with the moment of inertia constant in the frame of the mantle. The time derivative of the angular velocity of the core expressed in the mantle frame is given by where I c is the lunar core moment of inertia matrix.

Lunar Moments of Inertia
In the mantle frame, the undistorted moment of inertia of the mantle and the moment of inertia of the core are diagonal. The undistorted total moment of inertia I T  is given by where m M is the mass of the Moon, R M is the reference radius of the Moon, J M

2,
is the second-degree zonal harmonic of the undistorted Moon, and β L and γ L are ratios of the undistorted moments of inertia given by The undistorted total moment of inertia and the second-degree zonal harmonic of the undistorted Moon are not the same as the mean values since the tidal distortions have non-zero averages. We assume that the lunar rotation aligns an oblate coremantle boundary with the equator. Then the moment of inertia of the core I c is given by where α c = C c /C T is the ratio of the core polar moment of inertia to the undistorted total polar moment of inertia and f c is the core oblateness. Distortion of the core moment of inertia is not considered. The undistorted moment of inertia of the mantle is the difference between the undistorted total moment of inertia and the core moment of inertia: The moment of inertia of the mantle varies with time due to tidal distortion by Earth and spin distortion, where the position vector of Earth relative to Moon, r, and the angular velocity of the mantle, ω m , are evaluated at time  t − τ m ; k 2,M is the lunar potential Love number; m E is the mass of Earth; R M is the reference radius of the Moon; r is the Earth-Moon distance; x, y, and z are the components of the position of Earth relative to the Moon referred to the mantle frame; ω m,x , ω m,y , ω m,z are the components of ω m in the mantle frame; and n is the lunar mean motion. The rate of change of the mantleʼs moment of inertia is given by where p E  is the direction vector of Earthʼs pole and r EM  is the direction vector of Earth from the Moon, I M is the lunar moment of inertia tensor, R E is the reference radius of Earth, and j M,E is defined by . The torque on the mantle due to the interaction between the core and mantle is evaluated in the mantle frame and is given by where z m  is a unit vector in the mantle frame aligned with the polar axis. Parameter k v is for the viscous interaction at the core-mantle boundary. The torque on the core is the negative of the torque on the mantle.

Observational Data Used for Computing DE440 and DE441
The observations that have been used to compute DE440 and DE441 are summarized in Tables 3-5 for each body. An LLR observation measures the round-trip light time from an LLR station on Earth to a retroreflector on the Moon. There are five retroreflectors on the Moon: the Apollo 11, 14, and 15 landing sites and the Lunokhod 2 and 1 rovers. The LLR measurements started in 1970 following the first landings of astronauts and continue to the present. The LLR residuals can be expressed as one-way range residuals, i.e., one-way residual = t t c 2 measured computed ( ) . The LLR measurement accuracy has improved with time as technology for producing short-duration high-energy laser pulses and timing measurements has advanced. Spacecraft measurements are based on the Deep Space Network (DSN) radio range, Doppler, and VLBI measurements. For spacecraft in orbit about the planet, the Doppler measurements are typically used to estimate the position of the spacecraft with respect to the planetʼs center of mass and then range and VLBI measurements are used to estimate the orbit of the planet. For spacecraft flying by a planet, the range, Doppler, and VLBI data, as available, are used to estimate both the trajectory of the spacecraft and a 3D position of the planet, given as range, R.A., and decl.
Range measurements to spacecraft are usually made at regular intervals during a tracking pass, typically every 10 minutes, while Doppler measurements are made more frequently, typically every minute. Both range and Doppler measurements are based on the measurement of the phase of a radio signal, with the carrier signal used for Doppler and a ranging modulation signal used for range. Since the carrier signal is at a much higher frequency and usually has much higher signal strength, Doppler measurements change in range much more accurately than the range measurements. Because of the shorter wavelength associated with the higher frequency, the integer number of carrier wavelengths cannot be resolved, so Doppler measurements do not allow for an estimation of absolute range. Range measurements are more accurate measurements of round-trip light time. For plotting residuals, one-way residuals are used, which are essentially the same as the LLR residuals, i.e., one-way range = t t c 2 measured computed ( ) . The range measurement accuracy is often limited by a calibration of the signal path delay in the tracking station prior to each tracking pass. Since this calibration error is common to all range measurements in the tracking pass, there is only one  1985-1993165 La Palma 1984-19971030Tokyo 1986-198944 Washington 1926-19931783 statistically independent range point per pass. Therefore, only one range point per tracking pass was used in the data reduction, and the number of range measurements per spacecraft in Tables 3-5 reflect this.
Spacecraft VLBI measurements are usually made using two widely separated tracking stations. The measurements are made using a modulation on the carrier signal (delta-differential oneway range) and give one component of the direction to the    Astrometric measurements record the direction to the planet, namely, R.A. and decl., based on imaging relative to a star field. The accuracy of the star catalog is often the largest source of measurement error. The CCD type indicates more modern observations using electronic detectors, generally referred to star catalogs based on the Hipparcos mission launched in 1991 that are referred to the ICRF2 through the estimation of the positions of radio stars using VLBI. Older measurements were taken using photographic plates or transit methods, often referred to older star catalogs, though corrected to the Hipparcos catalog in some fashion. Barnard measured the angular separation between the outer planets and some of their satellites relative to angularly nearby stars at Yerkes Observatory. The positions of those stars are taken from modern star   catalogs, with accuracies limited by the knowledge of stellar proper motion. Transit observations cover a longer time span than the more modern spacecraft and astrometric measurements. Since the measurement noise is relatively large for the transit measurements, they do not contribute significantly to the ephemeris solution. The transit measurements are included mainly for historical comparison.
Occultation measurements of Pluto are included here, where the R.A. and decl. are determined from the timed disappearance and reappearance of a star occulted by Pluto (Desmars et al. 2019).

Residuals
The orbit of the Moon is determined from the LLR data, and Figure 4 shows the one-way residuals of the LLR data over about 50 years. In the last few years, data were available from the Observatoire de la Côte d'Azur, Apache Point, Matera, and Wettzell sites. The rms residual of the early LLR data (∼1970-1980) is about 20 cm while the rms residual of the recent LLR data is about 1.3 cm.
The orbit of Mercury is mainly determined from the MESSENGER radio range data, and Figure 5 shows the oneway range residuals over about 4 years. The rms residual of the MESSENGER ranges is about 0.7 m.
The orbit of Venus is mainly determined from the Venus Express radio range data, and Figure 6 shows the one-way range residuals over about 6.5 years. The rms residual of the Venus Express ranges is about 8 m.
The shape of Mars' orbit is mainly constrained by the radio range data of Mars orbiters, and Figure 7 shows the one-way range residuals of Mars Express (MEX), Mars Global Surveyor (MGS), Mars Odyssey (ODY), and Mars Reconnaissance Orbiter (MRO). The rms residual of the MEX ranges is about 2 m and the rms residuals of the MGS, ODY, and MRO ranges are about 0.7 m. The orientation of Mars' orbit is mainly constrained by VLBI measurements of the Mars orbiters, and Figure 8 shows the VLBI residuals of the Goldstone-Madrid and Goldstone-Canberra baselines. The rms residual of the Goldstone-Madrid baseline is about 0.25 mas and the rms residual of the Goldstone-Canberra baseline is about 0.18 mas.
The shape of Jupiterʼs orbit is mainly constrained by the Juno radio range data, and Figure 9 shows the one-way range residuals over about 4 years. The rms residual of the Juno ranges is about 13 m. The orientation of Jupiterʼs orbit is mainly constrained by the VLBA data of the Juno spacecraft, and Figure 10 shows the Juno VLBA residuals over about 2 years. The rms residual of the Juno VLBA data in R.A. is about 0.4 mas and the rms residual of the Juno VLBA data in decl. is about 0.6 mas.
The shape of Saturnʼs orbit is mainly constrained by the Cassini radio range data, and Figure 11 shows the one-way range residuals over about 13 years. The rms residual of the Cassini ranges is about 3 m. The orientation of Saturnʼs orbit is mainly constrained by the VLBA data of the Cassini spacecraft, and Figure 12 shows the Cassini VLBA residuals over about 13 years. The rms residual of the Cassini VLBA data in R.A. is about 0.6 mas and the rms residual of the Cassini VLBA data in decl. is about 0.8 mas if the two obvious outliers are included, and 0.35 and 0.36 mas if they are excluded.
The orbits of Uranus and Neptune are mainly constrained by astrometry and radio range measurements to the Voyager flybys and they are statistically consistent with DE430 (see Folkner et al. 2014).
The orbit of Pluto is mainly determined by astrometry, and Figure 13 shows the residuals of the stellar occultation measurements from Desmars et al. (2019). The rms residuals of the stellar occultations are about 8 mas in R.A. and about 11 mas in decl.

Conclusions
This paper presents JPLʼs new general-purpose ephemerides DE440 and DE441 created by fitting numerically integrated orbits to ground-based and space-based observations. Compared to DE430, the previous general-purpose ephemerides released in 2014, new data spanning over about 7 years have been added to compute DE440. The shapes of the orbits of Mercury, Venus, and Mars are constrained mainly by the radio range data of the MESSENGER, Venus Express, and Marsorbiting spacecraft, respectively. The orientations of inner planet orbits are tied through the ICRF via VLBI of Marsorbiting spacecraft. The orbit of the Moon is primarily determined by laser ranging to lunar retroreflectors with the data arc extended through 2020 March. The orbit accuracy of Jupiter has improved substantially by fitting the orbit to the new Juno radio ranges and VLBA data. The orbit of Saturn is determined by radio ranges and VLBA data of the Cassini spacecraft, with improved calibration of the radio range data. The orbits of Uranus and Neptune are constrained by astrometry and radio range measurements to the Voyager flybys. The orbit of Pluto is constrained by stellar occultation data reduced against the Gaia star catalog. DE440 spans years 1550-2650 and DE441 spans years −13,200 to +17,191. The ephemerides of DE440 are recommended for analyzing modern data while DE441 is recommended for analyzing historic data before the DE440 time span.