# Radio Occultation Data Analysis With Analytical Ray-Tracing

A. Caruso, A. Bourgoin, A. Togni, M. Zannoni, P. Tortora

► To cite this version:

A. Caruso, A. Bourgoin, A. Togni, M. Zannoni, P. Tortora. Radio Occultation Data Analysis With Analytical Ray-Tracing. Radio Science, 2023, 58, 10.1029/2023RS007740. insu-04475257

HAL Id: insu-04475257

https://insu.hal.science/insu-04475257

Submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Radio Science®

RESEARCH ARTICLE

10.1029/2023RS007740

**Key Points:**

* Novel approach to analyze radio occultation data by spherically symmetric atmospheres

* Analytical ray-tracing method is used to investigate one-way and two-way occultations

* Results are consistent with those obtained with Abel inversion and numerical ray-tracing methods while requiring low computational burden

**Correspondence to:**

A. Caruso,
andrea.caruso15@unibo.it

**Citation:**

Caruso, A., Bourgoin, A., Togni, A., Zannoni, M., & Tortora, P. (2023). Radio occultation data analysis with analytical ray-tracing. *Radio Science*, 58, e2023RS007740. https://doi.org/10.1029/2023RS007740

Received 24 APR 2023

Accepted 29 AUG 2023

© 2023. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_1.png" /></figure>

Radio Occultation Data Analysis With Analytical Ray-Tracing

A. Caruso¹ ID, A. Bourgoin², A. Togni¹, M. Zannoni¹,³ ID, and P. Tortora¹,³ ID

¹Department of Industrial Engineering, Alma Mater Studiorum—Università di Bologna, Forlì, Italy, ²SYRTE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC P6, LNE, Paris, France, ³Centro Interdipartimentale di Ricerca Industriale Aerospaziale, Alma Mater Studiorum—Università di Bologna, Forlì, Italy

