L E Bonnefoy

A Le Gall

E Lellouch

C Leyrat

M Janssen

R Sultana

Rhea's subsurface probed by the Cassini radiometer: insights into its thermal, structural, and compositional properties

Highlights

 Cassini radiometry and radar observations of Rhea at 2.2 cm-wavelength are reduced and analyzed.  Antenna temperatures are simulated by combining thermal, radiative, and emissivity models.  Large probed depths (~10 m) imply high purity of the icy regolith, especially around Inktomi crater.  High backscatter, low emissivity and dielectric constant suggest buried scattering structures.  South pole thermal inertia is higher than found in IR pointing to increasing compaction with depth.

Introduction

During its 13.5 years in the Saturn's system, the Cassini spacecraft has progressively revealed the variety of the surfaces of Saturn's mid-sized inner icy satellites, namely, by order of distance to the planet: Mimas, Enceladus, Tethys, Dione, and Rhea. Though probably all originating from Saturn's rings [START_REF] Charnoz | Accretion of Saturn's mid-sized moons during the viscous spreading of young massive rings: Solving the paradox of silicatepoor rings versus silicate-rich moons[END_REF], these surfaces have evolved differently due to their respective thermal histories and the dominant exogenic processes at play at their location in the Saturn's system. In particular, they have been processed by interactions with Saturn's dust rings, magnetosphere, and micro-meteoritic bombardment; their surfaces and near subsurfaces contain a record of these interactions and processes.

In particular, previous studies have found that Saturn's inner moons are among the brightest objects in the Solar System at centimetric wavelengths with anomalously high radar albedos (disk-integrated same-sense linear Radar albedo of 𝐴 𝑆𝐿 𝑑𝑖𝑠𝑘 = 1.2 -3.5 at 2.2 cm wavelength; Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF], hard to reconcile with commonly invoked scattering processes and their microwave emissivities (𝑒 = 0.45 -0.7, depending on satellite, Bond albedo, and temperature assumptions; [START_REF] Ostro | Cassini RADAR observations of Enceladus, Tethys, Dione, Rhea, Iapetus, Hyperion, and Phoebe[END_REF][START_REF] Ries | A large-scale anomaly in Enceladus' microwave emission[END_REF][START_REF] Gall | Thermally anomalous features in the subsurface of Enceladus's south polar terrain[END_REF]. While the reasons for these extreme radar albedos are not yet fully understood, they likely stem from significant volume scattering in the subsurface. Indeed, in a weakly absorbing medium such as a water ice regolith, voids and inhomogeneities cause internal reflections, which could enhance the reflectivity in the backscattering direction, and, concurrently decrease the probability of a wave emitted from the subsurface to escape (that is, decrease the surface emissivity). Most specifically, the coherent backscattering effect (CBE; [START_REF] Hapke | Coherent backscatter and the radar characteristics of outer planet satellites[END_REF] is so far the best explanation for the high radar albedos and polarization ratios of the Galilean satellites (Black et al., 2001). Though it is insufficient to explain the radar albedos of Enceladus, Mimas and Tethys (Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF] and even of the most radar-bright regions of Titan (Janssen et al., 2011), this effect may be invoked also for Dione and Rhea which at 2.2 cm are as radarbright as Europa (Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF].

The results previously mentioned are mainly based on data collected in the active mode of the RADAR on board the Cassini spacecraft, which operated at 2.2 cm (13.8 GHz, Ku band) [START_REF] Elachi | Radar: The Cassini Titan Radar Mapper[END_REF]. This instrument also included a passive (or radiometry) mode designed to record the thermal emission emanating from a surface at a wavelength of 2.2 cm. Through the measurement of antenna temperatures, the microwave radiometer data provide constraints on the chemical and thermophysical properties of a surface down to a depth of up to several meters for Saturn's icy satellites [START_REF] Ostro | Cassini RADAR observations of Enceladus, Tethys, Dione, Rhea, Iapetus, Hyperion, and Phoebe[END_REF][START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF]. On Iapetus and Enceladus, Cassini radiometry observations have indeed led to the identification of compositional [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF], structural [START_REF] Ries | A large-scale anomaly in Enceladus' microwave emission[END_REF], and thermal (Le [START_REF] Gall | Thermally anomalous features in the subsurface of Enceladus's south polar terrain[END_REF] anomalies. On Titan, they were used to map the surface dielectric constant and emissivity [START_REF] Janssen | Titan's surface at 2.2-cm wavelength imaged by the Cassini RADAR radiometer: Calibration and first results[END_REF]2016) and proved to be key to investigate the composition of landforms, such as dunes (e.g., [START_REF] Gall | Cassini SAR, radiometry, scatterometry, and altimetry observations of Titan's dune fields[END_REF][START_REF] Bonnefoy | Compositional and spatial variations in Titan dune and interdune regions from Cassini VIMS and RADAR[END_REF], craters (e.g., [START_REF] Werynski | Compositional variations of Titan's impact craters indicates active surface erosion[END_REF], and lakes (e.g., [START_REF] Gall | Composition, seasonal change, and bathymetry of Ligeia Mare, Titan, derived from its microwave thermal emission[END_REF]. After Titan, Rhea is the satellite which the Cassini radiometer observed the most; yet this dataset has so far not been exploited.

Rhea is the largest of Saturn's icy satellites (with a radius 𝑅 𝑅ℎ𝑒𝑎 =1527 km) and the furthest from Saturn (9𝑅 𝑆𝑎𝑡𝑢𝑟𝑛 ) among those hypothesized to have formed from the rings [START_REF] Charnoz | Accretion of Saturn's mid-sized moons during the viscous spreading of young massive rings: Solving the paradox of silicatepoor rings versus silicate-rich moons[END_REF]. It shows a heavily cratered surface, with bright fractures across its trailing hemisphere constituting a wispy terrain very similar to that also observed on Dione. Its leading hemisphere is marked by a fresh crater (~280 Ma old) called Inktomi (14°S, -112°E), and its bright icy ejecta [START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF][START_REF] Stephan | The Saturnian satellite Rhea as seen by Cassini VIMS[END_REF][START_REF] Scipioni | Spectroscopic classification of icy satellites of Saturn II: Identification of terrain units on Rhea[END_REF][START_REF] Dalle Ore | Impact craters: An ice study on Rhea[END_REF]. Both leading and trailing hemispheres show a distinct reddening, while the leading hemisphere is brighter, interpreted as a consequence of the infall of bright icy E-ring particles [START_REF] Verbiscer | Enceladus: Cosmic Graffiti Artist Caught in the Act[END_REF][START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF][START_REF] Stephan | The Saturnian satellite Rhea as seen by Cassini VIMS[END_REF]. Even though it is further away from the E-ring and seemingly older (i.e., more heavily cratered), Rhea is on average more radar-bright than Dione [START_REF] Ostro | Cassini RADAR observations of Enceladus, Tethys, Dione, Rhea, Iapetus, Hyperion, and Phoebe[END_REF]2010;[START_REF] Black | Arecibo radar observations of Rhea, Dione, Tethys, and Enceladus[END_REF]Le Gall et al., 2019). Rhea's enhanced radar brightness suggests a less contaminated water ice regolith, which may be explained by weaker space weathering further from Saturn and its intense radiation belts (Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF]. In this paper, we further investigate Rhea's subsurface properties and regional variations by analyzing all passive radiometry observations acquired on this satellite during the full Cassini mission.

Section 2 reviews the characteristics of the Cassini radiometry dataset, and of one single resolved active scatterometry observation which we included in our analysis. Section 3 details the combination of thermal, radiative transfer, and emissivity models which simulates the radiometry data. Section 4 presents the results of model fitting to the data and the derived new constraints on the thermophysical and compositional properties of the Rhea's subsurface. Lastly, in Section 5, we present our interpretation for the composition and structure of the nearsubsurface in different regions of Rhea.

Observations and data reduction

Rhea was observed by the Cassini Radar/radiometer in the course of 9 flybys, with resolutions ranging from 0.85𝑅 𝑅ℎ𝑒𝑎 to 0.1𝑅 𝑅ℎ𝑒𝑎 . The geometry and calibration of these data are described herein.

Observations of Rhea by the Cassini Radar/radiometer

Three different types of Cassini radiometry observations were used: "distant scans", "resolved scans", and the so-called "stares". All three types are shown in Fig. 1 on the example of the RH127 Rhea flyby that occurred on March 2, 2010.

Distant scans were designed for radiometry: the antenna scanned the target, regularly moving off-target for baseline calibration. Despite their primary objective being the measurement of the disk-integrated brightness temperatures, the distance to Rhea can significantly vary from one scan to another (from 65 000 to 200 000 km), and several distant scans do resolve features at the surface, with beam radii of at best 0.3𝑅 𝑅ℎ𝑒𝑎 (see Table 1 and Fig. 2). In addition, the integration time is longer (about 1 s) during distant observations than during stares, implying low photon noise.

A single resolved scan was performed on Rhea (during the aforementioned RH127 flyby), during which active scatterometry and passive radiometry data were simultaneously collected (with the same viewing geometry and spatial resolution) from a distance of 25 000 to 50 000 km. The resolved scatterometry data were initially calibrated by [START_REF] Wye | Radar scattering from Titan and Saturn's icy satellites using the Cassini spacecraft[END_REF]. The associated resolved radiometry observation is particularly useful as it has the highest available resolution (down to 0.1𝑅 𝑅ℎ𝑒𝑎 ) and can be analyzed jointly with concurrent scatterometry. However, it is also especially difficult to calibrate because i) the radar transmitter was ON during the scan, leading to instrument warming, ii) due to the lack of a "cold sky" observation either during or after the scan.

Stares were primarily intended for calculation of the average radar albedo of the observed hemisphere: the radar antenna targeted a single point on the surface for a long period (5-10 minutes) while many short bursts were transmitted to the surface. Radiometry measurements were performed in-between burst cycles, with a shorter integration time (<0.3 s) than for distant scans, leading to a higher noise. Analyses of the derived disk-integrated radar albedos have been published in [START_REF] Ostro | Cassini RADAR observations of Enceladus, Tethys, Dione, Rhea, Iapetus, Hyperion, and Phoebe[END_REF][START_REF] Ostro | New Cassini RADAR results for Saturn's icy satellites[END_REF] and Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF]. Calibration of the radiometry data acquired during stares is particularly challenging because i) the antenna does not move off-target, therefore a baseline calibration during the stare is impossible and ii) the fact that the transmitter is ON affects the receiver temperatures, especially if the integration time is short (long integration times let the receiver cool down in-between bursts). We therefore use caution whenever interpreting these data. Furthermore, they are only moderately useful as coverage is restricted to a few points on Rhea.

During the 9 different flybys of Rhea with microwave radiometry observations, there were 10 distant scans, 10 stares, and 1 resolved scan, as reported in Table 1. Note that the usable radiometry scans of Rhea were all taken either at night (local times at the sub-spacecraft point between 8.30 pm and midnight) or at the South pole during the southern summer and fall. As a consequence, there is no equatorial daytime data available. There is also little to no data on the North pole and on the trailing side of Rhea, while the South pole and sub-Saturn regions are the best sampled. (PIA 18438), projected to spacecraft view coordinates. Bottom: Calibrated antenna temperatures during the RH127 observation. The baseline is set to the Cosmic Microwave Background (CMB) temperature, i.e. 2.7 K with the Rayleigh-Jeans approximation (valid at 2.2 cm). The moment the radar transmitter turns ON is indicated with a red arrow. Periods when the radar transmitter was ON and the integration time was shorter than 0.3 s are indicated in red: during these periods, the background temperature is higher due to the heating of the receiver. We therefore only use the resolved scan data acquired after minute 366, when the receiver temperature is stable (see Section 2.2.2).

Radiometry data calibration

We follow the calibration method developed for Titan and described in [START_REF] Janssen | Titan's surface at 2.2-cm wavelength imaged by the Cassini RADAR radiometer: Calibration and first results[END_REF]2016), and adapt it to airless icy satellites such as Rhea (i.e., removing absorption in Titan's atmosphere). After using the antenna gain to convert from radiometer counts into Kelvins, the antenna temperatures measured by the Cassini radiometer correspond to the antenna beam pattern convolved with the brightness temperature map of Rhea, as viewed by the instrument, with an added temperature baseline 𝑇 𝑧𝑒𝑟𝑜 :

𝑇 𝑎,0 (𝑡) = ∫ ∫ 𝑇 𝑏 (𝜃 ′ , 𝜑 ′ )𝐺(𝜃 ′ -𝜃(𝑡), 𝜑 ′ -𝜑(𝑡))𝑠𝑖𝑛(𝜃 ′ )𝑑𝜃′𝑑𝜑′ 𝜋 0°2 𝜋 0°+ 𝑇 𝑧𝑒𝑟𝑜 (1)
where 𝜃(𝑡) and 𝜑(𝑡) correspond to the antenna boresight pointing direction at time 𝑡 and 𝑇 𝑏 (𝜃, 𝜑) is the brightness temperature distribution with 𝐺, the normalized gain of the antenna radiation pattern, acting as a spatial weighting function. The antenna boresight direction is calculated using the SPICE toolkit by NASA's Navigation and Ancillary Information Facility (NAIF) [START_REF] Acton | Ancillary data services of NASA's Navigation and Ancillary Information Facility[END_REF], while talking into account a pointing offset (see Appendix A for details).

The calibration then includes removal of the far sidelobe (FSL) contribution 𝑇 𝐹𝑆𝐿 and subtraction of the baseline offset 𝑇 𝑧𝑒𝑟𝑜 (including both the instrumental temperature background and the CMB (Cosmic Microwave Background) temperature of 2.7 K) obtained from observing the empty cold sky: 𝑇 𝑎 = 𝑇 𝑎,0 -𝑇 𝐹𝑆𝐿 -𝑇 𝑧𝑒𝑟𝑜 .

The uncertainty 𝜎 associated with our final value of the antenna temperature 𝑇 𝑎 is composed of i) random Gaussian noise 𝜎 𝐺 including instrumental noise and photon noise decreasing as the inverse square root of integration time, and ii) a 1% calibration uncertainty that globally takes into account uncertainties in the gain and the calibration factors [START_REF] Janssen | Titan's surface at 2.2-cm wavelength imaged by the Cassini RADAR radiometer: Calibration and first results[END_REF]2016). We determine 𝜎 𝐺 empirically by examining the standard deviation of the sky observations for each flyby (see Appendix B). The uncertainty on the antenna temperatures 𝑇 𝑎 is then given by the following equation:

𝜎 2 = 𝜎 𝐺 2 + (0.01𝑇 𝑎 ) 2
(2)

Removal of the far sidelobe contribution

