Assessment of Callisto Gravity-field Determination Using the Inter-satellite Range-rate Link

China will launch the “Tianwen-IV” mission around 2030, focusing on the orbiting exploration of Jupiter and Callisto, a moon of Jupiter. As part of this ambitious mission, a main satellite will carry another satellite that will be released in the Jupiter system to continue its journey toward Uranus. Considering the current mission planning, we propose an inter-satellite radio-observation mode that differs from the conventional observation mode of tracking from Earth to precisely determine the orbit of the satellites. Given the significance of the Callisto gravity field model in both science objectives and satellite navigation, we have conducted a series of simulation experiments to evaluate the potential of this inter-satellite range-rate data for accurately estimating the Callisto gravity field. The results obtained from the analysis demonstrate that by utilizing 40 days of ground station observations, it is possible to estimate the gravity field model of Callisto up to a degree of 70. Remarkably, when combining these ground station observations with inter-satellite observations, a comparable level of accuracy can be achieved with just 10 days of observations. Furthermore, with reduced inter-satellite observation noise, accuracy improves, enabling estimation up to 80 degrees or higher. Initial inter-satellite distance selection impacts estimation accuracy. These findings serve as a valuable test bed for the future “Tianwen-IV” mission to perform precise orbit determination and gravity field model estimation to reduce reliance on deep space stations.


