The Determination of the Rotational State and Interior Structure of Venus with VERITAS

Understanding the processes that led Venus to its current state and will drive its future evolution is a major objective of the next generation of orbiters. In this work we analyze the retrieval of the spin vector, the tidal response, and the moment of inertia of Venus with VERITAS, a NASA Discovery-class mission. By simulating a systematic joint analysis of Doppler tracking data and tie points provided by the onboard synthetic aperture radar, we show that VERITAS will provide accuracies (3σ) in the estimates of the tidal Love number k 2 to 4.6 × 10−4, its tidal phase lag to 0.°05, and the moment of inertia factor to 9.8 × 10−4 (0.3% of the expected value). Applying these results to recent models of the Venus interior, we show that VERITAS will provide much-improved constraints on the interior structure of the planet.


Introduction
The most comprehensive mapping of Venus was done by the Magellan mission in the early 1990s (Saunders et al. 1992). Magellan employed a combination of data from Doppler tracking and S-band Synthetic Aperture Radar (SAR), an altimeter and radiometer to make nearly global observations of the surface of Venus (Ford & Pettengill 1992). Magellan observations led to the most accurate in situ estimate of the planet's spin axis orientation, sidereal rotation period (Davies et al. 1992; see also Campbell et al. 2019 for a summary of several other observation campaigns), gravity field, and tidal response (Konopliv et al. 1999;Konopliv & Yoder 1996). The Magellan estimates, however, proved not sufficiently precise to constrain the structure of the mantle and core. As shown in Dumoulin et al. (2017), current estimates of the tidal response do not distinguish between a liquid and solid core, and the absence of a measurement of the tidal phase lag prevents us from inferring the viscous response of the interior. Until recently, models of Venus's interior relied solely on scaling Earth's interior structure to Venus's radius (e.g., Yoder 1995; Aitta 2012). A recent direct (ground-based) measurement of the moment of inertia factor (MOIF = C/MR 2 , where C is the polar moment of inertia and M and R are the planetary mass and radius, respectively) yields 7% fractional uncertainty and provides weak constraints on the internal density profile and core size (Margot et al. 2021). Improved measurements are needed to quantify the interior structure of Venus with precision.
The Venus Emissivity, Radio science, INSAR, Topography And Spectroscopy (VERITAS) mission (Freeman & Smrekar 2015; Smrekar & the VERITAS science team 2021) is a partnership led by NASA/JPL between US scientists and engineers, with strong collaborations and contributions of the German, Italian, and French Space Agencies. On 2021 June 2 NASA definitively selected VERITAS as one of the two winners of the Discovery 2019 competition. The launch is expected in the 2028-2030 time frame.
A key scientific objective of VERITAS is understanding the links between the interior, surface, and atmospheric evolution. The determination of the tidal response, tidal phase lag, and MOIF is specifically focused on pushing forward our understanding of the Venus interior. VERITAS will carry two science instruments: VISAR (Venus Interferometric Synthetic Aperture Radar), the X-band interferometric radar (Hensley et al. 2020); and VEM (Venus Emissivity Mapper), an infrared spectroscopic mapper (Helbert et al. 2020). Data from VISAR will be combined with two-way dual X-and Ka-band Doppler tracking data provided by the onboard telecom subsystem collected for the gravity science investigation and used to improve estimates of the Love number k 2 , the tidal phase lag d k 2 , and the MOIF in order to constrain the structure of the Venus interior.
Arriving at Venus after a 6-month cruise, VERITAS will begin an 11-month aerobraking phase, paused after 5 months for 5 months of VEM science observations, before continuing to its final nearly circular polar orbit (180 × 255 km in altitude, ∼85°.4 inclination, period ∼1.5 hr). VERITAS plans to operate for 4 Venus sidereal days (or four cycles, 243 Earth days each), providing a nearly global coverage of the planet for all its investigations (gravity science, VISAR, and VEM).
The goal of this work is to simulate the operational scenario of VERITAS's gravity experiment to assess the accuracy in the estimate of k 2 , d k 2 , and MOIF. In our work, alongside a typical orbit determination solution employing Earth Doppler tracking data, we explore a novel approach based on the systematic inclusion of VISAR landmark features observations, or tie points, to tighten the determination of the rotational state of the planet. The simulations presented here, moreover, represent the first assessment of the impact of recent advancements in the understanding of Venus's atmospheric dynamics, namely, the gravitational signature of atmospheric tides (Bills et al. 2020) and short-term sidereal period oscillations of the solid planet due to the transfer of atmospheric angular momentum (Margot et al. 2021).
In Section 2 we describe the concept and the assumptions used in our simulations for both Doppler and radar measurements (Sections 2.1 and 2.2, respectively) and their combination (Section 2.3). In Section 3 we discuss the simulation setup and observational scenario. In Section 4 we present and discuss the results of the simulations. Section 5 follows with concluding remarks. This manuscript is complemented by six appendices.