The beam pattern is composed of a Gaussian main beam with a half-power width of 0.37° and near sidelobes up to 2° that were both mapped scanning the Sun during the Cassini cruise, and far sidelobes (FSL) extending from 2° to 90°, all of which are given in [START_REF] Janssen | Titan's surface at 2.2-cm wavelength imaged by the Cassini RADAR radiometer: Calibration and first results[END_REF]. We estimate the FSL contribution 𝑇 𝐹𝑆𝐿 (𝑡) as a function of time 𝑡 throughout each scan by integrating Eq. ( 1) with 𝜃′ from 2° to 90°, while assuming for 𝑇 𝑏 (𝜃, 𝜑) a disk of constant temperature 𝑇 𝑑𝑖𝑠𝑘 on Rhea, against a 2.7 K sky background. If Saturn is present in the FSL, we also include it as a disk with the latitudinal 2.2-cm temperature profile derived by [START_REF] Janssen | Saturn's thermal emission at 2.2-cm wavelength as imaged by the Cassini RADAR radiometer[END_REF].

The disk temperature 𝑇 𝑑𝑖𝑠𝑘 is initially assumed to be 50 K to calculate a first estimate of 𝑇 𝐹𝑆𝐿 . We remove 𝑇 𝐹𝑆𝐿 and the baseline offset 𝑇 𝑧𝑒𝑟𝑜 (see Section 2.2.2) from the measured 𝑇 𝑎,0 to obtain a calibrated 𝑇 𝑎 . We then simulate the antenna temperatures by convolving a limbdarkened disk of uniform temperature 𝑇 𝑑𝑖𝑠𝑘 with the beam pattern along the scan path. The limb-darkening consists in a simple cos 0.04 (𝜃) brightness taper, where 𝜃 is the emission angle, as described in Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF]. We fit these simple simulated antenna temperatures to the observed calibrated 𝑇 𝑎 to obtain the best value of 𝑇 𝑑𝑖𝑠𝑘 for each distant scan. The value for 𝑇 𝑑𝑖𝑠𝑘 , which ranges from 42.5 K for the South pole flyby to 54.3 K in equatorial regions flybys, is then used to recalculate the FSL. For the resolved scan and stare observations, we assume the same disk temperature as the distant scan(s) undergone during the same flyby (because they are close enough in time that they observed the same region at the same local time). Overall, the Rhea observation resolution is low enough that the FSL always contribute less than 5% to the signal; it is therefore not necessary to include surface temperature variations when calculating 𝑇 𝐹𝑆𝐿 .

Determination of the baseline offset (zero level)

The temperature baseline offset 𝑇 𝑧𝑒𝑟𝑜 is determined after having subtracted the FSL contribution. During a distant scan, the antenna regularly moves away from the satellite disk and observes the sky for baseline calibration, leading to a comb-shaped time sequence (Fig. 1). Each time the radiometer points away from Rhea, we average the antenna temperatures of the 20 points with the least contribution from Rhea, and interpolate between observations. This zero level is mainly influenced by the instrument physical temperature, but if the antenna does not move far enough from Rhea, it may also include a contribution from the near sidelobes. This is taken into account by subtracting a zero level in the exact same way whenever we compute a simulated 𝑇 𝑎 (see Section 3.4).

The stares and the resolved scan, which were acquired while the Cassini radar transmitter was ON, are more challenging to calibrate. Indeed, the transmitter causes an increase in the temperature of the receiver, which plateaus at 2-3 K after 10-20 minutes (depending on the integration time); for example, the temperature increases during the RH127_3s stare in Fig. 1. In addition, in some instances, the transmitter was turned successively ON and OFF during the stares, causing jumps in the radiometry temperature measurements. Similarly, for longer integration times, the receiver cools down in-between bursts, leading to lower temperatures. Thus, if the integration time changes, the observation is too short, or the transmitter turns OFF before the temperature is stable, we simply cannot use the radiometry data. This is unfortunately the case for RH018_1s, RH022_1s, RH027_1s, RH047_1s, RH049_2s, and RH177_3s.

Several stares include radiometry data collected before the transmitter is turned ON: these data can be used without specific caution, but are only moderately useful as they only cover a single observation point. This is the case of RH022_4s, RH045_2s, and RH127_3s. The baseline 𝑇 𝑧𝑒𝑟𝑜 in this case is simply a linear interpolation of the baseline before the stare began (the baseline after the stare cannot be used as it is affected by the radar transmitter).

Lastly, the temperature is stable during part of the resolved scan RH127_4r. The cold sky observation directly before the resolved scan is only partially affected by the instrument warming because while the transmitter was ON, the integration time was 0.9 s, which is long enough for the receiver to cool down in-between bursts. The integration time then dropped at the beginning of the resolved scan (at minute 358 in Fig. 1), leading to a warming of the receiver over the first few minutes of resolved data. In order to assess the amplitude of the receiver heating and estimate after how long the system temperature stabilized, we applied the following steps. We first deconvolved (see Section 2.2.3 for the deconvolution method) the last half of the resolved data (minutes 380-415 in Fig. 1) which, >20 min after the integration time drop, have a stable receiver temperature. The temperature map thus obtained was then convolved with the beam pattern to obtain a simulated time-sequence. We find that the residual between the simulated time-sequence and the resolved observation stabilizes after 8 minutes; we therefore only kept the data acquired after this time (minute 366 in Fig. 1). We then calculated the disk-integrated temperature for these data and the two distant scans performed during the same flyby (i.e., on the same disk between minutes 230 and 325, see Fig. 1). Lastly, the difference in disk-integrated brightness temperature was subtracted from the resolved data, ensuring that this dataset is consistent with the unresolved and less noisy one.

Deconvolution

The goal of deconvolution is to map the brightness temperature distribution on Rhea at the resolution of the projected main beam size (given in Table 1). Once the FSL contribution and the baseline offset are removed, the calibrated antenna temperatures correspond to the antenna beam pattern (up to 2°) convolved with the brightness temperature map of Rhea, as viewed the instrument: 𝑇 𝑎 (𝑡) = ∫ ∫ 𝑇 𝑏 (𝜃 ′ , 𝜑 ′ )𝐺(𝜃 ′ -𝜃(𝑡), 𝜑 ′ -𝜑(𝑡))𝑠𝑖𝑛(𝜃 ′ )𝑑𝜃′𝑑𝜑′ 2°0 °2𝜋 0°

(3) where, again, 𝑇 𝑏 (𝜃, 𝜑) is the spatial distribution of the brightness temperature on Rhea as seen from the Cassini spacecraft and taking into account foreshortening effects.

We follow the iterative deconvolution approach applied by [START_REF] Zhang | Cassini microwave observations provide clues to the origin of Saturn's C ring[END_REF] on the Cassini radiometry data acquired on Saturn's C ring, with a small modification to adapt for the Rhea observations. We initially assume that Rhea is a disk of uniform temperature of 50 K (the initial value has no impact on the final result), to build a uniform surface temperature distribution 𝑇 𝑚𝑎𝑝 as seen from the instrument. We then convolve the antenna beam pattern with this map using Eq. (3) in order to compute a simulated antenna temperature 𝑇 𝑎,𝑠𝑖𝑚 . The residual Δ𝑇 = 𝑇 𝑎 -𝑇 𝑎,𝑠𝑖𝑚 is mapped by computing an average, weighted by the beam pattern, of all residuals. The residual map is then added to the initial temperature map 𝑇 𝑚𝑎𝑝 . We reiterate until the standard deviation of the residual is smaller than the mean theoretical noise level (Eq. ( 2)). The resolution of the map thus obtained is equal to the antenna main beam size, whose halfpower width is 0.37°. Any details below this resolution are artefacts of the deconvolution, originating from spatially overlapping data with different temperatures.

The results of the iterative deconvolution method are shown in Fig. 2 for the resolved RH127_4r observation and all distant scans of Rhea. We recall that all these observations were acquired on the night side or on the South pole. In all scans that do not include the poles, we notice a significant limb darkening, caused by the variations of the emissivity with the emission angle, averaged over the beam size. In spite of this limb darkening, temperature changes with latitude are clearly visible, as well as a cold area spatially correlated with the ejecta blanket of the young impact crater Inktomi (best visible in the higher-resolution data RH127 and RH177). Using CIRS (Cassini InfraRed Spectrometer) data which probes at most a few millimeters into the subsurface, [START_REF] Howett | Thermophysical property variations across Dione and Rhea[END_REF] found that the Inktomi ejecta blanket has a thermal inertia higher than its surroundings. However, at 2.2 cm with the Cassini radiometer, this area stands out (especially in the resolved data, RH127_4r; see Fig. 2) as a cold spot during the night; the exact opposite behavior would be expected for a high thermal inertia region with equal albedo and emissivity properties. This strongly suggests that the Inktomi ejecta blanket is associated with a locally low emissivity and/or high albedo, that overwhelms the thermal inertia effect.

Reduction of the active scatterometry data

During the resolved observation RH127_4r, passive radiometry observations were collected in between active scatterometry bursts. The scatterometry data have been reduced and published by [START_REF] Wye | Radar scattering from Titan and Saturn's icy satellites using the Cassini spacecraft[END_REF]; calibrated 𝜎 𝑎𝑣𝑔 0 were derived inverting the following form of the radar equation:

𝑃 𝑠 = 𝑃 𝑡 𝐺 𝑡 2 𝜆 2 (4𝜋) 3 𝜎 𝑎𝑣𝑔 0 ∫ 𝐺 2 (𝑑𝐴)𝑑𝐴 𝑅 4 (𝑑𝐴) 𝐴 𝑟𝑒𝑠 (4)
where 𝑃 𝑠 is the received echo signal power, 𝑃 𝑡 is the total transmitted power, 𝐺 𝑡 is the antenna peak gain, 𝜆 is the wavelength, 𝐺(𝑑𝐴) is the normalized gain of the antenna pattern at surface element 𝑑𝐴, 𝑅(𝑑𝐴) is the distance from the spacecraft to that surface element, and 𝜎 𝑎𝑣𝑔 0 is the beam-averaged normalized radar cross-section.

Icy satellite surfaces are primarily diffuse scattering surfaces, with little to no specular component [START_REF] Ostro | Cassini RADAR observations of Enceladus, Tethys, Dione, Rhea, Iapetus, Hyperion, and Phoebe[END_REF][START_REF] Ostro | New Cassini RADAR results for Saturn's icy satellites[END_REF][START_REF] Wye | Radar scattering from Titan and Saturn's icy satellites using the Cassini spacecraft[END_REF][START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF]; 𝜎 𝑎𝑣𝑔 0 therefore primarily varies as a function of the effective incidence angle 𝜃 𝑒𝑓𝑓 in the radar footprint (as defined in Wye, 2011) following a cosine power law of the form 𝐴𝑐𝑜𝑠 𝑛 𝜃 𝑒𝑓𝑓 where 𝐴 is a constant. This is the basis on which [START_REF] Wye | Radar scattering from Titan and Saturn's icy satellites using the Cassini spacecraft[END_REF] initially calibrated the RH127_4r scatterometry dataset, computing 𝜃 𝑒𝑓𝑓 for each observational time. However, since Wye (2011), a pointing offset was discovered in the antenna boresight coordinates of the Cassini radar. As reported by [START_REF] Zhang | Cassini microwave observations provide clues to the origin of Saturn's C ring[END_REF] and by [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF], the actual pointing direction is shifted by 0.036° (i.e., 10% of the main beam size) compared to the z-axis of the spacecraft. The offset in pointing affects both the effective incidence angles 𝜃 𝑒𝑓𝑓 and the distances to the surface in the boresight direction 𝑅(𝑑𝐴). We therefore re-analyzed the data while taking this offset into account, following the same approach as described in [START_REF] Wye | Radar scattering from Titan and Saturn's icy satellites using the Cassini spacecraft[END_REF]; the results are shown in Fig. 3. We derived the values of 𝐴 and 𝑛 from a least-square cosine power law fit through the data in three different regions: the Inktomi ejecta blanket (25° in longitude and latitude from the center of the crater), the trailing hemisphere (longitude 90 ± 70°E), and the leading hemisphere (longitude -90 ± 70°E, excluding the Inktomi ejecta blanket). The normalized disk-integrated radar albedo can then be estimated as follows [START_REF] Ostro | Cassini RADAR observations of Enceladus, Tethys, Dione, Rhea, Iapetus, Hyperion, and Phoebe[END_REF][START_REF] Wye | Radar scattering from Titan and Saturn's icy satellites using the Cassini spacecraft[END_REF][START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF]:

𝐴 𝑎𝑣𝑔 = 2𝐴 𝑛 + 1 (5)
Results are reported in Table 2.

Using a cosine power law fit through all the data to correct to a reference incidence angle of 32°, [START_REF] Wye | Radar scattering from Titan and Saturn's icy satellites using the Cassini spacecraft[END_REF] found an enhanced radar brightness around the Inktomi crater and on the trailing hemisphere of Rhea, and a decrease towards the northern high latitudes. Adopting the same method while taking the pointing offset into account, we now observe a mostly uniform backscattering surface except near the Inktomi crater, which remains much brighter than its surroundings (Fig. 3B). The leading hemisphere is slightly more radar-bright than the trailing one, which may indicate less contamination by non-ice impurities, consistent with the observations in the visible domain (e.g., [START_REF] Verbiscer | Enceladus: Cosmic Graffiti Artist Caught in the Act[END_REF]. However, the leading hemisphere enhancement in radar albedo is small enough that it can also be explained by contribution from the Inktomi crater ejecta blanket, which lies either in the antenna main beam or in the near sidelobes when acquiring data on this hemisphere. Similarly, the analysis of all disk-integrated radar albedos acquired during the Cassini mission did not show an unambiguous leading/trailing dichotomy on Rhea, due to the contribution to the signal of the Inktomi ejecta blanket (Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF].

Modeling and derivation of Rhea's subsurface thermophysical properties

In order to extract new constraints on the chemical and thermophysical properties of Rhea's surface, observations are compared to a model simulating the microwave thermal emission from an icy surface. The model closely follows the one presented in Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF], which we summarize here while detailing some modifications. As in Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF], our modeling approach is divided into three steps: a thermal model providing the vertical temperature profile within the subsurface 𝑇(𝑧), a radiative transfer model calculating the effective temperature 𝑇 𝑒𝑓𝑓 sensed by the Cassini radiometer (which can be very different from the surface temperature if the probed depth is large), and an emissivity model. In the microwave domain, the Rayleigh-Jeans approximation applies and the surface brightness temperature 𝑇 𝑏 is simply the product of the emissivity 𝑒 and the effective temperature 𝑇 𝑒𝑓𝑓 . Convolution of 𝑇 𝑏 with the beam pattern yields the simulated antenna temperature 𝑇 𝑎 , which can be directly compared to the Cassini radiometer observations.