Introduction
Callisto is the second largest moon of Jupiter and the outermost of the four Galilean moons.It has a diameter of about 4280 km and orbits about 1.88 million km from Jupiter.Callisto has a very thin atmosphere with a surface pressure of 7.5 × 10 −9 mbar, consisting mainly of carbon dioxide, with a composition that may include oxygen (Carlson 1999).It is outside the radiation belt of Jupiter, with no obvious tidally induced thermal heating (Khurana et al. 1998).This body formed slowly, and therefore retains information about the early evolution of the solar system.It likely hosts a subsurface ocean composed of liquid water (Khurana et al. 1998).
To gain comprehensive understanding of Callisto, China has planned a mission named "Tianwen-IV," currently in the engineering development and mission planning stage.Blanc et al. (2020) proposed two ambitious mission scenarios for China's Mission to the Jupiter System, named the Jupiter Callisto Orbiter and the Jupiter System Observer, to answer key scientific questions about the Jupiter System.Both use the combination of a main spacecraft and one or several specialized small platforms.Additionally, another proposal involves the main satellite exploring the Jupiter system while carrying a smaller satellite destined for Uranus (XU et al. 2022).While detailed background information has not been publicly disclosed, various sources outline the scientific objectives and potential mission scenarios for China's forthcoming exploration of the Jupiter system.This includes studies on the magnetic fields, plasmas, composition, and structure of Jupiter's atmosphere, as well as investigations into the surface characteristics and internal structure of Callisto.
The primary objective of our simulation experiments is to obtain a high-resolution gravity-field model to investigate the internal structure of Callisto.By precisely measuring the spacecraft's orbit and observing the perturbations induced by Callisto's gravity field, we can gain insights into its subsurface ocean and ice crust.Smith et al. (2021) obtained the density of the crust by analyzing the relationship between gravity and admittance, highlighting the necessity for a gravity resolution greater than l = 80.Genova et al. (2022) computed gravity/ topography admittance profiles using theoretical models of Callisto's interior, indicating that an accurate knowledge of the gravity field to a degree and order l > 80 is sufficient to fully characterize the ice shell density.Moreover, the gravity-field model of Callisto plays a critical role in precisely computing the orbits of Callisto's space satellites (Di Benedetto et al. 2021).Given the significance of the Callisto gravity-field model in both science objectives and satellite navigation, we propose an inter-satellite observation mode to improve the accuracy of gravity-field estimation for the "Tianwen-IV" mission.Our simulation experiments aim to evaluate the potential of this inter-satellite range-rate data for accurately estimating the Callisto gravity field.
The determination of the gravity field of Callisto has been studied in several past analyzes.Anderson et al. (2001) solved the gravitational parameter GM = 7179.292±0.009km 3 /s 2 and a second-degree gravity-field model for Callisto using Doppler data from the Galileo satellite.Additionally, they determined the mean radius of Callisto to be 2410.3±1.5 km by using four optical images.Desprats et al. (2021) evaluated the recovery of a gravity-field model from orbiter tracking for different orbital configurations.The scientific motivation for the JUICE mission is to explore Jupiter and its moons, including flyby observations of Callisto (Grasset et al. 2013).Cappuccio et al. (2022) evaluated the expected results of the Callisto gravity field from the JUICE's mission through simulation experiments.Di Benedetto et al. (2021) simulated the effect of different trajectories on the estimation of the gravity-field model and the Love number k 2 for Callisto.Genova et al. (2022) conceived a MAGIC mission concept for Callisto and designed comprehensive numerical simulations.In each of the abovementioned cases, the actual or simulated observation data used to solve for the physical parameters of Callisto was in the conventional two-range observation mode from Earth.
Due to the long distances and the influence of the Earth's ionosphere, troposphere, interplanetary plasma on the signal propagation, and the plasma environment around Jupiter, as well as the inherent near ill-posedness due to the observation geometry (Bonanno & Milani 2002), the ground-based observation mode has limitations, making it difficult to perform autonomous navigation and gravity-field solution (Leonard et al. 2012).Inter-satellite tracking between the two satellites orbiting the same celestial body can effectively mitigate these problems and its successful implementation can be seen in the notable GRACE and GRAIL missions.In the GRACE mission (Tapley et al. 2004), a series of Earth gravity-field models were obtained by using the GRACE data, such as GGM02s model and EIGEN-GRACE02S model (Tapley et al. 2005).Because of the large number of scientific results achieved by the GRACE mission, the GRACE-FO mission was launched in 2018 to continue its scientific goals.The GRACE-FO mission employs inter-satellite laser ranging interferometer to provide more precise and comprehensive information about the Earth's gravity field and its variations (Kornfeld et al. 2019).Similarly, in the GRAIL mission, a lunar gravity-field model up to degree and order of 1500 was obtained by employing inter-satellite observation data (Park et al. 2015).Therefore, in the future deep space exploration missions, the primary satellite could release one or more secondary satellites to form a constellation of satellites for inter-satellite observations, which not only reduce launch costs, but also greatly increases the scientific objectives (Radhakrishnan et al. 2014(Radhakrishnan et al. , 2016;;Murugan & Agrawal 2020).
In addition, there are many advantages of inter-satellite observations.Genova & Petricca (2021) proposed an alternative radio science system architecture, sharing certain equipment or functionalities, resulting in significant mass and power savings.Di Benedetto et al. (2019) described a mission concept based on a mothercraft-daughtercraft configuration, with small satellite going into polar orbit around Europa to probe the internal salty ocean and establishing a two-way relay link when the mother spacecraft during Europa flyby.Such a distributed mission in the harsh radiation environment of Jupiter reduces overall mission risk.Moreover, multiple satellites simultaneously observing the same planet, resulting in multisource data that can be fused for obtaining more scientific results, such as high-resolution gravity-field model (Tapley et al. 2005;Park et al. 2015).Furthermore, using satellite-to-satellite observations, the absolute position and orientation of the two satellites with respect to inertial space can be determined (Hill & Born 2007).This capability not only reduces the observational burden on ground stations, but also provides many purposes, such as data transfer, increasing autonomy, reducing costs and higher accuracy in orbit determination (Schlicht et al. 2020).In this work we will make use of this inter-satellite observation to demonstrate its contribution in Callisto gravity-field recovery.This paper is divided into the following sections: in Section 2, we describe the dynamical model and the inter-satellite observation mode in detail; in Section 3, we analyze the effects of different observation modes, observation noise levels, observation times, initial inter-satellite distances, orbital altitudes, and inclinations on estimating the gravity-field model for Callisto; in Section 4, we present our overall conclusions.

Methodology and Observation Mode
In this section, we describe in detail the dynamical model and key parameters for the two satellite orbit models, and introduce the inter-satellite observation mode.

Simulation Experiment Design
For this study, we used the Small body Precise Orbit determination Toolkit, a precision orbit determination program developed by Wuhan University.This software has been used in Rosetta tracking data processing and various simulation experiments for multiple missions (Gao et al. 2021a(Gao et al. , 2021b)).We adopted the DOP853 algorithm with relative and absolute tolerances of 1e-11 for multiarc integration.To determine the orbit of satellites and the gravity-field coefficients of Callisto, we utilized nonlinear least squares with constraints as an iterative method.The equation of motion for the satellite is defined by where r is the position vector of the satellite at time t in the J2000 frame; a Callisto denotes the acceleration due to gravity field of Callisto, a NB denotes N-body perturbation from other bodies in the solar system, a SRP denotes the solar radiation pressure, and a SOT denotes gravitational effect of the solid tide of Callisto.The gravity-field potential of Callisto is expanded into: where G is the gravitational constant.r, θ, and λ are the radial distance, latitude, and longitude of the satellite in a Callisto fixed frame.R is the reference radius of Callisto, and N is the maximum degree of the Callisto gravity-field model.P nm is the fully normalized associated Legendre function with degree n and order m.C nm , and S nm are the fully normalized spherical harmonic coefficients.The estimation of Callisto's gravity field is the joint estimation of the GM and C nm , S nm value.The acceleration a Callisto from the gravity field of Callisto is the gradient of Equation (2).
Callisto is not a strictly rigid body, and it is also affected by the tides generated by other celestial bodies.The corresponding tidal potential function can be expressed as The periodic variations of the Callisto's quadrupole field in response to tidal potential of Jupiter are controlled by the Love number k 2 (Di Benedetto et al. 2021), i.e., n = 2, therefore, this formula neglects Love numbers of degree > 2. m i denotes the mass of Jupiter.R i is the position vectors of Jupiter in the Callisto fixed coordinates.r is the satellite in the Callisto fixed coordinates.R is the reference radius of Callisto.γ i is the satellite-Callisto -Jupiter angle, and k 2 is the Love number of Callisto.The acceleration a SOT from the solid tide of Callisto is the gradient of Equation (3).The dynamic models used in this simulation experiment are described in detail in Table 1.
For this simulation, we used a synthetic gravity-field model of Callisto from Desprats et at. (2021), where its second degree and order coefficients are from the measured values determined by Anderson et al. (2001) and the other coefficients are scaled following the lunar gravity-field model.Ephemerides of the celestial bodies in N-body perturbation are from DE431 and JPL small body database browser.A cannonball model was used for the solar radiation pressure (Montenbruck & Gill 2000).The satellites for the "Tianwen-IV" mission are still under development, and precise details regarding their area and mass are unavailable.Consequently, for the purposes of simulation experiments, typical values of the surface-to-mass ratio (0.02 m 2 kg −1 ) for both satellites have been employed (Gao et al. 2021b).The solar pressure coefficient is set to 1.5 and estimated in each arc.Other key parameters are shown in Table 2.
The parameters in Table 2 were used as the nominal case in this simulation experiment.To mitigate the impact of unequal data quantities on the solutions, different observation times were set for the two observation modes.An arc length of 0.5 days was adopted to minimize the accumulation error of the integrator and enhance parallel computational efficiency.This approach involved conducting a single-arc estimation over 20 arcs of the initial state, enabling the global estimation of the gravity-field coefficients.At the initial propagation, we set the inter-satellite distance between the two satellites at the initial time to about 100 km.Both orbits are (almost) polar orbit with 200 km altitude.The elevation cutoff angle of the Earth ground station was set to 10°.In the following sections, the noise level, initial inter-satellite distance, orbital altitude, and inclination were changed to assess their impact on the solutions.

Inter-satellite Observation Mode
Two observation modes were designed, as shown in Figure 1.The white line represents the two-way range-rate measurement between the ground station on Earth and the primary satellite orbiting Callisto, hereinafter referred to as 2W.The red line represents the range-rate measurement between the primary satellite and the secondary satellite, hereinafter referred to as SST.
Like the GRACE and GRAIL missions, we used the highprecision range-rate date between the two satellites to estimate the gravity-field model.The simplified measurement equations for the SST observation mode are as follows (Leonard et al. 2012) In Equation (4), r P is the position vector of the primary satellite in the Barycentric Celestial Reference System (BCRS).r S is the position vector of the secondary satellite in BCRS.ρ = r P − r S is the position vector of the primary satellite with respect to the secondary satellite and    r r P S r =is the relative velocity between the two satellites. r represents the range-rate data of the primary satellite and the secondary satellite.Notably, due to the close proximity of the spacecraft, the primary error source originating from the ultrastable oscillators is effectively mitigated (Genova & Petriccal 2021).In practical scenarios, accounting for light-time delay becomes essential.This delay reflects the time it takes for light to travel from the primary satellite to the secondary satellite.Calculating this duration onboard is necessary to access the ephemerides precisely when needed.However, for the purpose of our study, the inclusion of this effect would not yield additional insights, and thus, it has been disregarded (Casini et al. 2023).

Results and Discussion
The accuracy of the estimated gravity field is assessed by several parameters.The "maximum degree" of the gravity field is typically determined from the point where the power spectra of the a posteriori uncertainty and nominal values gravity-field coefficients cross.To make this determination, the coefficient degree root mean square (rms), σ n , and the coefficient error rms, δ n , are used.We choose the same criterion for our study.The formulas to calculated the values are as in where C nm s and S nm s are formal errors of the gravity-field coefficient C nm and S nm , respectively.The σ n shows the "frequency" intensity of the gravity-field model.The δ n are retrieved from the a posteriori covariance matrix of the model.In the following section, we evaluate the effects of different observation modes, initial inter-satellite distances and orbital inclinations on estimating the gravity-field model of Callisto.

Observation Data Descriptions
The distance between the two satellites undergoes constant changes due to the differential perturbations experienced by their respective orbits.These perturbations arise from the slightly uneven mass distribution of Callisto (Kim & Tapley 2002) and the relative states of other perturbing bodies.Figure 2 depicts the variation in distance between the primary satellite and the secondary satellite over a 10 days period.
As depicted in Figure 2, after ten days of orbital motion, there is a noticeable difference in the orbital altitudes of the two satellites, exceeding 2 km.This variation in inter-satellite altitude consequently leads to fluctuations in the inter-satellite distance.Specifically, after ten days of orbiting, the distance between the two satellites has deviated by more than 2 km compared to the initial inter-satellite distance of 100 km.A similar pattern of inter-satellite distance variations was observed in a simulated SST observation mode for Mercury, Venus, and Mars, as demonstrated by Genova & Petricca (2021).This phenomenon serves as a significant indicator of the non-spherical nature of the gravity field.The variations of the range-rate data over time for both observation modes are given in Figure 3.
From Figure 3, we can see that after ten days of orbiting, the inter-satellite range-rate increased to 1 m s −1 , which is already much larger than inter-satellite noise level.The range-rate data from ground station on Earth are not continuous due to the relative motion/rotations between Earth and Callisto.Because of the relative motion of Earth/Jupiter around the Sun, and the motion of Callisto around Jupiter, and Earth rotation, the rangerate data from the station on Earth vary significantly, with a maximum value of approximately 30 km s −1 and a minimum value of about 10 km s −1 .

Comparing the Effects of Observation Modes
In this section, both types of range-rate data, namely the inter-satellite range-rate and the range-rate data from the Earth station, were combined to estimate the global gravity-field model of Callisto.The impact of the inter-satellite range-rate data on the estimation results was evaluated by analyzing the error spectra, as presented in Figure 4.
From Figure 4, we can see that the gravity-field coefficients up to degree 70 can be estimated using 40 days of observations from ground station in the 2W observation mode.After incorporating inter-satellite observations between the primary and secondary satellites, it becomes feasible to estimate with approximately the same level of accuracy as 40 days of ground station observations using only 10 days of observations.However, there is a notable distinction in the precision of the lower degree terms.
Drawing on the analysis of relevant prior research (Shan et al. 2018;Yan et al. 2018;Ye et al. 2018), we set the noise level of ground range-rate measurement to 0.1 mm s −1 , corresponding to the assumption of X-band frequency utilization exclusively in this case.Furthermore, if Ka-band intersatellite observations are employed, like the GRAIL mission, the data accuracy has the potential to achieve a few μm s −1 levels.Therefore, the inter-satellite range-rate measurement was conservatively set from 0.1 mm s −1 to 0.01 mm s −1 to verify the influence of noise levels on estimating the Callisto's gravity-field model.As shown in Figure 4, reducing the noise level by one order of magnitude improves the accuracy of the gravity-field coefficients by one order of magnitude, and the gravity-field coefficients of degree 80 or more than can be estimated.Therefore, in the SST observation mode, the accuracy of the gravity-field estimation is largely related to the measurement noise.

Comparing the Effects of Observation Time Length
Furthermore, the extent of subsatellite point coverage directly influences the acquisition of gravity-field signals from  Callisto, thereby playing a crucial role in the accurate determination of gravity-field models.
From Figure 5, by extending the observation time from 10 days to 35 days, the coverage density of the subsatellite points can be enhanced.When the observation time is 35 days, the satellite's subsatellite point almost fully covers the surface of Callisto.The greater coverage of subsatellite points enhances the observational opportunities, enabling a more comprehensive and precise characterization of Callisto's gravity field.Consequently, we evaluated the accuracy of the estimation based on the extended observation times.
From Figure 6, it becomes evident that extending the observation time while maintaining the inter-satellite observation noise at 0.01 mm s −1 yields an improved accuracy in estimating the gravity-field model.Notably, Genova et al. (2022) achieved a higher precision in estimating the Callisto gravity field utilizing the 2W observation mode, which was based on one year of observations.Notably, they assumed Titan's gravity coefficients are proportional to the ratio of the Callisto's radius for the degree 3, while higher-degree coefficients are following by topographic relief.The gravity-field model used in this paper is from Desprats et at. (2021).These two gravity-field models have different gravity-field signal strengths, with Genova et al. (2022)ʼs being two orders of magnitude lower, thus requiring longer observation times and lower orbital altitudes to obtain gravity-field information.Moreover, the inter-satellite observations are fused in this paper, which is more sensitive to the gravity-field signal.The further analysis of the inter-satellite observations is detailed in Figure 16 in Appendix A. The purpose of this simulation experiment is not to obtain a more accurate gravity-field model, but to compare the effect of observation time on estimating the high-degree gravity field for subsequent exploration missions.The relationship between the degree of the gravity-field model and admittance is described in Figures 17 and 18 in Appendix B.

Comparing the Effects of Initial Inter-satellite Distances
In the subsequent section, meticulous attention was given to designing the orbits of the two satellites operating in the 2W +SST observation mode.An essential consideration in this  design process was the establishment of an optimal intersatellite distance between the primary and secondary satellites.This parameter plays a pivotal role in ensuring the effectiveness and accuracy of the inter-satellite observations, thereby facilitating the acquisition of precise and reliable gravity-field data.The effect of the inter-satellite distance on estimates of the gravity-field model of Callisto was evaluated by using initial inter-satellite distances of 100, 200, and 300 km.
From Figure 7, we can see that increasing the initial intersatellite distance from 100 to 300 km resulted in a slight improvement in the medium degrees, but led to a deterioration of the high-degree ( > 54) coefficients.This is in agreement with the result in Yan et al. (2013).Resonance occurs when the inter-satellite angular in-plane separation coincides with several wavelengths of a particular harmonic degree, which follows where l denotes the resonant degree, and α is the angular separation between the two satellites.When the inter-satellite distance is 300 km, the separation angle is about 6.6°, the resonant harmonic degree should be around 54, which is consistent with the results in Figure 7.In addition, the intersatellite distance should not be too large because of the reduced  sensitivity of the inter-satellite observation data to a shallow mass distribution inside Callisto.

Comparing the Effects of Orbital Altitudes
The strong signal attenuation of the Callisto's gravitational field with the orbital altitude of the satellites of course plays a role in the retrieval of the gravity field.Higher-precision gravity information can be gained by adopting a lower orbital altitude for the two satellites.However, when the satellites orbits are too low, the atmosphere of Callisto will have adverse effect on the spacecraft.On the one hand, the atmospheric drag can cause the deterioration or even reduce its lifetime; on the other hand, the atmospheric perturbation may degrade the accuracy of the spacecraft orbit determination.Therefore, in order to evaluate the effect of different orbital altitudes on estimating the gravity-field model, orbital altitudes of 200, 300, and 400 km for the both satellites were used in this simulation experiment.The other parameters are same in the computation.
As shown in Figure 8, when the orbit altitude is 200 km, gravity-field models up to degree 80 can be computed.However, when increasing the orbit altitude to 400 km, only gravity-field models up to degree 50 can be computed, accompanied by a decrease in accuracy.When the two satellites are in lower orbit, gravity information in the medium-short-wavelength band is obtained.However, for sake of the compromise between the recovery accuracy of the gravity field and the satellite system life, the choice of orbital height for "Tianwen-IV" mission will require a trade-off analysis.

Comparing the Effects of Orbital Inclinations
The gravity-field coefficients of Callisto can be divided into 3 families: the zonal harmonic coefficients (l ≠ 0, m = 0), the sectorial harmonic coefficients (l = m ≠ 0) and the tesseral harmonic coefficients (l ≠ m ≠ 0).Different degrees of orbital inclination are sensitive to different types of Callisto's gravity-field coefficients of degree l and order m.Consequently, in this section, we focus on investigating the scenario where the two spacecraft orbits possess different inclinations.Specifically, we alter only the orbital inclination of the secondary satellite from a polar orbit to 80°or 70°while maintaining a polar orbit for the primary satellite, as illustrated in Figure 9.
In Figure 9, the primary satellite's orbit is represented by the orange line, indicating a polar orbit.On the other hand, the blue lines depict the orbit of the secondary satellite.When the two satellites are in different orbits, noticeable disparities arise in the orbital altitudes of the two satellites.These differences can be attributed to the uneven mass distribution within our Callisto model, as visually depicted in Figure 10.
Figure 10 demonstrates that when the two satellites are in a polar orbit, the orbital altitude difference is minimal.This characteristic ensures a stable formation over the long term, with no significant relative drifts between the orbit planes.However, by reducing the orbital inclination of the secondary  satellite to 80°, the orbital altitude difference between the two satellites becomes the largest, reaching approximately 25 km.This indicates that the two satellites experience more pronounced variations in gravitational force.Figure 11 showcases the gravity-field model solutions obtained from different orbital inclinations.
As depicted in Figure 11, the estimation accuracy shows a slight improvement when the two satellites are in a polar orbit.However, altering the orbital inclination of the secondary satellite has minimal impact on the estimation process.The advantage of a polar orbit lies in its ability to provide nearglobal coverage of satellite observations.Conversely, when the two satellites have different inclinations, the distinct orbital altitudes result in a more pronounced vertical gravity gradient, thereby offering additional information in that direction.Decreasing the orbital inclination of the secondary satellite further enhances coverage at low latitudes, which proves valuable for estimating the local gravity field of the study area.
Furthermore, Figure 11 reveals that the accuracy of lowdegree coefficients, up to degree 15, in the gravity-field estimation is significantly higher in cases where the satellites have different orbital inclinations compared to the polar orbit solution.This emphasizes the importance of incorporating various inclinations to improve the accuracy of low-degree  gravity-field coefficients.Therefore, for providing insight into the level of correlation between low-degree gravity-field coefficients, the normalized correlation coefficients for the gravity-field coefficients up to the 8th degree are selected to show in Figures 12-14.
Figures 12-14 illustrate that when tracking data is solved using polar orbits, the gravity-field coefficients exhibit significant off-diagonal cross-correlations.However, by reducing the orbital inclination of the secondary satellite from 90°t o 80°or 70°, the correlation between the gravity-field  coefficients is greatly diminished.This reduction in correlation can be attributed to the inclusion of cross-track information in the inter-satellite observation data.In the case where two satellites are in the same polar orbit, trailing each other, the lack of cross-track information limits the ability to estimate the gravity-field coefficients accurately.

Orbit Evolution
The analysis of the inter-satellite observation mode's orbit design has been conducted, with the estimation of gravitational field coefficients.Notably, the accuracy of gravity-field estimation is predominantly influenced by observation time, observation noise, and orbit altitude.The gravity-field model, being a critical factor in determining satellite orbits, is subjected to further investigation in this section.Specifically, the 80th degree gravity-field model is truncated to both the 60th degree and 40th degree.Subsequently, the impact of these different-degree gravity-field models on orbital evolution is meticulously examined.For this analysis, the initial orbital altitude is set at 200 km, the initial orbital inclination is fixed at 90°, and the integration time spans 10 days.By comparing the orbital evolution results of these three gravity-field models, the differences in satellite position and velocity are revealed, as illustrated in Figure 15.
As depicted in Figure 15, the positional difference between the orbit calculations based on an 80th degree gravity-field model and the model truncated to the 40th degree is 3 km, with a velocity difference of 2 m s −1 .Similarly, when comparing the results of the orbit calculations using the 80th degree gravityfield model and the model truncated to the 60th degree, the positional difference is reduced to 0.3 km, and the velocity difference is 0.2 m s −1 .While a 0.3 km position difference may appear small, it holds substantial significance in precision orbit determination.In addition, we also evaluated the performance  of the two observation modes in precise orbit determination and presented it in Table 3 in Appendix C.

Discussion and Conclusion
In this paper, we proposed a primary-secondary tracking mode for the future "Tianwen-IV" exploration mission, and compare the ability of this tracking data mode to estimate the gravity-field coefficients of Callisto w.r.t to Earth-only radio tracking mode.To assess the accuracy of the estimation, we studied the effect on the estimation of different inter-satellite observation noise levels, observation times, orbital altitudes, and orbital inclinations.Finally, we evolved the satellite's orbit by truncating gravity-field models of different degrees.
Our study demonstrates the efficacy of inter-satellite observations in enhancing the accuracy and resolution of gravity-field estimation for Callisto.The gravity-field coefficients up to degree 70 can be estimated through 40 days of ground station observations in the 2W observation mode.With the incorporation of inter-satellite observations between the primary and secondary satellites, achieving a comparable level of accuracy as 40 days of ground station observations is feasible within just 10 days of observations.In the SST observation mode, the noise of the inter-satellite data largely determines the accuracy of the gravity-field model estimate.By reducing the noise level of inter-satellite observation by an order of magnitude, the accuracy and maximum degree of the gravity-field coefficients is improved by an order of magnitude.Lowering the orbital altitude allows for greater sensitivity of a higher-degree gravity-field estimate.When the orbit altitude is 200 km, gravity-field models up to degree 80 can be computed.When the inter-satellite distance increases, there is little difference in the resolution of the estimated gravity-field model, but resonance occurs for some harmonic degrees of the gravity field, affecting the accuracy of the estimation.The gravity-field model plays the most important role in the evolution of satellite orbit.Conducting a thorough orbit analysis employing both the 60th and 80th degree gravity-field models, we unearth a noteworthy position difference of 0.3 km.
There are some ways of improvement in this simulation experiment.Accurate modeling of nonconservative forces is a challenge, such as atmospheric drag and solar radiation pressure.Hubble Telescope observations confirmed the presence of an oxygen-dominated collisional atmosphere on Callisto (Cunningham et al. 2015), which affects satellites.Moreover, solar radiation pressure also affects the satellites and may also influence the final covariance estimate if the uncertainty of the thermos-optical properties of the satellite surface is not properly considered.In the future exploration mission, the modeling of the gravity field may not be limited to the inter-satellite observation data, data fused from multiple sources however, such as onboard optical image data and radar altimetry data, would generate more types of observations.If these data can be processed by the onboard system in real time, then the accuracy of navigation, positioning, and landing could be greatly improved.This will be the focus of a subsequent experiment; we will carry out the fusion of multisource data from onboard sensors.
This work provides a reference for radio science to be performed by the "Tianwen-IV" mission, and demonstrates the value of implementing the satellite-to-satellite observation mode to improve the accuracy of the gravity-field model, consequently advancing the precision of satellite orbit determination.In addition, it would be worth to study the contribution of highresolution gravity-field model in interior structure analysis.Especially, for the future work we will study the relationship between the gravity field and Callisto's subsurface ocean structure, including the ocean thickness, depth and coverage.

Figure 1 .
Figure 1.Schematic diagram of the two observation modes.

Figure 2 .
Figure 2. Inter-satellite distance change and orbital altitude difference of the two satellites.

Figure 3 .
Figure 3. Range-rate data in the two observation modes.

Figure 4 .
Figure 4. Power spectra corresponding to the two observation modes (in black: the spectrum of the reference gravity field).

Figure 5 .
Figure 5.The subsatellite point coverage for 10 days (left) and 35 days (right) observation times.The sidereal period of Callisto is about 16.69 days.

Figure 6 .
Figure 6.Power spectra corresponding to the different observation times.

Figure 7 .
Figure 7. Power spectra corresponding to the initial satellite distances.

Figure 8 .
Figure 8. Power spectra corresponding to the orbital altitudes.

Figure 9 .
Figure 9. Orbital chart of the two satellites.

Figure 10 .
Figure 10.Orbital altitude difference between two satellites with different orbital inclinations.

Figure 11 .
Figure 11.Power spectra corresponding to the orbital inclinations.

Figure 12 .
Figure 12.Normalized correlation coefficients between the gravity-field coefficients for the orbital inclination of 90°.

Figure 14 .
Figure 14.Normalized correlation coefficients between the gravity-field coefficients for orbital inclinations of 70°and 90°.

Figure 15 .
Figure 15.Position and velocity difference for d/o 40 and d/o 80, d/o60 and d/o80.

Table 1
Full Descriptions of Dynamic Models on Spacecraft.order synthetic model from Desprats et al. (2021) N-body perturbation Sun, eight planets, Pluto, Io, Europa, Ganymede, Ceres, Vesta, and Pallas Solar radiation pressure area-to-mass ratio 0.02 m 2 /kg, solar radiation coefficient C r =1.5 for the two satellites Solid tide of Callisto k 2 : 0.5 (Di Benedetto et al. 2021)

Table 3
Formal Error of the Position Vector of Two Satellites