Data and Methods
It is well known that the sole knowledge of the gravitational field is not enough to infer the moments of inertia of a planet, which provide crucial constraints on its interior structure. To constrain the inertia tensor of a body, the gravity field information must be complemented by measurements of the rotational state. Precise Doppler tracking data, the primary observable quantity for gravity field recovery, are quite sensitive to the rotational state of the planet, but the attainable accuracy can be improved by augmenting the analysis with surface feature tracking. The latter provides direct observations of the rotational motion of the planet by measuring the inertial displacement of physical features located on the planet's surface. In this work, we make use of a novel approach (building on the technique proposed by Davies et al. 1992;Chodas et al. 1992) to combine Earth-spacecraft Doppler tracking data and repeated surface landmark observations (tie points) provided by the onboard interferometric SAR.

Spacecraft Doppler Tracking
Doppler measurements are the primary observables for reconstructing the orbit of the spacecraft and recovering the gravity field of a planet. These measurements are collected by recording the Doppler shift of a radio signal sent from the ground station to the spacecraft, which then coherently retransmits it back to Earth by means of an onboard transponder (two-way configuration). VERITAS's Doppler tracking system, with its heritage from ESA's BepiColombo (Iess et al. 2009(Iess et al. , 2021, is able to establish two simultaneous coherent radio links in the X (7.1-8.5 GHz) and  GHz) bands and to provide measurements of the range rate of the probe with an average accuracy of 0.01 mm s −1 (Ka band, 60 s integration time) under nominal operational conditions (Cappuccio et al. 2020). The duallink configuration can be used near superior solar conjunctions to suppress about 75% of the noise due to charged particles in the solar corona (Bertotti et al. 1993). In addition, the tracking system is capable of range measurements at the level of 2-3 cm at Ka band (Cappuccio et al. 2020).
The operational scenario of VERITAS consists of five to seven Doppler tracking passes a week, collected by NASA's Deep Space Network (DSN) ground stations. The VERITAS observation schedule entails approximately a daily contact to ground for 8 hr, as well as 16 hr of VISAR observations. For the gravity experiment we simulate 8 hr passes for 5 days a week collected by DSN station DSS 25 (Goldstone, CA). The integration time of the Doppler observables is set to 10 s, corresponding to a displacement of the spacecraft of ∼70 km on the surface, which is sufficient to resolve gravity field features as small as 190 km after four cycles or l > 100 globally (the high-resolution gravity mapping results are beyond the scope of this work and will be published separately; for preliminary results refer to Mazarico et al. 2019).

Radar Observations and Tie Points
VISAR is an X-band interferometric radar operating at 7.9 GHz (3.8 cm wavelength) and has a 20 MHz bandwidth from which radar imagery with 30 m ground resolution pixels and topographic data with 250 m spatial resolution and 5 m elevation accuracy is produced. The radar acquires data with a look angle of 30°(angle between the antenna boresight and spacecraft nadir) and images a swath width of 14.4 km (for the VISAR flight configuration and observing geometry see Figure  1 in Hensley et al. 2020). After an orbital period, planetary rotation shifts the VERITAS ground tracks by ∼10 km; thus, the swath width provides more than 2 km of overlap between swaths acquired on adjacent orbits, enabling coherent mapping of Venus's surface. Radar data are collected on 11 out of 16 orbits per day and downlinked to Earth on the remaining 5 orbits, when two-way X-and Ka-band tracking data are acquired. Therefore, radar and Doppler data are not collected simultaneously.
VISAR transmits pulses and records the received echoes to generate images of the backscatter signal from the surface. To achieve fine resolution in the radar along-track or azimuth direction, SAR image formation combines echoes from multiple pulses when a point is illuminated by the radar antenna beam. The pixel location in a radar image is determined by the range, i.e., distance from the platform to the pixel, and the Doppler frequency, i.e., projection of spacecraft velocity on the line of sight. For Venus, the range, derived from the delay between pulse transmit time and echo return time, must be corrected for the delay induced by the thick Venus atmosphere. Atmospheric contribution amounts to 200-400 m of additional range, depending on the pixel elevation and imaging geometry. Since VISAR is an interferometer, it solves for the threedimensional position of each pixel using the range, Doppler, and interferometric phase from two spatially separated antennas.
Surface features (landmarks) imaged on multiple orbits can be identified using automated matching software. The relative range and Doppler measurement errors depend on how accurately imagery acquired from different orbits can be matched. SAR image matching is hindered by speckle that results in a grainy appearance due to the coherent nature of imaging and from differences in imaging geometry, either incidence angle or look direction. Matching accuracy is a function of the image signal-to-noise ratio, the number of looks used to reduce both speckle and thermal noise, imaging geometry differences, and the amount of scene contrast (see Appendix A.1).
Identification of radar tie points will use an automated scene matching algorithm. The automated matching algorithm computes the cross-correlation for a search window that covers the largest expected offset due to ephemeris errors. To account for the spatially variable nature of the matching accuracy and the consequent range and Doppler measurement error, we adopt the match covariance matrix used in the automated matching algorithm to estimate the matching accuracy (Frankot et al. 1994). We tune the matching metric based on match accuracy statistics from Magellan stereo data that covered approximately 20% of the surface (see Appendix A.1).
The average accuracy of the range and Doppler observations of each radar tie point is 3 m and 10 Hz, respectively, derived from an average 0.2 pixel matching accuracy using a 32 × 32 matching window of 30 m resolution pixel imagery, where each pixel corresponds to 15 m of range and 40 Hz of Doppler.
To include radar tie points in our simulations, we generated a simulated data set of radar observations. Two types of radar tie points were simulated. The first type (local tie points) is observed in the swath overlap region of adjacent orbits. These measurements permit better orbiter trajectory determination by providing constraints between adjacent orbits when VERITAS is not tracked by the DSN. The second type of radar tie points are the so-called global tie points. A point on the surface can in principle be imaged up to eight times (excluding swath overlaps, i.e., local tie points) during the four-cycle mission: one time each on the descending and ascending passes, for each of the four cycles. Each observation is separated by half a Venus sidereal period, thus enabling us to place tight constraints on the inertial motion of surface features, directly related to the rotational state of the planet. The landmarks are defined in the Venus body fixed frame (see Appendix A.2).
For the simulation, we placed landmarks on a latitude/ longitude grid with approximately 150 km spacing separating points in both directions. We exclude orbits that are used for data downlink, in solar conjunction or in power-restricted orbits where data are not collected. A total of 967,605 tie points could be obtained from a set of 387,382 unique landmarks, but we apply downsampling in the simulations (Section 3).

Joint Inversion
Our approach consists of a systematic joint inversion of both Earth Doppler tracking and VISAR tie point data sets for the simultaneous retrieval of gravity, rotation, tidal response of Venus, and the location of the geodetic control network composed by the tracked landmarks. The simultaneous solution for the spacecraft orbit and the landmark positions allows us to place tight constraints between the planetary body fixed frame and the inertial frame, increasing by an order of magnitude the sensitivity (as shown in Section 4) to the rotational state of the planet. Differently from what was done for Magellan (e.g., Davies et al. 1992), in our work we do not solve separately for the spacecraft ephemerides and the landmark positions, or apply the joint inversion to only a limited subset of orbits. For the first time, we implement the joint inversion of the two data sets for the full gravity and rotation solution in a systematic way.
A two-step solution (i.e., the preliminary inversion of the tracking data and the subsequent inversion of the tie points) suffers the unescapable problem of the propagation of the orbital reconstruction errors in the geodetic control network solution, causing systematic errors that are difficult to evaluate and mitigate. The single-step solution adopted in this work, where Doppler data and tie points are jointly processed, overcomes this problem. An important aspect that needs to be emphasized is that with this approach no a priori information regarding the landmark registration accuracy is required. The accuracy of landmark position recovery is directly estimated in the inversion process.
The substantial increase in sensitivity becomes an efficient way to overcome the limitations for future high-precision gravity and rotation experiments at Venus due to atmospheric tides (Bills et al. 2020) and high variability of the sidereal period, recently observed by Margot et al. (2021). Both these aspects require a precise tying of the spacecraft orbit to the rotational motion of the planet. The tie point inclusion gives robustness to the solution by providing a direct observation of the rotational motion.

Numerical Simulations
To assess the capabilities of VERITAS to retrieve the rotational state, Love numbers, and MOIF, we conducted an extensive set of numerical simulations replicating the nominal operational scenario of VERITAS.
We assessed the capabilities of VERITAS through a covariance analysis. Using the JPL orbit determination software MONTE (Evans et al. 2018), we integrate the trajectory of the probe, generate synthetic Doppler and VISAR data according to the assumptions outlined in Sections 2.1 and 2.2, and superimpose white Gaussian noise. To account for the noise variability observed by Cappuccio et al. (2020), we draw the noise level for each arc from a uniform distribution ranging between 0.015 and 0.038 mm s −1 for Earth-spacecraft Doppler tracking. The noise assigned to the radar tie points is 3 m and 10 Hz in range and Doppler, respectively. We then combine all the data in a least-squares filter (ORACLE) developed at Sapienza University and validated with several space missions (e.g., Iess et al. 2018). The filter implements a multiarc approach that is best suited for data analysis of long-duration gravity experiments (e.g., Durante et al. 2020;Konopliv et al. 2013).
We randomly downsampled the full set of simulated landmarks to ∼12,000 and constructed both global and local tie points. The choice of simulating only a subset of landmarks is supported by two arguments. First, observations of a landmark-dense area might be highly correlated. Selecting only well-spaced points supports the assumption that the observations are statistically independent, therefore simplifying the analysis. Second, the outcome of the simulation can be considered a conservative estimate of what would be possible if the entire data set is processed (for a discussion on the influence of the number of measured landmarks, refer to Appendix A.3).
The dynamical model used to propagate the spacecraft trajectory includes the monopole gravitational acceleration of all main solar system bodies, a degree and order 50 static gravity field of Venus (derived from Konopliv et al. 1999, we limit the spherical expansion to degree 50 since higher degrees have negligible effects on the parameters of interest), the tidal response to the Sun, the nongravitational accelerations due to solar radiation pressure and atmospheric drag, and wheel desaturation maneuvers. To account for possible mismodeling of the nongravitational forces, we employ a large set of empirical accelerations with conservative a priori uncertainties (see Appendix A.4).
Our model also includes atmospheric tides, as the spacecraft tracking system will be sensitive to their effect (Goossens et al. 2018;Bills et al. 2020). The numerical results we report in the next paragraphs are based on the assumption of a knowledge of the atmospheric tidal model with 10% accuracy. Atmospheric tide modeling and the effect of the assumed a priori knowledge on the final results are discussed in Appendix A.5.
According to recent observation campaigns of the Venus rotation rate (Margot et al. 2021), the complex coupling between the atmosphere and the planet results in sidereal period variations significantly larger than what was predicted by general circulation models (GCMs; Lebonnois et al. 2010;Cottereau et al. 2011), leading to variations of the sidereal period up to ∼3 and ∼20 minutes over timescales of 1 Earth day and 117 Earth days, respectively. If not correctly modeled, these variations induce an error in the longitude positioning of surface features that grows in time. We accounted for this perturbation in our simulations by estimating a sidereal period every 2.5 days, setting the conservative a priori uncertainty value of 20 minutes over one arc. For the full set of estimated parameters and the detailed filter setup and assumptions, refer to Appendix A.4. Table 1 reports the uncertainties (all results in tables and text are given as three times the formal uncertainty, or 3σ) attainable for the Venus rotational parameters, the Love number, and the MOIF in the nominal VERITAS mission configuration for two cases: Doppler tracking data only, and Doppler tracking data combined with VISAR observations. The inclusion of VISAR tie point measurements in the orbit determination enables a large improvement in the determination of the rotational state of Venus, not attainable with Doppler data alone. The tie points increase the accuracy in the pole location and MOIF by about a factor of 10, while a smaller improvement (∼3) is found on k 2 and its tidal phase lag.

Results and Discussion
The current estimate of Venus Love number (0.295 ± 0.066, 2σ; Konopliv & Yoder 1996), coupled with the lack of a magnetic field, does not resolve a liquid or solid core (Dumoulin et al. 2017; see Figure 1). The analysis by Dumoulin et al. (2017) indicates that the state of the core and its size, the mantle composition, and the viscous response of the interior, or rather different classes of interior conditions, can be well constrained with a knowledge of k 2 to an accuracy smaller than 3% (0.01, the 1σ VERITAS requirement) and a precise measurement of the phase lag (VERITAS has a 1σ requirement of s =  d 0 . 25 k2 ). In this work we assume the classical definition of the tidal potential phase lag (Murray & Dermott 2000), adopted also in the models by Dumoulin et al. (2017).
Thanks to the augmentation provided by VISAR tie points, our simulations show that VERITAS will be able to determine these tidal quantities with an accuracy substantially better than these threshold values (see Table 1). With the joint processing of radio tracking data and radar tie points, the right ascension and declination of the pole (α 0 and δ 0 ) can be determined with an accuracy increased by an average factor of 10, improving the results obtained by Magellan by more than 100 times and the ground-based observations (Margot et al. 2021) by more than an order of magnitude. A comparable improvement is found for the obliquity ò (σ ò = 0.12 arcsec). The considerable improvement with respect to Magellan is mainly due to the more favorable VERITAS orbital geometry, the longer time span of the gravity observations, and the substantial improvements in the end-to-end radio tracking performance (the use of Ka-band, dedicated instrumentation for media calibration, for both charged particles and tropospheric water vapor, and open loop ground receivers). The use of radar tie points also leads to substantial improvements in orbit determination, and hence in the gravity and rotational state recovery (see Table 1). VERITAS will also measure the pole precession rate Ω and derive the MOIF. Assuming that the spin axis of Venus precedes in a conical motion about the orbit normal (as detailed in Appendix A.6), the precession rate can be determined to a level of 3.6 × 10 −3 deg century −1 . The corresponding relative uncertainty in the MOIF is 0.3% of its predicted central value (0.336; Cottereau & Souchay 2009). The accurate measurement of the MOIF provides an additional, strong constraint to models of the Venus interior, by reducing the uncertainty in the density of the core and the mantle, considering the core size as constrained by k 2 (see discussion below). As is well known (see, e.g., Bills & Rubincam 1995), MOIF alone cannot uniquely determine the interior structure, even for a two-layer model of the interior. Nonetheless, it is a crucial constraint to every geophysical model of Venus interior.
Based on current models (Dumoulin et al. 2017), the constraints on k 2 and d k 2 can determine the core state and distinguish different classes of interior conditions. For example, if the core is fluid with an Earth-like composition, the core size can be obtained to within 100 km and average mantle viscosity to within an order of magnitude (see Figure 1). The latter value strongly depends on the temperature distribution and the volatile content in the mantle and therefore provides information about the heat and volatile loss of Note. We report the uncertainties in the two cases of Doppler only and Doppler + tie point analyses (three times the formal uncertainty). The tie point improvement factor is the ratio between the uncertainties obtained without and with the inclusion of VISAR data. the planet. For example, a warm and wet mantle, representative of a planetary interior that has not cooled much and has lost little of its original water, has a low viscosity, while a cold and dry mantle, representative of an efficiently cooled and outgassed interior, has a high viscosity. These two extreme models would differ in viscosity by several orders of magnitude and could be distinguished by the measurement of the phase lag. Different formation scenarios lead to different compositional models based on cosmochemical assumptions and trends among Earth-like planets to model the interior of Venus. A major difference in the models is the FeO content of the mantle, which can vary between 0.42 and 18.7 wt.%. This results in different values of MOIF ranging between 0.33 and 0.342 (∼3% variation), with otherwise the same assumption about the thermal state and the core composition (Dumoulin et al. 2017). In addition to k 2 , knowledge of the MOIF with an accuracy of 0.3% will therefore further help to distinguish the mantle composition models.
The amount of light elements in the core, particularly important for a better understanding of Venus's magnetic field evolution and also informative about Venus's conditions during core formation, is not known. The two parameters together, k 2 and MOIF, help to better distinguish the models as has already been shown, for example, for Mars (Rivoldini et al. 2011; and recently confirmed with InSight seismic data by Stähler et al. 2021), and thus better than in the models of Dumoulin et al. (2017), for which MOIFs were assumed to be unknown. The information about the density distribution from the MOIF is not unique, i.e., for the same MOIF the core can be small and dense or relatively larger and lighter. If the core of Venus is liquid, the core size can be constrained independently with k 2 and core density can be constrained in combination with the MOIF. The inverse problem of interior structure determination is degenerate, and thus different models can lead to the same values of MOIF, k 2 , and d k 2 . Different models, however, are not equally physically likely, and their selection will rely on additional constraints arising from planetary geology and geophysics as has been done, e.g., on Mars All these fundamental quantities, such as core state, size and composition, and mantle composition and viscosity, are necessary to understand the formation of Venus and its thermal and magnetic evolution. They serve, for example, as inputs (core radius and core and mantle composition) or constraints (core state and present effective mantle viscosity) for modeling core and mantle processes and the thermal and magnetic evolution (e.g., O'Rourke et al. 2018).
Additional simulations (Appendix A.5) show that even under the very conservative assumption of an atmospheric tidal model with an uncertainty of 100%, the joint solution approach guarantees the insensitivity of the rotational state solution (and thus the MOIF) to the atmospheric tide. The tidal quantities (k 2 , d k 2 ) are still determined within the mission requirements with significant margin (see Table A2).
Based on the recent findings of Margot et al. (2021), and considering the high accuracy of the VERITAS tracking system, a constant sidereal period over the mission duration can no longer be assumed, as was done in the Magellan data analysis (Davies et al. 1992). We tackle the problem of the irregular rotation rate by estimating an average sidereal period every 2.5 days. The uncertainty in the estimated rotation period by processing data acquired in 2.5 days is approximately 9 minutes. Although this uncertainty is larger than the maximum observed variation over 1 Earth day by Margot et al. (2021; about 3 minutes), VERITAS will be able to measure variations of the sidereal period over longer timescales. Indeed, the uncertainty in the sidereal period is with good approximation inversely proportional to the time span of the observations. With the aforementioned retrieved uncertainty, VERITAS would be able to distinguish 3-minute variations over 8 days. Variations of 20 minutes over 117 days (Margot et al. 2021) will be measured to good accuracy by VERITAS. The joint inversion approach proves to be extremely valuable here as well: with Doppler data only, the achievable accuracy of the 2.5-dayaveraged sidereal period increases to 43 minutes. We point out that although short-term variations of the sidereal period (length of day) result in variations of the longitudinal position of the landmarks, the spin axis solution is robust. If r i (t) is the position of a landmark in the Venus-fixed, Venus-centered frame at time t, then the spin axis w  is determined by the orthogonality condition . Any small error in the longitude of the landmarks does not affect the orthogonality condition, in the presence of a slow precession. Thus, it is only the latitudinal position of the landmarks that affects the determination of the pole position.
As a by-product of the estimation process, we can refine the location of all observed landmarks and in particular the ones contributing to the global tie points, thus providing the backbone of an accurate geodetic control network. The median values of the recovered global landmark position accuracy in altitude, latitude, and longitude (mapped on the reference surface of Venus) are M alt = 3 m, M lat = 7 m, and M lon = 6 m, respectively. The determination of the landmark positions also enables the retrieval of the radial displacement associated with the tidal forcing, parameterized by the Love number h 2 . This tidal Love number has so far only been measured from orbit only for the Moon ) and Mercury (Bertone et al. 2021 , equivalent to ∼10 cm of maximum radial displacement) corresponds to a relative uncertainty on h 2 (predicted value of 0.45-0.75; Dumoulin et al. 2017) of 20%-33%. While the relative uncertainty is not sufficient to constrain interior models, it will serve as an additional check.
The inclusion of tie points improves the orbital solution by providing observability during periods in which the spacecraft is not tracked from Earth. This aspect is particularly important for providing a uniform positional coverage along the orbit and giving robustness to the determination of physical effects such as atmospheric drag, whose variability cannot be predicted with enough accuracy by a deterministic a priori model. The denser availability of orbit-related data is also the reason for a slight improvement in the retrieval of the low-degree gravity field, since the nearly continuous orbital coverage allows a better resolution of the large spatial scales (i.e., low degree) of the planetary gravitational field.

Conclusions
By simulating the nominal mission scenario, we show that the VERITAS mission to Venus, planned for launch in 2027-2028, has the capability to determine crucial parameters (tidal Love number k 2 and MOIF) needed to substantially improve models of the planet's interior structure. The precise characterization of the tidal response of Venus via the measurement of its complex Love number will allow us to place improved constraints on the state and size of the core and on the viscous response of the planet to tidal stresses, such that different classes of interior conditions can be obtained. We show that our data analysis approach, where Doppler tracking data and radar tie points are jointly processed, is extremely effective in the determination of the rotational state of the planet and the moment of inertia factor. It pushes further the possibility of understanding of the dynamical evolution of Earth's neighboring planet. In particular, the determination of k 2 and MOIF together will help better constrain the core size, as well as the core and mantle composition. These constraints are the first-order information necessary to address the question of why Venus lacks a dynamo. Models of interior structure, temperature, and composition compatible with the measurements (k 2 , d k 2 , and MOIF) would provide the present-day boundary condition for thermal evolution models of Venus. The current dearth of information on these fundamental characteristics of Venus's interior precludes meaningful comparisons between the different evolutionary paths of Venus and the other bodies of the inner solar system. Understanding the present-day state of Venus's interior and its past evolution will offer valuable clues as to how and why Venus evolved into an uninhabitable planet.
A portion of this research was conducted at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

A.1. Tie Points Simulation
The match covariance matrix is given by where k c is an empirical constant inferred from Magellan match statistics, I is the identity matrix, and H is the Hessian of the match correlation function ( ) c x y , . For a given image offset ( ) x y , the Hessian is given by where k s is an empirical value derived from Magellan match statistics and s v is the mean X-band backscatter value for Venus, roughly -10.5 dB.

A.2. Uncertainty on Venus's MOIF from Pole Precession Measurements
The torque of the Sun on Venus determines the precession of its spin axis in a conical motion about the orbit normal. The precession rate Ω is where J 2 is the unnormalized degree 2 zonal coefficient of the gravity field of Venus, k is the MOIF, ω is the sidereal spin rate, n is the mean motion, and ò is the obliquity (angle between spin axis and orbit normal As the precession period is much longer than the VERITAS observations, the precessional motion of Venus can be described by three first-degree polynomials for the spin vector right ascension, , with t 0 corresponding to the J2000 epoch and T sid being the sidereal period (see, e.g., Archinal et al. 2011). Due to the recent findings of Margot et al. (2021), a constant T sid cannot be employed, as its short-term variations would induce significant latitudinal variations of the position of the observed landmarks. Indeed, we estimate one T sid for each arc and compute ( ) w t as follows: is the prime meridian expression valid for the ith arc, t 0 i is the starting epoch of the ith arc, and T i sid is the sidereal period estimated at the ith arc.
The precession constant Ω can be directly associated (under the assumption that nutations are negligible and small deviation from the reference position, as shown and justified in Appendix A.6) with a  and d  as with the coefficients c 1 and c 2 determined by the orbital inclination of Venus, the reference position of the pole (α 0 , δ 0 ), and the longitude of the orbital node at the reference epoch. In our simulations we therefore estimate the pole polynomial coefficients and exploit the aforementioned procedure to assess the uncertainty on Ω and thus on the MOIF.

A.3. Effect of the Number of Landmarks
To explore the effectiveness of the inclusion of radar tie points, we performed a sensitivity analysis of the results to the number of included landmarks.
We analyzed the formal uncertainty improvement factor as a function of the number of observed landmarks. We ran simulations covering the range of 1000-8000 landmarks. Not surprisingly, the improvement factor P depends on the number of landmarks n as( ) P n n , 1 2 a consequence of the assumption that the measurements are statistically independent. The results that we report can be easily scaled to an arbitrarily higher number of landmarks. The increase in the accuracy of the rotational parameters, MOIF, and k 2 shows that, while the bulk of the information matrix comes from radio tracking data, tie points, being a largely independent data set, increase the overall information content by a quite significant amount.

A.4. Filter Setup
We report here the detailed setup of the orbit determination filter used in this work.
The complete list of parameters estimated in the filter includes the following: State of the orbiter (position and velocity), degree and order 50 gravity field spherical harmonics coefficients, complex Love number k 2 , atmospheric tide parameters as described in Appendix A.5, position (latitude, longitude, radius) of all observed landmarks, Venus pole location (right ascension, declination) and its precession rate, and Venus sidereal period.
Atmospheric density variations, depending on local time and solar activity (Müller-Wodarg et al. 2016;Kliore et al. 1992), induce accelerations on the spacecraft over typical timescales ranging from half to a quarter of the orbital period T. We account for this possible mismodeling by estimating cosinusoidal along-track accelerations with period T and T/2 and a time update of 2h with an a priori uncertainty corresponding to a 25% error on Cd. To account for possible misrepresentation of the solar radiation pressure, we estimate one scale factor per arc with an a priori uncertainty set to 25%. Moreover, we include the estimation of daily momentum desaturation maneuvers with an a priori uncertainty set to 6 mm s -1 in compliance with navigation assumptions (Wallace et al. 2019).
We run the described filter solving for the set of parameters of interest. In our multiarc approach, the tracking data are subdivided into 2.5-day arcs and the parameters are divided into two sets: local parameters (those affecting a single arc, e.g., position and velocity of the orbiter) and global parameters (parameters affecting all the arcs, e.g., the gravity field of the planet). The total number of global parameters amounts to 27,320.
In Table A1 we report the a priori uncertainty assumptions of the filter.

A.5. Atmospheric Tides Modeling
In the dynamical model used in our simulations we included the effect of atmospheric tides. As shown by Bills et al. (2020), the mass transport induced by solar heating of the atmosphere is not a negligible factor for high-precision radio science experiments at Venus. For a realistic assessment of the attainable accuracies of VERITAS, we modeled the timevariable gravity field induced by solar -heating-driven pressure variations of the atmosphere.
The spherical harmonics expansion of the total (static plus atmosphere) gravity field can be written as a function of time t as C lm is the total C coefficient of degree l and order m of the gravity field, C S lm is the static coefficient, and D ( ) C t lm is the correction due to the time-variable mass transport (the same formulation applies for S lm coefficients, here omitted for brevity).
To determine the time-variable atmospheric contribution, we employed the model developed by Garate-Lopez & Lebonnois (2018) for retrieving surface pressure variations and then converted these perturbations in the associated gravity field coefficient with a technique that includes the atmospheric loading contribution on the solid planet, originally developed for Earth (Petrov 2004) and applied also on Mars (Genova et al. 2016). This procedure produces the time series of spherical harmonics expansions of the atmospheric gravity field. The gravity field perturbation induced by solar heating is a periodic signal of fundamental frequency f 1 , equal to the main forcing  where f n = nf 0 with n = 1, 2, 3, 4 and A, B are coefficients derived from the Fourier analysis specific for each coefficient, degree, and order. In our simulations we assessed the necessity of including these effects in the dynamical model of VERITAS, as its extremely precise tracking system is sensitive to the main components of the thermal tide perturbation. In particular, we assessed that if thermal tides are not accounted for, significant biases might arise in the gravity field and rotational state solution, in particular affecting the Love number k 2 . The most recent analysis of the Venus gravity field accounted for the atmospheric contribution by forwardmodeling its effect (Goossens et al. 2017(Goossens et al. , 2018. We have chosen to adopt a conservative approach and account for the intrinsic uncertainty of the atmospheric tide modeling. We model the thermal tide field up to the degree and order that guarantees that the higher degrees produce no residual signal in the Doppler residuals (i.e., degree and order 18 for f 1 , 13 for f 2 , 7 for f 3 , and 10 for f 4 ) and considered the uncertainty associated with the correction coefficients A A B B , , , lm for the frequencies f 1 through f 4 .
We evaluated the effect of the assumed a priori knowledge of the atmospheric tide model, without delving into a detailed analysis of atmospheric-dynamics-related sources of uncertainty. We have chosen, then, to assume a certain level of uncertainty on the output of the model, i.e., the correction coefficients. In particular, we explored three cases by setting different a priori uncertainties. We considered an accurate model (model uncertainty equal to 10%), a medium-accuracy model (50% uncertainty), and a coarse-accuracy model (100% uncertainty). In Table A2 we report the results relative to each of the three assumptions. It is important to note how the results, when combining Doppler data and tie points, become significantly less sensitive to the accuracy of the model, for all the parameters except k 2 and d k 2 , which, not surprisingly, have significant sensitivity to the atmospheric tides. This indicates that even a coarse a priori knowledge of the model is sufficient to meet the scientific objectives of VERITAS.

A.6. Relating Precession and Pole Coordinates
In this section we will obtain the equations to express the motion of the pole as a function of the equatorial coordinates and their time derivatives. This relation has been used in the simulations as a constraint in the determination of the precession rate and the MOIF. Finally, we will show that the errors committed by neglecting the nutations of the pole have negligible consequences in the determination of the precession rate.
The Venus ecliptic (V E ) and the (usual) Earth ecliptic (E E ) reference frames are represented by the unit vectors {u V,x , u V,y , u V,z } and {u E,x , u E,y , u E,z }, respectively. The equatorial frame is represented by {u eq,x , u eq,y , u eq,z }.
We will use the following coordinates: 1. a( ) t , d ( ) t are the right ascension and declination (equatorial J2000 coordinates); 2. l ( ) t , b ( ) t are ecliptic coordinates referred to the E E reference frame at J2000.0; are ecliptic coordinates referred to the V E reference frame at J2000.0 We define 1. the direction (as a unit vector) P V of Venus's pole; 2. the direction (as a unit vector) P 0V of the normal to the Venus orbital plane (hereafter the "orbital pole").
From Equations (A6.2) and (A6.4) we obtain   A6.17 The a( ) t and d ( ) t coordinates as functions of b l ( ) t , V V are given by the following relations: The pole motion around the orbit pole is (t = 0 corresponds to J2000.0) where δβ V , δλ V are the nutations in obliquity and in longitude, respectively, and Ω is the precession rate of the Venus pole. The precession rate is the sum of the solar precession (∼44 74 yr −1 ) and the planetary precession (−10″ yr −1 ; Simon et al. 1994    Finally, we obtain the components of the initial velocity of the pole. The evolution of Venus's orbital elements due to planetary effects is already included in our setup, so here we will consider the solar precession rate only. From Equations (A6.27) and (A6.28) one can see that while the nutations in latitude (whose amplitudes are smaller than 0 1; Cottereau & Souchay 2009) when projected in α and δ coordinates are about unchanged (rescaling factors are +2.11 and −0.57, respectively), the nutations in longitude (the largest ones, smaller than 3″ in amplitude) are strongly reduced (rescaled by factors 0.07 and 0.04).
For this reason, the short-term oscillations do not affect the estimation of a  and d  .
Therefore, neglecting the nutations, the relations between a ḋ, and Ω are as follows: Finally, the ratio a ḋ ( )˙( ) / 0 0 (Equation (A6.22)) can be used as an a priori constraint between the two quantities.