Thermal model

The thermal model considers heat transfer only by conduction, not by radiation for simplicity. The temperature profile is derived at a given latitude/longitude point by solving the 1D heat equation down to six seasonal thermal skin depths (𝛿 𝑡ℎ 𝑠𝑒𝑎𝑠𝑜𝑛 , see equation ( 6) below), at which there are no more expected temperature variations with depth. The boundary condition at the surface imposes radiative equilibrium, taking into account the variations of the incident flux from the Sun along Rhea's orbit. Solar incident flux is computed using SPICE/NAIF libraries [START_REF] Acton | Ancillary data services of NASA's Navigation and Ancillary Information Facility[END_REF] for each time step, including solar eclipses by Saturn. Rhea is far enough from Saturn that the solar flux reflected off Saturn and the infrared thermal emission from Saturn are negligible as compared to the direct solar illumination. The boundary condition at the bottom of the computational layer assumes no internal source of heating and therefore a constant temperature at a depth of 6𝛿 𝑡ℎ 𝑠𝑒𝑎𝑠𝑜𝑛 and below (i.e., 3 -276 m for 𝐼 = 10 -1000 MKS). The number of spatial layers is estimated to ensure computational stability. The thermal model was run for every 10° in latitude and longitude and with a time step of 1.08 hoursi.e., 238113 steps over a Saturn year (29.46 Earth years) and 100 points over a Rhea day (4.518 Earth days). The details of the thermal model are given in [START_REF] Leyrat | Thermal properties of Saturn's rings: from CAMIRAS to the CASSINI mission[END_REF], [START_REF] Leyrat | Modelling the thermal emission from airless planetary surfaces and sub-surfaces[END_REF], and Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF], and its parameters are reported in Table 3. Unfortunately, there is no bolometric Bond albedo map for Rhea: instead, we assume four different values of albedo (constant over the disk), inspired from literature (Table 3). Spatial variations of the visible albedo, although present on optical images [START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF], cannot be explicitly taken into account. Models were run for seven values of thermal inertia 𝐼 from 10 to 1000 MKS, and results were interpolated for intermediate thermal inertia values. Note that the thermal inertia that enters the model does not in general represent the thermal inertia of the surface/subsurface bulk material. As we will find, best fit thermal inertia are much smaller than expected values for compact ices, which we interpret as due to a large porosity of the medium. The inferred values are therefore "effective" thermal inertia. The diurnal and seasonal thermal skin depths are related to the effective thermal inertia 𝐼 as follows:

𝛿 𝑡ℎ = 𝐼 𝜌𝑐 √ 𝑃 𝜋 (6)
where 𝑃 is the thermal wave period (a Rhea day for 𝛿 𝑡ℎ 𝑑𝑎𝑦 or year for 𝛿 𝑡ℎ 𝑠𝑒𝑎𝑠𝑜𝑛 ), 𝜌 is the bulk density of the medium (992 kg/m 3 for water ice), and 𝑐 is the bulk heat capacity. In the case of Rhea, 𝛿 𝑡ℎ 𝑠𝑒𝑎𝑠𝑜𝑛 = 49𝛿 𝑡ℎ 𝑑𝑎𝑦 . In the thermal model, we assume that variations in effective thermal conductivity 𝐾 𝐸 , and therefore in effective thermal inertia (𝐼 = √𝐾 𝐸 𝜌𝑐), are caused by the porosity 𝑝: 𝐾 𝐸 = (1 -𝑝)𝐾 𝐵 . The bulk heat capacity 𝑐 = 7.49 × 𝑇 0 + 90 Jkg -1 K -1 and the bulk thermal conductivity 𝐾 𝐵 = 567/𝑇 0 Wm -1 K are calculated for a temperature of 𝑇 0 = 75 K using the formulas for crystalline water ice presented in [START_REF] Klinger | Some consequences of a phase transition of water ice on the heat balance of comet nuclei[END_REF]. Rhea surface temperatures are actually between 25 and 100 K (based on our thermal models and [START_REF] Howett | Thermophysical property variations across Dione and Rhea[END_REF][START_REF] Howett | Thermal properties of Rhea's poles: Evidence for a meter-deep unconsolidated subsurface layer[END_REF]; assuming a temperature smaller than 75 K would yield larger thermal skin depths roughly inversely proportional to 𝑇 0 .

The output of the thermal model is the vertical temperature profile 𝑇(𝑧) for different assumed thermal inertias and albedos, over the whole surface. The model assumes homogeneous properties with depth: the parameters should therefore be regarded as an average over the sensed depth (that is, the electrical skin depth 𝛿 𝑒𝑙 at the instrument wavelength, see Table 3 and Section 3.2). In particular, the thermal inertia is expected to increase with depth as the medium becomes more compact.

Radiative transfer model

Radiative transfer models have been proposed to simulate the effects of scattering within both atmospheres (e.g., [START_REF] Kokhanovsky | Benchmark results in vector atmospheric radiative transfer[END_REF][START_REF] Ilyushin | Influence of a spatial structure of precipitates on polarization characteristics of the outgoing microwave radiation of the atmosphere[END_REF] and solid surfaces (e.g., [START_REF] Kuga | Millimeter-wave radar scattering from snow 1. Radiative transfer model[END_REF][START_REF] Wiesmann | Microwave Emission Model of Layered Snowpacks[END_REF][START_REF] Mishchenko | Bidirectional reflectance of flat, optically thick particulate layers: an efficient radiative transfer solution and applications to snow and soil surfaces[END_REF][START_REF] Moroz | Improvement of Mishchenko's T-matrix code for absorbing particles[END_REF][START_REF] Liang | The Effects of Layers in Dry Snow on Its Passive Microwave Emissions Using Dense Media Radiative Transfer Theory Based on the Quasicrystalline Approximation (QCA/DMRT)[END_REF][START_REF] Ulaby | Microwave radar and radiometric remote sensing[END_REF] on the outgoing radiation. Nonetheless, in order to limit both the complexity of the model and the number of parameters we have made the choice of a scatter-free radiative transfer model, like Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF] (see also Keihm et al., 2013 for asteroids). The degree of scattering and its regional variations are estimated afterward from the derived emissivity. The backscatter measured during the resolved scatterometry (Section 2.3) also provides an independent insight into the degree of scattering in the subsurface.

For a scatter-free uniform medium of varying temperature with depth, the radiative transfer equation has the following solution [START_REF] Ulaby | Microwave radar and radiometric remote sensing[END_REF]:

𝑇 𝑒𝑓𝑓 (𝜃) = ∫ 𝑇(𝑧)𝑒 -𝑧/(𝛿 𝑒𝑙 cos 𝜃 𝑡 ) 𝑑𝑧 +∞ 0 ∫ 𝑒 -𝑧/(𝛿 𝑒𝑙 cos 𝜃 𝑡 ) 𝑑𝑧 +∞ 0 (7)
The effective temperature sensed by the Cassini radiometer is thus a weighted function of the physical temperature profile, along the line of sight within the subsurface. The angle of transmission 𝜃 𝑡 in the subsurface is defined as follows:

𝑠𝑒𝑐𝜃 𝑡 = 1 𝑐𝑜𝑠𝜃 𝑡 = 1 √ 1 - sin 2 𝜃 𝜀′ (8)
where 𝜃 is the angle of emission at which the Cassini radiometer observes the surface locally, and 𝜀′ is the effective relative dielectric constant, which is the real part of the complex permittivity of the subsurface normalized by the permittivity of vacuum. Its value in the radiative transfer model is fixed at 𝜀′ = 1.15 (the value that will be determined from the emissivity model; see Section 3.3) because it only slightly influences the effective temperature, and only near the limbs. Even changing it from 𝜀′ = 1.15 to 𝜀′ = 2 would not affect the final results. The dielectric constant plays a much more important role in the emissivity, as detailed in Section 3.3.1. The electrical skin depth 𝛿 𝑒𝑙 expresses the characteristic depth down to which the radiometer senses the subsurface temperature. As explained in Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF], the signal sensed by the radiometer is in fact a function of the ratio between the electrical skin depth and the thermal skin depth, and we therefore consider the parameter 𝑟 = 𝛿 𝑒𝑙 𝛿 𝑡ℎ 𝑑𝑎𝑦 . A low ratio 𝑟 (𝑟 ≤ 1) indicates that the radiometer only probes the diurnal variations in temperature of the subsurface, whereas a high ratio 𝑟 (𝑟 ≫ 1) implies that the radiometer also captures seasonal variations in temperature. In other words, a high ratio 𝑟 is indicative of a medium very transparent to microwaves, likely porous and with very little contamination by non-ice impurities, whereas a low value of 𝑟 points to the presence of absorbing materials. The expression of the electrical skin depth 𝛿 𝑒𝑙 when assuming only losses by absorption in the subsurface is given by equation ( 13).

Effective dielectric constant 𝜀′ (normalized by that of vacuum)

1 to 5
The effective dielectric constant is the real part of the effective permittivity. For solid (zero-porosity) water ice at 2.2 cm, 𝜀 ′ = 3.13 [START_REF] Paillou | Microwave dielectric constant of Titan-relevant materials[END_REF].

Emissivity-backscatter model

As the last step, we use the emissivity model initially proposed by [START_REF] Janssen | Anomalous radar backscatter from Titan's surface?[END_REF] for Titan and since then applied on Iapetus and Enceladus (Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF]2017). This model relies on Kirchhoff's law of thermal radiation (emissivity =1 -reflectivity) using the formulation proposed by [START_REF] Peake | Interaction of electromagnetic waves with some natural surfaces[END_REF], which relates the emissivity to radar backscatter. The benefit of such a model is to exploit the 𝜎 0 measured during RH127 active scatterometry data (see Section 2.3) to provide constraints on the emissivity level and regional variations.

The emissivity modeled has two components: a quasi-specular component, which accounts for polarization and limb-darkening, and a diffuse component, which takes into account multiple scattering in the subsurface as well as the reflectivity/emissivity difference between Inktomi and the rest of Rhea:

𝑒(𝜃, 𝜑) = 1 -𝑅(𝜀′, 𝜃, 𝜑) ⏟ 𝑄𝑢𝑎𝑠𝑖-𝑠𝑝𝑒𝑐𝑢𝑙𝑎𝑟 𝑐𝑜𝑚𝑝𝑜𝑛𝑒𝑛𝑡 -( 1 + 𝜇 𝐿 2𝑛𝑓 𝐶𝐵𝐸 ) 𝜎 0 (𝜃, 𝜑) cos 𝑛 𝜃 ⏟ 𝐷𝑖𝑓𝑓𝑢𝑠𝑒 𝑐𝑜𝑚𝑝𝑜𝑛𝑒𝑛𝑡 (10)

Diffuse component:

Here, 𝑛 is found by fitting the function 𝜎 0 (𝜃 𝑒𝑓𝑓 ) = 𝐴𝑐𝑜𝑠 𝑛 𝜃 𝑒𝑓𝑓 through the scatterometry data while ignoring all data <30° in latitude or longitude away from Inktomi, where 𝜃 𝑒𝑓𝑓 is the effective incidence angle of each observation burst (see Section 2.3). Because the RH127 scatterometry data are not spatially complete (they only cover the anti-Saturn side of Rhea), we create a simplified map 𝜎 𝑚𝑎𝑝 0 (𝜃, 𝜑) = 𝜎 0 (𝜃,𝜑) 2𝑛𝑐𝑜𝑠 𝑛 𝜃 𝑒𝑓𝑓 by averaging all values far (>30° in latitude/longitude away) from Inktomi crater and assigning the resulting value to the entire disk, except for the region near Inktomi where a 2-dimensional Gaussian fit is applied to the data. The diffuse component of the model considers the coherent backscattering effect (CBE) which can boost the radar return in the backscattering direction by a factor of up to 2 hence: 1 ≤ 𝑓 𝐶𝐵𝐸 ≤ 2. Invoking such an effect proved necessary to explain the centimetric radar albedos of Galilean moons and, in particular, of the brightest of them, Europa (Black et al., 2001). Rhea is slightly more radar-bright than Europa (Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF]. 𝜇 𝐿 is the linear polarization ratio i.e., the ratio of the backscattering cross-section in orthogonal polarization to the backscattering cross-section in same polarization when transmission is linearly polarized. By definition 0 ≤ 𝜇 𝐿 ≤ 1. 𝜇 𝐿 is null for a smooth homogeneous surface while it equals 1 when the surface is very rough or when multiple/volume scattering occurs in the subsurface. Here 𝜇 𝐿 is assumed to be independent of the viewing geometry (𝜃, 𝜑). The parameter 𝑓 = 1+𝜇 𝐿 𝑓 𝐶𝐵𝐸 is derived by finding the value which yields the best fit for each combination of thermal inertia 𝐼, skin depth ratio 𝑟, and effective dielectric constant 𝜀′. Small values of 𝑓 indicate important scattering in the subsurface volume. Given the boundaries of 𝜇 𝐿 and 𝑓 𝐶𝐵𝐸 , 𝑓 should be bounded between 0.5 (𝑓 𝐶𝐵𝐸 = 2 and 𝜇 𝐿 = 0) and 2 (𝑓 𝐶𝐵𝐸 = 1 and 𝜇 𝐿 = 1). However, we always find 𝑓 < 0.5, which implies that the CBE is not sufficient to explain the high recorded backscattering cross sections; this is further discussed in Section 5.

Quasi-specular component:

Although Rhea's near subsurface is primarily a diffusely scattering medium like that of most of Saturn's icy satellites (e.g. [START_REF] Ostro | Cassini RADAR observations of Enceladus, Tethys, Dione, Rhea, Iapetus, Hyperion, and Phoebe[END_REF][START_REF] Ostro | New Cassini RADAR results for Saturn's icy satellites[END_REF], the emissivity model also includes a quasi-specular component: 1 -𝑅(𝜀′, 𝜃, 𝜑) where 𝑅(𝜀′, 𝜃, 𝜑) is the reflectivity of a smooth surface as given by the Fresnel equations. 𝑅(𝜀′, 𝜃, 𝜑) depends on the dielectric constant 𝜀′ and on the incidence and polarization angles (e.g., Ulaby, 1982). In particular, the quasi-specular component of the emissivity is expected to reach a peak at the Brewster angle 𝜃 𝐵 (see Fig. 4) which relates to 𝜀′ as follows: 𝜃 𝐵 = tan -1 (√𝜀′). 𝜀′ can therefore be estimated using observations acquired over all incidences and azimuthal angles.

Summary of the model and convolution with the beam

At this point, we produce a simulated brightness temperature map 𝑇 𝑏 (𝐼, 𝑟, 𝜀′, 𝑓, 𝜃, 𝜑) = 𝑒(𝜀′, 𝑓, 𝜃, 𝜑) × 𝑇 𝑒𝑓𝑓 (𝐼, 𝑟, 𝜃, 𝜑), which depends on the observation geometry (𝜃, 𝜑), two thermal/electrical parameters (𝐼, 𝑟), and two emissivity parameters (𝜀′, 𝑓). An example of simulated emissivity, effective temperature and brightness temperature maps is shown in Fig. 4.

The next step consists in convolving the brightness temperature map with the beam pattern in order to obtain simulated antenna temperatures. To limit the computation time, we separate the diffuse and quasi-specular components of the emissivity before convolution, yielding the following equation for the antenna temperature 𝑇 𝑎,0 𝑚𝑜𝑑𝑒𝑙 :

𝑇 𝑎,0 𝑚𝑜𝑑𝑒𝑙 (𝐼, 𝑟, 𝜀 ′ , 𝑓) = ∬ (1 -𝑅(𝜀′))𝑇 𝑒𝑓𝑓 (𝐼, 𝑟)𝐺𝑑𝛺 4𝜋 -𝑓 ∬ 𝜎 𝑚𝑎𝑝 0 𝑇 𝑒𝑓𝑓 (𝐼, 𝑟)𝐺𝑑𝛺 4𝜋 (11)
where 𝑑𝛺 = sin(𝜃) 𝑑𝜃𝑑𝜑 is the solid angle in the viewing direction (𝜃, 𝜑); all variables, except 𝑓, are functions of the viewing direction. We apply this convolution with the beam pattern for the geometry of each time step, thus obtaining a time sequence which can be directly compared to the calibrated antenna temperatures measured by the Cassini radiometer (e.g., Fig. 1). Before being able to compare the modeled antenna temperature to the observations, there is one last required step: the baseline must be subtracted from the simulated data. Indeed, during the distant scans, Rhea is almost always present in the near sidelobes during the cold sky observations, and the modeled temperature thus does not reach zero. Since this contribution was removed from the observed 𝑇 𝑎 𝑜𝑏𝑠 when subtracting the radiometer temperature baseline (see Section 2.2.2), the exact same method must be applied to the modeled 𝑇 𝑎,0 𝑚𝑜𝑑𝑒𝑙 in order to obtain final modeled antenna temperature time sequence 𝑇 𝑎 𝑚𝑜𝑑𝑒𝑙 .

Fitting method

The free parameters of the model are the thermal inertia 𝐼, the ratio of diurnal thermal skin depth to electrical skin depth 𝑟, and the effective dielectric constant 𝜀′. Although Eq. ( 11) may suggest that 𝑓 could be derived directly by comparison to the data for each combination of (𝐼, 𝑟, 𝜀′), this is in fact not possible because the subtraction of the baseline is also slightly dependent of 𝑓. Instead, for each combination of the parameters 𝐼, 𝑟, and 𝜀′ and for 15 values of 𝑓 from 0.1 to 0.8, we compute 𝑇 𝑎 𝑚𝑜𝑑𝑒𝑙 and the reduced weighted chi-squared statistic 𝜒 𝑟 2 , as follows:

𝜒 𝑟 2 = ∑ 𝑤 × ( 𝑇 𝑎 𝑜𝑏𝑠 -𝑇 𝑎 𝑚𝑜𝑑𝑒𝑙 𝜎 ) 2 𝑡 𝑓 𝑡 0 ∑ 𝑤 𝑡 𝑓 𝑡 0 (12)
where 𝜎 is the noise defined in Section 2.2 (Eq. ( 2)), 𝑡 0 and 𝑡 𝑓 are the start and end time of the observation, and 𝑤 is a weight calculated according to the angular distance to Inktomi and the distance between consecutive observations (described below). Plotting 𝜒 𝑟 2 versus 𝑓 yields a parabolic curve, whose minimum indicates the best value for 𝑓.

The rationale behind the weights and their calculation are described in the following. First, when we include the stares or the resolved scan in the fit, the observations are not uniformly distributed in time and space. Indeed, the stares observe a unique point at Rhea's surface many times during tens of minutes. These consecutive observations of the same region are thus not independent, and should not be given more weight than the observation of the same region during a much faster but independent distant scan. This also poses a problem for the resolved dataset (RH127_4r), which has repetitive coverage near the limbs (and not near nadir). In order to take into account the varying resolution and spatial coverage, each observation is weighted by the angular distance ∆𝑑 separating it from the previous observation.

Second, looking at the RH127 scatterometry data (Fig. 3) and the concurrent deconvolved brightness temperature map (RH127_4r in Fig. 2), it is obvious that areas near or within the Inktomi ejecta blanket region (IEBR) are more radar-bright and radiometrically colder than their surroundings, strongly suggesting different thermal properties and/or emissivity in this region. The thermal emission model we have developed only accounts for differences in emissivity (using the scatterometry observations, see Section 3.3), while the local thermal inertia, ratio, and albedo may also be different. We further note that the bolometric Bond albedo of the IEBR is not straightforward to compute from ISS and VIMS data at different wavelengths and phase angles, and such a derivation is beyond the scope of the present paper. The brightness of the IEBR varies considerably with wavelength [START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF], and even its appearance in 3-color composite maps changes with the color ratios chosen. Indeed, the IEBR is much brighter than its surroundings in the map of [START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF], but is indistinguishable in the updated map of Schenk et al. (2018) (PIA18438). From CIRS thermal infrared data, [START_REF] Howett | Thermophysical property variations across Dione and Rhea[END_REF] do not find a clear albedo change near Inktomi either. Furthermore, because the IEBR is a fairly small area, there is not sufficient radiometry coverage in local hours or seasons to constrain its properties independently from those of its surroundings.

Instead, we treated the IEBR in two different ways: i) using the backscattering crosssection map 𝜎 𝑚𝑎𝑝 0 described in Section 3.3 (see also Fig. 4), or ii) fitting only the data with no Inktomi contribution and weighting the rest of the data to account for the contribution of the IEBR in the near sidelobes. To calculate the percent contribution of the IEBR to the antenna temperature, we first delineated this region using the enhanced three color global map of Rhea from [START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF]. We then convolved this map with the beam pattern, yielding the percent filling factor 𝜈 𝐼𝐸𝐵𝑅 of Inktomi for each particular observation. Data "with no Inktomi contribution" are defined as those where the Inktomi areal filling factor is 𝜈 𝐼𝐸𝐵𝑅 < 0.02 (2%). Because even these data are not entirely devoid of Inktomi contribution, we weighed them with the factor (0.02 -𝜈 𝐼𝐸𝐵𝑅 ). The final weights applied when calculating 𝜒 𝑟 2 are therefore 𝑤 = ∆𝑑 × (0.02 -𝜈 𝐼𝐸𝐵𝑅 ).

The best fitting value for each parameter (𝐼, 𝑟, 𝜀′) is the one giving the smallest 𝜒 𝑟 2 (𝜒 𝑟,𝑚𝑖𝑛 2 ); they are given in Table 4. The 2-σ error bars on the parameters are found using the 95% confidence intervals from the 𝜒 2 statistics, taking into account the number of DOF (degrees-of-freedom). We ran all models for 4 different values of the bolometric Bond albedo 𝐴 𝐵 from [START_REF] Howett | Thermal inertia and bolometric Bond albedo values for Mimas, Enceladus, Tethys, Dione, Rhea and Iapetus as derived from Cassini/CIRS measurements[END_REF][START_REF] Howett | Thermal properties of Rhea's poles: Evidence for a meter-deep unconsolidated subsurface layer[END_REF] and [START_REF] Pitman | Disk-integrated bolometric Bond albedos and rotational light curves of saturnian satellites from Cassini Visual and Infrared Mapping Spectrometer[END_REF] (see Table 3).

Results

As a first approach to the data, we calculate disk-integrated emissivities for each distant scan of Rhea, using the thermal and radiative transfer but not the emissivity-backscatter model. This preliminary analysis yields some constraints on the probed depths of the South pole and on emissivity variations over Rhea. However, it is very limited: both an emissivity model and a comprehensive analysis of all overlapping observations is necessary. The results of applying the complete model and fitting method presented in Section 3 are then described, providing new constraints on the thermal and dielectric properties of the subsurface in different regions of Rhea.

Preliminary analysis of disk-integrated emissivities

Before searching for the best-fit thermal and dielectric parameters, we conduct a preliminary analysis of the disk-integrated emissivity for each distant scan ("u" observations in Table 1). For each combination of 𝐼 and 𝑟, we average the effective temperature map, obtained from Eq. ( 7), over the disk, obtaining 𝑇 𝑒𝑓𝑓 𝑑𝑖𝑠𝑘 (𝐼, 𝑟). In parallel, we determine a disk-integrated brightness temperature 𝑇 𝑏 𝑑𝑖𝑠𝑘 by modeling Rhea as a limb-darkened disk of uniform temperature, and finding the temperature that best fits the antenna temperatures 𝑇 𝑎 , as described in Section 2.2.1. The disk-integrated emissivity for each set of parameters is then given by 𝑒 𝑑𝑖𝑠𝑘 (𝐼, 𝑟) = 𝑇 𝑏 𝑑𝑖𝑠𝑘 𝑇 𝑒𝑓𝑓 𝑑𝑖𝑠𝑘 ⁄ (𝐼, 𝑟). This method is the same as that used by Le Gall et al. ( 2014) on Iapetus. Fig. 5 displays the disk-integrated emissivity for each combination of thermal inertia 𝐼 and skin depth ratio 𝑟, for each distant scan. A uniform bolometric Bond albedo of 0.55 was assumed; decreasing the albedo would correspondingly decrease the emissivities. Note that, while we explore a wide parameter space (𝑟 from 0.1 to 1000 and 𝐼 from 10 to 1000 MKS), electric skin depths below 10 cm (<5 times the wavelength) are very unlikely in a transparent medium such as a water-ice regolith. The thermal inertia / ratio combinations corresponding to such short skin depths, which generally yield very high emissivities, should therefore be discarded. For the same parameters (𝐼,𝑟), we find higher emissivities for flybys RH018 and RH049, which are on the sub-Saturn side, than on the rest of Rhea. This suggests the presence of a high-emissivity or low-albedo anomaly in this region. For most of Rhea, this simple method cannot further constrain the thermal or dielectric parameters. In the South pole, however, data were collected during the Southern summer (RH011 in July 2005) and during the southern fall (RH177 in December 2012), both centered at a latitude of 77°S. This temporal sampling offers the opportunity for especially reliable parameter inversion in the South polar regions (see Fig. 10 and Section 5.2).

The disk-integrated emissivities of individual flybys provide initial but limited insights into the properties of the regolith in the South polar region. In order to further constrain the thermal and dielectric parameters of Rhea's near subsurface, the emissivity-backscatter model and fitting method described in Section 3.3 -3.5 are required.

Fitting results

To investigate the thermal and electrical properties of different regions of interest of Rhea, we apply the fitting method described above to the time-sequence data observed in three different regions: the South pole (defined as latitudes poleward of 50°S), the sub-Saturn side (0±45°E; <50° latitude), and the anti-Saturn side (180±45°E; <50° latitude). There is too little coverage on the North pole and the trailing hemisphere for their analysis, and the leading hemisphere is contaminated by the Inktomi crater ejecta blanket, which clearly has very specific chemical and/or thermophysical properties. We also extract best-fit parameters for all the equatorial (<50° latitude) data, and for all available Rhea data. Regional differences in albedo and emissivity are taken into account together as described in Section 4.4. The IEBR is either masked or modeled as an emissivity anomaly (see Section 3.4); the results obtained in both cases are reported in Table 4 when Inktomi is within the region of interest (i.e., not for the sub-Saturn side). Fig. 6 illustrates the 2-σ confidence regions in the (𝐼, 𝑟) space, for best fit values of 𝜀′ and 𝑓. ) along with their 2-σ error bars as obtained assuming a bolometric Bond albedo of 𝐴 𝐵 = 0.55. These parameters are then converted into diurnal, seasonal (Eq. 6), and electrical skin depths, assuming a temperature of 𝑇 0 = 75 K when calculating the specific heat capacity. Using a lower 𝑇 0 would yield higher thermal and electric skin depths and a lower loss tangent. The emissivity 𝑒 is derived from the dielectric constant 𝜀′ and the factor 𝑓 (Eq. 10). For a low-loss medium like porous water ice and assuming all losses are due to absorption, the loss tangent is derived from the wavelength 𝜆 (2.2 cm), the electric skin depth 𝛿 𝑒𝑙 , and the dielectric constant 𝜀′ with the following equation:

Figure 5: Disk-integrated emissivities for different combinations of thermal inertia and skin depth ratio obtained from the 10 distant scans of Rhea (see Table 1). A uniform bolometric

tan 𝛿 = 𝜀 ′′ 𝜀 ′ = 𝜆 2𝜋√𝜀′𝛿 𝑒𝑙 (13)
The same equation allows us to find the imaginary component of the permittivity, 𝜀′′. A low value of 𝜀′′ (and of tan 𝛿) corresponds to a low-loss, transparent medium, through which the signal can travel easily: it is therefore associated with large electric skin depths 𝛿 𝑒𝑙 . However, the low dielectric constants and high radar brightness suggest important losses due to scattering (in addition to absorption). Thus the derived values for the loss tangent are strict upper bounds.

For the emissivity 𝑒, we find values in the range 0.62 -0.76 assuming 𝐴 𝐵 = 0.55. This range is extended to 0.59 -0.84 if the albedo varies from 0.42 to 0.72. This is largely consistent with the value of 0.59 derived by [START_REF] Ostro | Cassini RADAR observations of Enceladus, Tethys, Dione, Rhea, Iapetus, Hyperion, and Phoebe[END_REF] using fewer observations, a 0.45 Bond albedo, and a very simple temperature model. For the ranges of (𝐼,𝑟) shown in Fig. 6 and for the flybys in each considered regions, the disk-integrated emissivity (shown in Fig. 5) is within error of the emissivity given in Table 4. For example, on the South pole, the diskintegrated emissivity within the 2-σ best fitting ellipses in the (𝐼, 𝑟) domain extends from 0.62 to 0.66, within error of the value shown in Table 4 (0.65 -0.03 +0.03 ), both for a bolometric Bond albedo of 0.55. Slightly smaller disk-integrated emissivities (by ~0.02) are generally obtained in regions with Inktomi and its ejecta blanket, because this region is included in the diskintegrated data.