**Abstract** Radio occultation experiments are a sensing technique dedicated to the remote sounding of planetary atmospheres. The technique exploits the frequency shift of a radio signal due to refraction in a planetary atmosphere. The aim is to infer the physical properties of the neutral atmosphere (e.g., pressure and temperature) and ionosphere (e.g., the electron number density). For one-way occultations, the data processing usually relies on Abel transform algorithms when the atmosphere is spherically symmetric. For two-way occultations, such techniques require the introduction of approximate relationships for the bending experienced by the signal to be obtained. In this context, we introduce a new method to process two-way occultations data by spherically symmetric atmospheres using a ray-tracing approach. However, the numerical integration of the geometrical optics equation through the atmosphere requires a significant computational time due to initial pointing issues. For this reason, our novel algorithm exploits a closed-form solution to the equations of geometrical optics (Bourgoin et al., A&A, 624, A41, 2019, https://doi.org/10.1051/0004-6361/201834962) applied to a spherically symmetric atmosphere. Within this approach, the bending is directly provided by the analytical solution and no numerical integration is required. In addition, we develop a procedure enabling us to disentangle the contributions from dispersive and neutral media in the frequency shift. This procedure is validated by comparing our vertical profiles to those obtained using Abel inversion or numerical ray-tracing for Mars and Titan occultation experiments. We show that our algorithm provides similar results to purely numerical ray-tracing algorithms while significantly decreasing the computational time.

**Plain Language Summary** Radio occultation experiments are a remote sensing technique used to sound the atmospheres of planets and moons to infer their physical characteristics. For this purpose, radio occultations take advantage of the frequency shift brought on by refraction when a radio signal passes through the atmospheric medium. A fundamental parameter required to derive atmospheric properties is the bending experienced by the radio ray while traversing the atmosphere. For one-way experiments, the bending angle can be calculated using the frequency shift once the occultation geometry is known. When dealing with two-way occultations, it cannot be directly inferred from the frequency shift. In these cases, numerical ray-tracing algorithms are required, which, however, are computationally expensive. In this context, we develop a novel procedure based on a closed-form solution to the equations of geometrical optics. This new method can produce results consistent with that obtained by fully numerical algorithms but require a lower computational burden.

# 1. Introduction

A radio occultation event occurs when a radio link is established between a spacecraft and a ground-based antenna once the spacecraft is passing behind a planetary atmosphere as seen from the Earth. Then, the vertical physical properties of the occulting atmosphere can be inferred by exploiting the frequency changes of the radio signal as it progressively sinks into the optical medium (Kliore et al., 1965). Hence, the processing of radio occultation data is based on an accurate modeling of refraction.

At the geometrical optics approximation, the effect of refraction on light propagation is twofold. First, it causes the phase of the signal to speed up or slow down when it crosses through a neutral or a dispersive medium. Second, it bends the light trajectory toward regions of higher refractivity. All in all, refraction thus generates both an excess path delay and a geometrical delay with respect to a signal that would have been transmitted in a neat vacuum (Bourgoin et al., 2021). The variations of these atmospheric delays with respect to the positions of the emitter and receiver eventually describe how refraction may disrupt the frequency of the electromagnetic signal (see e.g., Bourgoin, 2020; Bourgoin et al., 2021, 2022; Synge, 1960). These are the mechanisms that one should properly account for in order to provide an accurate processing of radio occultation data.
Radio Science®

RESEARCH ARTICLE

10.1029/2023RS007740

**Key Points:**

* Novel approach to analyze radio occultation data by spherically symmetric atmospheres

* Analytical ray-tracing method is used to investigate one-way and two-way occultations

* Results are consistent with those obtained with Abel inversion and numerical ray-tracing methods while requiring low computational burden

**Correspondence to:**

A. Caruso,
andrea.caruso15@unibo.it

**Citation:**

Caruso, A., Bourgoin, A., Togni, A., Zannoni, M., & Tortora, P. (2023). Radio occultation data analysis with analytical ray-tracing. *Radio Science*, 58, e2023RS007740. https://doi.org/10.1029/2023RS007740

Received 24 APR 2023

Accepted 29 AUG 2023

© 2023. The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

Radio Occultation Data Analysis With Analytical Ray-Tracing

A. Caruso¹ ID, A. Bourgoin², A. Togni¹, M. Zannoni¹,³ ID, and P. Tortora¹,³ ID

¹Department of Industrial Engineering, Alma Mater Studiorum—Università di Bologna, Forlì, Italy, ²SYRTE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC P6, LNE, Paris, France, ³Centro Interdipartimentale di Ricerca Industriale Aerospaziale, Alma Mater Studiorum—Università di Bologna, Forlì, Italy

**Abstract** Radio occultation experiments are a sensing technique dedicated to the remote sounding of planetary atmospheres. The technique exploits the frequency shift of a radio signal due to refraction in a planetary atmosphere. The aim is to infer the physical properties of the neutral atmosphere (e.g., pressure and temperature) and ionosphere (e.g., the electron number density). For one-way occultations, the data processing usually relies on Abel transform algorithms when the atmosphere is spherically symmetric. For two-way occultations, such techniques require the introduction of approximate relationships for the bending experienced by the signal to be obtained. In this context, we introduce a new method to process two-way occultations data by spherically symmetric atmospheres using a ray-tracing approach. However, the numerical integration of the geometrical optics equation through the atmosphere requires a significant computational time due to initial pointing issues. For this reason, our novel algorithm exploits a closed-form solution to the equations of geometrical optics (Bourgoin et al., A&A, 624, A41, 2019, https://doi.org/10.1051/0004-6361/201834962) applied to a spherically symmetric atmosphere. Within this approach, the bending is directly provided by the analytical solution and no numerical integration is required. In addition, we develop a procedure enabling us to disentangle the contributions from dispersive and neutral media in the frequency shift. This procedure is validated by comparing our vertical profiles to those obtained using Abel inversion or numerical ray-tracing for Mars and Titan occultation experiments. We show that our algorithm provides similar results to purely numerical ray-tracing algorithms while significantly decreasing the computational time.

**Plain Language Summary** Radio occultation experiments are a remote sensing technique used to sound the atmospheres of planets and moons to infer their physical characteristics. For this purpose, radio occultations take advantage of the frequency shift brought on by refraction when a radio signal passes through the atmospheric medium. A fundamental parameter required to derive atmospheric properties is the bending experienced by the radio ray while traversing the atmosphere. For one-way experiments, the bending angle can be calculated using the frequency shift once the occultation geometry is known. When dealing with two-way occultations, it cannot be directly inferred from the frequency shift. In these cases, numerical ray-tracing algorithms are required, which, however, are computationally expensive. In this context, we develop a novel procedure based on a closed-form solution to the equations of geometrical optics. This new method can produce results consistent with that obtained by fully numerical algorithms but require a lower computational burden.

**1. Introduction**

A radio occultation event occurs when a radio link is established between a spacecraft and a ground-based antenna once the spacecraft is passing behind a planetary atmosphere as seen from the Earth. Then, the vertical physical properties of the occulting atmosphere can be inferred by exploiting the frequency changes of the radio signal as it progressively sinks into the optical medium (Kliore et al., 1965). Hence, the processing of radio occultation data is based on an accurate modeling of refraction.

At the geometrical optics approximation, the effect of refraction on light propagation is twofold. First, it causes the phase of the signal to speed up or slow down when it crosses through a neutral or a dispersive medium. Second, it bends the light trajectory toward regions of higher refractivity. All in all, refraction thus generates both an excess path delay and a geometrical delay with respect to a signal that would have been transmitted in a neat vacuum (Bourgoin et al., 2021). The variations of these atmospheric delays with respect to the positions of the emitter and receiver eventually describe how refraction may disrupt the frequency of the electromagnetic signal (see e.g., Bourgoin, 2020; Bourgoin et al., 2021, 2022; Synge, 1960). These are the mechanisms that one should properly account for in order to provide an accurate processing of radio occultation data.
Up to now, numerous occultation experiments have been successfully conducted, from the occultation of Mariner IV by Mars on 15 July 1965 (Kliore et al., 1965) to the more recent Ganymede-Juno occultation on 7 June 2021 (Buccino et al., 2022). Comprehensive reviews of objectives, methods, and results of some of the radio occultation experiments carried out during the last 50 years can be found in Fjeldbo and Eshleman (1968), Fjeldbo et al. (1971), Lindal (1992), Jenkins et al. (1994), Pätzold et al. (2005, 2007), Schinder et al. (2011a, 2011b), Hinson et al. (2017), Withers et al. (2018), and Dalba and Withers (2019).

Distinctions must be made between one-way and two-way occultation experiments. In one-way downlink mode, the onboard frequency source, for example, an Ultra Stable Oscillator (USO), provides a stable frequency reference for the transmitted signal to the Earth. In a one-way uplink configuration, an Earth-based station transmits a reference signal to the spacecraft, which records and re-transmits it back to Earth via telemetry. Uplink occultations benefit from the fact that they can provide a higher signal-to-noise ratio than typical downlink experiments. Instead, for two-way occultations, the frequency of an uplink signal (from the ground station to the spacecraft's antenna) is used as a value of reference for the downlink signal (from the spacecraft back to the Earth).

One of the most widely used methods to process one-way radio occultation data for spherically symmetric atmospheres is based on Abel transform algorithms (Fjeldbo et al., 1971; Phinney & Anderson, 1968). This method was applied to the atmospheres of many planets and moons, such as Mars (Pätzold et al., 2016; Petricca et al., 2021; Withers & Moore, 2020; Withers et al., 2014), Venus (Ando et al., 2020; Bocanegra-Bahamón et al., 2019; Fjeldbo et al., 1971; Gramigna et al., 2023; Häusler et al., 2006), and Titan (Schinder et al., 2011b, 2020). It directly returns the refractivity profile from the bending angle and the impact parameter of the radio signal. In one-way occultation experiments, both the bending angle and the impact parameter can be deduced from the frequency shift and the geometry of the occultation at every time step (Withers et al., 2014). For two-way occultation experiments, two bending angles and two impact parameters—one bending angle and one impact parameter per uplink and downlink—must be derived from a single frequency shift and from the geometry of the occultation. All in all, these represent three independent equations for four unknowns (Withers & Moore, 2020) and a fourth relationship is thus needed to find a solution. However, an exact fourth condition does not exist, so approximate relationships are usually employed to close the system (Jenkins et al., 1994; Lipa & Tyler, 1979; Withers & Moore, 2020).

This issue can be overcome by using ray-tracing algorithms, which do not suffer this problem (Grandin et al., 2014; Schinder et al., 2015, 2020). In fact, this method allows one to determine at the same time and sequentially the (uplink and downlink) bending angles and the refractivity gradient inside each shell of a layered atmosphere (starting from the top of the atmosphere). This is done by solving a square system of nonlinear equations enforcing that (a) each individual (uplink and downlink) ray must point in the direction of the receiving antenna after leaving the atmosphere, (b) the computed frequency at reception must match exactly the observed one. Hence, within this approach, no approximation is required. Moreover, the ray-tracing method is more general than the Abel transform since it can also be used to study atmospheric properties of non-spherical bodies. Indeed, it was successfully applied to the analysis of oblate atmospheres, such as for the Cassini-Saturn occultation experiments (Schinder et al., 2011a, 2015), or the occultations of Voyager spacecraft by Jupiter (Lindal et al., 1981), Saturn (Lindal et al., 1985), Uranus (Lindal et al., 1987), and Neptune (Lindal, 1992). Numerical ray-tracing approach usually requires very significant computational time, to such an extent that the resolution of the atmospheric properties' profiles must be reduced to speed up the computation (Schinder et al., 2020). Moreover, for two-way occultation data processing, the numerical ray-tracing usually focuses on the neutral part of the atmospheres (Schinder et al., 2015, 2020). To our knowledge, a ray-tracing method for two-way ionospheric occultations does not exist, even though, in general, numerical ray-tracing could be applied for two-way ionospheric occultations. Instead, approaches commonly used in this case are based on the Abel transform (Dalba & Withers, 2019; Withers & Moore, 2020).

Bourgoin et al. (2019) have recently found an analytical solution to the differential equations describing the path of the signal when it crosses a spherically symmetric atmosphere consisting of a single layer. Assuming that the refractivity varies linearly with the gravitational potential of the planet/satellite hosting the occulting atmosphere, the authors showed that the light ray trajectory is exactly a hyperbola within the approximation of geometrical optics. In addition, they re-expressed the geometrical optics equations as a set of osculating equations describing the way first integrals of the light path are modified when considering more complete refractivity profiles. In the current work, we exploit the hyperbola solution to establish a new ray-tracing algorithm able to process two-way
Up to now, numerous occultation experiments have been successfully conducted, from the occultation of Mariner IV by Mars on 15 July 1965 (Kliore et al., 1965) to the more recent Ganymede-Juno occultation on 7 June 2021 (Buccino et al., 2022). Comprehensive reviews of objectives, methods, and results of some of the radio occultation experiments carried out during the last 50 years can be found in Fjeldbo and Eshleman (1968), Fjeldbo et al. (1971), Lindal (1992), Jenkins et al. (1994), Pätzold et al. (2005, 2007), Schinder et al. (2011a, 2011b), Hinson et al. (2017), Withers et al. (2018), and Dalba and Withers (2019).

Distinctions must be made between one-way and two-way occultation experiments. In one-way downlink mode, the onboard frequency source, for example, an Ultra Stable Oscillator (USO), provides a stable frequency reference for the transmitted signal to the Earth. In a one-way uplink configuration, an Earth-based station transmits a reference signal to the spacecraft, which records and re-transmits it back to Earth via telemetry. Uplink occultations benefit from the fact that they can provide a higher signal-to-noise ratio than typical downlink experiments. Instead, for two-way occultations, the frequency of an uplink signal (from the ground station to the spacecraft's antenna) is used as a value of reference for the downlink signal (from the spacecraft back to the Earth).

One of the most widely used methods to process one-way radio occultation data for spherically symmetric atmospheres is based on Abel transform algorithms (Fjeldbo et al., 1971; Phinney & Anderson, 1968). This method was applied to the atmospheres of many planets and moons, such as Mars (Pätzold et al., 2016; Petricca et al., 2021; Withers & Moore, 2020; Withers et al., 2014), Venus (Ando et al., 2020; Bocanegra-Bahamón et al., 2019; Fjeldbo et al., 1971; Gramigna et al., 2023; Häusler et al., 2006), and Titan (Schinder et al., 2011b, 2020). It directly returns the refractivity profile from the bending angle and the impact parameter of the radio signal. In one-way occultation experiments, both the bending angle and the impact parameter can be deduced from the frequency shift and the geometry of the occultation at every time step (Withers et al., 2014). For two-way occultation experiments, two bending angles and two impact parameters—one bending angle and one impact parameter per uplink and downlink—must be derived from a single frequency shift and from the geometry of the occultation. All in all, these represent three independent equations for four unknowns (Withers & Moore, 2020) and a fourth relationship is thus needed to find a solution. However, an exact fourth condition does not exist, so approximate relationships are usually employed to close the system (Jenkins et al., 1994; Lipa & Tyler, 1979; Withers & Moore, 2020).

This issue can be overcome by using ray-tracing algorithms, which do not suffer this problem (Grandin et al., 2014; Schinder et al., 2015, 2020). In fact, this method allows one to determine at the same time and sequentially the (uplink and downlink) bending angles and the refractivity gradient inside each shell of a layered atmosphere (starting from the top of the atmosphere). This is done by solving a square system of nonlinear equations enforcing that (a) each individual (uplink and downlink) ray must point in the direction of the receiving antenna after leaving the atmosphere, (b) the computed frequency at reception must match exactly the observed one. Hence, within this approach, no approximation is required. Moreover, the ray-tracing method is more general than the Abel transform since it can also be used to study atmospheric properties of non-spherical bodies. Indeed, it was successfully applied to the analysis of oblate atmospheres, such as for the Cassini-Saturn occultation experiments (Schinder et al., 2011a, 2015), or the occultations of Voyager spacecraft by Jupiter (Lindal et al., 1981), Saturn (Lindal et al., 1985), Uranus (Lindal et al., 1987), and Neptune (Lindal, 1992). Numerical ray-tracing approach usually requires very significant computational time, to such an extent that the resolution of the atmospheric properties' profiles must be reduced to speed up the computation (Schinder et al., 2020). Moreover, for two-way occultation data processing, the numerical ray-tracing usually focuses on the neutral part of the atmospheres (Schinder et al., 2015, 2020). To our knowledge, a ray-tracing method for two-way ionospheric occultations does not exist, even though, in general, numerical ray-tracing could be applied for two-way ionospheric occultations. Instead, approaches commonly used in this case are based on the Abel transform (Dalba & Withers, 2019; Withers & Moore, 2020).

Bourgoin et al. (2019) have recently found an analytical solution to the differential equations describing the path of the signal when it crosses a spherically symmetric atmosphere consisting of a single layer. Assuming that the refractivity varies linearly with the gravitational potential of the planet/satellite hosting the occulting atmosphere, the authors showed that the light ray trajectory is exactly a hyperbola within the approximation of geometrical optics. In addition, they re-expressed the geometrical optics equations as a set of osculating equations describing the way first integrals of the light path are modified when considering more complete refractivity profiles. In the current work, we exploit the hyperbola solution to establish a new ray-tracing algorithm able to process two-way
occultation data by spherically symmetric atmospheres. We extend the idealized scenario of a single-layer atmosphere presented by Bourgoin et al. (2019) to operational capabilities. In particular, we provide a model which better represents a realistic scenario. We model the occulting atmosphere as a stack of spherical layers where we assume that the refractivity in each shell evolves linearly with the gravitational potential. This allows us to chain hyperbolic trajectories inside each layer and to obtain a closed-form solution for the total bending experienced by the radio ray. The benefit of this approach is twofold, first, because of its simplicity, it lets us study the physical properties of the ionosphere and neutral atmosphere simultaneously for both one-way and two-way experiments. Then, because the analytical solutions avoid the hassle of integrating numerically the equations for the radio ray at each time step, the method provides better performances than numerical ray-tracing in terms of computational time. It is, in fact, hundreds of times faster than a numerical approach based on the method described in Schinder et al. (2015). However, in its current form, our algorithm is limited to the data processing of radio occultations by spherically symmetric atmospheres. Therefore, it cannot be applied to the analysis of oblate atmospheres such as those of gas giants, that is, Jupiter and Saturn. Nevertheless, this could be changed by complementing the method with first-order solutions derived by Bourgoin et al. (2019) that go beyond the linear approximation in the refractivity profile of each layer.

The paper is organized as follows. In Section 2, we introduce the mathematical modeling that is used hereafter in Section 4 to process radio occultations data. We first start by recalling some results of Bourgoin et al. (2019) that are important for the present study. For instance, we present the hyperbolic solution which is derived from the linear evolution of the refractivity profile with the gravitational potential of the occulting planet. The solution is then adapted to the case of an atmosphere made with a stack of concentric layers. In Section 3, we introduce the main observable quantity in radio occultation experiments, namely the Doppler frequency shift. We specialize the discussion to one-way, and then, two-way experiments. We outline the procedure to process radio occultations' Doppler data up to the inference of the vertical profile of physical properties. In Section 4, we assess the precision of our new approach by comparing the retrieved vertical profiles to those already published in the literature from numerical ray-tracing and Abel inversion algorithms. Finally, we give our conclusions in Section 5.

## 2. Mathematical Modeling of a Radio Signal Path

This paper presents a ray-tracing approach for analyzing data from a radio occultation experiment. The main idea of such a method is to consider the atmosphere as composed of a series of concentric shells, and the properties of the gas filling each shell are computed layer by layer with an iterative procedure. This method is based on the geometric optics approximation, which assumes that each sky frequency is related to a ray path traversing the atmosphere. For one-way experiments, it is assumed that the spacecraft position at the time of signal transmission and the ground antenna's position at the time of reception are known with respect to the occulting body. In the case of two-way experiments, we need to know the positions of the ground antenna at the time of uplink signal transmission and downlink signal reception, as well as the spacecraft position at the time of signal reception and retransmission. The calculation of the ray path through the atmosphere is performed by integrating the eikonal equation. Finally, the problem is translated into solving a system of nonlinear equations: we look for the ray's initial direction and the refractivity gradient inside the shells to ensure that the ray direction points toward the receiving antenna upon leaving the atmosphere (how its boundary is defined will be discussed in Section 3) and that the signal reaches the ground antenna with the recorded sky frequency.

All these steps, including the methodology used to build each layer inside the atmosphere, will be discussed in detail throughout the paper, while this section focuses on the analytical determination of the ray trajectory through the atmosphere. In particular, analytical expressions for the position and direction of the ray at the exit from the atmosphere are derived as a function of the initial position of the transmitting antenna, the ray's initial direction, and the refractivity gradient inside each shell. These results will be used in Section 3 for the analysis of both one-way and two-way occultation data.

Let us assume that an antenna, labeled A, is transmitting a radio signal at time $t_A$. The signal is then received at time $t_B$ (with $t_B > t_A$) by another antenna, labeled B (see Figure 1). Moreover, we assume that the receiving antenna is being occulted by a planetary atmosphere as viewed from the transmitting antenna. Therefore, the signal passes through the occulting atmosphere and experiences refraction at two different levels, through the change of the phase velocity and the bending of the trajectory toward regions of higher refractivity.
occultation data by spherically symmetric atmospheres. We extend the idealized scenario of a single-layer atmosphere presented by Bourgoin et al. (2019) to operational capabilities. In particular, we provide a model which better represents a realistic scenario. We model the occulting atmosphere as a stack of spherical layers where we assume that the refractivity in each shell evolves linearly with the gravitational potential. This allows us to chain hyperbolic trajectories inside each layer and to obtain a closed-form solution for the total bending experienced by the radio ray. The benefit of this approach is twofold, first, because of its simplicity, it lets us study the physical properties of the ionosphere and neutral atmosphere simultaneously for both one-way and two-way experiments. Then, because the analytical solutions avoid the hassle of integrating numerically the equations for the radio ray at each time step, the method provides better performances than numerical ray-tracing in terms of computational time. It is, in fact, hundreds of times faster than a numerical approach based on the method described in Schinder et al. (2015). However, in its current form, our algorithm is limited to the data processing of radio occultations by spherically symmetric atmospheres. Therefore, it cannot be applied to the analysis of oblate atmospheres such as those of gas giants, that is, Jupiter and Saturn. Nevertheless, this could be changed by complementing the method with first-order solutions derived by Bourgoin et al. (2019) that go beyond the linear approximation in the refractivity profile of each layer.

The paper is organized as follows. In Section 2, we introduce the mathematical modeling that is used hereafter in Section 4 to process radio occultations data. We first start by recalling some results of Bourgoin et al. (2019) that are important for the present study. For instance, we present the hyperbolic solution which is derived from the linear evolution of the refractivity profile with the gravitational potential of the occulting planet. The solution is then adapted to the case of an atmosphere made with a stack of concentric layers. In Section 3, we introduce the main observable quantity in radio occultation experiments, namely the Doppler frequency shift. We specialize the discussion to one-way, and then, two-way experiments. We outline the procedure to process radio occultations' Doppler data up to the inference of the vertical profile of physical properties. In Section 4, we assess the precision of our new approach by comparing the retrieved vertical profiles to those already published in the literature from numerical ray-tracing and Abel inversion algorithms. Finally, we give our conclusions in Section 5.

## 2. Mathematical Modeling of a Radio Signal Path

This paper presents a ray-tracing approach for analyzing data from a radio occultation experiment. The main idea of such a method is to consider the atmosphere as composed of a series of concentric shells, and the properties of the gas filling each shell are computed layer by layer with an iterative procedure. This method is based on the geometric optics approximation, which assumes that each sky frequency is related to a ray path traversing the atmosphere. For one-way experiments, it is assumed that the spacecraft position at the time of signal transmission and the ground antenna's position at the time of reception are known with respect to the occulting body. In the case of two-way experiments, we need to know the positions of the ground antenna at the time of uplink signal transmission and downlink signal reception, as well as the spacecraft position at the time of signal reception and retransmission. The calculation of the ray path through the atmosphere is performed by integrating the eikonal equation. Finally, the problem is translated into solving a system of nonlinear equations: we look for the ray's initial direction and the refractivity gradient inside the shells to ensure that the ray direction points toward the receiving antenna upon leaving the atmosphere (how its boundary is defined will be discussed in Section 3) and that the signal reaches the ground antenna with the recorded sky frequency.

All these steps, including the methodology used to build each layer inside the atmosphere, will be discussed in detail throughout the paper, while this section focuses on the analytical determination of the ray trajectory through the atmosphere. In particular, analytical expressions for the position and direction of the ray at the exit from the atmosphere are derived as a function of the initial position of the transmitting antenna, the ray's initial direction, and the refractivity gradient inside each shell. These results will be used in Section 3 for the analysis of both one-way and two-way occultation data.

Let us assume that an antenna, labeled A, is transmitting a radio signal at time $t_A$. The signal is then received at time $t_B$ (with $t_B > t_A$) by another antenna, labeled B (see Figure 1). Moreover, we assume that the receiving antenna is being occulted by a planetary atmosphere as viewed from the transmitting antenna. Therefore, the signal passes through the occulting atmosphere and experiences refraction at two different levels, through the change of the phase velocity and the bending of the trajectory toward regions of higher refractivity.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic representation of a radio signal path between the emitter A and the receiver B. The signal trajectory is depicted in red in the reference frame $R(O; x, y, z)$ whose vector basis is $(\hat{e}_x, \hat{e}_y, \hat{e}_z)$. The angle $\omega$ represents the argument of the pericenter for the light trajectory. The angles $\nu$ and $\theta$ are the true anomaly and polar angle, respectively. The angles $\psi$ and $\epsilon$ are the argument of the bending and the total bending of the light path, respectively. For a spherically symmetric atmosphere, the light ray propagates into a fixed plane so that only two components are actually needed to locate a point in the $(x, y)$-plane, namely $(r, \theta)$.</figcaption></figure>

In this work, we place ourselves within the geometric optics approximation, where the wavelength of the radio signal is assumed to be small when compared to the characteristic size of the medium. In this framework, light rays are the curves whose tangents at each point coincide with the direction of propagation of the electromagnetic waves. The fundamental equations governing the propagation of light rays are given by the vector eikonal equations (Born & Wolf, 1980)

$$ \frac{dt}{ds} = \frac{n}{c}, \quad \frac{d\mathbf{r}}{ds} = \hat{\mathbf{s}}, \quad \frac{d\hat{\mathbf{s}}}{ds} = \frac{\nabla n - \hat{\mathbf{s}}(\hat{\mathbf{s}} \cdot \nabla n)}{n}, \qquad (1) $$

where $n$ is the refractive index, $\mathbf{r}$ is the position vector of a point P along the ray path, $\hat{\mathbf{s}}$ is the unit vector representing the direction of the ray at P (see Figure 1), $s$ represents the curvilinear abscissa that labels each point along the ray, and $t$ is the corresponding elapsed time.

In addition to the differential equations, the trajectory of the light path is constrained by the following boundary conditions:

$$ t(s = s_A) = t_A, \quad \mathbf{r}(t_A) = \mathbf{r}_A, \quad \hat{\mathbf{s}}(t_A) = \hat{\mathbf{s}}_A, \qquad (2) $$

$$ t(s = s_B) = t_B, \quad \mathbf{r}(t_B) = \mathbf{r}_B, \quad \hat{\mathbf{s}}(t_B) = \hat{\mathbf{s}}_B, \qquad (3) $$

where $\hat{\mathbf{s}}_A$ and $\hat{\mathbf{s}}_B$ are the unit vectors representing the direction of the ray at the level of the emitter and receiver, respectively. The quantities $\mathbf{r}_A$ and $\mathbf{r}_B$ are the position of the emitter and receiver, respectively. Finally, as mentioned previously, $t_A$ and $t_B$ represent the time of emission and reception at the curvilinear abscissa $s_A$ and $s_B$, respectively (hereafter, we set $s_A = 0$ for convenience).

In the following discussion we will assume that the planetary atmosphere extends up to an upper boundary, above which is vacuum. For the trivial case of a ray propagating in a vacuum, the index of refraction is $n = 1$ all along the light trajectory so Equation 1 immediately lead to $d\hat{\mathbf{s}}/ds = 0$ and, therefore, $\hat{\mathbf{s}}_A = \hat{\mathbf{s}}_B = \hat{\mathbf{s}}_{vac}$ with

$$ c(t_B - t_A) = s_B = \|\mathbf{r}_B - \mathbf{r}_A\|, \quad \hat{\mathbf{s}}_{vac} = \frac{\mathbf{r}_B - \mathbf{r}_A}{s_B}. \qquad (4) $$

For a light propagating in an optical medium (i.e., with $n \neq 0$), the presence of a refractive index gradient causes the signal to depart from a rectilinear trajectory, so $\hat{\mathbf{s}}_A \neq \hat{\mathbf{s}}_B$. Consequently, for the radio signal to reach the receiving antenna B, the initial direction of the ray must be different from the vacuum direction $\hat{\mathbf{s}}_{vac}$.

Moreover, for radio occultation experiments where the receiver and emitter are considered to be outside of the optical medium, we may consider that $n = 1$ outside the planetary atmosphere, meaning that $\hat{\mathbf{s}}_A$ coincides with the direction of the ray when it enters into the atmosphere whereas $\hat{\mathbf{s}}_B$ coincides with the direction of the ray when it leaves the atmosphere on his way to the receiver. We thus have the additional geometrical condition on the direction of the exiting ray:

$$ \hat{\mathbf{s}}_B - \frac{\mathbf{r}_B - \mathbf{r}_E}{\|\mathbf{r}_B - \mathbf{r}_E\|} = 0, \qquad (5) $$

where $\mathbf{r}_E$ is the position of the ray when it leaves the atmosphere on his way to the receiver B.

In general, the ray path can be determined by solving numerically Equation 1 with conditions 2, 3, and 5. However, in the particular case of a signal passing through a single-shell spherically symmetric atmosphere, Bourgoin et al. (2019) have found an exact solution whose formulation is recalled in the next subsection.

## 2.1. Ray Path Through a Single Shell

This subsection reports all the relationships required to get the ray path through the atmosphere starting from the knowledge of the occultation geometry and the properties of the atmospheric medium. In particular, this subsection leverages the analytical solution to Equation 1 obtained by Bourgoin et al. (2019), which lets us avoid
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic representation of a radio signal path between the emitter A and the receiver B. The signal trajectory is depicted in red in the reference frame $R(O; x, y, z)$ whose vector basis is $(\hat{e}_x, \hat{e}_y, \hat{e}_z)$. The angle $\omega$ represents the argument of the pericenter for the light trajectory. The angles $\nu$ and $\theta$ are the true anomaly and polar angle, respectively. The angles $\psi$ and $\epsilon$ are the argument of the bending and the total bending of the light path, respectively. For a spherically symmetric atmosphere, the light ray propagates into a fixed plane so that only two components are actually needed to locate a point in the $(x, y)$-plane, namely $(r, \theta)$.</figcaption></figure>

In this work, we place ourselves within the geometric optics approximation, where the wavelength of the radio signal is assumed to be small when compared to the characteristic size of the medium. In this framework, light rays are the curves whose tangents at each point coincide with the direction of propagation of the electromagnetic waves. The fundamental equations governing the propagation of light rays are given by the vector eikonal equations (Born & Wolf, 1980)

$$ \frac{dt}{ds} = \frac{n}{c}, \quad \frac{d\mathbf{r}}{ds} = \hat{\mathbf{s}}, \quad \frac{d\hat{\mathbf{s}}}{ds} = \frac{\nabla n - \hat{\mathbf{s}}(\hat{\mathbf{s}} \cdot \nabla n)}{n}, \qquad (1) $$

where $n$ is the refractive index, $\mathbf{r}$ is the position vector of a point P along the ray path, $\hat{\mathbf{s}}$ is the unit vector representing the direction of the ray at P (see Figure 1), $s$ represents the curvilinear abscissa that labels each point along the ray, and $t$ is the corresponding elapsed time.

In addition to the differential equations, the trajectory of the light path is constrained by the following boundary conditions:

$$ t(s = s_A) = t_A, \quad \mathbf{r}(t_A) = \mathbf{r}_A, \quad \hat{\mathbf{s}}(t_A) = \hat{\mathbf{s}}_A, \qquad (2) $$

$$ t(s = s_B) = t_B, \quad \mathbf{r}(t_B) = \mathbf{r}_B, \quad \hat{\mathbf{s}}(t_B) = \hat{\mathbf{s}}_B, \qquad (3) $$

where $\hat{\mathbf{s}}_A$ and $\hat{\mathbf{s}}_B$ are the unit vectors representing the direction of the ray at the level of the emitter and receiver, respectively. The quantities $\mathbf{r}_A$ and $\mathbf{r}_B$ are the position of the emitter and receiver, respectively. Finally, as mentioned previously, $t_A$ and $t_B$ represent the time of emission and reception at the curvilinear abscissa $s_A$ and $s_B$, respectively (hereafter, we set $s_A = 0$ for convenience).

In the following discussion we will assume that the planetary atmosphere extends up to an upper boundary, above which is vacuum. For the trivial case of a ray propagating in a vacuum, the index of refraction is $n = 1$ all along the light trajectory so Equation 1 immediately lead to $d\hat{\mathbf{s}}/ds = 0$ and, therefore, $\hat{\mathbf{s}}_A = \hat{\mathbf{s}}_B = \hat{\mathbf{s}}_{\text{vac}}$ with

$$ c(t_B - t_A) = s_B = \|\mathbf{r}_B - \mathbf{r}_A\|, \quad \hat{\mathbf{s}}_{\text{vac}} = \frac{\mathbf{r}_B - \mathbf{r}_A}{s_B}. \qquad (4) $$

For a light propagating in an optical medium (i.e., with $n \neq 0$), the presence of a refractive index gradient causes the signal to depart from a rectilinear trajectory, so $\hat{\mathbf{s}}_A \neq \hat{\mathbf{s}}_B$. Consequently, for the radio signal to reach the receiving antenna B, the initial direction of the ray must be different from the vacuum direction $\hat{\mathbf{s}}_{\text{vac}}$.

Moreover, for radio occultation experiments where the receiver and emitter are considered to be outside of the optical medium, we may consider that $n = 1$ outside the planetary atmosphere, meaning that $\hat{\mathbf{s}}_A$ coincides with the direction of the ray when it enters into the atmosphere whereas $\hat{\mathbf{s}}_B$ coincides with the direction of the ray when it leaves the atmosphere on his way to the receiver. We thus have the additional geometrical condition on the direction of the exiting ray:

$$ \hat{\mathbf{s}}_B - \frac{\mathbf{r}_B - \mathbf{r}_E}{\|\mathbf{r}_B - \mathbf{r}_E\|} = 0, \qquad (5) $$

where $\mathbf{r}_E$ is the position of the ray when it leaves the atmosphere on his way to the receiver B.

In general, the ray path can be determined by solving numerically Equation 1 with conditions 2, 3, and 5. However, in the particular case of a signal passing through a single-shell spherically symmetric atmosphere, Bourgoin et al. (2019) have found an exact solution whose formulation is recalled in the next subsection.

## 2.1. Ray Path Through a Single Shell

This subsection reports all the relationships required to get the ray path through the atmosphere starting from the knowledge of the occultation geometry and the properties of the atmospheric medium. In particular, this subsection leverages the analytical solution to Equation 1 obtained by Bourgoin et al. (2019), which lets us avoid
the numerical integration of the eikonal equation, but it is based on the assumption of a single-shell atmosphere. In what follows, the position of a generic point along the ray trajectory is referred to the center-of-mass *O* of the occulting body.

For a spherically symmetric atmosphere the refractive index is a function of the radial distance only, namely $n = n(r)$ with $r = \|\mathbf{r}\|$. This assumption can be valid both for planetary ionospheres (Withers et al., 2014) and neutral atmospheres (Schinder et al., 2020). In particular, in the latter case, within the hydrostatic equilibrium assumption and considering a non-rotating fluid body, the surface of constant index of refraction $n$ coincides with the surface of constant gravitational potential $\Phi$ (Bourgoin et al., 2019; Poisson & Will, 2014), so $n$ may eventually be expressed as a function of $\Phi$ (Schinder et al., 2015). Moreover, in the spherically symmetric case, $\Phi$ is the Newtonian potential, that is, $\Phi(r) = -\mu/r$ where $\mu$ is the gravitational parameter of the occulting body. In principle, extended versions of this method may be developed for oblate objects, but here we focus exclusively on spherically symmetric atmospheres.

Within the spherical symmetry assumption, the impact parameter vector $\mathbf{K} = \mathbf{r} \times n\hat{\mathbf{s}}$ is constant all along the light path (Born & Wolf, 1980). This means the signal propagates into a plane perpendicular to $\mathbf{K}$. It is thus convenient to introduce a reference frame $\mathcal{R}(O; x, y, z)$ whose vector basis $(\hat{e}_x, \hat{e}_y, \hat{e}_z)$ is defined such as

$$ \hat{e}_x = -\hat{s}_{\text{vac}}, \quad \hat{e}_y = \hat{e}_z \times \hat{e}_x, \quad \hat{e}_z = \frac{\mathbf{K}}{K}, \qquad (6) $$

where $K = \|\mathbf{K}\|$. Thus, the ray path lies in the (x, y)-plane; see Figure 1.

If we assume that the refractive index is a linear function of the gravitational potential $\Phi(r) = -\mu/r$, viz.

$$ n(r) = \eta + \frac{\alpha\mu}{r}, \qquad (7) $$

then the signal path is a hyperbola (Bourgoin et al., 2019). In Equation 7, $\eta$ is the value of $n$ when $r$ goes to infinity, and $\alpha = -dn/d\Phi$ is the (constant) refractive index gradient. Therefore, the ray trajectory is mathematically given by

$$ r(v) = \frac{p}{1 + e \cos \kappa v} \quad \text{with} \quad \kappa = \sqrt{1 - \frac{1}{e^2}}, \qquad (8) $$

where $e$ is the eccentricity and $p$ is the semi-latus rectum of the hyperbolic path (we consider non-periodic solutions, i.e., $e > 1$, and therefore $0 < \kappa < 1$) and whose expressions are given by

$$ e = \frac{K}{\alpha\mu}, \quad p = \frac{\alpha\mu}{\eta}(e^2 - 1). \qquad (9) $$

The variable $\nu$ in Equation 8 is the true anomaly, that is, the angle measured counterclockwise from the direction of ray pericenter to $\mathbf{r}$ (see Figure 1). In the context of the radio occultation experiments, the argument of the refractive bending $\psi$, which is defined as the angle between the y-axis and the tangent to the ray, is an important angle for the determination of the Doppler effect. It can be computed from the true anomaly such as

$$ \psi(\nu) = \omega + \nu - 2 \arctan \left[ \tanh \left( \frac{F(\nu)}{2} \right) \right], \qquad (10) $$

where $F$ is the hyperbolic anomaly:

$$ \cosh F(\nu) = \frac{e + \cos \kappa \nu}{1 + e \cos \kappa \nu}, \quad \sinh F(\nu) = \frac{\sqrt{e^2 - 1} \sin \kappa \nu}{1 + e \cos \kappa \nu}, \qquad (11) $$

and $\omega$ is the argument of pericenter, namely the constant angle between the x-axis and the direction of the ray pericenter. The total refractive bending $\varepsilon$ can be obtained by evaluating the change in $\psi$ between the position of the transmitting and receiving antenna, $\varepsilon = \psi(\nu_B) - \psi(\nu_A)$, with $\nu_A$ and $\nu_B$ the true anomalies at the location of the transmitter and receiver, respectively. Note that, within this manuscript, $\varepsilon$ is defined as positive when the direction vector $\hat{s}$ rotates in the counterclockwise direction. It is negative otherwise. Also, in the case of no bending, $\hat{s}_A$ and $\hat{s}_B$ coincide with $\hat{s}_{\text{vac}}$, and the value of $\psi$ remains unchanged between the position of A and B and equal to 90 deg.
the numerical integration of the eikonal equation, but it is based on the assumption of a single-shell atmosphere. In what follows, the position of a generic point along the ray trajectory is referred to the center-of-mass *O* of the occulting body.

For a spherically symmetric atmosphere the refractive index is a function of the radial distance only, namely $n = n(r)$ with $r = \|\mathbf{r}\|$. This assumption can be valid both for planetary ionospheres (Withers et al., 2014) and neutral atmospheres (Schinder et al., 2020). In particular, in the latter case, within the hydrostatic equilibrium assumption and considering a non-rotating fluid body, the surface of constant index of refraction $n$ coincides with the surface of constant gravitational potential $\Phi$ (Bourgoin et al., 2019; Poisson & Will, 2014), so $n$ may eventually be expressed as a function of $\Phi$ (Schinder et al., 2015). Moreover, in the spherically symmetric case, $\Phi$ is the Newtonian potential, that is, $\Phi(r) = -\mu/r$ where $\mu$ is the gravitational parameter of the occulting body. In principle, extended versions of this method may be developed for oblate objects, but here we focus exclusively on spherically symmetric atmospheres.

Within the spherical symmetry assumption, the impact parameter vector $\mathbf{K} = \mathbf{r} \times n\hat{\mathbf{s}}$ is constant all along the light path (Born & Wolf, 1980). This means the signal propagates into a plane perpendicular to $\mathbf{K}$. It is thus convenient to introduce a reference frame $\mathcal{R}(O; x, y, z)$ whose vector basis $(\hat{e}_x, \hat{e}_y, \hat{e}_z)$ is defined such as

$$ \hat{e}_x = -\hat{s}_{\text{vac}}, \quad \hat{e}_y = \hat{e}_z \times \hat{e}_x, \quad \hat{e}_z = \frac{\mathbf{K}}{K}, \qquad (6) $$

where $K = \|\mathbf{K}\|$. Thus, the ray path lies in the (x, y)-plane; see Figure 1.

If we assume that the refractive index is a linear function of the gravitational potential $\Phi(r) = -\mu/r$, viz.

$$ n(r) = \eta + \frac{\alpha\mu}{r}, \qquad (7) $$

then the signal path is a hyperbola (Bourgoin et al., 2019). In Equation 7, $\eta$ is the value of $n$ when $r$ goes to infinity, and $\alpha = -dn/d\Phi$ is the (constant) refractive index gradient. Therefore, the ray trajectory is mathematically given by

$$ r(v) = \frac{p}{1 + e \cos \kappa v} \quad \text{with} \quad \kappa = \sqrt{1 - \frac{1}{e^2}}, \qquad (8) $$

where $e$ is the eccentricity and $p$ is the semi-latus rectum of the hyperbolic path (we consider non-periodic solutions, i.e., $e > 1$, and therefore $0 < \kappa < 1$) and whose expressions are given by

$$ e = \frac{K}{\alpha\mu}, \quad p = \frac{\alpha\mu}{\eta}(e^2 - 1). \qquad (9) $$

The variable $\nu$ in Equation 8 is the true anomaly, that is, the angle measured counterclockwise from the direction of ray pericenter to $\mathbf{r}$ (see Figure 1). In the context of the radio occultation experiments, the argument of the refractive bending $\psi$, which is defined as the angle between the y-axis and the tangent to the ray, is an important angle for the determination of the Doppler effect. It can be computed from the true anomaly such as

$$ \psi(\nu) = \omega + \nu - 2 \arctan \left[ \tanh \left( \frac{F(\nu)}{2} \right) \right], \qquad (10) $$

where $F$ is the hyperbolic anomaly:

$$ \cosh F(\nu) = \frac{e + \cos \kappa \nu}{1 + e \cos \kappa \nu}, \quad \sinh F(\nu) = \frac{\sqrt{e^2 - 1} \sin \kappa \nu}{1 + e \cos \kappa \nu}, \qquad (11) $$

and $\omega$ is the argument of pericenter, namely the constant angle between the x-axis and the direction of the ray pericenter. The total refractive bending $\varepsilon$ can be obtained by evaluating the change in $\psi$ between the position of the transmitting and receiving antenna, $\varepsilon = \psi(\nu_B) - \psi(\nu_A)$, with $\nu_A$ and $\nu_B$ the true anomalies at the location of the transmitter and receiver, respectively. Note that, within this manuscript, $\varepsilon$ is defined as positive when the direction vector $\hat{s}$ rotates in the counterclockwise direction. It is negative otherwise. Also, in the case of no bending, $\hat{s}_A$ and $\hat{s}_B$ coincide with $\hat{s}_{\text{vac}}$, and the value of $\psi$ remains unchanged between the position of A and B and equal to 90 deg.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left panel: Schematic representation of ray path (red line) through a multilayered atmosphere with $m = 4$. We denote by $\alpha_i$ the gradient of the refractive index in the $i$th shell whose outer radius is $r_i$. Right panel: Within the spherical symmetry assumption, the ray path lies on the (x, y)-plane of $\mathcal{R}$. The atmosphere boundary has a radius $r_1$. The signal enters the atmosphere at $P_1$ whose polar angle is denoted by $\theta_1$. The polar angle at the position of the transmitting antenna is $\theta_A$. The variable $\beta_A$ is the angle between the initial ray direction at $P_1$ and $\hat{e}_x$. Also, the distance between the ray asymptote from $O$, namely the impact parameter, is denoted by $K$. Finally, the polar angle at the point $\vec{E}$ where the signal leaves the atmosphere is $\vartheta_E$.</figcaption></figure>

In Figure 1, we also depicted the polar angle $\vartheta = \omega + \nu$ which is measured counterclockwise from the direction of the x-axis to $\mathbf{r}$.

Despite the availability of an analytical expression for the ray path, a single shell atmospheric model with constant refractive index gradient $\alpha$ is not well suited to describe realistic planetary atmospheres. To achieve better precision, many methods (see e.g., Lindal, 1992; Schinder et al., 2015 and also Section 3.5 of Bourgoin et al. (2019)) proposed to model planetary atmospheres as a series of concentric shells. Indeed, when the thickness of each layer is much smaller than the refractivity scale height, the gradient of the refractive index with respect to the gravitational potential can be considered constant in each shell. Accordingly, the analytical solution reported in this section still holds inside each spherical shell, and the total ray path is composed of a series of hyperbolic arcs patched at the shells' boundaries.

We recall that for the analysis of radio occultation data, we need to compute the state of a ray (position and direction) at the exit point from the atmosphere. A detailed description of the procedure used to find the final state as a function of the initial position of the transmitting antenna, the ray's initial direction, and the refractivity gradient inside each shell is given in the next subsection assuming the ray is traversing a multilayered atmosphere.

## 2.2. Ray Path Through a Multilayered Atmosphere

This subsection provides the analytical expressions used to compute the position and direction vectors of a ray when it leaves the atmosphere boundary. In particular, this subsection extends the results of the previous subsection to the case of a multilayered atmosphere.

Let us assume a spherically symmetric planetary atmosphere that extends up to a maximum radial distance $r_1$ from the occulting body's center-of-mass (for $r \ge r_1$, we impose the vacuum condition, namely $n = 1$). Moreover, the entire atmosphere is assumed composed of $m$ concentric spherical shells, with the $i$th shell filled with a fluid of constant refractive index gradient. Let us denote by $\alpha_i$ the gradient of the refractive index in the $i$th shell. We denote the outer radius of the $i$th shell by $r_i$, so each spherical shell is bounded between two spheres of radius $r_i$ and $r_{i+1}$. Also, if the occulting body has a solid surface, the $m$th shell extends down to the planet radius $R_p$, with $R_p = r_{m+1}$. In the case of a planet without a solid surface or for Venus, in which critical refraction prevents the ray from penetrating inside the lower atmosphere below 40 km (Gramigna et al., 2023), it may be possible to guess an arbitrary minimum radius $r_{m+1}$ below which we know the ray cannot penetrate or simply putting $r_{m+1} = 0$. In Figure 2 to the left, we depicted a light ray trajectory across a layered atmosphere.

Note that, if the transmitting antenna is located outside the planetary atmosphere, the signal initially travels through a vacuum, and the ray path is a straight line to the point $P_1$ where the ray enters into the atmosphere. Therefore, the ray direction at the point $P_1$ coincides with the direction $\hat{s}_A$ when the signal is emitted. In the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Left panel: Schematic representation of ray path (red line) through a multilayered atmosphere with $m = 4$. We denote by $\alpha_i$ the gradient of the refractive index in the $i$th shell whose outer radius is $r_i$. Right panel: Within the spherical symmetry assumption, the ray path lies on the (x, y)-plane of $\mathcal{R}$. The atmosphere boundary has a radius $r_1$. The signal enters the atmosphere at $P_1$ whose polar angle is denoted by $\theta_1$. The polar angle at the position of the transmitting antenna is $\theta_A$. The variable $\beta_A$ is the angle between the initial ray direction at $P_1$ and $\hat{e}_x$. Also, the distance between the ray asymptote from $O$, namely the impact parameter, is denoted by $K$. Finally, the polar angle at the point $\vec{E}$ where the signal leaves the atmosphere is $\vartheta_E$.</figcaption></figure>

In Figure 1, we also depicted the polar angle $\vartheta = \omega + \nu$ which is measured counterclockwise from the direction of the x-axis to $\mathbf{r}$.

Despite the availability of an analytical expression for the ray path, a single shell atmospheric model with constant refractive index gradient $\alpha$ is not well suited to describe realistic planetary atmospheres. To achieve better precision, many methods (see e.g., Lindal, 1992; Schinder et al., 2015 and also Section 3.5 of Bourgoin et al. (2019)) proposed to model planetary atmospheres as a series of concentric shells. Indeed, when the thickness of each layer is much smaller than the refractivity scale height, the gradient of the refractive index with respect to the gravitational potential can be considered constant in each shell. Accordingly, the analytical solution reported in this section still holds inside each spherical shell, and the total ray path is composed of a series of hyperbolic arcs patched at the shells' boundaries.

We recall that for the analysis of radio occultation data, we need to compute the state of a ray (position and direction) at the exit point from the atmosphere. A detailed description of the procedure used to find the final state as a function of the initial position of the transmitting antenna, the ray's initial direction, and the refractivity gradient inside each shell is given in the next subsection assuming the ray is traversing a multilayered atmosphere.

## 2.2. Ray Path Through a Multilayered Atmosphere

This subsection provides the analytical expressions used to compute the position and direction vectors of a ray when it leaves the atmosphere boundary. In particular, this subsection extends the results of the previous subsection to the case of a multilayered atmosphere.

Let us assume a spherically symmetric planetary atmosphere that extends up to a maximum radial distance $r_1$ from the occulting body's center-of-mass (for $r \ge r_1$, we impose the vacuum condition, namely $n = 1$). Moreover, the entire atmosphere is assumed composed of $m$ concentric spherical shells, with the $i$th shell filled with a fluid of constant refractive index gradient. Let us denote by $\alpha_i$ the gradient of the refractive index in the $i$th shell. We denote the outer radius of the $i$th shell by $r_i$, so each spherical shell is bounded between two spheres of radius $r_i$ and $r_{i+1}$. Also, if the occulting body has a solid surface, the $m$th shell extends down to the planet radius $R_p$, with $R_p = r_{m+1}$. In the case of a planet without a solid surface or for Venus, in which critical refraction prevents the ray from penetrating inside the lower atmosphere below 40 km (Gramigna et al., 2023), it may be possible to guess an arbitrary minimum radius $r_{m+1}$ below which we know the ray cannot penetrate or simply putting $r_{m+1} = 0$. In Figure 2 to the left, we depicted a light ray trajectory across a layered atmosphere.

Note that, if the transmitting antenna is located outside the planetary atmosphere, the signal initially travels through a vacuum, and the ray path is a straight line to the point $P_1$ where the ray enters into the atmosphere. Therefore, the ray direction at the point $P_1$ coincides with the direction $\hat{s}_A$ when the signal is emitted. In the
context of radio occultation data processing, we are particularly interested in computing (a) $\mathbf{r}_E = \overrightarrow{OE}$, the position vector of the point E where the ith signal is exiting the atmosphere and (b) $\hat{\mathbf{s}}_E$, the direction of the ith light ray at the same point E. Both are functions of the initial pointing direction $\hat{\mathbf{s}}_A$ and of coefficients $\alpha_i$ with $i = \{1, ..., m\}$.

### 2.2.1. Position of the Exit Point

First, let us determine the vector position of E. To this end, we first introduce the angle $\beta_A$ from $\hat{\mathbf{e}}_x$ to $\hat{\mathbf{s}}_A$ (see Figure 2) such that the components of $\hat{\mathbf{s}}_A$ in $\mathcal{R}$ are defined by

$$ \hat{\mathbf{s}}_A = (\cos\beta_A, \sin\beta_A). \qquad (12) $$

As seen from Figure 2, the polar angle $\theta_1$ from $\hat{\mathbf{e}}_x$ to $\overrightarrow{OP_1}$ can be expressed such as (Fjeldbo et al., 1971)

$$ \theta_1 = \beta_A - \gamma - \frac{\pi}{2}, \qquad (13) $$

where the angle $\gamma$ is given by the geometry of the occultation according to

$$ \gamma = \arctan \left[ \frac{\sqrt{r_1^2 - K^2}}{K} \right]. \qquad (14) $$

The norm of the impact parameter can be expressed as a function of $\beta_A$ and $\theta_A$, the polar angle of the transmitter, according to Fjeldbo et al. (1971)

$$ K = r_A \sin(\beta_A - \theta_A) \qquad (15) $$

when assuming that the transmitting antenna is outside the atmosphere, namely $n(r_A) = 1$.

Considering that the index of refraction is a linear function of the gravitational potential inside the ith shell, we have the relation:

$$ n_i(r) = \eta_i + \frac{\alpha_i \mu}{r}, \quad \forall r \in [r_{i+1}, r_i] \quad \text{with} \quad i = \{1, \dots, m\}, \qquad (16) $$

where for $i > 1$, $\eta_i$ is chosen to ensure a smooth continuation of the index of refraction at the interface between two layers. Under this assumption, Equation 8 holds inside each shell, so the true anomaly $\nu_i$ along the ith hyperbolic arc of the ray path can be calculated as

$$ \nu_i = \begin{cases} +\frac{1}{\kappa_i}\arccos\left[\left(\frac{p_i}{r} - 1\right)\frac{1}{e_i}\right], & \text{if } \mathbf{r} \cdot \hat{\mathbf{s}} > 0, \\ 0, & \text{if } \mathbf{r} \cdot \hat{\mathbf{s}} = 0, \\ -\frac{1}{\kappa_i}\arccos\left[\left(\frac{p_i}{r} - 1\right)\frac{1}{e_i}\right], & \text{if } \mathbf{r} \cdot \hat{\mathbf{s}} < 0, \end{cases} \qquad (17) $$

where $\kappa_i$ has the same form as in Equation 8 and where $e_i$ and $p_i$ are the eccentricity and semi-latus rectum of the light path in the ith layer as shown in Equation 9. Their expressions are given by

$$ \kappa_i = \sqrt{1 - \frac{1}{e_i^2}}, \quad e_i = \frac{K}{\alpha_i \mu}, \quad p_i = \frac{\alpha_i \mu}{\eta_i} (e_i^2 - 1). \qquad (18) $$

In order to determine the components of the vector position of the point E in $\mathcal{R}$, we must determine the polar angle $\theta_E$. Note that, for $i < m$, a ray traverses each shell twice (except for the deepest shell) as shown in Figure 2. Due to the spherical symmetry assumption, the variation in the polar angle $\Delta\theta_i$ in the first half of a shell (where $\mathbf{r} \cdot \hat{\mathbf{s}} < 0$) is the same as that computed in the second half (where $\mathbf{r} \cdot \hat{\mathbf{s}} > 0$). Therefore, we can only focus on the first half of the ray path, up to the ray pericenter, which lies on the axis of symmetry of the ray path. Note that the variation of the polar angle $\Delta\theta_i$ is equal to the variation in the true anomaly $\Delta\nu_i$. As for the first $m-1$ shells, this can be simply computed as

$$ \Delta\nu_i = \nu_i^{(1)} - \nu_i^{(0)} \quad \text{with} \quad i = \{1, \dots, m-1\}, \qquad (19) $$
context of radio occultation data processing, we are particularly interested in computing (a) $\mathbf{r}_E = \overrightarrow{OE}$, the position vector of the point E where the ith signal is exiting the atmosphere and (b) $\hat{\mathbf{s}}_E$, the direction of the ith light ray at the same point E. Both are functions of the initial pointing direction $\hat{\mathbf{s}}_A$ and of coefficients $\alpha_i$ with $i = \{1, ..., m\}$.

### 2.2.1. Position of the Exit Point

First, let us determine the vector position of E. To this end, we first introduce the angle $\beta_A$ from $\hat{\mathbf{e}}_x$ to $\hat{\mathbf{s}}_A$ (see Figure 2) such that the components of $\hat{\mathbf{s}}_A$ in $\mathcal{R}$ are defined by

$$ \hat{\mathbf{s}}_A = (\cos\beta_A, \sin\beta_A). \qquad (12) $$

As seen from Figure 2, the polar angle $\vartheta_1$ from $\hat{\mathbf{e}}_x$ to $\overrightarrow{OP_1}$ can be expressed such as (Fjeldbo et al., 1971)

$$ \vartheta_1 = \beta_A - \gamma - \frac{\pi}{2}, \qquad (13) $$

where the angle $\gamma$ is given by the geometry of the occultation according to

$$ \gamma = \arctan \left[ \frac{\sqrt{r_1^2 - K^2}}{K} \right]. \qquad (14) $$

The norm of the impact parameter can be expressed as a function of $\beta_A$ and $\vartheta_A$, the polar angle of the transmitter, according to Fjeldbo et al. (1971)

$$ K = r_A \sin(\beta_A - \vartheta_A) \qquad (15) $$

when assuming that the transmitting antenna is outside the atmosphere, namely $n(r_A) = 1$.

Considering that the index of refraction is a linear function of the gravitational potential inside the ith shell, we have the relation:

$$ n_i(r) = \eta_i + \frac{\alpha_i \mu}{r}, \quad \forall r \in [r_{i+1}, r_i] \quad \text{with} \quad i = \{1, \dots, m\}, \qquad (16) $$

where for $i > 1$, $\eta_i$ is chosen to ensure a smooth continuation of the index of refraction at the interface between two layers. Under this assumption, Equation 8 holds inside each shell, so the true anomaly $\nu_i$ along the ith hyperbolic arc of the ray path can be calculated as

$$ \nu_i = \begin{cases} +\frac{1}{\kappa_i}\arccos\left[\left(\frac{p_i}{r} - 1\right)\frac{1}{e_i}\right], & \text{if } \mathbf{r} \cdot \hat{\mathbf{s}} > 0, \\ 0, & \text{if } \mathbf{r} \cdot \hat{\mathbf{s}} = 0, \\ -\frac{1}{\kappa_i}\arccos\left[\left(\frac{p_i}{r} - 1\right)\frac{1}{e_i}\right], & \text{if } \mathbf{r} \cdot \hat{\mathbf{s}} < 0, \end{cases} \qquad (17) $$

where $\kappa_i$ has the same form as in Equation 8 and where $e_i$ and $p_i$ are the eccentricity and semi-latus rectum of the light path in the ith layer as shown in Equation 9. Their expressions are given by

$$ \kappa_i = \sqrt{1 - \frac{1}{e_i^2}}, \quad e_i = \frac{K}{\alpha_i \mu}, \quad p_i = \frac{\alpha_i \mu}{\eta_i} (e_i^2 - 1). \qquad (18) $$

In order to determine the components of the vector position of the point E in $\mathcal{R}$, we must determine the polar angle $\vartheta_E$. Note that, for $i < m$, a ray traverses each shell twice (except for the deepest shell) as shown in Figure 2. Due to the spherical symmetry assumption, the variation in the polar angle $\Delta\vartheta_i$ in the first half of a shell (where $\mathbf{r} \cdot \hat{\mathbf{s}} < 0$) is the same as that computed in the second half (where $\mathbf{r} \cdot \hat{\mathbf{s}} > 0$). Therefore, we can only focus on the first half of the ray path, up to the ray pericenter, which lies on the axis of symmetry of the ray path. Note that the variation of the polar angle $\Delta\vartheta_i$ is equal to the variation in the true anomaly $\Delta\nu_i$. As for the first $m-1$ shells, this can be simply computed as

$$ \Delta\nu_i = \nu_i^{(1)} - \nu_i^{(0)} \quad \text{with} \quad i = \{1, \dots, m-1\}, \qquad (19) $$
where $v_i^{(0)} = v_i(r_i)$ and $v_i^{(1)} = v_i(r_{i+1})$ are computed using Equation 17 with negative sign since we are only considering the first half of the ray where $\mathbf{r} \cdot \hat{\mathbf{s}} < 0$.

As for the mth shell (the deepest one), since the true anomaly is zero at the ray pericenter, the change in the true anomaly inside the first half of the deepest shell is simply

$$ \Delta v_m = \frac{1}{\kappa_m} \arccos\left[\left(\frac{p_m}{r_m} - 1\right) \frac{1}{e_m}\right]. \qquad (20) $$

Finally, the value of $\vartheta_E$ is

$$ \vartheta_E = \vartheta_1 + 2 \sum_{i=1}^{m} \Delta v_i, \qquad (21) $$

where $\vartheta_1$ is given by Equation 13. Therefore, the position vector $\mathbf{r}_E$ is given by

$$ \mathbf{r}_E = (r_1 \cos\vartheta_E, r_1 \sin\vartheta_E) \qquad (22) $$

Another important variable that will be useful in the following discussion is the ray pericenter radius $r_K$. It can be easily obtained observing that the closest approach point is located inside the deepest shell, at the point where the true anomaly is equal to zero. Therefore, from Equation 8 it follows that

$$ r_K = \frac{p_m}{1 + e_m}. \qquad (23) $$

### 2.2.2. Direction of the Ray at the Exit Point

In order to compute the total bending angle, which is needed for the determination of the Doppler effect, we need to infer an expression for the direction of the light ray at E, the point where the ray is exiting the atmosphere. As previously mentioned, within the spherical symmetry assumption, we can only focus on the first half of the light ray trajectory. Note that, the bending produced inside the ith shell can be obtained (recalling the definition of $\psi$) by

$$ \varepsilon_i = \psi_i(v_i^{(1)}) - \psi_i(v_i^{(0)}) \quad \text{with} \quad i = \{1, \dots, m-1\}, \qquad (24) $$

where

$$ \psi_i(\nu) = \omega_i + \nu - 2 \arctan \left[ \tanh \left( \frac{F_i(\nu)}{2} \right) \right], \qquad (25) $$

$\omega_i$ is the argument of pericenter in the ith shell, and $F_i$ is the hyperbolic anomaly evaluated at $\nu$ in the ith shell, viz.

$$ \cosh F_i(\nu) = \frac{e_i + \cos \kappa_i \nu}{1 + e_i \cos \kappa_i \nu}, \qquad \sinh F_i(\nu) = \frac{\sqrt{e_i^2 - 1} \sin \kappa_i \nu}{1 + e_i \cos \kappa_i \nu}. \qquad (26) $$

Evaluating Equation 25 at $\nu = v_i^{(0)}$ and $\nu = v_i^{(1)}$, then substituting the results into Equation 24, we infer

$$ \varepsilon_i = \Delta v_i - 2 \left\{ \arctan \left[ \tanh \left( \frac{F_i^{(1)}}{2} \right) \right] - \arctan \left[ \tanh \left( \frac{F_i^{(0)}}{2} \right) \right] \right\}, \qquad (27) $$

where $\Delta v_i$ is given in Equation 19, and where we recall that $i = \{1, \dots, m-1\}$. $F_i^{(1)}$ (resp. $F_i^{(0)}$) is the value of hyperbolic anomaly computed at $\nu = v_i^{(1)}$ (resp. $\nu = v_i^{(0)}$). For the deepest shell, since the hyperbolic anomaly is zero at ray pericenter, the expression of the bending angle is given by

$$ \varepsilon_m = \Delta v_m + 2 \arctan \left[ \tanh \left( \frac{F_m^{(0)}}{2} \right) \right], \qquad (28) $$

where the expression of $\Delta v_m$ is given in Equation 20 and $F_m^{(0)}$ is the hyperbolic anomaly evaluated at $\nu_m(r_m)$. Therefore, $\varepsilon$, the total bending angle, is finally expressed such as
where $v_i^{(0)} = v_i(r_i)$ and $v_i^{(1)} = v_i(r_{i+1})$ are computed using Equation 17 with negative sign since we are only considering the first half of the ray where $r \cdot \hat{s} < 0$.

As for the mth shell (the deepest one), since the true anomaly is zero at the ray pericenter, the change in the true anomaly inside the first half of the deepest shell is simply

$$ \Delta v_m = \frac{1}{\kappa_m} \arccos\left[\left(\frac{p_m}{r_m} - 1\right) \frac{1}{e_m}\right]. \qquad (20) $$

Finally, the value of $\vartheta_E$ is

$$ \vartheta_E = \vartheta_1 + 2 \sum_{i=1}^{m} \Delta v_i, \qquad (21) $$

where $\vartheta_1$ is given by Equation 13. Therefore, the position vector $\mathbf{r}_E$ is given by

$$ \mathbf{r}_E = (r_1 \cos\vartheta_E, r_1 \sin\vartheta_E) \qquad (22) $$

Another important variable that will be useful in the following discussion is the ray pericenter radius $r_K$. It can be easily obtained observing that the closest approach point is located inside the deepest shell, at the point where the true anomaly is equal to zero. Therefore, from Equation 8 it follows that

$$ r_K = \frac{p_m}{1 + e_m}. \qquad (23) $$

### 2.2.2. Direction of the Ray at the Exit Point

In order to compute the total bending angle, which is needed for the determination of the Doppler effect, we need to infer an expression for the direction of the light ray at E, the point where the ray is exiting the atmosphere. As previously mentioned, within the spherical symmetry assumption, we can only focus on the first half of the light ray trajectory. Note that, the bending produced inside the ith shell can be obtained (recalling the definition of $\psi$) by

$$ \varepsilon_i = \psi_i(v_i^{(1)}) - \psi_i(v_i^{(0)}) \quad \text{with} \quad i = \{1, \dots, m-1\}, \qquad (24) $$

where

$$ \psi_i(\nu) = \omega_i + \nu - 2 \arctan \left[ \tanh \left( \frac{F_i(\nu)}{2} \right) \right], \qquad (25) $$

$\omega_i$ is the argument of pericenter in the ith shell, and $F_i$ is the hyperbolic anomaly evaluated at $\nu$ in the ith shell, viz.

$$ \cosh F_i(\nu) = \frac{e_i + \cos\kappa_i\nu}{1 + e_i \cos\kappa_i\nu}, \qquad \sinh F_i(\nu) = \frac{\sqrt{e_i^2 - 1} \sin\kappa_i\nu}{1 + e_i \cos\kappa_i\nu}. \qquad (26) $$

Evaluating Equation 25 at $\nu = v_i^{(0)}$ and $\nu = v_i^{(1)}$, then substituting the results into Equation 24, we infer

$$ \varepsilon_i = \Delta v_i - 2 \left\{ \arctan \left[ \tanh \left( \frac{F_i^{(1)}}{2} \right) \right] - \arctan \left[ \tanh \left( \frac{F_i^{(0)}}{2} \right) \right] \right\}, \qquad (27) $$

where $\Delta v_i$ is given in Equation 19, and where we recall that $i = \{1, \dots, m-1\}$. $F_i^{(1)}$ (resp. $F_i^{(0)}$) is the value of hyperbolic anomaly computed at $\nu = v_i^{(1)}$ (resp. $\nu = v_i^{(0)}$). For the deepest shell, since the hyperbolic anomaly is zero at ray pericenter, the expression of the bending angle is given by

$$ \varepsilon_m = \Delta v_m + 2 \arctan \left[ \tanh \left( \frac{F_m^{(0)}}{2} \right) \right], \qquad (28) $$

where the expression of $\Delta v_m$ is given in Equation 20 and $F_m^{(0)}$ is the hyperbolic anomaly evaluated at $\nu_m(r_m)$. Therefore, $\varepsilon$, the total bending angle, is finally expressed such as
$$ \epsilon = 2 \sum_{i=1}^{m} \epsilon_i, \qquad (29) $$

and hence, the direction of the light ray at E is given by

$$ \hat{s}_E = (\cos\beta_E, \sin\beta_E) \quad \text{with} \quad \beta_E = \beta_A + \epsilon. \qquad (30) $$

Note that, when $m = 1$, the proposed atmospheric modeling describes a single-layered atmosphere. In this case, the polar angle $\theta_E$ reduces to $\theta_E = \theta_1 + 2\Delta\nu_1$, where $\Delta\nu_1 = \Delta\nu_m$ is given by Equation 20. The total bending angle is $\epsilon = 2\epsilon_1$, with $\epsilon_1 = \epsilon_m$ given by Equation 28. Thus, $\mathbf{r}_E$ and $\hat{s}_E$ are still given by Equations 22 and 30.

At this point, it is worth noting that both the position and direction vectors $\mathbf{r}_E$ and $\hat{s}_E$ are only functions of $\beta_A$ and $\alpha_i$ for $i = \{1, \dots, m\}$. $\eta_i$ are not free variables since they may be computed using Equation 16 and enforcing the continuity conditions on the value of $n$ at the shells' boundaries (once the series of values $\alpha_i$ is given).

### 3. Radio Occultation Data Processing

In this section, we present a new method for processing radio occultation data based on the analytical ray-tracing method. Contrary to the method presented by Schinder et al. (2015), our technique is specialized to the case of spherically symmetric atmospheres and does not require any numerical integration, since, as explained in the previous section, we are able to compute analytically the ray path through the atmosphere. On the other hand, Schinder's method is more general, and, as such, can be used both for spherical and oblate atmospheres, but at the cost of a much longer computational time (as will be shown in the next section).

Radio occultation data processing starts with the analysis of the signal recorded at the ground station antenna to generate a time series of received frequencies. To this end, we may use, for example, a Fourier transform-based technique, as the one presented by Togni et al. (2021). In a ray-tracing method, each frequency value in the time series corresponds to a ray that passes through the atmosphere. In Section 2 we gave a description of the method used to propagate a radio signal path through the atmosphere, and this can be used both for one-way and two-way occultations. The only difference is that for one-way occultations we apply this methodology just once (from the spacecraft position to the Earth's antenna), whereas for two-way occultation this method is applied both to the uplink (from Earth to spacecraft) and downlink (from spacecraft to Earth) rays separately. The determination of the atmospheric refractivity profile is performed sequentially, namely, the atmosphere is built layer by layer, and the refractivity gradient in the deepest shell reached by a given light ray depends on the value of the refractivity gradient of previous shells. The next subsections explain how to process radio occultation data for one-way and two-way occultation experiments.

#### 3.1. Evaluation of Doppler Shift

In addition to the modeling of the direction and position of the ray at the exit from the atmosphere, another quantity that must be computed for processing radio occultation data is the Doppler shift. The Doppler frequency shift that a signal experiences when traveling entirely in a vacuum can be evaluated as (Schinder et al., 2015)

$$ \left( \frac{f_B}{f_A} \right)_{\text{vac}} = \frac{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_B/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_A/c} \sqrt{\frac{1 - 2U(\mathbf{r}_A)/c^2 - \nu_A^2/c^2}{1 - 2U(\mathbf{r}_B)/c^2 - \nu_B^2/c^2}} \qquad (31) $$

with $\nu_A$ and $\nu_B$ the norm of $\mathbf{v}_A$ and $\mathbf{v}_B$ being the velocities of the emitter and receiver, respectively, computed with respect to the occulting body. The Equation 31 is accurate to the order $c^{-2}$ with $c$ being the speed of light in a vacuum, and $U(\mathbf{r})$ is the total gravitational potential computed at the position $\mathbf{r}$.

According to Bourgoin et al. (2019), the deviation of the Doppler frequency shift $f_B/f_A$ (when atmospheric refraction occurs) from the vacuum contribution $(f_B/f_A)_{\text{vac}}$ can be expressed as

$$ \frac{f_B/f_A}{(f_B/f_A)_{\text{vac}}} = \frac{1 - n(\mathbf{r}_B) \frac{(\hat{s}_B - \hat{s}_{\text{vac}}) \cdot \mathbf{v}_B/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_B/c} - [n(\mathbf{r}_B) - 1] \frac{\hat{s}_{\text{vac}} \cdot \mathbf{v}_B/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_B/c}}{1 - n(\mathbf{r}_A) \frac{(\hat{s}_A - \hat{s}_{\text{vac}}) \cdot \mathbf{v}_A/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_A/c} - [n(\mathbf{r}_A) - 1] \frac{\hat{s}_{\text{vac}} \cdot \mathbf{v}_A/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_A/c}} \qquad (32) $$
$$ \epsilon = 2 \sum_{i=1}^{m} \epsilon_i, \qquad (29) $$

and hence, the direction of the light ray at E is given by

$$ \hat{s}_E = (\cos\beta_E, \sin\beta_E) \quad \text{with} \quad \beta_E = \beta_A + \epsilon. \qquad (30) $$

Note that, when $m = 1$, the proposed atmospheric modeling describes a single-layered atmosphere. In this case, the polar angle $\theta_E$ reduces to $\theta_E = \theta_1 + 2\Delta\nu_1$, where $\Delta\nu_1 = \Delta\nu_m$ is given by Equation 20. The total bending angle is $\epsilon = 2\epsilon_1$, with $\epsilon_1 = \epsilon_m$ given by Equation 28. Thus, $\mathbf{r}_E$ and $\hat{s}_E$ are still given by Equations 22 and 30.

At this point, it is worth noting that both the position and direction vectors $\mathbf{r}_E$ and $\hat{s}_E$ are only functions of $\beta_A$ and $\alpha_i$ for $i = \{1, \dots, m\}$. $\eta_i$ are not free variables since they may be computed using Equation 16 and enforcing the continuity conditions on the value of $n$ at the shells' boundaries (once the series of values $\alpha_i$ is given).

### 3. Radio Occultation Data Processing

In this section, we present a new method for processing radio occultation data based on the analytical ray-tracing method. Contrary to the method presented by Schinder et al. (2015), our technique is specialized to the case of spherically symmetric atmospheres and does not require any numerical integration, since, as explained in the previous section, we are able to compute analytically the ray path through the atmosphere. On the other hand, Schinder's method is more general, and, as such, can be used both for spherical and oblate atmospheres, but at the cost of a much longer computational time (as will be shown in the next section).

Radio occultation data processing starts with the analysis of the signal recorded at the ground station antenna to generate a time series of received frequencies. To this end, we may use, for example, a Fourier transform-based technique, as the one presented by Togni et al. (2021). In a ray-tracing method, each frequency value in the time series corresponds to a ray that passes through the atmosphere. In Section 2 we gave a description of the method used to propagate a radio signal path through the atmosphere, and this can be used both for one-way and two-way occultations. The only difference is that for one-way occultations we apply this methodology just once (from the spacecraft position to the Earth's antenna), whereas for two-way occultation this method is applied both to the uplink (from Earth to spacecraft) and downlink (from spacecraft to Earth) rays separately. The determination of the atmospheric refractivity profile is performed sequentially, namely, the atmosphere is built layer by layer, and the refractivity gradient in the deepest shell reached by a given light ray depends on the value of the refractivity gradient of previous shells. The next subsections explain how to process radio occultation data for one-way and two-way occultation experiments.

#### 3.1. Evaluation of Doppler Shift

In addition to the modeling of the direction and position of the ray at the exit from the atmosphere, another quantity that must be computed for processing radio occultation data is the Doppler shift. The Doppler frequency shift that a signal experiences when traveling entirely in a vacuum can be evaluated as (Schinder et al., 2015)

$$ \left( \frac{f_B}{f_A} \right)_{\text{vac}} = \frac{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_B/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_A/c} \sqrt{\frac{1 - 2U(\mathbf{r}_A)/c^2 - \nu_A^2/c^2}{1 - 2U(\mathbf{r}_B)/c^2 - \nu_B^2/c^2}} \qquad (31) $$

with $\nu_A$ and $\nu_B$ the norm of $\mathbf{v}_A$ and $\mathbf{v}_B$ being the velocities of the emitter and receiver, respectively, computed with respect to the occulting body. The Equation 31 is accurate to the order $c^{-2}$ with $c$ being the speed of light in a vacuum, and $U(\mathbf{r})$ is the total gravitational potential computed at the position $\mathbf{r}$.

According to Bourgoin et al. (2019), the deviation of the Doppler frequency shift $f_B/f_A$ (when atmospheric refraction occurs) from the vacuum contribution $(f_B/f_A)_{\text{vac}}$ can be expressed as

$$ \frac{f_B/f_A}{(f_B/f_A)_{\text{vac}}} = \frac{1 - n(\mathbf{r}_B) \frac{(\hat{s}_B - \hat{s}_{\text{vac}}) \cdot \mathbf{v}_B/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_B/c} - [n(\mathbf{r}_B) - 1] \frac{\hat{s}_{\text{vac}} \cdot \mathbf{v}_B/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_B/c}}{1 - n(\mathbf{r}_A) \frac{(\hat{s}_A - \hat{s}_{\text{vac}}) \cdot \mathbf{v}_A/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_A/c} - [n(\mathbf{r}_A) - 1] \frac{\hat{s}_{\text{vac}} \cdot \mathbf{v}_A/c}{1 - \hat{s}_{\text{vac}} \cdot \mathbf{v}_A/c}} \qquad (32) $$
where $(f_B/f_A)_{vac}$ can be obtained from Equation 31, $f_A$ is the frequency of the radio signal transmitted by the emitter A, and $f_B$ is the frequency of the signal when it is received at the antenna B. Note that, if both the transmitter and receiver are assumed to be in a vacuum, that is $n(r_B) = n(r_A) = 1$, then Equation 32 reduces to

$$ \frac{f_B/f_A}{(f_B/f_A)_{vac}} = \frac{1 - \frac{(\hat{s}_B - \hat{s}_{vac})v_B/c}{1 - \hat{s}_{vac}v_B/c}}{1 - \frac{(\hat{s}_A - \hat{s}_{vac})v_A/c}{1 - \hat{s}_{vac}v_A/c}} \equiv \chi(\beta_A, \alpha_1, \dots, \alpha_m), \quad (33) $$

where the right-hand side is a function $\chi$ of the variables $\beta_A$ and $\alpha_i$ with $i = \{1, \dots, m\}$. Indeed, Equation 33 gives a relation between the ray asymptotes directions $\hat{s}_A$ and $\hat{s}_B$ and the observables $f_B$. In turn, $\hat{s}_A$ and $\hat{s}_B$ can be computed following the methodology given in Section 2. In particular, $\hat{s}_A$ can be directly obtained from Equation 12 once $\beta_A$ is known, and $\hat{s}_B$ is the direction from $r_E$ to $r_B$ (see Equation 5), and $r_E$ is obtained from Equation 22. Therefore, Equation 33 lets us define a relation between $\beta_A$ and $\alpha_i$ with $i = \{1, \dots, m\}$ and the observables $f_B$, which will be useful in the following analyses of one-way and two-way occultations. Note that, we assume that the refractivity is equal to 1 at the ground station because the contribution of Earth's atmosphere remains constant during all the occultation, unlike the contribution coming from the central occulting body, so its effect can be absorbed when fitting the frequency residuals during the baseline as will be discussed hereafter.

## 3.2. One-Way Radio Occultation

In a one-way occultation experiment, the transmitting antenna is mounted on a spacecraft in orbit around a celestial body, whereas the receiving antenna is located on a ground station. In this case, an USO is usually employed on-board the spacecraft to provide a stable frequency reference. In the following analysis we make the assumption that the position of the spacecraft at time $t_A$ and ground antenna at $t_B$ with respect to the occulting body are known. In particular, we assume $t_B$ is an input, whereas $t_A$ is computed assuming a photon travels on a straight line between the transmitting antenna A and the receiving antenna B at the light speed in a vacuum (Withers et al., 2014). Also, the position of the occulting body is computed at the spacecraft time $t_A$.

The first step is to isolate the contribution to the frequency shift due to atmospheric refraction. In this context, let us define the residual frequency $f_{res}$ by

$$ f_{res} = f_B - f_{vac}, \quad (34) $$

where $f_{vac}$ represents the frequency that would have been recorded if the signal were to travel entirely in a vacuum,

$$ f_{vac} = \left( \frac{f_B}{f_A} \right)_{vac} f_A, \quad (35) $$

and $(f_B/f_A)_{vac}$ is given by Equation 31.

Just before an ingress occultation (or just after an egress occultation), there exists a time interval, called the baseline, during which the signal travels through a vacuum. The residual frequency should be zero during the baseline, but it usually departs from zero due to the effect of interplanetary plasma, Earth's troposphere and ionosphere, the long-term drift of the USO, and errors in the spacecraft trajectory. For this reason, a polynomial least square fit is performed on the residuals along the baseline, and the transmitted frequency $f_A$ is computed such that any trend or bias is removed, obtaining residuals that have zero mean value (Schinder et al., 2015).

The atmospheric refractivity $N = n - 1$ is given by the sum of two contributions: the part due to the planetary ionosphere $N_{iono} < 0$, and the one due to the neutral atmosphere $N_{neut} > 0$, that is, $N = N_{iono} + N_{neut}$ (see Section 6 of Withers (2010)). In general, in the upper part of the planetary atmosphere, the effect of the ionosphere is dominant, so $N_{neut}$ can be neglected with respect to $N_{iono}$. On the contrary, at lower altitudes, $N_{iono}$ is small when compared to $N_{neut}$, that is, we may consider the effect of neutral molecules only. There will be an intermediate region where both the ionosphere and neutral atmosphere refractivities offset each other, so $N \approx 0$. Usually, in that region, both $N_{iono}$ and $N_{neut}$ are assumed equal to zero (Withers, 2010). Note that, the above assumptions are useful when analyzing single-frequency experiments. Indeed, in the case of dual-frequency experiments (in which the transmitted signals are coherently related), data collected at two different frequency bands may be combined in
where $(f_B/f_A)_{vac}$ can be obtained from Equation 31, $f_A$ is the frequency of the radio signal transmitted by the emitter A, and $f_B$ is the frequency of the signal when it is received at the antenna B. Note that, if both the transmitter and receiver are assumed to be in a vacuum, that is $n(r_B) = n(r_A) = 1$, then Equation 32 reduces to

$$ \frac{f_B/f_A}{(f_B/f_A)_{vac}} = \frac{1 - \frac{(\hat{s}_B - \hat{s}_{vac})v_B/c}{1 - \hat{s}_{vac}v_B/c}}{1 - \frac{(\hat{s}_A - \hat{s}_{vac})v_A/c}{1 - \hat{s}_{vac}v_A/c}} \equiv \chi(\beta_A, \alpha_1, \dots, \alpha_m), \quad (33) $$

where the right-hand side is a function $\chi$ of the variables $\beta_A$ and $\alpha_i$ with $i = \{1, \dots, m\}$. Indeed, Equation 33 gives a relation between the ray asymptotes directions $\hat{s}_A$ and $\hat{s}_B$ and the observables $f_B$. In turn, $\hat{s}_A$ and $\hat{s}_B$ can be computed following the methodology given in Section 2. In particular, $\hat{s}_A$ can be directly obtained from Equation 12 once $\beta_A$ is known, and $\hat{s}_B$ is the direction from $r_E$ to $r_B$ (see Equation 5), and $r_E$ is obtained from Equation 22. Therefore, Equation 33 lets us define a relation between $\beta_A$ and $\alpha_i$ with $i = \{1, \dots, m\}$ and the observables $f_B$, which will be useful in the following analyses of one-way and two-way occultations. Note that, we assume that the refractivity is equal to 1 at the ground station because the contribution of Earth's atmosphere remains constant during all the occultation, unlike the contribution coming from the central occulting body, so its effect can be absorbed when fitting the frequency residuals during the baseline as will be discussed hereafter.

## 3.2. One-Way Radio Occultation

In a one-way occultation experiment, the transmitting antenna is mounted on a spacecraft in orbit around a celestial body, whereas the receiving antenna is located on a ground station. In this case, an USO is usually employed on-board the spacecraft to provide a stable frequency reference. In the following analysis we make the assumption that the position of the spacecraft at time $t_A$ and ground antenna at $t_B$ with respect to the occulting body are known. In particular, we assume $t_B$ is an input, whereas $t_A$ is computed assuming a photon travels on a straight line between the transmitting antenna A and the receiving antenna B at the light speed in a vacuum (Withers et al., 2014). Also, the position of the occulting body is computed at the spacecraft time $t_A$.

The first step is to isolate the contribution to the frequency shift due to atmospheric refraction. In this context, let us define the residual frequency $f_{res}$ by

$$ f_{res} = f_B - f_{vac}, \quad (34) $$

where $f_{vac}$ represents the frequency that would have been recorded if the signal were to travel entirely in a vacuum,

$$ f_{vac} = \left( \frac{f_B}{f_A} \right)_{vac} f_A, \quad (35) $$

and $(f_B/f_A)_{vac}$ is given by Equation 31.

Just before an ingress occultation (or just after an egress occultation), there exists a time interval, called the baseline, during which the signal travels through a vacuum. The residual frequency should be zero during the baseline, but it usually departs from zero due to the effect of interplanetary plasma, Earth's troposphere and ionosphere, the long-term drift of the USO, and errors in the spacecraft trajectory. For this reason, a polynomial least square fit is performed on the residuals along the baseline, and the transmitted frequency $f_A$ is computed such that any trend or bias is removed, obtaining residuals that have zero mean value (Schinder et al., 2015).

The atmospheric refractivity $N = n - 1$ is given by the sum of two contributions: the part due to the planetary ionosphere $N_{iono} < 0$, and the one due to the neutral atmosphere $N_{neut} > 0$, that is, $N = N_{iono} + N_{neut}$ (see Section 6 of Withers (2010)). In general, in the upper part of the planetary atmosphere, the effect of the ionosphere is dominant, so $N_{neut}$ can be neglected with respect to $N_{iono}$. On the contrary, at lower altitudes, $N_{iono}$ is small when compared to $N_{neut}$, that is, we may consider the effect of neutral molecules only. There will be an intermediate region where both the ionosphere and neutral atmosphere refractivities offset each other, so $N \approx 0$. Usually, in that region, both $N_{iono}$ and $N_{neut}$ are assumed equal to zero (Withers, 2010). Note that, the above assumptions are useful when analyzing single-frequency experiments. Indeed, in the case of dual-frequency experiments (in which the transmitted signals are coherently related), data collected at two different frequency bands may be combined in
such a way as to isolate the ionospheric effect ($N_{\text{iono}}$) from the neutral atmosphere's contribution ($N_{\text{neut}}$) (Dalba & Withers, 2019; Schinder et al., 2011a).

For the computation of the refractive index profile, we first need to define an arbitrary boundary of the atmosphere. This can be done by tracing a first ray as a straight line and computing its minimum distance from the occulting body center-of-mass O. This distance defines the radius $r_1$ of the first shell, whose outer boundary is assigned a value of refractive index equal to 1, that is, $n(r_1) = 1$. Then, the entire atmosphere is modeled as a single spherical shell with constant refractive index gradient $\alpha_1$. At this point, a second ray is traced using the procedure described in Section 2.2 with $m = 1$. In particular, the ray final state ($r_E$, $\hat{s}_E$) (or, equivalently, ($\theta_E$, $\beta_E$) at the point where the ray leaves the atmosphere) is given by Equations 22 and 30 as a function of $\beta_A$ and $\alpha_1$. Thus, the unknown values of $\beta_A$ and $\alpha_1$ can be obtained by enforcing the following conditions. The first one requires that the ray, after leaving the atmosphere, must point toward the position of the receiving antenna at time $t_B$. Moreover, the signal must arrive at the receiving antenna with the recorded frequency $f_B$. In order for the ray to point to the receiving antenna at time $t_B$, the following equation should be satisfied

$$ \hat{s}_E = \hat{s}_B, \qquad (36) $$

where it must be reminded that $\hat{s}_B$ is the unit vector from the point where the ray leaves the atmosphere toward the receiver; see Equation 5. Since $\hat{s}_E$ is a unit vector, it is sufficient to enforce the equality of Equation 36 on the y-component only, that is,

$$ \sin(\beta_E) = \frac{y_B - r_1 \sin\theta_E}{\sqrt{(x_B - r_1 \cos\theta_E)^2 + (y_B - r_1 \sin\theta_E)^2}}, \qquad (37) $$

with $r_B = x_B \hat{e}_x + y_B \hat{e}_y$. Note that, since the bending angle is a small quantity, the x-component of $\hat{s}_E$ is always negative.

An expression for the frequency received at a ground station can be obtained from Equation 33 as

$$ f_B = f_A \left( \frac{f_B}{f_A} \right)_{\text{vac}} \chi \qquad (38) $$

where $\chi$ is defined in Equation 33. The set of two Equations 37 and 38 can be solved for the two unknowns $\beta_A$ and $\alpha_1$ using a numerical method. At the same time, the ray pericenter radius is evaluated using Equation 23. This value is used as the outer radius $r_2$ of a second shell, whereas the first shell is assigned a fixed value of the refractive index gradient equal to $\alpha_1$.

Then, the unknown value of refractive index gradient $\alpha_2$ of the second shell is computed (along with $\beta_A$) by applying the procedure just described. In particular, the final state ($r_E$, $\hat{s}_E$) of a new traced ray is evaluated using the procedure of Section 2.2 with a number of shells $m = 2$. Note that, the state at point E for the second light ray depends on $\beta_A$, $\alpha_1$, and $\alpha_2$. However, the value of $\alpha_1$ is assumed fixed and equal to that obtained for the first ray. Therefore, the system of Equations 37 and 38 is solved for $\beta_A$ and $\alpha_2$ only.

The procedure is then repeated to find the initial direction $\beta_A$ and the value of refractive index gradient $\alpha_m$ of the deepest shell for all the subsequent rays, that is, for all frequency values of the time series, building the atmosphere layer by layer. Finally, once all values of refractive index gradient have been determined, the function $n = n(r)$ within each shell can be calculated using Equation 16 with the initial condition $n(r_1) = 1$.

### 3.3. Two-Way Radio Occultation

In case of two-way radio occultations, a signal is transmitted from the ground station, labeled A, to the spacecraft, labeled B (uplink signal). The frequency of the uplink signal is taken as a reference value. Indeed, the onboard electronics multiplies this frequency by a factor $\Lambda$ (the turn-around ratio) and transmits this back to the ground station, labeled C (downlink signal). For a two-way experiment, the uplink emitter A and downlink receiver C are the same observer (i.e., $A \equiv C$) at two different times. Both the uplink and downlink signals traverse the planetary atmosphere, whose physical properties are unknown. In the following discussion, we will use the superscript “up” when referring to the uplink ray, and “dw” for the downlink ray.
such a way as to isolate the ionospheric effect ($N_{\text{iono}}$) from the neutral atmosphere's contribution ($N_{\text{neut}}$) (Dalba & Withers, 2019; Schinder et al., 2011a).

For the computation of the refractive index profile, we first need to define an arbitrary boundary of the atmosphere. This can be done by tracing a first ray as a straight line and computing its minimum distance from the occulting body center-of-mass O. This distance defines the radius $r_1$ of the first shell, whose outer boundary is assigned a value of refractive index equal to 1, that is, $n(r_1) = 1$. Then, the entire atmosphere is modeled as a single spherical shell with constant refractive index gradient $\alpha_1$. At this point, a second ray is traced using the procedure described in Section 2.2 with $m = 1$. In particular, the ray final state ($r_E$, $\hat{s}_E$) (or, equivalently, ($\theta_E$, $\beta_E$) at the point where the ray leaves the atmosphere) is given by Equations 22 and 30 as a function of $\beta_A$ and $\alpha_1$. Thus, the unknown values of $\beta_A$ and $\alpha_1$ can be obtained by enforcing the following conditions. The first one requires that the ray, after leaving the atmosphere, must point toward the position of the receiving antenna at time $t_B$. Moreover, the signal must arrive at the receiving antenna with the recorded frequency $f_B$. In order for the ray to point to the receiving antenna at time $t_B$, the following equation should be satisfied

$$ \hat{s}_E = \hat{s}_B, \qquad (36) $$

where it must be reminded that $\hat{s}_B$ is the unit vector from the point where the ray leaves the atmosphere toward the receiver; see Equation 5. Since $\hat{s}_E$ is a unit vector, it is sufficient to enforce the equality of Equation 36 on the y-component only, that is,

$$ \sin(\beta_E) = \frac{y_B - r_1 \sin\theta_E}{\sqrt{(x_B - r_1 \cos\theta_E)^2 + (y_B - r_1 \sin\theta_E)^2}}, \qquad (37) $$

with $r_B = x_B \hat{e}_x + y_B \hat{e}_y$. Note that, since the bending angle is a small quantity, the x-component of $\hat{s}_E$ is always negative.

An expression for the frequency received at a ground station can be obtained from Equation 33 as

$$ f_B = f_A \left( \frac{f_B}{f_A} \right)_{\text{vac}} \chi \qquad (38) $$

where $\chi$ is defined in Equation 33. The set of two Equations 37 and 38 can be solved for the two unknowns $\beta_A$ and $\alpha_1$ using a numerical method. At the same time, the ray pericenter radius is evaluated using Equation 23. This value is used as the outer radius $r_2$ of a second shell, whereas the first shell is assigned a fixed value of the refractive index gradient equal to $\alpha_1$.

Then, the unknown value of refractive index gradient $\alpha_2$ of the second shell is computed (along with $\beta_A$) by applying the procedure just described. In particular, the final state ($r_E$, $\hat{s}_E$) of a new traced ray is evaluated using the procedure of Section 2.2 with a number of shells $m = 2$. Note that, the state at point E for the second light ray depends on $\beta_A$, $\alpha_1$, and $\alpha_2$. However, the value of $\alpha_1$ is assumed fixed and equal to that obtained for the first ray. Therefore, the system of Equations 37 and 38 is solved for $\beta_A$ and $\alpha_2$ only.

The procedure is then repeated to find the initial direction $\beta_A$ and the value of refractive index gradient $\alpha_m$ of the deepest shell for all the subsequent rays, that is, for all frequency values of the time series, building the atmosphere layer by layer. Finally, once all values of refractive index gradient have been determined, the function $n = n(r)$ within each shell can be calculated using Equation 16 with the initial condition $n(r_1) = 1$.

### 3.3. Two-Way Radio Occultation

In case of two-way radio occultations, a signal is transmitted from the ground station, labeled A, to the spacecraft, labeled B (uplink signal). The frequency of the uplink signal is taken as a reference value. Indeed, the onboard electronics multiplies this frequency by a factor $\Lambda$ (the turn-around ratio) and transmits this back to the ground station, labeled C (downlink signal). For a two-way experiment, the uplink emitter A and downlink receiver C are the same observer (i.e., $A \equiv C$) at two different times. Both the uplink and downlink signals traverse the planetary atmosphere, whose physical properties are unknown. In the following discussion, we will use the superscript “up” when referring to the uplink ray, and “dw” for the downlink ray.
For the uplink leg, the ground station antenna transmits a signal with frequency $f_A^{up}$ in its own rest frame at time $t_A^{up}$. The spacecraft antenna then receives a signal whose frequency is $f_B^{up}$ at time $t_B^{up}$. For the downlink leg, the spacecraft antenna transmits a signal at time $t_B^{dw}$ with frequency $f_B^{dw} = \Lambda f_B^{up}$, which is finally received at the ground station at time $t_C^{dw}$ and has a frequency equal to $f_C^{dw}$. In the following we assume $t_B^{dw} = t_B^{up}$, so the position $\mathbf{r}_B^{up}$ of the spacecraft when it receives the uplink signal coincides with the position $\mathbf{r}_B^{dw}$ when it transmits the downlink signal. Note that the position $\mathbf{r}_A^{up}$ of the ground station at time $t_A^{up}$ when it transmits the uplink signal is different from its position $\mathbf{r}_C^{dw}$ at reception time $t_C^{dw}$. Moreover, the uplink transmitted frequency $f_A^{up}$ is assumed known as a function of time, and the time series of downlink received frequencies $f_C^{dw}$ can be retrieved from the signal recorded at the ground station (Togni et al., 2021).

Similarly to the one-way case, the residual frequency can be defined as

$$f_{\text{res}} = f_C^{\text{dw}} - f_{\text{vac}}^{\text{dw}} \quad (39)$$

where $f_{\text{vac}}^{dw}$ would be the frequency of the downlink signal received by the ground station if both uplink and downlink rays traveled in a vacuum, and it can be computed as

$$f_{\text{vac}}^{\text{dw}} = \Lambda \left( \frac{f_C}{f_B} \right)_{\text{vac}}^{\text{dw}} \left( \frac{f_B}{f_A} \right)_{\text{vac}}^{\text{up}} f_A^{\text{up}} \quad (40)$$

where the terms $(f_B/f_A)_{\text{vac}}^{up}$ and $(f_C/f_B)_{\text{vac}}^{dw}$ represent the Doppler shift in a vacuum that can be computed using Equation 31 once the trajectories of all the objects involved in the occultation experiment are known. Attention must be paid when using Equation 31 since the emitter and receiver switch their roles between the uplink and downlink rays. Also in this case, a polynomial least square fit of $f_{\text{res}}$ is performed along the baseline, and the transmitted frequency $f_A^{up}$ is corrected in such a way as to obtain a residual frequency that has zero mean value during the baseline (Schinder et al., 2015).

The procedure shown in the following is based on the one described in Schinder et al. (2015, 2020), but with a slight difference. In particular, Schinder et al. (2020) are only interested in the analysis of Titan's neutral atmosphere, so they assume the refractive index is independent of the frequency. Therefore, they neglect the effect of Titan's ionosphere as it is weak and well above the portion of the atmosphere they analyze. A more general treatment should nevertheless consider the dependence of $n$ on the frequency when the signal traverses the ionosphere.

Let us define $n^{dw}$ as the refractive index at the frequency of the downlink signal. As explained in the previous subsection, the refractivity in the upper atmosphere is usually due to the ionosphere only, so $n^{dw} = n_{\text{iono}}^{dw}$ at high altitudes. Instead, in the lower atmosphere, the refraction is mainly produced by the neutral molecules, so $n^{dw} = n_{\text{neut}}^{dw}$ at lower altitudes. For two-way experiments, since the downlink and uplink signals have slightly different frequencies (their ratio being equal to $\Lambda$), the dependence of $n$ on the frequency must be taken into account when tracing the downlink and uplink rays through the ionosphere. In particular, the refractivity $N_{\text{iono}}$ due to the ionosphere is proportional to $f^{-2}$ (see Section 4 of Withers (2010)), so it follows that

$$\frac{n_{\text{iono}}^{\text{up}} - 1}{n_{\text{iono}}^{\text{dw}} - 1} = \Lambda^2 \quad (41)$$

Thus, the downlink ray is traced through the ionosphere assuming the refractive index is $n_{\text{iono}}^{dw}$. Instead, the uplink ray must be traced, enforcing that the refractive index is equal to $n_{\text{iono}}^{up}$, which can be obtained from Equation 41 once $n_{\text{iono}}^{dw}$ is known. Below a certain altitude $r = \tilde{r}$, that is, within the neutral atmosphere, the refractive index is independent of the frequency, so we may assume $n_{\text{neut}}^{dw} = n_{\text{neut}}^{up}$. The value of $\tilde{r}$ may be known from past literature results, or, otherwise, may be obtained with a trial-and-error procedure. Note that, for dual-frequency experiments, since the ionosphere's effect can be removed from the recorded data, the whole atmosphere can be thought of as consisting only of the neutral part as in Schinder et al. (2020).

At this point, the procedure to retrieve the refractivity profile $n^{dw} = n^{dw}(r)$ at the frequency of the downlink signal is similar to that described for the one-way case. In particular, each value of received frequency $f_C^{dw}$ in the time series is associated with a pair of uplink/downlink rays. First, the atmosphere boundary is to be defined. To this
For the uplink leg, the ground station antenna transmits a signal with frequency $f_A^{up}$ in its own rest frame at time $t_A^{up}$. The spacecraft antenna then receives a signal whose frequency is $f_B^{up}$ at time $t_B^{up}$. For the downlink leg, the spacecraft antenna transmits a signal at time $t_B^{dw}$ with frequency $f_B^{dw} = \Lambda f_B^{up}$, which is finally received at the ground station at time $t_C^{dw}$ and has a frequency equal to $f_C^{dw}$. In the following we assume $t_B^{dw} = t_B^{up}$, so the position $\mathbf{r}_B^{up}$ of the spacecraft when it receives the uplink signal coincides with the position $\mathbf{r}_B^{dw}$ when it transmits the downlink signal. Note that the position $\mathbf{r}_A^{up}$ of the ground station at time $t_A^{up}$ when it transmits the uplink signal is different from its position $\mathbf{r}_C^{dw}$ at reception time $t_C^{dw}$. Moreover, the uplink transmitted frequency $f_A^{up}$ is assumed known as a function of time, and the time series of downlink received frequencies $f_C^{dw}$ can be retrieved from the signal recorded at the ground station (Togni et al., 2021).

Similarly to the one-way case, the residual frequency can be defined as

$$f_{\text{res}} = f_C^{\text{dw}} - f_{\text{vac}}^{\text{dw}} \qquad (39)$$

where $f_{\text{vac}}^{dw}$ would be the frequency of the downlink signal received by the ground station if both uplink and downlink rays traveled in a vacuum, and it can be computed as

$$f_{\text{vac}}^{\text{dw}} = \Lambda \left( \frac{f_C}{f_B} \right)_{\text{vac}}^{\text{dw}} \left( \frac{f_B}{f_A} \right)_{\text{vac}}^{\text{up}} f_A^{\text{up}} \qquad (40)$$

where the terms $(f_B/f_A)_{\text{vac}}^{up}$ and $(f_C/f_B)_{\text{vac}}^{dw}$ represent the Doppler shift in a vacuum that can be computed using Equation 31 once the trajectories of all the objects involved in the occultation experiment are known. Attention must be paid when using Equation 31 since the emitter and receiver switch their roles between the uplink and downlink rays. Also in this case, a polynomial least square fit of $f_{\text{res}}$ is performed along the baseline, and the transmitted frequency $f_A^{up}$ is corrected in such a way as to obtain a residual frequency that has zero mean value during the baseline (Schinder et al., 2015).

The procedure shown in the following is based on the one described in Schinder et al. (2015, 2020), but with a slight difference. In particular, Schinder et al. (2020) are only interested in the analysis of Titan's neutral atmosphere, so they assume the refractive index is independent of the frequency. Therefore, they neglect the effect of Titan's ionosphere as it is weak and well above the portion of the atmosphere they analyze. A more general treatment should nevertheless consider the dependence of $n$ on the frequency when the signal traverses the ionosphere.

Let us define $n^{dw}$ as the refractive index at the frequency of the downlink signal. As explained in the previous subsection, the refractivity in the upper atmosphere is usually due to the ionosphere only, so $n^{dw} = n_{\text{iono}}^{dw}$ at high altitudes. Instead, in the lower atmosphere, the refraction is mainly produced by the neutral molecules, so $n^{dw} = n_{\text{neut}}^{dw}$ at lower altitudes. For two-way experiments, since the downlink and uplink signals have slightly different frequencies (their ratio being equal to $\Lambda$), the dependence of $n$ on the frequency must be taken into account when tracing the downlink and uplink rays through the ionosphere. In particular, the refractivity $N_{\text{iono}}$ due to the ionosphere is proportional to $f^{-2}$ (see Section 4 of Withers (2010)), so it follows that

$$\frac{n_{\text{iono}}^{\text{up}} - 1}{n_{\text{iono}}^{\text{dw}} - 1} = \Lambda^2 \qquad (41)$$

Thus, the downlink ray is traced through the ionosphere assuming the refractive index is $n_{\text{iono}}^{dw}$. Instead, the uplink ray must be traced, enforcing that the refractive index is equal to $n_{\text{iono}}^{up}$, which can be obtained from Equation 41 once $n_{\text{iono}}^{dw}$ is known. Below a certain altitude $r = \tilde{r}$, that is, within the neutral atmosphere, the refractive index is independent of the frequency, so we may assume $n_{\text{neut}}^{dw} = n_{\text{neut}}^{up}$. The value of $\tilde{r}$ may be known from past literature results, or, otherwise, may be obtained with a trial-and-error procedure. Note that, for dual-frequency experiments, since the ionosphere's effect can be removed from the recorded data, the whole atmosphere can be thought of as consisting only of the neutral part as in Schinder et al. (2020).

At this point, the procedure to retrieve the refractivity profile $n^{dw} = n^{dw}(r)$ at the frequency of the downlink signal is similar to that described for the one-way case. In particular, each value of received frequency $f_C^{dw}$ in the time series is associated with a pair of uplink/downlink rays. First, the atmosphere boundary is to be defined. To this
end, we trace the first pair of uplink and downlink rays as straight lines and compute their minimum distances $r_1^{up}$ and $r_1^{dw}$ from O. The atmosphere radius is chosen as $r_1 = \min\{r_1^{up}, r_1^{dw}\}$. The atmosphere boundary is assigned a value of refractive index $n^{dw}(r_1) = 1$, and the whole atmosphere is initially modeled as a single shell filled with a fluid of unknown constant refractive index gradient $\alpha_1$.

Thus, a second pair of uplink and downlink signals is considered, and they are separately traced through the atmosphere using the procedure outlined in Section 2.2 with $m = 1$. Due to the possible onset of numerical errors, we suggest tracing the uplink ray from the spacecraft antenna (the receiver) to the past position of the ground antenna (the transmitter), and, then, reversing the direction vector $\hat{s}$ (Schinder et al., 2015). The final state (at the exit point from the atmosphere) of both uplink and downlink rays is thus given as a function of the initial direction $\beta_B^{up}$ and $\beta_B^{dw}$ (both at the location of the spacecraft), and the refractive index gradient of the deepest shell $\alpha_1$. Their values may be obtained by enforcing that both uplink and downlink rays point to the position of the ground antenna once they leave the atmosphere, and the frequency of the (downlink) signal received at the ground station antenna must coincide with the recorded one. This latter condition can be obtained as in Withers and Moore (2020). In particular, by applying Equation 33 to the downlink signal, the received frequency can be written as

$$f_C^{dw} = \chi^{dw} \left( \frac{f_C}{f_B} \right)_{\text{vac}}^{dw} f_B^{dw} = \Lambda \chi^{dw} \left( \frac{f_C}{f_B} \right)_{\text{vac}}^{dw} f_B^{up} \quad (42)$$

where we used $f_B^{dw} = \Lambda f_B^{up}$ for the second equality. If Equation 33 is also applied to the uplink signal, we can find an expression for $f_B^{up}$, viz.

$$f_B^{up} = \chi^{up} \left( \frac{f_B}{f_A} \right)_{\text{vac}}^{up} f_A^{up}. \quad (43)$$

By substituting Equation 43 into 42, we finally get

$$f_C^{dw} = \Lambda \chi^{dw} \chi^{up} \left( \frac{f_C}{f_B} \right)_{\text{vac}}^{dw} \left( \frac{f_B}{f_A} \right)_{\text{vac}}^{up} f_A^{up}. \quad (44)$$

Therefore, the two equations for the pointing (which are analogous to Equation 37, but written both for the uplink and downlink rays) along with Equation 44 constitute a set of three equations in three unknowns. Once the values of $\beta_B^{up}$, $\beta_B^{dw}$ and $\alpha_1$ have been obtained, the pericenter radii $r_K^{up}$ and $r_K^{dw}$ of both uplink and downlink rays are computed using Equation 23, and the outer radius of a new (second) shell is defined as $r_2 = \min\{r_K^{up}, r_K^{dw}\}$. From now on, the first shell has a fixed value of refractive index gradient equal to $\alpha_1$, and the atmosphere is now constituted by two spherical shells, with the second shell assumed filled with a substance of constant refractive index gradient $\alpha_2$.

The procedure above is therefore repeated for a third pair of uplink/downlink rays, which are traced through the atmosphere as explained in Section 2.2 with $m = 2$. In particular, the state at the exit point of uplink and downlink rays is a function of $\beta_B^{up}$, $\beta_B^{dw}$, $\alpha_1$ and $\alpha_2$. However, as explained before, the value of $\alpha_1$ is assumed to be fixed and equal to that obtained when analyzing the second pair of rays. Therefore, a new system of three equations in three unknowns is numerically solved to get the initial (uplink and downlink) ray direction and the refractive index gradient $\alpha_2$ of the deepest shell. Then, the outer radius of a new shell is computed, and the whole procedure is repeated again for all the successive rays, that is, for all the frequencies $f_C^{dw}$ in the time series.

Finally, once all values $\alpha_i$ have been determined, the refractive index at the frequency of the downlink signal can be calculated using Equation 16 inside each shell. Thus, the uplink refractive index inside the ionosphere can be obtained from Equation 41, whereas it equals the one of the downlink signal inside the neutral atmosphere. Note that the approach described in this subsection can be simply extended to three-way occultations assuming the downlink and uplink ground-based antennas are not the same.

We want to point out another important issue that may arise in this analysis. Let us assume that at a certain iteration, the total number of shells is $m$. If, for example, the (downlink and/or uplink) ray's closest approach point is within the $m$th layer, then the procedure described in Section 2 can be used as is to compute the final state of the
end, we trace the first pair of uplink and downlink rays as straight lines and compute their minimum distances $r_1^{up}$ and $r_1^{dw}$ from O. The atmosphere radius is chosen as $r_1 = \min\{r_1^{up}, r_1^{dw}\}$. The atmosphere boundary is assigned a value of refractive index $n^{dw}(r_1) = 1$, and the whole atmosphere is initially modeled as a single shell filled with a fluid of unknown constant refractive index gradient $\alpha_1$.

Thus, a second pair of uplink and downlink signals is considered, and they are separately traced through the atmosphere using the procedure outlined in Section 2.2 with $m = 1$. Due to the possible onset of numerical errors, we suggest tracing the uplink ray from the spacecraft antenna (the receiver) to the past position of the ground antenna (the transmitter), and, then, reversing the direction vector $\hat{s}$ (Schinder et al., 2015). The final state (at the exit point from the atmosphere) of both uplink and downlink rays is thus given as a function of the initial direction $\beta_B^{up}$ and $\beta_B^{dw}$ (both at the location of the spacecraft), and the refractive index gradient of the deepest shell $\alpha_1$. Their values may be obtained by enforcing that both uplink and downlink rays point to the position of the ground antenna once they leave the atmosphere, and the frequency of the (downlink) signal received at the ground station antenna must coincide with the recorded one. This latter condition can be obtained as in Withers and Moore (2020). In particular, by applying Equation 33 to the downlink signal, the received frequency can be written as

$$f_C^{dw} = \chi^{dw} \left( \frac{f_C}{f_B} \right)_{\text{vac}}^{dw} f_B^{dw} = \Lambda \chi^{dw} \left( \frac{f_C}{f_B} \right)_{\text{vac}}^{dw} f_B^{up} \quad (42)$$

where we used $f_B^{dw} = \Lambda f_B^{up}$ for the second equality. If Equation 33 is also applied to the uplink signal, we can find an expression for $f_B^{up}$, viz.

$$f_B^{up} = \chi^{up} \left( \frac{f_B}{f_A} \right)_{\text{vac}}^{up} f_A^{up}. \quad (43)$$

By substituting Equation 43 into 42, we finally get

$$f_C^{dw} = \Lambda \chi^{dw} \chi^{up} \left( \frac{f_C}{f_B} \right)_{\text{vac}}^{dw} \left( \frac{f_B}{f_A} \right)_{\text{vac}}^{up} f_A^{up}. \quad (44)$$

Therefore, the two equations for the pointing (which are analogous to Equation 37, but written both for the uplink and downlink rays) along with Equation 44 constitute a set of three equations in three unknowns. Once the values of $\beta_B^{up}$, $\beta_B^{dw}$ and $\alpha_1$ have been obtained, the pericenter radii $r_K^{up}$ and $r_K^{dw}$ of both uplink and downlink rays are computed using Equation 23, and the outer radius of a new (second) shell is defined as $r_2 = \min\{r_K^{up}, r_K^{dw}\}$. From now on, the first shell has a fixed value of refractive index gradient equal to $\alpha_1$, and the atmosphere is now constituted by two spherical shells, with the second shell assumed filled with a substance of constant refractive index gradient $\alpha_2$.

The procedure above is therefore repeated for a third pair of uplink/downlink rays, which are traced through the atmosphere as explained in Section 2.2 with $m = 2$. In particular, the state at the exit point of uplink and downlink rays is a function of $\beta_B^{up}$, $\beta_B^{dw}$, $\alpha_1$ and $\alpha_2$. However, as explained before, the value of $\alpha_1$ is assumed to be fixed and equal to that obtained when analyzing the second pair of rays. Therefore, a new system of three equations in three unknowns is numerically solved to get the initial (uplink and downlink) ray direction and the refractive index gradient $\alpha_2$ of the deepest shell. Then, the outer radius of a new shell is computed, and the whole procedure is repeated again for all the successive rays, that is, for all the frequencies $f_C^{dw}$ in the time series.

Finally, once all values $\alpha_i$ have been determined, the refractive index at the frequency of the downlink signal can be calculated using Equation 16 inside each shell. Thus, the uplink refractive index inside the ionosphere can be obtained from Equation 41, whereas it equals the one of the downlink signal inside the neutral atmosphere. Note that the approach described in this subsection can be simply extended to three-way occultations assuming the downlink and uplink ground-based antennas are not the same.

We want to point out another important issue that may arise in this analysis. Let us assume that at a certain iteration, the total number of shells is $m$. If, for example, the (downlink and/or uplink) ray's closest approach point is within the $m$th layer, then the procedure described in Section 2 can be used as is to compute the final state of the
ray at the exit point. Instead, if the ray path does not reach the deepest shell, but, for example, the closest approach point is inside the (m - 1)-th shell, then the value of Δν$_{m-1}$ computed using Equation 19 is not a real number. This is because the ray does not intersect the shell with radius r$_{m}$. In this case, the problem can still be solved simply assuming that the atmosphere is composed of m - 1 shells and then using the procedure described in Section 2 to derive the final state of the ray.

3.4. From Refractivity to Physical Properties of Atmospheres

The refractivity profile obtained through the ray-tracing method introduced above can be used to infer the phys-
ical properties of atmospheres. First, in the case in which magnetic fields and electron-neutral collisions can
be neglected, the electron number density n_e in the ionosphere can be directly derived from the refractivity as
(Withers, 2010)

$$
n_e(f) = - \frac{8\pi^2 m_e \epsilon_0 f^2}{e^2} N_{\text{iono}} \qquad (45)
$$

where $m_e$ is the electron mass, $\epsilon_0$ is the vacuum permittivity, $e$ is the electron charge, and $f$ is the signal frequency.

Instead, for the pressure and temperature to be determined, it is first necessary to make some assumption on the neutral atmosphere composition, and, therefore, on the value of mean molecular mass $\bar{m}$ and mean refractive volume $N_v$ (Fjeldbo & Eshleman, 1968). Indeed, the total number density $n_n$ can be retrieved from the neutral refractivity as

$$
n_n = \frac{N_{\text{neut}}}{N_v} \tag{46}
$$

Also, for a fluid in hydrostatic equilibrium the following equality holds

$$
dP = -\bar{m} n_n \frac{\mu}{r^2} dr \qquad (47)
$$

where *P* is the pressure, which (within the assumption of ideal gas) can also be written as

$$
P = k_B n_n T \tag{48}
$$

where k_B is the Boltzmann constant and T is the temperature. An expression for the temperature can be obtained combining Equations 47 and 48 as follows

$$
T = T(r_0) \frac{n_n(r_0)}{n_n} - \frac{1}{k_B n_n} \int_{r_0}^{r} \bar{m} \frac{\mu}{r^2} dr \quad (49)
$$

where r₀ is a reference radius at which a suitable upper boundary condition *T(r₀)* is enforced. Therefore, once the temperature profile has been computed with Equation 49, the pressure can be simply obtained from Equation 48.

4. Test Cases

In this section, the novel method presented in this paper is tested and validated over a set of radio occultation experiments. The next subsections report some results of radio occultation data analyses in terms of temperature and pressure profiles for neutral atmospheres, and electron number density profiles of ionospheres. Both one-way and two-way experiments are examined, and the results are compared to those presented in the literature.

All the information concerning the position and velocity of the celestial bodies and objects involved in the radio
occultation experiments are provided by kernels used as input to the SPICE toolkit (Acton, 1996).

**4.1. One-Way Radio Occultation**

In order to validate our novel method, in this subsection we analyze a one-way occultation experiment. In particu-
lar, a comparison is performed between the results obtained using our ray-tracing approach and those obtained
ray at the exit point. Instead, if the ray path does not reach the deepest shell, but, for example, the closest approach point is inside the (m - 1)-th shell, then the value of Δν<sub>m-1</sub> computed using Equation 19 is not a real number. This is because the ray does not intersect the shell with radius r<sub>m</sub>. In this case, the problem can still be solved simply assuming that the atmosphere is composed of m - 1 shells and then using the procedure described in Section 2 to derive the final state of the ray.

3.4. From Refractivity to Physical Properties of Atmospheres

The refractivity profile obtained through the ray-tracing method introduced above can be used to infer the phys-
ical properties of atmospheres. First, in the case in which magnetic fields and electron-neutral collisions can
be neglected, the electron number density n_e in the ionosphere can be directly derived from the refractivity as
(Withers, 2010)

$$
n_e(f) = - \frac{8\pi^2 m_e \epsilon_0 f^2}{e^2} N_{\text{iono}} \qquad (45)
$$

where *m*<sub>e</sub> is the electron mass, *ε*<sub>0</sub> is the vacuum permittivity, *e* is the electron charge, and *f* is the signal frequency.

Instead, for the pressure and temperature to be determined, it is first necessary to make some assumption on the neutral atmosphere composition, and, therefore, on the value of mean molecular mass m and mean refractive volume N<sub>v</sub> (Fjeldbo & Eshleman, 1968). Indeed, the total number density n<sub>n</sub> can be retrieved from the neutral refractivity as

$$
n_n = \frac{N_{\text{neut}}}{N_v} \tag{46}
$$

Also, for a fluid in hydrostatic equilibrium the following equality holds

$$
dP = - \bar{m} n_n \frac{\mu}{r^2} dr \qquad (47)
$$

where *P* is the pressure, which (within the assumption of ideal gas) can also be written as

$$
P = k_B n_n T \tag{48}
$$

where k<sub>B</sub> is the Boltzmann constant and T is the temperature. An expression for the temperature can be obtained combining Equations 47 and 48 as follows

$$
T = T(r_0) \frac{n_n(r_0)}{n_n} - \frac{1}{k_B n_n} \int_{r_0}^{r} \bar{m} \frac{\mu}{r^2} dr \quad (49)
$$

where r₀ is a reference radius at which a suitable upper boundary condition *T(r₀*) is enforced. Therefore, once the temperature profile has been computed with Equation 49, the pressure can be simply obtained from Equation 48.

4. Test Cases

In this section, the novel method presented in this paper is tested and validated over a set of radio occultation experiments. The next subsections report some results of radio occultation data analyses in terms of temperature and pressure profiles for neutral atmospheres, and electron number density profiles of ionospheres. Both one-way and two-way experiments are examined, and the results are compared to those presented in the literature.

All the information concerning the position and velocity of the celestial bodies and objects involved in the radio
occultation experiments are provided by kernels used as input to the SPICE toolkit (Acton, 1996).

**4.1. One-Way Radio Occultation**

In order to validate our novel method, in this subsection we analyze a one-way occultation experiment. In particu-
lar, a comparison is performed between the results obtained using our ray-tracing approach and those obtained
with a classical Abel transform. We examined Mars Global Surveyor (MGS) ingress occultation by Mars that occurred on 27 December 1998. This one-way occultation experiment was already investigated by Withers et al. (2014), which used X-band data recorded at the DSN-25 antenna. Note that, to perform this analysis, we use the same time series of post-fit residuals as the one used by Withers et al. (2014), with a frequency estimation every 0.4 s (Withers, 2014).

Using the procedure outlined in Section 3.2 and 3.4, the pressure and temperature profiles can be computed. Note that, for each value of the received frequency in the time series, the set of non-linear Equations 37 and 38 has been solved with the built-in MATLAB function fsolve (Powell, 1970) and using the same kernels as those listed in Withers et al. (2014). In particular, we assume Mars' atmosphere is characterized by a constant value of mean molecular mass $\bar{m} = 7.221 \times 10^{-26}$ kg and a mean refractive volume $N_v = 1.804 \times 10^{-29}$ m³ (Hinson et al., 1999). The resulting pressure and temperature profiles are reported in the top and middle panels of Figure 3, and they are also compared to the results reported in Withers et al. (2014). Also, the electron density profile for this occultation can be calculated using Equation 45 and is shown in the bottom panel of Figure 3. As shown in the right panels of Figure 3, the atmospheric properties obtained with our approach are fully consistent with those reported in the literature. Indeed, the difference between our results and those reported by Withers et al. (2014) is always smaller than the uncertainty values reported in the Planetary Data System (PDS) for this occultation (Hinson, 1999, 2000).

## 4.2. Two-Way Radio Occultation

In this subsection, two test cases are shown in which we use the procedure outlined in Section 3.3 to process time series of radio frequencies collected during two-way occultation experiments. In particular, we examine an occultation of the MAVEN spacecraft by Mars and another of the Cassini spacecraft by Titan.

### 4.2.1. MAVEN-Mars Occultation

The MAVEN Radio Occultation Science Experiment is aimed at measuring the electron density profiles of Mars' ionosphere (Withers et al., 2018, 2020). It consists of two-way, single-frequency (X-band) radio occultation experiments with a turn-around ratio $\Lambda = 880/749$. Here we examine the ingress occultation that occurred on 14 February 2017. All the input data required by the ray-tracing method (time series of received frequency and uplink transmitted frequency) have been downloaded from the Planetary Plasma Interactions Node of PDS (Withers, 2023). In particular, our analysis uses data collected at the DSS-55 antenna and the same residual frequencies reported in the PDS, with a frequency estimation given every 0.5 s.

Figures 4 and 5 depict the atmospheric properties obtained for this experiment. In particular, the upper panels of Figure 4 depict a comparison between the refractivity (at the frequency of the downlink signal) computed with the procedure introduced in this paper and that reported in the PDS, which has been obtained using the procedure described by Withers and Moore (2020). Note that, in this example, we have assumed the neutral atmosphere extends from the surface up to $\tilde{r} = 3460$ km, where we have observed that $N \approx 0$.

The refractive index profile can be used to infer the electron number density inside Mars' ionosphere using Equation 45, and this is shown in the bottom panels of Figure 4. Note that, the root mean square of the difference between our result and the one archived in the PDS is about $260 \text{ cm}^{-3}$ and the maximum difference is about $870 \text{ cm}^{-3}$, a value smaller than the $1\sigma$ uncertainty of $1,044 \text{ cm}^{-3}$ (as shown in the right bottom panel of Figure 4).

The radio data can be also exploited to derive the properties of the neutral atmosphere at lower altitudes (where we have observed that $N > 0$). Note that Withers' approach is only able to derive the refractivity inside the ionosphere (Withers & Moore, 2020), so it can not be used to infer the physical properties of the neutral atmosphere. Indeed, the right top panel of Figure 4 depicts a significant difference (greater than the $3\sigma$ uncertainty) between the refractivity values computed at low altitudes with our approach and those reported in the PDS. Pressure and temperature profiles can be retrieved using Equations 48 and 49. For this computation we used a mean molecular mass of $7.221 \times 10^{-26}$ kg and a mean refractive volume $N_v = 1.804 \times 10^{-29}$ m³ for Mars' atmosphere (Hinson et al., 1999). A boundary condition for the computation of pressure and temperature may be generated using the approximate procedure described by Withers et al. (2003, 2014). In particular, at a given radial distance $r_0 = 3,440$ km we may obtain an initial value for the pressure as
with a classical Abel transform. We examined Mars Global Surveyor (MGS) ingress occultation by Mars that occurred on 27 December 1998. This one-way occultation experiment was already investigated by Withers et al. (2014), which used X-band data recorded at the DSN-25 antenna. Note that, to perform this analysis, we use the same time series of post-fit residuals as the one used by Withers et al. (2014), with a frequency estimation every 0.4 s (Withers, 2014).

Using the procedure outlined in Section 3.2 and 3.4, the pressure and temperature profiles can be computed. Note that, for each value of the received frequency in the time series, the set of non-linear Equations 37 and 38 has been solved with the built-in MATLAB function fsolve (Powell, 1970) and using the same kernels as those listed in Withers et al. (2014). In particular, we assume Mars' atmosphere is characterized by a constant value of mean molecular mass $\bar{m} = 7.221 \times 10^{-26}$ kg and a mean refractive volume $N_v = 1.804 \times 10^{-29}$ m³ (Hinson et al., 1999). The resulting pressure and temperature profiles are reported in the top and middle panels of Figure 3, and they are also compared to the results reported in Withers et al. (2014). Also, the electron density profile for this occultation can be calculated using Equation 45 and is shown in the bottom panel of Figure 3. As shown in the right panels of Figure 3, the atmospheric properties obtained with our approach are fully consistent with those reported in the literature. Indeed, the difference between our results and those reported by Withers et al. (2014) is always smaller than the uncertainty values reported in the Planetary Data System (PDS) for this occultation (Hinson, 1999, 2000).

## 4.2. Two-Way Radio Occultation

In this subsection, two test cases are shown in which we use the procedure outlined in Section 3.3 to process time series of radio frequencies collected during two-way occultation experiments. In particular, we examine an occultation of the MAVEN spacecraft by Mars and another of the Cassini spacecraft by Titan.

### 4.2.1. MAVEN-Mars Occultation

The MAVEN Radio Occultation Science Experiment is aimed at measuring the electron density profiles of Mars' ionosphere (Withers et al., 2018, 2020). It consists of two-way, single-frequency (X-band) radio occultation experiments with a turn-around ratio $\Lambda = 880/749$. Here we examine the ingress occultation that occurred on 14 February 2017. All the input data required by the ray-tracing method (time series of received frequency and uplink transmitted frequency) have been downloaded from the Planetary Plasma Interactions Node of PDS (Withers, 2023). In particular, our analysis uses data collected at the DSS-55 antenna and the same residual frequencies reported in the PDS, with a frequency estimation given every 0.5 s.

Figures 4 and 5 depict the atmospheric properties obtained for this experiment. In particular, the upper panels of Figure 4 depict a comparison between the refractivity (at the frequency of the downlink signal) computed with the procedure introduced in this paper and that reported in the PDS, which has been obtained using the procedure described by Withers and Moore (2020). Note that, in this example, we have assumed the neutral atmosphere extends from the surface up to $\tilde{r} = 3460$ km, where we have observed that $N \approx 0$.

The refractive index profile can be used to infer the electron number density inside Mars' ionosphere using Equation 45, and this is shown in the bottom panels of Figure 4. Note that, the root mean square of the difference between our result and the one archived in the PDS is about $260 \text{ cm}^{-3}$ and the maximum difference is about $870 \text{ cm}^{-3}$, a value smaller than the $1\sigma$ uncertainty of $1,044 \text{ cm}^{-3}$ (as shown in the right bottom panel of Figure 4).

The radio data can be also exploited to derive the properties of the neutral atmosphere at lower altitudes (where we have observed that $N > 0$). Note that Withers' approach is only able to derive the refractivity inside the ionosphere (Withers & Moore, 2020), so it can not be used to infer the physical properties of the neutral atmosphere. Indeed, the right top panel of Figure 4 depicts a significant difference (greater than the $3\sigma$ uncertainty) between the refractivity values computed at low altitudes with our approach and those reported in the PDS. Pressure and temperature profiles can be retrieved using Equations 48 and 49. For this computation we used a mean molecular mass of $7.221 \times 10^{-26}$ kg and a mean refractive volume $N_v = 1.804 \times 10^{-29}$ m³ for Mars' atmosphere (Hinson et al., 1999). A boundary condition for the computation of pressure and temperature may be generated using the approximate procedure described by Withers et al. (2003, 2014). In particular, at a given radial distance $r_0 = 3,440$ km we may obtain an initial value for the pressure as
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Profiles for Mars' atmospheric properties retrieved during Mars Global Surveyor one-way occultation: the left panels depict the profiles retrieved with the analytical ray-tracing method (black dotted lines) and with the procedure outlined in Withers et al. (2014) (red solid lines). The right panels depict the difference between the two profiles and the $1\sigma$ uncertainty reported in the Planetary Data System.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Profiles for Mars' atmospheric properties retrieved during Mars Global Surveyor one-way occultation: the left panels depict the profiles retrieved with the analytical ray-tracing method (black dotted lines) and with the procedure outlined in Withers et al. (2014) (red solid lines). The right panels depict the difference between the two profiles and the $1\sigma$ uncertainty reported in the Planetary Data System.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4. Profiles for Mars' atmospheric properties retrieved during MAVEN two-way occultation: the upper panels depict a comparison between the refractivity profiles retrieved with the analytical ray-tracing method and that stored in the Planetary Data System. The bottom panels depict a comparison of the electron density profiles.</figcaption></figure>

$$P(r_0) = \bar{m} n_n(r_0) \frac{\mu}{r_0^2} H \qquad (50)$$

where *H* is the scale height of neutral refractivity profile evaluated between 3,430 and 3,440 km. Thus, the initial value of temperature can be stated from the ideal gas Equation 48, viz.

$$T(r_0) = \frac{\mathcal{N}_v P(r_0)}{k_B N_{\text{neut}}(r_0)} \qquad (51)$$

With this procedure we obtain $P(r_0) \approx 3$ Pa and $T(r_0) \approx 165$ K. The whole pressure and temperature profiles are shown in Figure 5 along with error bars that represent an estimate of the uncertainties in the pressure and
<figure><img src="image_6.png" /><figcaption>Figure 4. Profiles for Mars' atmospheric properties retrieved during MAVEN two-way occultation: the upper panels depict a comparison between the refractivity profiles retrieved with the analytical ray-tracing method and that stored in the Planetary Data System. The bottom panels depict a comparison of the electron density profiles.</figcaption></figure>

$$P(r_0) = \bar{m} n_n(r_0) \frac{\mu}{r_0^2} H \qquad (50)$$

where *H* is the scale height of neutral refractivity profile evaluated between 3,430 and 3,440 km. Thus, the initial value of temperature can be stated from the ideal gas Equation 48, viz.

$$T(r_0) = \frac{\mathcal{N}_v P(r_0)}{k_B N_{\text{neut}}(r_0)} \qquad (51)$$

With this procedure we obtain $P(r_0) \approx 3$ Pa and $T(r_0) \approx 165$ K. The whole pressure and temperature profiles are shown in Figure 5 along with error bars that represent an estimate of the uncertainties in the pressure and
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Profiles for Mars' atmospheric properties retrieved during MAVEN two-way occultation: comparison between the temperature and pressure profiles obtained using the ray-tracing method and those retrieved from the Mars Climate Database. The error bars show an estimate of the uncertainty due to thermal noise; see Bourgoin et al. (2022).</figcaption></figure>

temperature values due to thermal noise, obtained with the approximate expressions introduced by Bourgoin et al. (2022). Figure 5 depicts that the temperature and pressure values are consistent with those expected by the numerical simulation of Mars' atmosphere presented in the Mars Climate Database; see Forget et al. (1999) and Millour et al. (2021).

### 4.2.2. Cassini-Titan Occultation

This subsection reports the results of our analysis of Cassini-Titan two-way ingress occultation T102 occurred on 18 June 2014. For this test, we used X-band data collected at the DSS-43 antenna. The same antenna was also used to transmit the uplink signal. Note that, for this analysis a frequency estimation is given every 0.128 s and we perform the baseline correction using a linear least square fit. Also, we consider a value of $\bar{m} \approx 4.624 \times 10^{-26}$ kg and $N_v \approx 1.107 \times 10^{-29} \text{ m}^3$ for Titan's atmosphere. For the retrieval of pressure and temperature profiles we use the hydrostatic equilibrium equation and the Bender equation of state (Bender, 1973; Buhner et al., 1981). The pressure and temperature profiles are shown in Figure 6.

Note that the results obtained with the method proposed in this paper are consistent with those archived in the PDS (Schinder, 2020) and obtained with the method outlined by Schinder et al. (2015), except in the lower part of the atmosphere where the difference between the two pressure profiles is slightly larger than the uncertainty due to thermal noise reported in the PDS, as shown in the right top panel of Figure 6. Anyway, this difference is still very small (~2 × 10⁻³ bar) if compared to the value of pressure at the surface (~1.4 bar) and does not significantly affect the scientific interpretation of the results for Titan's atmosphere.

We have also retrieved the electron density profile of Titan's ingress occultation T102 and compared the result with that published in the manuscript by Dalba and Withers (2019) and archived in the PDS (Withers & Huber, 2020). In particular, for this analysis, we have used the S and X band data collected at the DSS-43 antenna. Figure 7 depicts that the difference between the PDS profile and the result of the analytical ray-tracing algorithm is always below the uncertainty value except for a few points where the difference is just slightly higher than the 1σ uncertainty reported in the PDS (Withers & Huber, 2020). Therefore, we conclude that the results obtained with the two methodologies are fully consistent.

Moreover, a method based on analytical ray-tracing gives better performance in terms of computational time if compared to the fully numerical method. Indeed, the time needed by the latter is about 220 times longer than that required by our new approach when using a count time of 2 s between frequency samples, 300 times longer when using a count time of 1 s, and 400 times longer if the count time is 0.128 s. The lower computational effort required by the analytical approach may enable an increase in the resolution of the atmospheric properties' profiles (still taking into account limitations imposed by diffraction phenomena and given by the value of the first Fresnel zone (Bocanegra-Bahamón et al., 2023)) while still requiring a reasonable computational burden.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Profiles for Mars' atmospheric properties retrieved during MAVEN two-way occultation: comparison between the temperature and pressure profiles obtained using the ray-tracing method and those retrieved from the Mars Climate Database. The error bars show an estimate of the uncertainty due to thermal noise; see Bourgoin et al. (2022).</figcaption></figure>

temperature values due to thermal noise, obtained with the approximate expressions introduced by Bourgoin et al. (2022). Figure 5 depicts that the temperature and pressure values are consistent with those expected by the numerical simulation of Mars' atmosphere presented in the Mars Climate Database; see Forget et al. (1999) and Millour et al. (2021).

### 4.2.2. Cassini-Titan Occultation

This subsection reports the results of our analysis of Cassini-Titan two-way ingress occultation T102 occurred on 18 June 2014. For this test, we used X-band data collected at the DSS-43 antenna. The same antenna was also used to transmit the uplink signal. Note that, for this analysis a frequency estimation is given every 0.128 s and we perform the baseline correction using a linear least square fit. Also, we consider a value of $\bar{m} \approx 4.624 \times 10^{-26}$ kg and $N_v \approx 1.107 \times 10^{-29} \text{ m}^3$ for Titan's atmosphere. For the retrieval of pressure and temperature profiles we use the hydrostatic equilibrium equation and the Bender equation of state (Bender, 1973; Buhner et al., 1981). The pressure and temperature profiles are shown in Figure 6.

Note that the results obtained with the method proposed in this paper are consistent with those archived in the PDS (Schinder, 2020) and obtained with the method outlined by Schinder et al. (2015), except in the lower part of the atmosphere where the difference between the two pressure profiles is slightly larger than the uncertainty due to thermal noise reported in the PDS, as shown in the right top panel of Figure 6. Anyway, this difference is still very small (~2 × 10⁻³ bar) if compared to the value of pressure at the surface (~1.4 bar) and does not significantly affect the scientific interpretation of the results for Titan's atmosphere.

We have also retrieved the electron density profile of Titan's ingress occultation T102 and compared the result with that published in the manuscript by Dalba and Withers (2019) and archived in the PDS (Withers & Huber, 2020). In particular, for this analysis, we have used the S and X band data collected at the DSS-43 antenna. Figure 7 depicts that the difference between the PDS profile and the result of the analytical ray-tracing algorithm is always below the uncertainty value except for a few points where the difference is just slightly higher than the 1σ uncertainty reported in the PDS (Withers & Huber, 2020). Therefore, we conclude that the results obtained with the two methodologies are fully consistent.

Moreover, a method based on analytical ray-tracing gives better performance in terms of computational time if compared to the fully numerical method. Indeed, the time needed by the latter is about 220 times longer than that required by our new approach when using a count time of 2 s between frequency samples, 300 times longer when using a count time of 1 s, and 400 times longer if the count time is 0.128 s. The lower computational effort required by the analytical approach may enable an increase in the resolution of the atmospheric properties' profiles (still taking into account limitations imposed by diffraction phenomena and given by the value of the first Fresnel zone (Bocanegra-Bahamón et al., 2023)) while still requiring a reasonable computational burden.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Pressure and temperature profiles for T102 two-way occultation: the left panels depict the profiles retrieved with the analytical ray-tracing method (black dotted line) and that reported in the Planetary Data System (PDS) (red solid line). The right panels depict the difference between the two profiles and the uncertainty values reported in the PDS.</figcaption></figure>

**5. Conclusion**

This paper presented a novel approach for processing data collected during radio occultation experiments to observe planetary atmospheres. Based on recent literature results, a mathematical model has been developed to analytically trace the signal path through a multilayered and spherically symmetric atmosphere. This model lets an analytical expression be obtained for the ray bending angle, so no approximate additional assumptions are required for its computation as in Abel transform-based approaches for two-way occultations.

We have shown that the proposed method can provide solutions that are fully consistent with past literature results or with existing models of atmospheric properties. Also, a procedure has been proposed to deal with single-frequency experiments, such that the properties of a celestial body's ionosphere and neutral atmosphere can be retrieved simultaneously and independently. Finally, the analytical ray-tracing-based method has been
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Pressure and temperature profiles for T102 two-way occultation: the left panels depict the profiles retrieved with the analytical ray-tracing method (black dotted line) and that reported in the Planetary Data System (PDS) (red solid line). The right panels depict the difference between the two profiles and the uncertainty values reported in the PDS.</figcaption></figure>

**5. Conclusion**

This paper presented a novel approach for processing data collected during radio occultation experiments to observe planetary atmospheres. Based on recent literature results, a mathematical model has been developed to analytically trace the signal path through a multilayered and spherically symmetric atmosphere. This model lets an analytical expression be obtained for the ray bending angle, so no approximate additional assumptions are required for its computation as in Abel transform-based approaches for two-way occultations.

We have shown that the proposed method can provide solutions that are fully consistent with past literature results or with existing models of atmospheric properties. Also, a procedure has been proposed to deal with single-frequency experiments, such that the properties of a celestial body's ionosphere and neutral atmosphere can be retrieved simultaneously and independently. Finally, the analytical ray-tracing-based method has been
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Electron density profile for T102 two-way occultation: the left panel depicts the profiles retrieved with the analytical ray-tracing method (black dotted line) and that reported in the Planetary Data System (PDS) (red solid line). The right panel depicts the difference between the two profiles and the uncertainty values reported in the PDS.</figcaption></figure>

shown to be hundreds of times faster than the existing numerical approach, so that the resolution of the atmospheric properties' profiles may be improved (within limits imposed by diffraction effects) without requiring an excessive increase of the computational time.

## Data Availability Statement

Mars' atmospheric properties profiles and time series of residual frequency for MGS occultation can be retrieved in Withers (2014). Uncertainties for this occultation can be found in Hinson (1999, 2000). Data for the MAVEN ROSE experiment can be downloaded from Withers (2023). Finally, atmospheric profiles derived from radio science data acquired during Cassini occultations by Titan are available in Schinder (2020), whereas the ionospheric profiles can be found in Withers and Huber (2020).

## References

* Acton, C. H. (1996). Ancillary data services of NASA's navigation and ancillary information facility. *Planetary and Space Science*, 44(1), 65-70. https://doi.org/10.1016/0032-0633(95)00107-7

* Ando, H., Imamura, T., Tellmann, S., Pätzold, M., Häusler, B., Sugimoto, N., et al. (2020). Thermal structure of the Venusian atmosphere from the sub-cloud region to the mesosphere as observed by radio occultation. *Scientific Reports*, 10(1), 3448. https://doi.org/10.1038/s41598-020-59278-8

* Bender, E. (1973). An equation of state for predicting vapour-liquid equilibria of the system N2-Ar-O2. *Cryogenics*, 13(1), 11-18. https://doi.org/10.1016/0011-2275(73)90258-0

* Bocanegra-Bahamón, T. M., Ao, C. O., Wang, K. N., & Vergados, P. (2023). Phase matching method for inversion of Venus radio occultation signals. *Radio Science*, 58(3), e2022RS007579. https://doi.org/10.1029/2022RS007579

* Bocanegra-Bahamón, T. M., Molera Calvés, G., Gurvits, L. I., Cimò, G., Dirkx, D., Duev, D. A., et al. (2019). Venus Express radio occultation observed by PRIDE. *Astronomy and Astrophysics*, 624, A59. https://doi.org/10.1051/0004-6361/201833160

* Born, M., & Wolf, E. (1980). In *Principles of optics* (6th ed.). Pergamon Press Ltd. Oxford, United Kingdom. (chap. 3).

* Bourgoin, A. (2020). General expansion of time transfer functions in optical spacetime. *Physical Review D*, 101(6), 064035. https://doi.org/10.1103/PhysRevD.101.064035

* Bourgoin, A., Gramigna, E., Zannoni, M., Gomez Casajus, L., & Tortora, P. (2022). Determination of uncertainty profiles in neutral atmospheric properties measured by radio occultation experiments. *Advances in Space Research*, 70(8), 2555-2570. https://doi.org/10.1016/j.asr.2022.07.015

* Bourgoin, A., Zannoni, M., Gomez Casajus, L., Tortora, P., & Teyssandier, P. (2021). Relativistic modeling of atmospheric occultations with time transfer functions. *Astronomy and Astrophysics*, 648, A46. https://doi.org/10.1051/0004-6361/202040269

* Bourgoin, A., Zannoni, M., & Tortora, P. (2019). Analytical ray-tracing in planetary atmospheres. *Astronomy and Astrophysics*, 624, A41. https://doi.org/10.1051/0004-6361/201834962

* Buccino, D., Parisi, M., Gramigna, E., Gomez-Casajus, L., Tortora, P., Zannoni, M., et al. (2022). Ganymede's ionosphere observed by a dual-frequency radio occultation with Juno. *Geophysical Research Letters*, 49(23), e2022GL098420. https://doi.org/10.1029/2022GL098420

* Buhner, K., Maurer, G., & Bender, E. (1981). Pressure-enthalpy diagrams for methane, ethane, propane, ethylene and propylene. *Cryogenics*, 21(3), 157-164. https://doi.org/10.1016/0011-2275(81)90267-8

* Dalba, P. A., & Withers, P. (2019). Cassini radio occultation observations of Titan's ionosphere: The complete set of electron density profiles. *Journal of Geophysical Research*, 124(1), 643-660. https://doi.org/10.1029/2018JA025693

## Acknowledgments

The authors are grateful to the Italian Space Agency (ASI) for financial support through ASI-INAF Agreement No. 2018-25-HH.0 and 2023-6-HH.0 in the context of ESA's JUICE mission, and ASI-CRAS Agreement No. 2022.16.HH.0 in the context of NASA's Juno mission.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Electron density profile for T102 two-way occultation: the left panel depicts the profiles retrieved with the analytical ray-tracing method (black dotted line) and that reported in the Planetary Data System (PDS) (red solid line). The right panel depicts the difference between the two profiles and the uncertainty values reported in the PDS.</figcaption></figure>

shown to be hundreds of times faster than the existing numerical approach, so that the resolution of the atmospheric properties' profiles may be improved (within limits imposed by diffraction effects) without requiring an excessive increase of the computational time.

## Data Availability Statement

Mars' atmospheric properties profiles and time series of residual frequency for MGS occultation can be retrieved in Withers (2014). Uncertainties for this occultation can be found in Hinson (1999, 2000). Data for the MAVEN ROSE experiment can be downloaded from Withers (2023). Finally, atmospheric profiles derived from radio science data acquired during Cassini occultations by Titan are available in Schinder (2020), whereas the ionospheric profiles can be found in Withers and Huber (2020).

## References

* Acton, C. H. (1996). Ancillary data services of NASA's navigation and ancillary information facility. *Planetary and Space Science*, 44(1), 65-70. https://doi.org/10.1016/0032-0633(95)00107-7

* Ando, H., Imamura, T., Tellmann, S., Pätzold, M., Häusler, B., Sugimoto, N., et al. (2020). Thermal structure of the Venusian atmosphere from the sub-cloud region to the mesosphere as observed by radio occultation. *Scientific Reports*, 10(1), 3448. https://doi.org/10.1038/s41598-020-59278-8

* Bender, E. (1973). An equation of state for predicting vapour-liquid equilibria of the system N2-Ar-O2. *Cryogenics*, 13(1), 11-18. https://doi.org/10.1016/0011-2275(73)90258-0

* Bocanegra-Bahamón, T. M., Ao, C. O., Wang, K. N., & Vergados, P. (2023). Phase matching method for inversion of Venus radio occultation signals. *Radio Science*, 58(3), e2022RS007579. https://doi.org/10.1029/2022RS007579

* Bocanegra-Bahamón, T. M., Molera Calvés, G., Gurvits, L. I., Cimò, G., Dirkx, D., Duev, D. A., et al. (2019). Venus Express radio occultation observed by PRIDE. *Astronomy and Astrophysics*, 624, A59. https://doi.org/10.1051/0004-6361/201833160

* Born, M., & Wolf, E. (1980). In *Principles of optics* (6th ed.). Pergamon Press Ltd. Oxford, United Kingdom. (chap. 3).

* Bourgoin, A. (2020). General expansion of time transfer functions in optical spacetime. *Physical Review D*, 101(6), 064035. https://doi.org/10.1103/PhysRevD.101.064035

* Bourgoin, A., Gramigna, E., Zannoni, M., Gomez Casajus, L., & Tortora, P. (2022). Determination of uncertainty profiles in neutral atmospheric properties measured by radio occultation experiments. *Advances in Space Research*, 70(8), 2555-2570. https://doi.org/10.1016/j.asr.2022.07.015

* Bourgoin, A., Zannoni, M., Gomez Casajus, L., Tortora, P., & Teyssandier, P. (2021). Relativistic modeling of atmospheric occultations with time transfer functions. *Astronomy and Astrophysics*, 648, A46. https://doi.org/10.1051/0004-6361/202040269

* Bourgoin, A., Zannoni, M., & Tortora, P. (2019). Analytical ray-tracing in planetary atmospheres. *Astronomy and Astrophysics*, 624, A41. https://doi.org/10.1051/0004-6361/201834962

* Buccino, D., Parisi, M., Gramigna, E., Gomez-Casajus, L., Tortora, P., Zannoni, M., et al. (2022). Ganymede's ionosphere observed by a dual-frequency radio occultation with Juno. *Geophysical Research Letters*, 49(23), e2022GL098420. https://doi.org/10.1029/2022GL098420

* Buhner, K., Maurer, G., & Bender, E. (1981). Pressure-enthalpy diagrams for methane, ethane, propane, ethylene and propylene. *Cryogenics*, 21(3), 157-164. https://doi.org/10.1016/0011-2275(81)90267-8

* Dalba, P. A., & Withers, P. (2019). Cassini radio occultation observations of Titan's ionosphere: The complete set of electron density profiles. *Journal of Geophysical Research*, 124(1), 643-660. https://doi.org/10.1029/2018JA025693

## Acknowledgments

The authors are grateful to the Italian Space Agency (ASI) for financial support through ASI-INAF Agreement No. 2018-25-HH.0 and 2023-6-HH.0 in the context of ESA's JUICE mission, and ASI-CRAS Agreement No. 2022.16.HH.0 in the context of NASA's Juno mission.
Fjeldbo, G., & Eshleman, V. R. (1968). The atmosphere of Mars analyzed by integral inversion of the Mariner IV occultation data. *Planetary and Space Science*, **16**(8), 1035–1059. https://doi.org/10.1016/0032-0633(68)90020-2

Fjeldbo, G., Kliore, A. J., & Eshleman, V. R. (1971). The neutral atmosphere of Venus as studied with the Mariner V radio occultation experiments. *The Astronomical Journal*, **76**(2), 123–140. https://doi.org/10.1086/111096

Forget, F., Hourdin, F., Fournier, R., Hourdin, C., Talagrand, O., Collins, M., et al. (1999). Improved general circulation models of the Martian atmosphere from the surface to above 80 km. *Journal of Geophysical Research*, **104**(E10), 24155–24175. https://doi.org/10.1029/1999JE001025

Gramigna, E., Parisi, M., Buccino, D., Casajus, L. G., Zannoni, M., Bourgoin, A., et al. (2023). Analysis of NASA's DSN Venus Express radio occultation data for year 2014. *Advances in Space Research*, **71**(1), 1198–1215. https://doi.org/10.1016/j.asr.2022.10.070

Grandin, M., Bleilly, P.-L., Witasse, O., & Marchaudon, A. (2014). Mars Express radio-occultation data: A novel analysis approach. *Journal of Geophysical Research: Space Physics*, **119**(12), 10621–10632. https://doi.org/10.1002/2014JA020698

Häusler, B., Pätzold, M., Tyler, G. L., Simpson, R., Bird, M., Dehant, V., et al. (2006). Radio science investigations by VeRa onboard the Venus Express spacecraft. *Planetary and Space Science*, **54**(13–14), 1315–1335. https://doi.org/10.1016/j.pss.2006.04.032

Hinson, D. P. (1999). Archived uncertainty data on pressure and temperature for MGS occultation occurred on 27 December 1998 (8361m48a.tps) [Dataset]. NASA Planetary Data System. Retrieved from https://atmos.nmsu.edu/PDS/data/mors_1101/tps/1998_358/

Hinson, D. P. (2000). Archived uncertainty data on electron density for MGS occultation occurred on 27 December 1998 (8361m48a.eds) [Dataset]. NASA Planetary Data System. Retrieved from https://atmos.nmsu.edu/PDS/data/mors_1102/eds/1998_358/

Hinson, D. P., Linscott, I. R., Young, L. A., Tyler, G., Stern, S., Beyer, R., et al. (2017). Radio occultation measurements of Pluto's neutral atmosphere with New Horizons. *Icarus*, **290**, 96–111. https://doi.org/10.1016/j. icarus.2017.02.031

Hinson, D. P., Simpson, R. A., Twicken, J. D., Tyler, G. L., & Flasar, F. M. (1999). Initial results from radio occultation measurements with Mars Global Surveyor. *Journal of Geophysical Research*, **104**(E11), 26997–27012. https://doi.org/10.1029/1999JE001069

Jenkins, J. M., Steffes, P. G., Hinson, D. P., Twicken, J. D., & Tyler, G. (1994). Radio occultation studies of the Venus atmosphere with the Magellan spacecraft: 2. Results from the October 1991 experiments. *Icarus*, **110**(1), 79–94. https://doi.org/10.1006/icar.1994.1108

Kliore, A., Cain, D. L., Levy, G. S., Eshleman, V. R., Fjeldbo, G., & Drake, F. D. (1965). Occultation experiment: Results of the first direct measurement of Mars's atmosphere and ionosphere. *Science*, **149**(3689), 1243–1248. https://doi.org/10.1126/science.149.3689.1243

Lindal, G. F. (1992). The atmosphere of Neptune: An analysis of radio occultation data acquired with Voyager 2. *The Astronomical Journal*, **103**(3), 967–982. https://doi.org/10.1086/116119

Lindal, G. F., Lyons, J. R., Sweetnam, D. N., Eshleman, V. R., Hinson, D. P., & Tyler, G. L. (1987). The atmosphere of Uranus: Results of radio occultation measurements with Voyager 2. *Journal of Geophysical Research*, **92**(A13), 14987–15001. https://doi.org/10.1029/JA092iA13p14987

Lindal, G. F., Sweetnam, D. N., & Eshleman, V. R. (1985). The atmosphere of Saturn: An analysis of the Voyager radio occultation measurements. *The Astronomical Journal*, **90**(6), 1136–1146. https://doi.org/10.1086/113820

Lindal, G. F., Wood, G. E., Levy, J. D., Anderson, J. D., Sweetnam, D. N., Hotz, H. B., et al. (1981). The atmosphere of Jupiter: An analysis of the Voyager radio occultation measurements. *Journal of Geophysical Research*, **86**(A10), 8721–8727. https://doi.org/10.1029/JA086iA10p08721

Lipa, B., & Tyler, L. (1979). Statistical and computational uncertainties in atmospheric profiles from radio occultation: Mariner 10 at Venus. *Icarus*, **39**(2), 192–208. https://doi.org/10.1016/0019-1035(79)90163-5

Millour, E., Forget, F., Spiga, A., Vals, M., Zakharov, V., Montabone, L., et al. (2021). The Mars climate Database (version 5.3). In *Scientific workshop: From Mars Express to ExoMars*. ESAC.

Pätzold, M., Häusler, B., Bird, M. K., Tellmann, S., Mattei, R., Asmar, S. W., et al. (2007). The structure of Venus' middle atmosphere and ionosphere. *Nature*, **450**(7170), 657–660. https://doi.org/10.1038/nature06239

Pätzold, M., Häusler, B., Tyler, G. L., Andert, T., Asmar, S., Bird, M., et al. (2016). Mars Express 10 years at Mars: Observations by the Mars Express radio science experiment (MaRS). *Planetary and Space Science*, **127**, 44–90. https://doi.org/10.1016/j.pss.2016.02.013

Pätzold, M., Tellmann, S., Häusler, B., Hinson, D., Schaa, R., & Tyler, G. L. (2005). A sporadic third layer in the ionosphere of Mars. *Science*, **310**(5749), 837–839. https://doi.org/10.1126/science.1117755

Petricca, F., Cascioli, G., & Genova, A. (2021). A technique for the analysis of radio occultation data to retrieve atmospheric properties and associated uncertainties. *Radio Science*, **56**(5), e2020RS007205. https://doi.org/10.1029/2020RS007205

Phinney, R. A., & Anderson, D. L. (1968). On the radio occultation method for studying planetary atmospheres. *Journal of Geophysical Research*, **73**(5), 1819–1827. https://doi.org/10.1029/JA073i005p01819

Poisson, E., & Will, C. M. (2014). *Gravity*. Cambridge University Press. (chap. 2). https://doi.org/10.1017/CBO9781139507486

Powell, M. J. D. (1970). In P. Rabinowitz (Ed.), *Numerical methods for nonlinear algebraic equations. A Fortran subroutine for solving systems of nonlinear algebraic equations*. Gordon & Breach Science Publishers. (chap. 7).

Schinder, P. J. (2020). Atmospheric profiles derived from Radio Science data acquired during Cassini occultations by Titan [Dataset]. NASA Planetary Data System. Retrieved from https://atmos.nmsu.edu/PDS/data/PDS4/titan_profiles Rochet/data/

Schinder, P. J., Flasar, F. M., Marouf, E. A., French, R. G., Anabtawi, A., Barbinis, E., et al. (2020). The structure of Titan's atmosphere from Cassini radio occultations: One- and two-way occultations. *Icarus*, **345**, 113720. https://doi.org/10.1016/j.icarus.2020.113720

Schinder, P. J., Flasar, F. M., Marouf, E. A., French, R. G., Anabtawi, A., Barbinis, E., & Kliore, A. J. (2015). A numerical technique for two-way radio occultations by oblate axisymmetric atmospheres with zonal winds. *Radio Science*, **50**(7), 712–727. https://doi.org/10.1002/2015RS005690

Schinder, P. J., Flasar, F. M., Marouf, E. A., French, R. G., McGhee, C. A., Kliore, A. J., et al. (2011a). Saturn's equatorial oscillation: Evidence of descending thermal structure from Cassini radio occultations. *Geophysical Research Letters*, **38**(8), L08205. https://doi.org/10.1029/2011GL047191

Schinder, P. J., Flasar, F. M., Marouf, E. A., French, R. G., McGhee, C. A., Kliore, A. J., et al. (2011b). The structure of Titan's atmosphere from Cassini radio occultations. *Icarus*, **215**(2), 460–474. https://doi.org/10.1016/j. icarus.2011.07.030

Synge, J. L. (1960). *Relativity: The general theory*. North-Holland Publishing Company. (Russian translation: IL (Foreign Literature), Moscow, 1963).

Togni, A., Zannoni, M., Gomez Casajus, L., & Tortora, P. (2021). An FFT-based method for Doppler observables estimation in deep space tracking. In *2021 IEEE 8th international workshop on metrology for AeroSpace* (pp. 294–299). https://doi.org/10.1109/MetroAeroSpace51421.2021.9511674

Withers, P. (2010). Prediction of uncertainties in atmospheric properties measured by radio occultation experiments. *Advances in Space Research*, **46**(1), 58–73. https://doi.org/10.1016/j.asr.2010.03.004

Withers, P. (2014). Archived data for MGS occultation occurred on 27 December 1998 [Dataset]. Boston University. Retrieved from http://sirius.bu.edu/withers/radioocc2014suppinfo/

Withers, P. (2023). MAVEN ROSE electron density profile data collection. [Dataset]. NASA Planetary Data System. Retrieved from https://pds-ppi.igpp.ucla.edu/mission/MAVEN/MAVEN/ROSE
Fjeldbo, G., & Eshleman, V. R. (1968). The atmosphere of Mars analyzed by integral inversion of the Mariner IV occultation data. *Planetary and Space Science*, **16**(8), 1035–1059. https://doi.org/10.1016/0032-0633(68)90020-2

Fjeldbo, G., Kliore, A. J., & Eshleman, V. R. (1971). The neutral atmosphere of Venus as studied with the Mariner V radio occultation experiments. *The Astronomical Journal*, **76**(2), 123–140. https://doi.org/10.1086/111096

Forget, F., Hourdin, F., Fournier, R., Hourdin, C., Talagrand, O., Collins, M., et al. (1999). Improved general circulation models of the Martian atmosphere from the surface to above 80 km. *Journal of Geophysical Research*, **104**(E10), 24155–24175. https://doi.org/10.1029/1999JE001025

Gramigna, E., Parisi, M., Buccino, D., Casajus, L. G., Zannoni, M., Bourgoin, A., et al. (2023). Analysis of NASA's DSN Venus Express radio occultation data for year 2014. *Advances in Space Research*, **71**(1), 1198–1215. https://doi.org/10.1016/j.asr.2022.10.070

Grandin, M., Bleilly, P.-L., Witasse, O., & Marchaudon, A. (2014). Mars Express radio-occultation data: A novel analysis approach. *Journal of Geophysical Research: Space Physics*, **119**(12), 10621–10632. https://doi.org/10.1002/2014JA020698

Häusler, B., Pätzold, M., Tyler, G. L., Simpson, R., Bird, M., Dehant, V., et al. (2006). Radio science investigations by VeRa onboard the Venus Express spacecraft. *Planetary and Space Science*, **54**(13–14), 1315–1335. https://doi.org/10.1016/j.pss.2006.04.032

Hinson, D. P. (1999). Archived uncertainty data on pressure and temperature for MGS occultation occurred on 27 December 1998 (8361m48a.tps) [Dataset]. NASA Planetary Data System. Retrieved from https://atmos.nmsu.edu/PDS/data/mors_1101/tps/1998_358/

Hinson, D. P. (2000). Archived uncertainty data on electron density for MGS occultation occurred on 27 December 1998 (8361m48a.eds) [Dataset]. NASA Planetary Data System. Retrieved from https://atmos.nmsu.edu/PDS/data/mors_1102/eds/1998_358/

Hinson, D. P., Linscott, I. R., Young, L. A., Tyler, G., Stern, S., Beyer, R., et al. (2017). Radio occultation measurements of Pluto's neutral atmosphere with New Horizons. *Icarus*, **290**, 96–111. https://doi.org/10.1016/j. icarus.2017.02.031

Hinson, D. P., Simpson, R. A., Twicken, J. D., Tyler, G. L., & Flasar, F. M. (1999). Initial results from radio occultation measurements with Mars Global Surveyor. *Journal of Geophysical Research*, **104**(E11), 26997–27012. https://doi.org/10.1029/1999JE001069

Jenkins, J. M., Steffes, P. G., Hinson, D. P., Twicken, J. D., & Tyler, G. (1994). Radio occultation studies of the Venus atmosphere with the Magellan spacecraft: 2. Results from the October 1991 experiments. *Icarus*, **110**(1), 79–94. https://doi.org/10.1006/icar.1994.1108

Kliore, A., Cain, D. L., Levy, G. S., Eshleman, V. R., Fjeldbo, G., & Drake, F. D. (1965). Occultation experiment: Results of the first direct measurement of Mars's atmosphere and ionosphere. *Science*, **149**(3689), 1243–1248. https://doi.org/10.1126/science.149.3689.1243

Lindal, G. F. (1992). The atmosphere of Neptune: An analysis of radio occultation data acquired with Voyager 2. *The Astronomical Journal*, **103**(3), 967–982. https://doi.org/10.1086/116119

Lindal, G. F., Lyons, J. R., Sweetnam, D. N., Eshleman, V. R., Hinson, D. P., & Tyler, G. L. (1987). The atmosphere of Uranus: Results of radio occultation measurements with Voyager 2. *Journal of Geophysical Research*, **92**(A13), 14987–15001. https://doi.org/10.1029/JA092iA13p14987

Lindal, G. F., Sweetnam, D. N., & Eshleman, V. R. (1985). The atmosphere of Saturn: An analysis of the Voyager radio occultation measurements. *The Astronomical Journal*, **90**(6), 1136–1146. https://doi.org/10.1086/113820

Lindal, G. F., Wood, G. E., Levy, J. D., Anderson, J. D., Sweetnam, D. N., Hotz, H. B., et al. (1981). The atmosphere of Jupiter: An analysis of the Voyager radio occultation measurements. *Journal of Geophysical Research*, **86**(A10), 8721–8727. https://doi.org/10.1029/JA086iA10p08721

Lipa, B., & Tyler, L. (1979). Statistical and computational uncertainties in atmospheric profiles from radio occultation: Mariner 10 at Venus. *Icarus*, **39**(2), 192–208. https://doi.org/10.1016/0019-1035(79)90163-5

Millour, E., Forget, F., Spiga, A., Vals, M., Zakharov, V., Montabone, L., et al. (2021). The Mars climate Database (version 5.3). In *Scientific workshop: From Mars Express to ExoMars*. ESAC.

Pätzold, M., Häusler, B., Bird, M. K., Tellmann, S., Mattei, R., Asmar, S. W., et al. (2007). The structure of Venus' middle atmosphere and ionosphere. *Nature*, **450**(7170), 657–660. https://doi.org/10.1038/nature06239

Pätzold, M., Häusler, B., Tyler, G. L., Andert, T., Asmar, S., Bird, M., et al. (2016). Mars Express 10 years at Mars: Observations by the Mars Express radio science experiment (MaRS). *Planetary and Space Science*, **127**, 44–90. https://doi.org/10.1016/j.pss.2016.02.013

Pätzold, M., Tellmann, S., Häusler, B., Hinson, D., Schaa, R., & Tyler, G. L. (2005). A sporadic third layer in the ionosphere of Mars. *Science*, **310**(5749), 837–839. https://doi.org/10.1126/science.1117755

Petricca, F., Cascioli, G., & Genova, A. (2021). A technique for the analysis of radio occultation data to retrieve atmospheric properties and associated uncertainties. *Radio Science*, **56**(5), e2020RS007205. https://doi.org/10.1029/2020RS007205

Phinney, R. A., & Anderson, D. L. (1968). On the radio occultation method for studying planetary atmospheres. *Journal of Geophysical Research*, **73**(5), 1819–1827. https://doi.org/10.1029/JA073i005p01819

Poisson, E., & Will, C. M. (2014). *Gravity*. Cambridge University Press. (chap. 2). https://doi.org/10.1017/CBO9781139507486

Powell, M. J. D. (1970). In P. Rabinowitz (Ed.), *Numerical methods for nonlinear algebraic equations. A Fortran subroutine for solving systems of nonlinear algebraic equations*. Gordon & Breach Science Publishers. (chap. 7).

Schinder, P. J. (2020). Atmospheric profiles derived from Radio Science data acquired during Cassini occultations by Titan [Dataset]. NASA Planetary Data System. Retrieved from https://atmos.nmsu.edu/PDS/data/PDS4/titan_profiles Rochet/data/

Schinder, P. J., Flasar, F. M., Marouf, E. A., French, R. G., Anabtawi, A., Barbinis, E., et al. (2020). The structure of Titan's atmosphere from Cassini radio occultations: One- and two-way occultations. *Icarus*, **345**, 113720. https://doi.org/10.1016/j.icarus.2020.113720

Schinder, P. J., Flasar, F. M., Marouf, E. A., French, R. G., Anabtawi, A., Barbinis, E., & Kliore, A. J. (2015). A numerical technique for two-way radio occultations by oblate axisymmetric atmospheres with zonal winds. *Radio Science*, **50**(7), 712–727. https://doi.org/10.1002/2015RS005690

Schinder, P. J., Flasar, F. M., Marouf, E. A., French, R. G., McGhee, C. A., Kliore, A. J., et al. (2011a). Saturn's equatorial oscillation: Evidence of descending thermal structure from Cassini radio occultations. *Geophysical Research Letters*, **38**(8), L08205. https://doi.org/10.1029/2011GL047191

Schinder, P. J., Flasar, F. M., Marouf, E. A., French, R. G., McGhee, C. A., Kliore, A. J., et al. (2011b). The structure of Titan's atmosphere from Cassini radio occultations. *Icarus*, **215**(2), 460–474. https://doi.org/10.1016/j. icarus.2011.07.030

Synge, J. L. (1960). *Relativity: The general theory*. North-Holland Publishing Company. (Russian translation: IL (Foreign Literature), Moscow, 1963).

Togni, A., Zannoni, M., Gomez Casajus, L., & Tortora, P. (2021). An FFT-based method for Doppler observables estimation in deep space tracking. In *2021 IEEE 8th international workshop on metrology for AeroSpace* (pp. 294–299). https://doi.org/10.1109/MetroAeroSpace51421.2021.9511674

Withers, P. (2010). Prediction of uncertainties in atmospheric properties measured by radio occultation experiments. *Advances in Space Research*, **46**(1), 58–73. https://doi.org/10.1016/j.asr.2010.03.004

Withers, P. (2014). Archived data for MGS occultation occurred on 27 December 1998 [Dataset]. Boston University. Retrieved from http://sirius.bu.edu/withers/radioocc2014suppinfo/

Withers, P. (2023). MAVEN ROSE electron density profile data collection. [Dataset]. NASA Planetary Data System. Retrieved from https://pds-ppi.igpp.ucla.edu/mission/MAVEN/MAVEN/ROSE
Withers, P., Felici, M., Mendillo, M., Moore, L., Narvaez, C., Vogt, M. F., et al. (2020). The MAVEN radio occultation science experiment (ROSE). *Space Science Reviews*, **126**(4), 61. https://doi.org/10.1007/s11214-020-00687-6

Withers, P., Felici, M., Mendillo, M., Moore, L., Narvaez, C., Vogt, M. F., & Jakosky, B. M. (2018). First ionospheric results from the MAVEN radio occultation science experiment (ROSE). *Journal of Geophysical Research: Space Physics*, **123**(5), 4171–4180. https://doi.org/10.1029/2018JA025182

Withers, P., & Huber, L. (2020). Cassini orbiter radio science subsystem occultation and electron density data archive [Dataset]. NASA Planetary Data System Atmospheres (ATM) Node. https://doi.org/10.17189/1518927

Withers, P., & Moore, L. (2020). How to process radio occultation data: 2. From time series of two-way, single-frequency frequency residuals to vertical profiles of ionospheric properties. *Radio Science*, **55**(8), e2019RS007046. https://doi.org/10.1029/2019RS007046

Withers, P., Moore, L., Cahoy, K., & Beerer, I. (2014). How to process radio occultation data: 1. From time series of frequency residuals to vertical profiles of atmospheric and ionospheric properties. *Planetary and Space Science*, **101**, 77–88. https://doi.org/10.1016/j.pss.2014.06.011

Withers, P., Towner, M. C., Hathi, B., & Zarnnecki, J. C. (2003). Analysis of entry accelerometer data: A case study of Mars Pathfinder. *Planetary and Space Science*, **51**(9–10), 541–561. https://doi.org/10.1016/S0032-0633(03)00077-1
Withers, P., Felici, M., Mendillo, M., Moore, L., Narvaez, C., Vogt, M. F., et al. (2020). The MAVEN radio occultation science experiment (ROSE). *Space Science Reviews*, 126(4), 61. https://doi.org/10.1007/s11214-020-00687-6

Withers, P., Felici, M., Mendillo, M., Moore, L., Narvaez, C., Vogt, M. F., & Jakosky, B. M. (2018). First ionospheric results from the MAVEN radio occultation science experiment (ROSE). *Journal of Geophysical Research: Space Physics*, 123(5), 4171–4180. https://doi.org/10.1029/2018JA025182

Withers, P., & Huber, L. (2020). Cassini orbiter radio science subsystem occultation and electron density data archive [Dataset]. NASA Planetary Data System Atmospheres (ATM) Node. https://doi.org/10.17189/1518927

Withers, P., & Moore, L. (2020). How to process radio occultation data: 2. From time series of two-way, single-frequency frequency residuals to vertical profiles of ionospheric properties. *Radio Science*, 55(8), e2019RS007046. https://doi.org/10.1029/2019RS007046

Withers, P., Moore, L., Cahoy, K., & Beerer, I. (2014). How to process radio occultation data: 1. From time series of frequency residuals to vertical profiles of atmospheric and ionospheric properties. *Planetary and Space Science*, 101, 77–88. https://doi.org/10.1016/j.pss.2014.06.011

Withers, P., Towner, M. C., Hathi, B., & Zarnnecki, J. C. (2003). Analysis of entry accelerometer data: A case study of Mars Pathfinder. *Planetary and Space Science*, 51(9–10), 541–561. https://doi.org/10.1016/S0032-0633(03)00077-1