We further note that these results are obtained while neglecting surface roughness. However, we have checked that this assumption does not change the results by also deriving the quasi-specular component of the emissivity according to the model proposed by [START_REF] White | Lunar polarization studies at 3.1 mm wavelength[END_REF], which includes large-scale surface roughness as a new parameter. For a representative sub-section of the data (only one scan for each flyby, thus ignoring repeat scans and the resolved observation), we always obtain the same best-fitting thermal inertias 𝐼 and skin depth ratios 𝑟 as with the smooth model (as given in Table 4). For the South pole and the anti-Saturn side the derived dielectric constant remains the same and the best-fitting roughness is minimal (1° rms slope). On the anti-Saturn side, the highest tested value of roughness (35° rms slope) yields the best fit, with a correspondingly higher dielectric constant (1.75 instead of 1.3-1.4). The large number of parameters of this model means error bars should be very high and interpretations on surface roughness would be speculative. Nonetheless, the fact that we always find the same best 𝐼 and 𝑟 while the dielectric constant 𝜀 ′ remains < 2 provides confidence in the robustness of the results obtained assuming a smooth surface.

Effect of the bolometric Bond albedo

Modifying the albedo uniformly across Rhea mostly affects the absolute values of the brightness temperature, rather than its spatial or temporal variations. A higher albedo leads to colder modeled physical temperatures (at first approximation they are proportional to (1 -𝐴 𝐵 ) 1 4 ) and therefore the emissivity required to fit the model must be higher (Fig. 7), which further implies a higher scattering factor 𝑓. Thus 𝑓 and the bolometric Bond albedo 𝐴 𝐵 are essentially degenerate parameters, while the remaining parameters change little with albedo. This is true only to first order: we note a small difference in the derived thermal inertias, and consequently in the electric skin depth and loss tangent, as a function of 𝐴 𝐵 (see Fig. 7 for the case of the South Pole). Including the influence of the global value of the albedo would increase the error bars on the derived parameters by about 10%, but would not change the overall results (except for the scattering factor 𝑓 and emissivity 𝑒) which is why we only show the case where 𝐴 𝐵 = 0.55 in Table 4.

4.4.Regional albedo/emissivity variations

The local value of the albedo, however, has a significant influence on the observed temperature map. Indeed, applying our fitting methods to the data from all flybys at once yields imperfect fits (𝜒 𝑟 2 > 2.8), with the residual clearly showing, whether Inktomi is masked or not, an anomaly on the sub-Saturn side (longitude of 0°E; Fig. 8b &c). More localized anomalies are also apparent on the leading hemisphere (longitudes of -180 to 0°E). At first order, these poor fits may be explained by albedo and/or emissivity variations at the surface of Rhea. A higher albedo yields lower scattering factors 𝑓 and higher emissivities 𝑒, but barely changes the other derived parameters. The Inktomi ejecta blanket region is handled with two different methods: either it is masked and the fit is applied to the rest of the region ("Masked"), or it is included but treated as an emissivity anomaly, with a Gaussian dependence with distance to Inktomi itself ("Included"); see Section 3.4 for details. Note that Inktomi being located on the anti-Saturn Southern hemisphere, it has no influence on the sub-Saturn data. When the error bars reach the boundaries of the range of tested values for the thermal inertia (10 < 𝐼 < 1000), the actual error bars cannot be known, and MIN or MAX are indicated instead. < 0.9

Inktomi

1 No upper bound can be extracted for tan 𝛿 if no lower bound was found for the thermal inertia. 4 for the other parameters) for each of the 5 regions where the fitting method was applied, and for each of the two different ways to treat the Inktomi ejecta blanket (masked or not). The minimum reduced 𝜒 2 value (𝜒 𝑟,𝑚𝑖𝑛 2) is indicated for each fit. The bolometric Bond albedo used is 𝐴 𝐵 = 0.55. Note that, while the thermal inertia and skin depth ratio vary, the best-fit electric skin depth 𝛿 𝑒𝑙 always remains within 5-15 m. Interestingly, these variations can be captured by taking advantage of the model degeneracy between the subsurface scattering parameter 𝑓 and the bolometric Bond albedo 𝐴 𝐵 . In this approach, Rhea's surface is divided into 45° latitude and longitude bins, a size large enough to contain enough data to constrain 𝑓 and small enough to take into account spatial variations and in particular the IEBR and the leading/trailing dichotomy. In each of these areas and for each set of parameters (𝐼, 𝑟, 𝜀′), we find the best-fit value of 𝑓 following the method described in Section 3.5. We then combine the data over all of Rhea and calculate a global 𝜒 𝑟 2 . The parameters (𝐼, 𝑟, 𝜀′) are assumed uniform over all of Rhea, while the subsurface scattering and albedo are allowed to vary spatially (but cannot be disentangled). With this method, much better fits are obtained; the residuals 𝜒 𝑟,𝑚𝑖𝑛 2 are indicated in Fig. 8e & f for both different ways of treating the IEBR (masked or included and modeled using the scatterometry data), with and without dividing 𝑓 into smaller regions. Calculating 𝑓 in 45 × 45° regions yields a much more uniform residual, and thus sufficiently takes into account the albedo and/or emissivity variations. This method is applied to all data at once, and to the equatorial data only (i.e, at latitude within ±50°), with and without Inktomi (Table 4); albedo/emissivity variations are less important within smaller regions where this method is not used (South pole and sub-and anti-Saturn sides) .

To examine the influence of the IEBR, Fig. 8 shows the residual of a global fit obtained assuming the emissivity is uniform over Inktomi, without (Fig. 8a) then with (Fig. 8d) dividing 𝑓 into 45° squares. The strong negative residual near Inktomi and the poor quality of the fit testify of the necessity to either i) mask the IEBR (Fig. 8b &e) or ii) model the IEBR as an emissivity anomaly (Fig. 8c &f). Once 𝑓 has been calculated in each 45° region, we can compute the emissivity in each region using Eq. ( 10), assuming a constant albedo; the resulting map is shown in Fig. 9. In practice, emissivity and albedo effects cannot be disentangled using this method. Thus regions with low emissivity in this map are more scattering and/or have a higher albedo than their surroundings (i.e., their brightness temperature is colder, due either to low emissivity and/or high albedo). We find a low-emissivity/high-albedo anomaly spatially correlated with the IEBR, and a high-emissivity/low-albedo anomaly at the northern midlatitudes of the sub-Saturn side (see Fig. 9). Conversely, the partial bolometric Bond albedo map derived from CIRS observations by [START_REF] Howett | Thermophysical property variations across Dione and Rhea[END_REF] shows an increase in albedo from -30° to 60°N. This contradiction may be explained either by a higher emissivity in high-latitude sub-Saturn regions, or by the lower reliability of high northern latitude observations which were acquired at high emission angles. On the other hand, our results are consistent with a higher albedo at the South pole than at the Equator (except near Inktomi), as also found by [START_REF] Howett | Thermal properties of Rhea's poles: Evidence for a meter-deep unconsolidated subsurface layer[END_REF]. Fig. 8c andf show a positive residual to the West of Inktomi: this is because the Inktomi crater ejecta is, in fact, butterfly-shaped rather than circular as was modeled. Lastly, Fig. 8 reveals another cold spot to the North-west of Inktomi (-180--100°E, 45°N), which is associated with a slight increase in the backscatter image (Fig. 3). This cold spot is not correlated to any optical or IR anomaly, but does coincide with a large part of Tirawa crater, the largest crater on Rhea. Further interpretation of this anomaly would require more highresolution data which cannot be obtained in the foreseeable future; it thus remains to be interpreted. Figure 9: Emissivity map assuming a uniform albedo 𝐴 𝐵 = 0.55, using all available data. The emissivity is assumed constant in every 45x45° latitude/longitude bin. There is no data available in the white regions; the emissivities between 45 and 65°N are based on a few limb observations and are less reliable. Because emissivity and albedo cannot be separated, apparent lower emissivity regions observed in this map may instead be caused by higher bolometric Bond albedos (and vice versa for high emissivity regions).

Regional variations in terms of thermal and dielectric properties

Because all of the observations were taken during the local night (sub-spacecraft point between 20:00 and 00:30 local time, see Table 1), the thermal inertia of the equatorial regions can only be approximately constrained. However, at high latitudes, seasonal variations in temperature are much more prevalent. In particular, at the South pole of Rhea, radiometry data were collected during the southern summer and fall, as detailed in Section 4.1. Combining summer and fall observations provides reliable constraints on the thermal properties of the nearsubsurface. Fits for the South pole are especially good (𝜒 𝑟,𝑚𝑖𝑛 2 ~1), and are shown in Fig. 10, where the effects of the different parameters are illustrated separately. Lower values of the thermal inertia yield warmer modeled temperatures in the summer (2005), and colder in fall (2012), as the subsurface retains less heat (Fig. 10f). The same is true for the ratio parameter: at smaller ratios, we only probe the top layers, where the temperature changes more quickly with incoming solar flux, leading to warmer summer temperatures and colder fall temperatures (Fig. 10g). The 𝜒 𝑟 2 maps (Fig. 10a, b, &c) show that the derived thermal inertia (𝐼 = 70 -850 MKS), ratio (𝑟 = 20 -100), and dielectric constant (𝜀 ′ = 1.09 -1.35) are still somewhat correlated, while the best determined quantity is the electric skin depth (4 -10 m).

While constraining thermal properties requires to have a sample of data collected at different times of the day or of the year, the determination of the dielectric constant requires observations at different emission and polarization angles. As expressed in Eq. ( 10) and described in Section 3.3, the dielectric constant mainly intervenes in the quasi-specular component of the emissivity model, described by Fresnel's coefficients. More specifically, the dielectric constant controls both limb darkening and the emissivity variations around the Brewster angle. However, observing either of these characteristics and deriving the effective dielectric constant is possible only with relatively high-resolution observations such as RH127 (sub-Saturn; resolution of 0.1-0.5𝑅 𝑅ℎ𝑒𝑎 ) and RH177 (South pole; resolution of 0.3-0.4𝑅 𝑅ℎ𝑒𝑎 ). As seen in Fig. 10e, only the last two flybys (RH127 and RH177) show changes of the modeled temperature with 𝜀′, and are therefore useful for the derivation of this parameter.

Interpretations and discussion

In this section, we examine the implications of the results of fitting the simulated antenna temperatures to Rhea radiometry observations. We look first at global compositional and structural characteristics derived for Rhea, then at regional variability. Finally, we compare our results to those obtained on Iapetus with the same method (Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF].

Global results

Composition

As a general result, we find that the Cassini radiometer is observing Rhea's subsurface down to depths of 𝛿 𝑒𝑙 = 6 -13 m, with variable error bars. At the South pole, excluding Inktomi and for all considered albedo values, the minimum electrical skin depth is 3.8 m, which is 170 times the radiometer wavelength (2.2 cm). Large penetration depths indicate a low loss tangent (tan 𝛿 < 1.1 × 10 -3 at the South pole; tan 𝛿 < 4.7 × 10 -3 elsewhere) and a correspondingly small imaginary part of the electrical permittivity, pointing to a very transparent, non-absorbing medium. This further implies that the near subsurface of Rhea is composed of very pure and probably also porous water ice. Impurities, such as organics or silicates originating from impacts or dust from the Saturn system, would tend to reduce the transparency of the medium. RH011) are lower even though they are in summer, because the resolution is not as high as for the later data (see Table 1). e) Modeled 𝑇 𝑎 for 𝐼 and 𝑟 constant, and 𝜀′ varying. f) Modeled 𝑇 𝑎 for 𝜀′ and 𝑟 constant, and 𝐼 varying. g) Modeled 𝑇 𝑎 for 𝜀′ and 𝐼 constant, and 𝑟 varying. Note that both the thermal inertia and the skin depth ratio are reasonably well constrained due to the observations at different seasons, whereas the dielectric constant is determined from the higher-resolution datasets (RH127 and RH177).

This interpretation is consistent with the high recorded backscattering coefficients (Fig. 3) and disk-integrated radar albedos (Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF] which are both indicative of ultraclean water ice in the top few meters. The upper crust of Rhea was already suspected to be composed of pure water ice down to a depth of a few km based on the strong H2O spectral signature identified by Cassini VIMS in the trailing hemisphere tectonic scarps and in the walls of recent craters [START_REF] Stephan | The Saturnian satellite Rhea as seen by Cassini VIMS[END_REF]. The trailing hemisphere darkening and the leading hemisphere reddening apparent in visible and near-infrared data [START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF] must then be due to a thin (i.e., at most a few centimeters) layer of contaminants, through which the radar/radiometer can easily see. Indeed, a non-icy layer tens centimeters thick would increase the emissivity and reduce the probed depth, as on the dark side of Iapetus (Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF].

Iapetus, the only other icy satellite of Saturn where the same method was applied, shows a lower radar brightness (𝐴 𝑎𝑣𝑔 = 0.3 -1.0), higher emissivities (𝑒 = 0.74 -0.87), and shallower depths probed by the radiometer (𝛿 𝑒𝑙 = 1.2 -4.1 m) than on Rhea (Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF]2019). These values are explained by the presence of absorbing contaminants, especially on Iapetus's leading side, introduced by the incoming Phoebe ring particles. Meanwhile, at Rhea, fluxes from Saturn's Phoebe ring, E ring, charged particles from the magnetosphere, and other solid impactors are expected to be very small [START_REF] Verbiscer | Saturn's largest ring[END_REF][START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF][START_REF] Howett | Ring and Magnetosphere Interactions with Satellite Surfaces[END_REF][START_REF] Hendrix | Surface Composition of Saturn's Icy Moons[END_REF], in agreement with a thin non-icy coating. An almost pure water ice upper crust suggests a relatively young surface; however, to accurately determine the age of the surface, a model of the micrometeoroid and large impactor flux and composition would be necessary, which is particularly difficult given the multiple sources of impactors (Kirchoff et al., 2018 and references therein). Moreover, these fluxes are expected to be different in different regions of Rhea (polar/equatorial, leading/trailing, sub-Saturn/anti-Saturn), so the well-constrained polar loss tangent is not applicable to all of Rhea.

Structure

For very diffuse scattering surfaces such as those of the icy satellites, the effective dielectric constant is mostly a measure of the degree of depolarization of the incident waves by the regolith. It thus provides insight into the subsurface structure (porosity, heterogeneity) rather than composition (although a pure water ice regolith will favor scattering in the subsurface volume by allowing large penetration depth). This is why very low effective dielectric constants are possible for a depolarized water ice medium, even though the bulk dielectric constant of water ice is 3.13 at cm wavelengths (e.g., Mätzler, 1996;[START_REF] Paillou | Microwave dielectric constant of Titan-relevant materials[END_REF]. The low derived 𝜀 ′ values (𝜀 ′ = 1.07 -1.5, close to that of vacuum) indicate that little to no increase of the emissivity is detected at the Brewster angles: the near-surface mainly consists in an unpolarized, likely porous, regolith.

We recall that diffuse scattering in the subsurface was modeled when calculating the emissivity as 𝑓 × 𝜎 0 (𝜃,𝜑) 2𝑛𝑐𝑜𝑠 𝑛 𝜃 𝑒𝑓𝑓 . Thus when 𝜎 0 is held constant (i.e., everywhere except near Inktomi), the parameter 𝑓 = 1+𝜇 𝐿 𝑓 𝐶𝐵𝐸 represents the amount of scattering. As discussed in section 3.4, 𝑓 should theoretically lie between 0.5 and 2, and more likely be > 1 since the derived very low 𝜀′ points to a highly depolarizing subsurface and therefore to a linear polarization ratio 𝜇 𝐿 close to 1. Instead, we find 𝑓 = 0.21 -0.42 for a bolometric Bond albedo of 𝐴 𝐵 = 0.72, and 𝑓 = 0.39 -0.60 for 𝐴 𝐵 = 0.42. For 𝜇 𝐿 = 1 and 𝐴 𝐵 = 0.55, this would imply a coherent backscattering factor 𝑓 𝐶𝐵𝐸 = 3.6 -6 while it is limited to the theoretical value of 2. This further means that the coherent backscattering effect is not sufficient to explain the radar backscatter recorded on Rhea. This result is consistent with the analysis of the disk-integrated radar albedo measured by the Cassini radar in the Saturnian inner system (Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF]. It suggests that some exotic scattering processes are at play in the subsurface of Rhea (as well as on Saturn's other moons including in some regions on Titan, [START_REF] Janssen | Anomalous radar backscatter from Titan's surface?[END_REF], likely triggered by the presence of structures of centimetric size (i.e., of the order of the wavelength; [START_REF] Hapke | Coherent backscatter and the radar characteristics of outer planet satellites[END_REF] very efficient at returning the radar waves in the backscattering direction.

The thermal inertias values derived from the Cassini microwave radiometry dataset also provide insights into the state of the regolith, especially by comparison to values inferred from measurements in the infrared which are sensitive to shallower depths (a few millimeters against few meters). The thermal inertias we derive on Rhea at 2.2 cm (best fits around 50 -300 MKS, Table 4) are higher than those measured in the infrared (1 -50 MKS; [START_REF] Howett | Thermophysical property variations across Dione and Rhea[END_REF][START_REF] Howett | Thermal properties of Rhea's poles: Evidence for a meter-deep unconsolidated subsurface layer[END_REF], implying compaction of the medium with depth, as also found on Iapetus by Le [START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF]. Yet these thermal inertia values remain low compared to the bulk thermal inertia of water ice (2000 MKS). Low thermal inertias have been measured at infrared and millimeter wavelengths on outer Solar System bodies, including Jupiter's satellites (𝐼 = 50 -70 MKS, [START_REF] Spencer | Temperatures on Europa from Galileo Photopolarimeter-Radiometer: Nighttime Thermal Anomalies[END_REF]), Saturn's icy satellites (𝐼 = 1.5 -70 MKS, [START_REF] Howett | Thermal inertia and bolometric Bond albedo values for Mimas, Enceladus, Tethys, Dione, Rhea and Iapetus as derived from Cassini/CIRS measurements[END_REF], and trans-Neptunian objects (𝐼 = 0.1 -10, [START_REF] Lellouch | TNOs are Cool": A survey of the transNeptunian region -IX. Thermal properties of Kuiper belt objects and Centaurs from combined Herschel and Spitzer observations[END_REF]. These values have been generally interpreted as primarily caused by high porosities in the top few mm-cm of the surface, though [START_REF] Ferrari | Low thermal inertias of icy planetary surfaces -Evidence for amorphous ice?[END_REF] showed a dependence on grain size, grain arrangement, the presence of amorphous water ice, and temperature (if heat transfer by radiation is important). Amorphous water ice has a lower bulk thermal inertia (dependent on temperature), a lower thermal conductivity, and a slightly higher dielectric constant than crystalline water ice [START_REF] Ferrari | Low thermal inertias of icy planetary surfaces -Evidence for amorphous ice?[END_REF]. Based on the heliocentric distance variation of the thermal inertia in outer Solar System objects, [START_REF] Ferrari | Low thermal inertias of icy planetary surfaces -Evidence for amorphous ice?[END_REF] favor the presence of amorphous ice at cm depths, under a thin crystalline coating, both for Mimas and TNOs. VIMS data indicate the coexistence of amorphous and crystalline ice on Rhea's surface [START_REF] Dalle Ore | Impact craters: An ice study on Rhea[END_REF]. Thus our observations are in general agreement with the interpretations of [START_REF] Ferrari | Low thermal inertias of icy planetary surfaces -Evidence for amorphous ice?[END_REF], additionally pointing to an increase of thermal inertia with depth, consistent with decreasing porosity and/or amorphous ice fraction with depth.

South Pole

Because the South pole was observed at different seasons, the thermal inertia 𝐼 and skin depth ratio 𝑟 are well determined, while also providing an excellent fit (𝜒 𝑟,𝑚𝑖𝑛 2 ~ 1). For a bolometric Bond albedo of 𝐴 𝐵 = 0.72 [START_REF] Howett | Thermal properties of Rhea's poles: Evidence for a meter-deep unconsolidated subsurface layer[END_REF], when Inktomi is masked, we find a thermal inertia of 𝐼 = 209 -140 +662 MKS, an electric skin depth of 𝛿 𝑒𝑙 = 7.9 -4.1 +4.0 m, and a dielectric constant of 𝜀′ = 1.2 -0.11 +0.15 . These values can be compared to the ones found by [START_REF] Howett | Thermal properties of Rhea's poles: Evidence for a meter-deep unconsolidated subsurface layer[END_REF]: assuming that both poles share the same thermal properties and combining CIRS northern winter and southern summer observations, they inferred a bolometric Bond albedo 𝐴 𝐵 = 0.70 -0.74 and a thermal inertia 𝐼 = 1 -46 MKS. They interpreted these low thermal inertias as indicative of a meter-deep unconsolidated regolith at Rhea's poles. The larger thermal inertias derived from Cassini radiometry imply a gradient in thermal inertia with depth: while the upper tens of cm of the polar surface may be a very fluffy, snow-like medium, at meter depths the subsurface is more compacted. Specifically, the thermal inertia may increase from 𝐼~25 MKS in the first meter to 𝐼~210 at ~5 -10 m depth.

As discussed in Section 5.1, 𝜀′, the real part of the permittivity, depends primarily on the structure of the subsurface, whereas the imaginary part 𝜀′′ is very sensitive to the amount of non-icy contaminants and can thus offer constraints on the composition. To deduce the allowed volume fraction of contaminants, we assume a ternary mixture composed of water ice, vacuum, and compacted tholins, which is a possible relatively low-loss component of the optically dark material seen on Rhea's leading and trailing hemisphere [START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF]. We use the bulk permittivity values measured by [START_REF] Paillou | Microwave dielectric constant of Titan-relevant materials[END_REF]: 𝜀 𝑡ℎ𝑜𝑙𝑖𝑛 = 2.33 + 𝑗20.6 × 10 -3 and 𝜀 𝑤𝑎𝑡𝑒𝑟 𝑖𝑐𝑒 = 3.13 + 𝑗1.3 × 10 -3 . The Hashin-Shtrikman bounds, derived from the Maxwell-Garnett mixing laws, allow us to constrain the volume fraction of each constituent [START_REF] Hashin | A Variational Approach to the Theory of the Effective Magnetic Permeability of Multiphase Materials[END_REF][START_REF] Sihvola | Mixing Rules with Complex Dielectric Coefficients[END_REF]. Since we assumed only losses through absorption when deriving the loss tangent and it is known that scattering is important on Rhea (Section 5.1), only a maximal volume fraction of tholins can be obtained, with the actual value likely being significantly smaller.

The compositional constraints inferred from the imaginary part of the permittivity at the South pole are shown in the ternary diagram of Fig. 11. The low loss tangent indicates that the average porosity over the probed depths is necessarily > 10%, which remains consistent with the thermal inertia variations described above. We also find that tholin compounds can only be present in small quantities (< 10%) at the South pole. While ferrous oxides or silicates are also candidate components for this dark material, their loss tangent is at least two orders of magnitude higher than that of tholins, implying that they can only be present in very small quantities (i.e., a thin layer at the surface).

Equatorial vs polar regions

The interpretation of the radiometry dataset collected in Rhea's equatorial regions suffers from several limitations: i) all data were collected during the local night, ii) while some data were acquired 5 years apart, seasonal temperature variations are small near the equator, iii) the equatorial region is not homogeneous: the residuals shown in Fig. 8 reveal cold and warm spots, at the +/-2 K level, that the model is unable to explain, even after masking or mapping the IEBR. These are the reasons why the fits obtained at the Equator (𝜒 𝑟,𝑚𝑖𝑛 2 > 1.7) are never as good as those obtained in the South pole (𝜒 𝑟,𝑚𝑖𝑛 2 ~1), and are associated with large error bars (Table 4, Fig. 6).

Nevertheless, it appears that the equatorial regions tend to have lower thermal inertias (and correspondingly higher skin depth ratios, keeping the electric skin depth of the same order) than the South Polar region, especially on the sub-Saturn side (Table 4 and Fig. 6). This could indicate higher equatorial porosities, larger amounts of amorphous water ice, larger grain sizes, and/or looser contacts between grains [START_REF] Ferrari | Low thermal inertias of icy planetary surfaces -Evidence for amorphous ice?[END_REF] near the Equator. Most of these explanations could result from the enhanced intensity of impact gardening at low latitudes.

The real part of the dielectric constant, while still low (1.20 -1.50 including 2-σ error bars), tends to be larger at the Equator than at the South pole (1.08 -1.30), which can be explained by a smaller porosity, a larger concentration of scattering and depolarizing structures in the subsurface, or more contamination by non-ice compounds at low latitudes. In addition, the resolved scatterometry data do not show any significant increase in scattering near the South pole which argues in favor of more contamination by impurities in the Equatorial region. Such impurities may have been introduced by impact gardening. Unfortunately, the loss tangent cannot provide complementary information on the relative polar and equatorial impurity content, as no lower bounds can be constrained. Another explanation for the decrease in dielectric constant with increasing latitude would be the deposition of CO2 frost at the poles as suggested by Teolis and Waite (2016) (CO2 has a lower dielectric constant than water ice, [START_REF] Paillou | Microwave dielectric constant of Titan-relevant materials[END_REF].

Sub-Saturn vs anti-Saturn

Rhea's trailing hemisphere was not well observed, and the data collected on its leading hemisphere are too affected by the IEBR for the reliable derivation of parameters. Instead, two distant scans centered on the sub-Saturn side and four distant scans, one resolved scan, and two stares on the anti-Saturn side are available. This dataset can be exploited with in mind the caveats that i) both of these regions are fairly heterogeneous (as they are located in a narrow brighter longitudinal band in-between the redder leading and trailing hemispheres [START_REF] Schenk | Plasma, plumes and rings: Saturn system dynamics as recorded in global color patterns on its midsize icy satellites[END_REF] and ii) the IEBR as well as the cold spot near 45°N and -180 to -100°E (Fig. 7) may affect the anti-Saturn results.

We find that, if the albedo is kept constant at 𝐴 𝐵 = 0.55, the sub-Saturn side has a slightly higher emissivity (0.70 vs 0.67 with an error bar of about 0.01) than the anti-Saturn side (Table 4), especially in the northern mid-latitudes (Fig. 10). The most likely explanation is that the trailing hemisphere terrain is affecting the sub-Saturn side; indeed, one of the two sub-Saturn observations is centered at 15°E, i.e., within the optically dark trailing region, whereas all of the anti-Saturn observations are centered between -136°E and -167°E, i.e., within the optically brighter leading hemisphere. Assuming that 𝑒 varies with bolometric Bond albedo 𝐴 𝐵 as (1 -𝐴 𝐵 ) -1 4 , we find that, in order to have the same emissivity on both sides, the albedos of the sub-Saturn/trailing and of the anti-Saturn/leading must be, respectively, 0.55 and 0.62, or 0.46 and 0.55, depending upon which side we impose an albedo of 0.55. This albedo contrast is of the same order as those measured between the trailing and leading by Howett et al. (2010) (0.57 -0.26 +0.20 and 0.63 -0.12 +0.11 ) and by [START_REF] Pitman | Disk-integrated bolometric Bond albedos and rotational light curves of saturnian satellites from Cassini Visual and Infrared Mapping Spectrometer[END_REF] (0.42 ± 0.10 and 0.55 ± 0.08) which suggests that the emissivity contrast we derived is most likely, in fact, an albedo difference. This would further imply that the optically dark layer that covers the trailing hemisphere, close to the sub-Saturn, is thin (at most a few cm).

Nevertheless, it remains puzzling that, in the emissivity map derived from dividing Rhea into 45° latitude/longitude regions, the positive emissivity anomaly is on the sub-Saturn side, and not strictly on the trailing side. While this could be attributed to poor constraints on limb data, it could also mean that a positive emissivity anomaly in the sub-Saturn region is really the cause of this dichotomy, rather than the albedo dichotomy discussed above. An alternative explanation is that Saturn shine is heating the sub-Saturn side. The effects of both the thermal IR flux from Saturn and the visible flux reflected by Saturn, which contribute of the order of 3\% of the total incident flux, were assumed negligible in the thermal model (Section~\ref{sec5-1}). It is possible that we are detecting a local increase in temperature at the sub-Saturn side due to these optical and IR fluxes from Saturn.

Inktomi

Based on the CIRS instrument, which probes the top millimeters of the surface, [START_REF] Howett | Thermophysical property variations across Dione and Rhea[END_REF] detected a locally higher thermal inertia in the IEBR (~19 MKS instead of the ~10 MKS equatorial thermal inertia), but no local variation in albedo. The authors attributed this detection to either larger grain sizes, or large blocks of icy ejecta among a small grain-sized surface. VIMS spectra of Inktomi and its surroundings also suggest larger grain sizes nearer to Inktomi, while the rest of Rhea's leading hemisphere is covered by a processed, fine-grained regolith [START_REF] Stephan | The Saturnian satellite Rhea as seen by Cassini VIMS[END_REF]. VIMS spectra also indicate a lower amorphous to crystalline ice ratio near recent craters, especially Inktomi [START_REF] Dalle Ore | Impact craters: An ice study on Rhea[END_REF], which would increase the thermal inertia in agreement with CIRS results.

While there is too little spatial and temporal coverage to derive parameters from the IEBR alone, the emissivity map clearly shows a low emissivity/high albedo anomaly in the IEBR (Fig. 9). This feature is actually already visible in the deconvolved data prior to any modeling (Fig. 2). In addition, we find that masking and including the IEBR lead to the same derived thermal and dielectric parameters for all regions (Table 4 and Fig. 6), implying that the IEBR does not have any highly unusual (beyond the error bars) thermal or dielectric properties at the depths sensed by the radiometer (6 -13 m). However, small amplitude (within the error bars) and/or localized (e.g., only very close to the crater) thermal and permittivity anomalies remain possible. Using the same relationship between emissivity and albedo as in Section 5.4, we find that the emissivity over the IEBR (Fig. 9) would be the same as over its surroundings if it had an albedo ~0.08 higher, which is within the range of albedo variations on Rhea. Nonetheless, the high radar brightness measured during the scatterometry experiment (Fig. 3) strongly suggests that the IEBR, at metric depths, is associated with a strong emissivity anomaly. Both the thermal inertia and the albedo may be different in this region, conceivably contributing to the behavior shown in Fig. 9, but the high scattering properties of the IEBR provide a more natural explanation to its apparently low radiometric emissivity.

Even assuming that the IEBR anomaly is only caused by a lower emissivity, the precise magnitude of this emissivity anomaly cannot be measured, as both radiometry and scatterometry have footprint sizes larger than the crater itself, and the emissivity difference decreases with distance to Inktomi. An estimate can nonetheless be obtained as follows. Averaged over a 45° region, the emissivity is about 6% lower in the IEBR (Fig. 9). The better resolved scatterometry map (Fig. 3) shows a peak increase in the backscatter of 37% at Inktomi compared to its surroundings, which would lead to a 20% peak decrease in emissivity (𝑒 = 1 -𝑅 -𝑓 × 𝜎 0 from Eq. ( 10)). Such high scattering, as well as the very high radar albedo of the IEBR (2.12 ± 0.06; Table 2) would be explained by the presence of reflecting subsurface structures back-scattering the signal to the Cassini spacecraft, more than in the rest of Rhea. Both a high 2.2-cm radar brightness and a low emissivity are likely caused by wavelength-scale inhomogeneities, such as multiple fractures in the ice blocks forming the ejecta, or large pores or grain sizes (~a few cm).

Summary and conclusion

We have presented, reduced, and analyzed all available resolved and unresolved Cassini radiometry observations of Saturn's largest icy inner satellite Rhea, as well as a resolved scatterometry observation. With a combination of thermal, radiative transfer, and emissivity models, we have simulated antenna temperatures for each observation, in order to derive global and regional constraints on the thermal, physical, and compositional properties of Rhea's regolith.

The Cassini radiometer is found to probe the subsurface down to 6 -13 m below the surface; it was the only instrument on board the Cassini spacecraft able to investigate such depths. Overall, the microwave dataset points to a very transparent (and therefore rich in pure water ice) and unpolarized (due to multiple/volume scattering) subsurface. The thermal inertia values obtained from Cassini radiometry are generally higher (> 69 MKS on the South pole, best fitting values of 50 -300 MKS everywhere) than those inferred from IR measurements (1 -46 MKS; [START_REF] Howett | Thermophysical property variations across Dione and Rhea[END_REF][START_REF] Howett | Thermal properties of Rhea's poles: Evidence for a meter-deep unconsolidated subsurface layer[END_REF], indicating an increasing degree of compaction in the regolith with depth. At the South pole, where it is better determined, the low loss tangent (< 1 × 10 -3 ) indicates at most 10% of contaminants and > 10% porosity averaged over the top 5-10 m. Very high backscatter is measured on Rhea during the resolved scatterometry observation, yielding high radar albedos of ~1.7, consistent with those found from distant diskintegrated observations by Le [START_REF] Gall | Dust and Snow Cover on Saturn's Icy Moons[END_REF], using the same instrument. As observed on Enceladus (Le [START_REF] Gall | Thermally anomalous features in the subsurface of Enceladus's south polar terrain[END_REF] and the radar-bright Xanadu region on Titan [START_REF] Janssen | Anomalous radar backscatter from Titan's surface?[END_REF], Rhea's inferred emissivities (𝑒 = 0.65 -0.70 for 𝐴 𝐵 = 0.55) seem to be not small enough to correspond to such high backscatter when invoking common random scattering mechanisms. In particular, the CBE (e.g., Black et al., 2001) predicts insufficient scattering. The presence of (maybe organized) scattering cm-scale structures in the subsurface of Rhea (and Saturn's other icy satellites) may be at the origin of both the high recorded radar backscattering coefficients and, by depolarizing the signal, of the very low inferred dielectric constants (𝜀 ′ = 1.1 -1.5). The structure of the subsurface remains to be fully investigated, as the nature and formation mechanism of these scattering subsurface structures is as yet unknown.

Regional variations suggest higher thermal inertias and smaller dielectric constants at the South pole, while there does not seem to be a significant difference in scattering properties between the high southern latitudes and the Equator. Our best interpretation is that exogenous processes such as (micro)meteoroid impacts, incoming dust from the E-ring and Phoebe ring, and electron and ionized particle impacts, which preferentially affect lower latitudes, are introducing small amounts of material with a slightly higher dielectric constant than porous water ice while also loosening the structure of the subsurface (hence a lower thermal inertia). The radiometry and scatterometry data have also shown that the Inktomi crater ejecta blanket is a low-emissivity and high-backscatter region, but these data are too limited in resolution and local time sampling to detect a thermal inertia or dielectric constant anomaly. It could be that increased scattering in the subsurface of this region is caused by numerous and organized cm-scale inhomogeneities and/or a very high density of randomly oriented fractures in the ice blocks of the ejecta blanket, consistent with the interpretations of VIMS and CIRS data [START_REF] Stephan | The Saturnian satellite Rhea as seen by Cassini VIMS[END_REF][START_REF] Howett | Thermophysical property variations across Dione and Rhea[END_REF]. High water ice purity near Inktomi (observed by VIMS: [START_REF] Stephan | The Saturnian satellite Rhea as seen by Cassini VIMS[END_REF][START_REF] Scipioni | Spectroscopic classification of icy satellites of Saturn II: Identification of terrain units on Rhea[END_REF] would also increase the probed depth, and thus increase the number of opportunities for the signal to be scattered, leading to lower emissivity and higher radar brightness. As over the rest of Rhea, it is unclear which subsurface structures and organization would lead to sufficiently strong backscatter while keeping the emissivity high.

In order to investigate the aspect and formation of these scattering structures, both theoretical and numerical modeling of the cm-wavelength scattering properties of porous, fractured, or otherwise inhomogeneous media is necessary. Laboratory studies of the thermal, dielectric, and scattering properties of ices of varying porosity and grain size at 2.2 cm would also provide the ground truth both to constrain the models and to interpret the observations. Further information on the structure and composition of Rhea's subsurface may be obtained from ground-based multi-wavelength radiometry observations. While it is challenging to observe Saturn's innermost satellites in the micro-waves due to the proximity of Saturn and its rings (affecting strongly the sidelobes), Rhea is far enough from Saturn that it could be resolved at millimetric wavelengths by radio interferometry with ALMA (the Atacama Large Millimeter Array) and centimetric wavelengths by the VLA (Very Large Array). Similar past and ongoing work on Iapetus has built microwave spectra for its leading and trailing hemisphere, allowing i) the detection of a possible absorption feature maybe related to the grain size on the trailing side (Ries, 2012), ii) similarities between Iapetus leading side and Phoebe, and iii) likely changes in thermophysical and chemical properties with depth on the leading side [START_REF] Bonnefoy | Probing the subsurface of the two faces of Iapetus[END_REF]. On Rhea, mm-cm radiometry observations would similarly help characterize the structure of the subsurface, especially by comparison with its 2.2 and 13-cm radar brightness [START_REF] Black | Arecibo radar observations of Rhea, Dione, Tethys, and Enceladus[END_REF][START_REF] Gall | Iapetus' near surface thermal emission modeled and constrained using Cassini RADAR Radiometer microwave observations[END_REF]. Moreover, Ku-band radiometry from the Earth would provide crucial daytime data to complement the Cassini equatorial night-time observations: VLA observations would help constrain the thermal inertia of Rhea's equatorial regions. The extended VLA configuration also permits resolutions very similar to those of the Cassini radiometer (beam radius of 0.2 𝑅 𝑅ℎ𝑒𝑎 for the Q band, i.e. 7 mm), and would detect features such as the IEBR or the sub-Saturn anomalous region.

While ground-based observations would provide essential information, spacecraftbased radiometry remains unique in both resolution and geometry; indeed, local night-time data are unreachable from the Earth. After the end of the Cassini-Huygens mission and its myriad discoveries, several future missions will explore icy moons. The Dragonfly mission [START_REF] Lorenz | Dragonfly: A rotorcraft lander concept for scientific exploration at Titan[END_REF], which will send a rotorcraft to Titan's equatorial regions by 2034, will conduct the in situ compositional, thermal, and physical characterization of a surface which, in spite of undeniable differences with Rhea due to the abundance of organics and erosional processes, is also an icy surface with high volume scattering [START_REF] Janssen | Titan's surface at 2.2-cm wavelength imaged by the Cassini RADAR radiometer: Calibration and first results[END_REF]2016). Jupiter's icy satellites will be visited by the JUICE (Jupiter Icy Moon Explorer) and Europa Clipper missions. In particular, the SWI (Submillimetre Wave Instrument; [START_REF] Hartogh | The Submillimetre Wave Instrument on JUICE[END_REF] aboard JUICE will lead the thermal characterization of the near-subsurface of Ganymede. The SWI observations will be especially valuable to compare to thermal observations of Jupiter's and Saturn's moons, and will help us understand the formation and evolution of icy surfaces in the Solar System.

Figure 1 :

 1 Figure 1: Observations of Rhea on March 2, 2010 (flyby RH127), showing all three types of radiometry observations: two distant scans, a stare, and a resolved scan. Top: image of the disk of Rhea as seen by the Cassini spacecraft during each type of observation. The boresight pointing directions are shown with blue dots for each observation. The half-power beamwidth at the beginning and end of each sequence is shown with a magenta circle. The direction of the sensed electric field (direction of parallel polarization) is indicated with a black double arrow. The background Rhea map is the global enhanced 3-color Rhea mosaic (PIA 18438), projected to spacecraft view coordinates. Bottom: Calibrated antenna temperatures during the RH127 observation. The baseline is set to the CosmicMicrowave Background (CMB) temperature, i.e. 2.7 K with the Rayleigh-Jeans approximation (valid at 2.2 cm). The moment the radar transmitter turns ON is indicated with a red arrow. Periods when the radar transmitter was ON and the integration time was shorter than 0.3 s are indicated in red: during these periods, the background temperature is higher due to the heating of the receiver. We therefore only use the resolved scan data acquired after minute 366, when the receiver temperature is stable (see Section 2.2.2).

Figure 2 :

 2 Figure 2: Deconvolved brightness temperatures measured on Rhea, during all 10 distant scans and the unique available well resolved observation RH127_4r. North is towards the top. The position of the Inktomi crater is marked with a black dot. The half-power beam width is shown with a magenta circle; this corresponds to the mean resolution of each image. Note that all observations were acquired during the local night-time, so the sub-solar point is never visible.

Figure 3 :

 3 Figure 3: Analysis of the RH127 scatterometry data. a) Normalized radar cross-section 𝜎 𝑎𝑣𝑔 0 versus effective incidence angle 𝜃 𝑒𝑓𝑓 . The data collected on the leading hemisphere (center longitude > -160°E), the trailing hemisphere (center longitude < 160°E) and the position of the Inktomi crater are indicated in different colors and are fitted separately with the function 𝑓(𝜃 𝑒𝑓𝑓 ) = 𝐴𝑐𝑜𝑠 𝑛 𝜃 𝑒𝑓𝑓 . The derived averaged radar albedo of each region is indicated. b) Values of 𝜎 𝑐𝑜𝑟𝑟 0 , the normalized radar cross-section corrected to an incidence angle of 32° using a fitting law through all the data, plotted at the center of each observation. The position of the leading hemisphere (LH), the trailing hemisphere (TH), and the ejecta blanket of crater Inktomi are indicated. On the background map, latitudes are plotted every 15°; longitudes are plotted every 30°. Note the increased backscatter near Inktomi. The background Rhea map is the global enhanced 3-color Rhea mosaic (PIA 18438), projected to spacecraft view coordinates.

Figure 4 :

 4 Figure 4 : Visualization of the thermal emission model for the following parameter values: 𝜀′ = 1.15, 𝐼 = 250 MKS and 𝑟 = 100. The quasi-specular component 1 -𝑅(𝜀′) and the simplified scatterometry map 𝜎 0 2𝑛𝑐𝑜𝑠 𝑛 𝜃 used in the diffuse component of the emissivity model are shown in the first two diagrams, followed by the effective temperature map 𝑇 𝑒𝑓𝑓 (𝐼, 𝑟) and the resulting brightness temperature distribution 𝑇 𝑏 (𝐼, 𝑟, 𝜀′, 𝑓) = 𝑒(𝜀′, 𝑓) × 𝑇 𝑒𝑓𝑓 (𝐼, 𝑟). The positions of the Brewster angles (where 𝑅(𝜀′) = 0) are shown with circles; Inktomi is marked with an x. All maps are shown as seen by the Cassini spacecraft at the time of the RH127 observation (the sub-spacecraft latitude was 0°N).

  Figure 5: Disk-integrated emissivities for different combinations of thermal inertia and skin depth ratio obtained from the 10 distant scans of Rhea (see Table 1). A uniform bolometric Bond albedo of 0.55 is assumed. Oblique lines indicate constant values of the electric skin depths. Note that electric skin depths below 10 cm (<5 times the wavelength), which yield the highest emissivities, are unlikely in a water-ice regolith. On the central map of Rhea, the position of the Inktomi crater ejecta blanket is indicated with a white circle. The sub-Saturn point is at 0°E and the anti-Saturn point at +/-180°. The background Rhea map is the global enhanced 3-color Rhea mosaic (PIA 18438).

Figure 6 :

 6 Figure 6: 2-σ best fitting ellipses in the (𝐼, 𝑟) domain, for best-fit values of 𝜀′ (see Table4for the other parameters) for each of the 5 regions where the fitting method was applied, and for each of the two different ways to treat the Inktomi ejecta blanket (masked or not). The minimum reduced 𝜒 2 value (𝜒 𝑟,𝑚𝑖𝑛

Figure 7 :

 7 Figure 7: 2-σ best fitting ellipses in the (𝐼, 𝑟) domain for the South pole (> 50°S), for four different values of the bolometric Bond albedo 𝐴 𝐵 and for 𝜀 ′ = 1.2. Note that the derived ratio remains the same, but the inferred thermal inertia and the electric skin depth 𝛿 𝑒𝑙 decrease with increasing albedo. The value of the emissivity at nadir for the best fit parameters is indicated; as expected it increases with albedo.

Figure 8 :

 8 Figure 8: Residual maps using all data together a) Assuming emissivity is uniform over the whole satellite b) Excluding data with an Inktomi contribution > 2% c) Modeling the Inktomi ejecta blanket as a change in emissivity, Gaussian with distance from Inktomi. d), e), and f) use the same treatment of Inktomi as a), b), and c), respectively, while dividing f into 45° latitude and longitude bins. In each map, the position of the Inktomi ejecta blanket is indicated by a black circle.

Figure 10 :

 10 Figure 10: Results of the fitting procedure in the South pole of Rhea (>50°S), assuming a uniform emissivity map and a bolometric Bond albedo of 𝐴 𝐵 = 0.72. a) 2-σ best fitting ellipses in the (𝐼, 𝑟) domain. b) 2-σ best fitting ellipses in the (𝐼, 𝜀′) domain. c) 2-σ best fitting ellipses in the ( 𝑟, 𝜀′) domain. d) Observed and best-fit modeled antenna temperatures 𝑇 𝑎 , given as a time sequence. Vertical dotted lines separate the data from different flybys, in order: RH011 (July 2005), RH045 (May 2007), RH049 (August 2007), RH127 distant then resolved (March 2010), and RH177 (December 2012). The southern fall equinox occurred in August 2009, between the RH049 and RH127 flybys, as indicated. Note that the 2005 temperatures (RH011) are lower even though they are in summer, because the resolution is not as high as for the later data (see Table1). e) Modeled 𝑇 𝑎 for 𝐼 and 𝑟 constant, and 𝜀′ varying. f) Modeled 𝑇 𝑎 for 𝜀′ and 𝑟 constant, and 𝐼 varying. g) Modeled 𝑇 𝑎 for 𝜀′ and 𝐼 constant, and 𝑟 varying. Note that both the thermal inertia and the skin depth ratio are reasonably well constrained due to the observations at different seasons, whereas the dielectric constant is determined from the higher-resolution datasets (RH127 and RH177).

Figure 11 :

 11 Figure 11: Ternary diagram showing the compositional constraints for Rhea's near-subsurface, assuming that it consists in porous water ice with tholin inclusions. The maximum values for 𝜀′′ include the effect of the albedo range (0.42 -0.72). The region corresponding to acceptable values of the imaginary components of the permittivity is shaded in green for the equatorial regions and in blue for the South pole region. Based on the constraints from the fit while masking Inktomi: 𝜀 ′′ < 1.1 × 10 -3 on the South pole region, and 𝜀 ′′ < 4.4 × 10 -3 for the equatorial regions. To help reading this diagram, an example composition is shown in red, corresponding to a mixture of 20% vacuum, 30% tholins, and 50% water ice.

Table 1 :

 1 Cassini radiometry observations of Rhea. Distant scans are identified with a "u" in the observation IDs, stares with an "s", and resolved scans with an "r". The main observed region is indicated next to the sub-spacecraft point coordinates. The symbols ☾ and ☼ indicate nighttime and daytime observations, respectively (neither is included for polar observations). Unusable data (due to temperature changes caused by the radar transmitter being ON or to lack of baseline observations for calibration) are indicated with the symbol X, whereas the data that we used are indicated with the symbol ✓. The beam radius is given for the beginning and end of each scan.

	Observation	Date	Beam	Sub-spacecraft	Sub-solar	Central
	ID	Start time -End time	radius (in	point	point	local time
			𝑅 𝑅ℎ𝑒𝑎 )	(°E, °N)	(°E, °N)	(hh:mm)
	RH11_1s	2005 JUL 14	0.85 -0.81	(-68, -77)	(177, -20)	19:42
	X	08:15:02 -09:00:03		South pole
	RH11_2u	2005 JUL 14	0.81 -0.78	(-48, -75)	(175, -20)	21:08
	✓	09:05:44 -09:41:44		South pole	☾/☼
	RH18_1s	2005 NOV 27	0.47 -0.69	(-18,0)	(161,-19)	00:02
	X	02:47:03 -04:33:24		Sub-Saturn	☾
	RH18_2u	2005 NOV 27	0.69 -0.79	(-22, 0)	(157, -19)	00:04
	✓	04:35:35 -05:21:55		Sub-Saturn	☾
	RH22_1s	2006 MAR 21	0.39 -0.45	(-114, 1)	(59, -18)	00:30
	X	09:40:22 -10:55:01		Leading	☾
	RH22_2u	2006 MAR 21	0.47 -0.57	(-136, 1)	(52, -18)	23:30
	✓	11:03:26 -12:57:06		Saturn Leading, Anti-	☾
	RH22_3u	2006 MAR 21	0.60 -0.67	(-153, 1)	(46, -18)	22:47
	✓	13:03:16 -14:57:06		Saturn Leading, Anti-	☾
	RH22_4s	2006 MAR 21	0.72 -0.84	(-166, 1)	(40, -18)	22:19
	✓	15:05:25 -16:28:45		Anti-Saturn	☾
	RH27_1s	2006 AUG 17	0.80 -0.75	(41, 24)	(84, -16)	09:07
	X	05:38:14 -06:24:47		Saturn Trailing, Sub-	☼
	RH45_1u	2007 MAY 27	0.85 -0.83	(-58, -44)	(-157, -12)	18:36
	✓	11:21:28 -11:49:44		pole Leading, South	☾/☼
	RH45_2s	2007 MAY 27	0.83 -0.79	(-64, -45)	(-158, -12)	18:19
	✓	11:49:43 -13:00:12		pole Leading, South	☾/☼
	RH47_1s	2007 JUN 28	0.61 -0.68	(104, -3)	(163, -11)	08:04
	X	14:40:19 -15:13:45		Trailing	☼
	RH49_1u	2007 AUG 29	0.55 -0.47	(15, 0)	(-114,-11)	20:39
	✓	20:05:06 -20:47:36		Sub-Saturn	☾
	RH49_2s	2007 AUG 29	0.47 -0.37	(12, 0)	(-117, -11)	20:37
	X	20:47:35 -21:40:45		Sub-Saturn	☾

Table 2 :

 2 Results of the re-analysis of the scatterometry data, after correcting the antenna pointing offset.

		𝐴	𝑛	𝐴 𝑎𝑣𝑔
	Western leading (LH)	2.14 ± 0.01	1.39 ± 0.02	1.79 ± 0.02
	Eastern trailing (TH)	2.08 ± 0.02	1.50 ± 0.02	1.67 ± 0.03
	Inktomi	2.40 ± 0.04	1.26 ± 0.03	2.12 ± 0.06

Table 3 :

 3 Parameters of the thermal and radiative transfer models

	Parameter	Values	Comments/references
	Bolometric Bond	0.42, 0.55,	Pitman et al. (2010) used VIMS data to derive bolometric
	albedo 𝐴 𝐵	0.63, 0.72	Bond albedos of 0.55 ± 0.08 and 0.42 ± 0.10, while Howett et al. (2010) found 0.63 -0.26 +0.20 and 0.57 -0.12 +0.11 , for the leading
			and trailing hemispheres respectively, and Howett et al.
			(2016) derived 0.72 ± 0.02 for the South pole.
	Thermal inertia	10, 20, 50,	Howett et al. (2014) found 𝐼 = 10 -20 MKS at the surface
	𝐼 (Jm -2 K -1 s -1/2 )	100, 250,	in the equatorial regions. Howett et al. (2016) found 𝐼 = 1 -
		500, 1000	46 MKS at the South pole down to about 1 m. The bulk
			thermal inertia of crystalline water ice is 2000 MKS.
	Skin depth ratio 𝑟 = 𝛿 𝑒𝑙 𝛿 𝑡ℎ 𝑑𝑎𝑦	0.01 to 1000 The diurnal thermal skin depth 𝛿 𝑡ℎ 𝑑𝑎𝑦 can be derived from the thermal inertia using equation (6) with 𝑃 = 4,518 Earth days.

Table 4

 4

	gives the best-fit parameters (𝐼, 𝑟 = 𝛿 𝑒𝑙 𝛿 𝑡ℎ 𝑑𝑎𝑦 ⁄	, 𝜀′, and 𝑓 =	1+µ 𝐿 𝑓 𝐶𝐵𝐸

Table 4 :

 4 Best-fit parameters and associated 2-sigma error bars in different regions of Rhea, assuming a uniform global bolometric Bond albedo of 𝐴 𝐵 = 0.55.

Acknowledgements

The authors acknowledge financial support from the Île-de-France region DIM-ACAV+ (Domain d'Intérêt Majeur en Astrophysique et Conditions d'Apparition de la Vie) project, as well as the IUF (Institut Universitaire de France). We also wish to thank the Cassini radar team for providing the unique dataset used herein, as well as two anonymous reviewers for their insightful comments, which improved the quality of this work.

Appendix A: Correction of the pointing offset

A constant pointing offset of 0.036° (i.e. 10% of the main beam width) has been reported by [START_REF] Zhang | Cassini microwave observations provide clues to the origin of Saturn's C ring[END_REF], which we have corrected for by modifying the pointing coordinates in the NAIF SPICE Cassini spacecraft frame definition kernel. In addition to this constant offset, there are also small variations in the pointing offset, of the order of 𝑥 𝑜𝑓𝑓𝑠𝑒𝑡 𝑦 𝑜𝑓𝑓𝑠𝑒𝑡 ± 0.01°, where 𝑥 and 𝑦 correspond to the 𝑥and 𝑦-axes of the Cassini Radar antenna coordinate system. While generally limited, these offsets have a significant impact on the interpretation of the antenna temperatures near the limbs; it is therefore crucial to correctly determine them. There is also an inaccuracy of 𝑡 𝑜𝑓𝑓𝑠𝑒𝑡 < 0.6 s in the observation time.

For "distant scans", we fit a simulated Rhea disk convolved with the known beam pattern to the observed data near the limbs. We only use the data collected near the limbs, since data with a large contribution from Rhea could include actual surface brightness temperature variations that are not due to a pointing offset, and data entirely off the limbs (in the cold sky) are only affected by noise and variations in the temperature baseline, independently of the pointing offset. At each observation time, the spacecraft position and orientation with respect to Rhea are known thanks to the SPICE kernels and can be used to compute Rhea's angular size, the antenna polarization and its theoretical pointing direction (𝑥, 𝑦). These latter parameters are interpolated to time 𝑡 + 𝑡 𝑜𝑓𝑓𝑠𝑒𝑡 , then the pointing direction is adjusted to (𝑥 + 𝑥 𝑜𝑓𝑓𝑠𝑒𝑡 , 𝑦 + 𝑦 𝑜𝑓𝑓𝑠𝑒𝑡 ). To find the offsets, we first run the fit with a simple model of Rhea's surface temperatures, which allows for a latitudinal temperature variations such that

. This simple model has five free parameters (𝑥 𝑜𝑓𝑓𝑠𝑒𝑡 , 𝑦 𝑜𝑓𝑓𝑠𝑒𝑡 , 𝑡 𝑜𝑓𝑓𝑠𝑒𝑡 , 𝑇 𝑒𝑞𝑢𝑎𝑡𝑜𝑟 , ∆𝑇). While it produces a reasonable fit to the data near the limbs, it appears unable to simulate longitudinal temperature variations on Rhea, which contribute even to the data near the limbs. A more accurate model of Rhea surface temperatures is described in Section 3. Using the offsets given by the latitudinal model, we look for the best fitting thermal inertia, skin depth ratio, and dielectric constant following to the method detailed in Section 3.5. We then select the effective surface Rhea temperature map for these best parameters: this temperature map is then used to fit for the offsets again, with (𝑥 𝑜𝑓𝑓𝑠𝑒𝑡 , 𝑦 𝑜𝑓𝑓𝑠𝑒𝑡 , 𝑡 𝑜𝑓𝑓𝑠𝑒𝑡 ) being the only parameters. If the offsets found this way are considerably different from those given by the latitudinal model, we reiterate. Two examples of the offset fitting results are given in Fig. A1. For the "stares" and "resolved scans", there is very little data near the limbs, and we cannot use the same method. Instead, we take the same offsets as those found for the "distant scans" performed during the same flyby. 2017) are applied to the SPICE kernels. c) After applying the initial offsets to the kernels, the data are fitted by the model described in Section 3 and a small additional offset is inferred. Note that an anomalously large offset was found for RH049_1u, even after fixing the spacecraft pointing direction in the SPICE kernels.

Appendix B: Gaussian component of the measurement noise

As described in Section 2.2, the uncertainty associated with the calibrated antenna temperature 𝑇 𝑎 measured by the Cassini radiometer is composed of i) random Gaussian noise 𝜎 𝐺 , and ii) a 1% calibration uncertainty that globally takes into account uncertainties in the gain and the calibration factors [START_REF] Janssen | Titan's surface at 2.2-cm wavelength imaged by the Cassini RADAR radiometer: Calibration and first results[END_REF]2016). We determine 𝜎 𝐺 empirically by examining the standard deviation of the data as follows.

For each flyby, observations with constant average temperatures are selected: these are either observations of the sky or of a fixed point on Rhea during stares. The calibrated antenna temperatures are smoothed over 20 consecutive bursts. The difference between the observed and smoothed antenna temperatures, plotted in Fig. B1a, provides an estimate of the noise level. For each flyby and for each integration time, the standard deviation of the noise is calculated. As expected, we find that the amplitude of the Gaussian noise is inversely proportional to the square root of the integration time, with an added offset taking into account the constant instrument read noise: 𝜎 𝐺 = 0.057/√𝑡 𝑖𝑛𝑡 + 0.046 (Fig. B1b). Because the noise level varies near 0.8 to 1 s integration times, we use the empirically determined value of 𝜎 𝐺 for each integration time when calculating the 𝜒 𝑟 2 (Eq. ( 12)). We note that, for an integration time of 𝑡 𝑖𝑛𝑡 = 1 s, we find a measurement noise of ~0.1 K, which is about four times higher than the theoretical value 0.026 K given in [START_REF] Janssen | Titan's surface at 2.2-cm wavelength imaged by the Cassini RADAR radiometer: Calibration and first results[END_REF]2016). After re-examination of this previous work, an error of a factor √10 has been found in the calculation of the measurement noise, which should therefore be 0.075 K for 𝑡 𝑖𝑛𝑡 = 1 s, that is, much closer to the value found here empirically (~0.1 K). The small error found here does not affect previously published results, for which the uncertainty was empirically determined.