Thermal properties of the leading hemisphere of Callisto inferred from ALMA observations

We present a thermal observation of Callisto's leading hemisphere obtained using the Atacama Large Millimeter/submillimeter Array (ALMA) at 0.87 mm (343 GHz). The angular resolution achieved for this observation was $\sim$$0.16^{\prime\prime}$, which for Callisto at the time of this observation ($D\sim 1.05^{\prime\prime}$) was equivalent to $\sim$6 elements across the surface. Our disk-integrated brightness temperature of 116 $\pm$ 5 K (8.03 $\pm$ 0.40 Jy) is consistent with prior disk-integrated observations. Global surface properties were derived from the observation using a thermophysical model (de Kleer et al. 2021) constrained by spacecraft data. We find that models parameterized by two thermal inertia components more accurately fit the data than single thermal inertia models. Our best-fit global parameters adopt a lower thermal inertia of 15-50 $\text{J}\:\text{m}^{-2}\:\text{K}^{-1}\:\text{s}^{-1/2}$ and a higher thermal inertia component of 1200-2000 $\text{J}\:\text{m}^{-2}\:\text{K}^{-1}\:\text{s}^{-1/2}$, with retrieved millimeter emissivities of 0.89-0.91. We identify several thermally anomalous regions, including spots $\sim$3 K colder than model predictions co-located with the Valhalla impact basin and a complex of craters in the southern hemisphere; this indicates the presence of materials possessing either a higher thermal inertia or a lower emissivity. A warm region confined to the mid-latitudes in these leading hemisphere data may be indicative of regolith property changes due to exogenic sculpting.


INTRODUCTION
The Galilean moons of Jupiter -Io, Europa, Ganymede, and Callisto -form a continuum of geologic activity for which the inactive endmember is Callisto. A satellite whose surface is dominated almost exclusively by impact craters, Callisto retains one of the oldest surfaces (∼4.5 Gyr) in the Solar System. Although it is the third largest planetary moon and may harbor a subsurface ocean (Zimmer 2000), Callisto bears little to no evidence for any tectonics, volcanism, or other manifestations of surface-interior interactions (Moore et al. 2004). By contrast, geologic activity is readily observed on the other Galilean satellites. For instance: Io is volcanically active in the present day and thermal emission from its colorful eruptions is routinely observed using ground-based telescopes (de Kleer et al. 2019) and spacecraft (Mura et al. 2020;Davies et al. 2006); Europa has a young, water-rich surface accompanied by evidence for surface-interior interactions (Carlson et al. 2009); and Ganymede is encircled by bright, ribbon-like, tectonically-deformed regions suggestive of major past activity (Pappalardo et al. 2004). Therefore, as the only Galilean moon with a surface dominated by exogenic sculpting, Callisto stands out as one of the best records of long-term impact modification the Solar System has to offer.
As a consequence of its inactivity, Callisto's surface at the global scale comprises just a few geologic units. These units include cratered plains, light plains, and smooth/dark plains, as well as units linked to multi-ring basins or other impacts (Moore et al. 2004). In particular, Callisto's multi-ring impact basins are the largest in the Solar System. Valhalla, the largest of these basins, has outer rings extending to a diameter of ∼3800 km and occupies ∼16% of Callisto's surface area. Structures like Valhalla are possible relics of a time in Callisto's history where the surface had a higher thermal gradient and was more ductile (Schenk 1995) and represent regions of interest to look for exposed subsurface materials. While the multi-ring basins are the most obvious large features on Callisto, the dominant landform is the heavily cratered plains. These swaths of impact-scarred terrain are often blanketed by a dark material for which the exact composition and origin remains uncertain, with endogenic explanations often focusing on surface erosion/sublimation processes (Moore et al. 2004) and exogenic theories usually invoking dust infall from the irregular Jovian satellites (Bottke et al. 2013). Although Callisto retains the largest share of dark material (its albedo is 0.2, the lowest of the Galileans, Moore et al. 2004), Europa and Ganymede also have similarly dark products on their surfaces (Greeley et al. 2004;Pappalardo et al. 2004) . Altogether, these features contribute to the fact that Callisto's quiescent surface remains poorly understood.
For many Solar System objects, unknowns pertaining to unique geologic features or the balance of exo-and endogenic processes on their surfaces have been investigated by inferring material properties from thermal observations. At wavelengths in the infrared and (sub)millimeter, radiation originates from the surface/near-subsurface and can be leveraged to infer properties such as emissivity, thermal inertia ( = J m −2 K −1 s −1∕2 ; henceforth thermal inertia will be given without units but always uses these SI units), porosity, etc. The variation of such properties across a surface encodes information as to how that surface evolved and provides clues regarding the nature of the external modification environment. For example, in the Saturnian system, thermal infrared Cassini Composite Infrared Spectrometer (CIRS) measurements of the small icy moons Tethys (Howett et al. 2012(Howett et al. , 2019 and Mimas (Howett et al. 2011) revealed both satellites have a high thermal inertia feature centered at mid-latitudes on their leading hemispheres, demonstrating the Saturnian particle environment is capable of altering surface texture (Howett et al. 2011;Schaible et al. 2017;Schenk et al. 2011). At the Moon, measurements across the ∼13-400 µm regime from the Diviner Lunar Radiometer onboard the Lunar Reconnaissance Orbiter showed that the Moon's regolith is uniform at global scales, suggesting processes like impact gardening work quickly to homogenize the upper 10 cm; additionally, maps of local thermal inertia variations helped assign plausible impact ejecta deposits to parent craters (Hayne et al. 2017). In the Jovian system, recent work with the Atacama Large Millimeter/submillimeter Array (ALMA) has proved useful for characterizing the surface of Io (de Pater et al. 2020), confirming the surficial origin of apparent nighttime hot spots on Europa (Trumbo et al. 2017(Trumbo et al. , 2018, and identifying plausible exogenic thermal trends on Ganymede (de Kleer et al. 2021a). With all of these studies, surface trends were identified with hemisphere-scale or global-scale coverage of the planetary body. In particular, some (sub)millimeter/centimeter observations may probe below the strongest time-of-day temperature variations (i.e., below the diurnal thermal skin depth) and below the most heavily processed upper surface layers.
At present, the thermal properties of Callisto's surface are poorly constrained because there are only a few disk-resolved studies at thermal wavelengths. In the infrared, the Voyager Infrared Interferometer Spectrometer and Radiometer (IRIS) instrument returned a featureless spectrum of Callisto, offering little information regarding composition (Spencer 1987). However, the eclipse cooling curves showed that Callisto's surface is not thermally uniform, as the best-fit model adopted a two-component thermal inertia surface with one higher component of = 300 and a lower component of = 15 (Spencer 1987). Moreover, the IRIS spectra exhibited a trend of increasing brightness temperature with decreasing wavelength, suggesting unresolved temperature contrasts on the surface; this trend also varied with solar incidence angle, hinting that Callisto's surface may be less smooth than Ganymede's, which did not exhibit this solar incidence angle dependence (Spencer 1987). The thermal IR temperature of Callisto's surface was measured to be ∼158 K at 20 µm (Morrison et al. 1972;Moore et al. 2004) and ∼140-150 K at ∼10 µm (de Pater et al. 2021), while the H 2 O ice surface temperature derived using H 2 O spectral features is somewhat lower at ∼115 K (Grundy et al. 1999). Past disk-integrated radio thermal wavelength measurements of Callisto's brightness temperature show a decrease from 135 K (de Pater et al. 1989) in the sub-mm regime to about ∼90-100 K at cm wavelengths (Pauliny-Toth et al. 1974;Muhleman et al. 1986;Butler 2012;Berge & Muhleman 1975;de Pater et al. 1984). Although other disk-integrated results for Callisto fill in the gaps between these wavelengths, there is a scarcity of spatially resolved datasets.
In this work, we present the first high-resolution thermal observation of Callisto using ALMA by mapping the leading hemisphere at 0.87 mm (343 GHz). Section 2 describes the observations, data analysis, and flux calibration. Section 3 details the thermophysical model used to interpret the data and describes how surface properties are derived. Section 4 presents the final, calibrated leading hemisphere image of Callisto as obtained with ALMA and includes the results and interpretation of the thermophysical modeling analysis. Section 5 summarizes the conclusions of this work.

Observations
Our observations of Callisto were acquired using ALMA. Located on the Chajnantor plateau in the Atacama desert, Chile, the main ALMA array is composed of 54 12-meter radio antennas that are linked via a correlator to function as an interferometer. As an interferometer, ALMA can achieve higher spatial resolution on planetary objects than is possible using single-dish facilities at the same wavelength. The spatial resolution ALMA can achieve, combined with its frequency coverage, make this telescope facility ideal for conducting thermal studies of the Galilean satellites from Earth.
We observed the leading hemisphere of Callisto using ALMA on 2016 November 01. The data were taken using the Band 7 receiver at a central frequency of 343 GHz (0.87 mm) with four spectral windows collectively spanning ∼8 GHz and an on-source integration time of 121 s. During observations, the array was in configuration C-5 with 39 12-m antennae. The baselines ranged from 18.6 m to 1.1 km, affording an angular resolution of ∼0.16 ′′ . The angular diameter of Callisto during data acquisition was ∼1.05 ′′ , so ∼6 resolution elements across the disk were achieved, equivalent to a spatial footprint of ∼790 km at the satellite's distance. At the time of the observation, the sub-observer longitude was 50 • W. Calibration for array pointing, flux density, bandpass response, and phase were obtained via observations of quasars J1256-0547 and J1232-0224.

Data Analysis
After observations, the data were provided as a pipeline-calibrated measurement set containing the visibilities, which are the fundamental quantities measured by an interferometer and are equal to the Fourier transform of the sky brightness temperature distribution at discrete spatial frequencies (Muders et al. 2014). For a review of the interferometric observation and imaging of Solar System objects, see Butler & Bastian (1999). Because continuum observations of the Galilean satellites typically have very high signal-to-noise ratios (SNR) and easily modeled shapes (disks), we improved the SNR of the final pipeline-generated images by applying an iterative self-calibration procedure (for an in-depth explanation of self-calibration, see Brogan et al. (2018)). To conduct this processing, we used the Common Astronomy Software Application (CASA) package (McMullin et al. 2007). As ALMA observations at 0.87 mm (343 GHz) probe continuum emission from the subsurface, we integrated the signal over the entire 8 GHz of bandwidth using multi-frequency synthesis (Sault & Wieringa 1994). Our phase-only self-calibration using CASA followed as such: for the first round of imaging, a Lambertian disk scaled to the size and brightness temperature of Callisto was used as a startmodel for the tclean task (Rau & Cornwell 2011); no additional clean components were added. Next, the complex antenna-based gains were calculated using the CASA gaincal task on a interval spanning the full integration time (i.e., > 121 s) and then applied using applycal. Then, the corrected visibilities were imaged again adopting a shallow clean, using the previous tclean output image as the startmodel. This cycle of cleaning and computing gains continued down to an interval of 2 s, at which the image SNR ceased to substantially improve. A dynamic range of ∼580 (peak flux/rms in units of Jy/beam) was obtained for the final image. For tclean parameters, we adopted a Briggs weighting scheme with a robust parameter of 0.0, and used the "Clark" deconvolver (Clark 1980). A primary beam correction was applied to the self-calibrated image product. The dimensions of the final clean beam were 0.16 ′′ × 0.19 ′′ . The rms (root-mean-square) noise of the final image (measured from a non-source region of the non-primary beam corrected product) was 0.585 mJy/beam, or about 0.2 K.
After self-calibrating the data, we derived the disk-integrated flux density directly from the visibilities. We report disk-integrated flux densities from these fits rather than from the images because the visibilities are the more direct data product. Reporting the disk-integrated flux is useful for placing our results in context with previous radio observations of Callisto that are largely spatially unresolved. The visibilities were fit with Bessel functions using the CASA task uvmodelfit assuming a uniform disk model, fitting to all spectral windows, and by excluding baselines >200 m for an optimal fit, as the longer baselines are sensitive to the variations across the disk, rather than just the disk-integrated flux. The final disk-integrated brightness temperature ( ) was calculated by inverting the following equation 1 : where is the disk-integrated flux density in Janskys, is the radius of Callisto in arcseconds and, in SI units, ℎ is Planck's constant, is the speed of light, is the observation frequency, is the Boltzmann constant, and is the cosmic microwave background (2.7 K). The uncertainties in the final incorporate flux density scale calibration errors, which was parameterized by a 5% error of the total flux density, and visibility fitting errors on the order of 0.4 mJy. Our final calibrated image is presented in Fig. 1. Result of subtracting a sample Lambertian disk from (a) to highlight difference in across the disk. In both images, the ellipse in the lower left corner represents the FWHM (full-width at half-maximum) of the synthesized ALMA beam. The latitude/longitude grid is spaced at 30 • increments, and the 0 • longitude line is marked with dashes. Callisto's north pole is aligned with the image axis. Notably, the Valhalla impact crater appears slightly colder than other disk regions viewed at similar incident and emission angles.

ALMA Flux Density Calibration Check
The ALMA pipeline applies a flux density calibration to the target data during initial processing. Because the ALMA calibrator sources are quasars (which are variable) and the ALMA pipeline occasionally does not use the most complete calibrator information, manual flux calibration correction is sometimes required (e.g., Trumbo et al. 2018;de Kleer et al. 2021a). We checked whether such a correction was required for these data by comparing the pipeline flux model to the ALMA Calibrator Catalogue 2 of all observations of our calibration source in the period surrounding our observations. The quasar used for flux calibration was J1256-0547, which the calibrator catalogue recorded observations of at Bands 3 and 7. A check of flux density measurements in these bands between 2016-10-27 and 2016-11-01 (i.e., before and after our data were obtained) indicated it is not necessary to correct the pipeline calibration of these data.

THERMOPHYSICAL MODEL
The final product of calibrating and imaging the ALMA Callisto visibilities is a map of brightness temperature across the satellite's disk ( Fig. 1). In order to interpret this observed map, we fit the data with a thermophysical model that treats the transport of heat by conduction and radiation through Callisto's surface. For each latitude and longitude on the satellite, the model constructs a temperature profile that evolves with time t and depth z into the surface. Practically, the thermophysical model numerically solves the 1D heat diffusion equation, which is given by: where is the density, is the heat capacity, is the temperature, and is the thermal conductivity. For this differential equation, the boundary conditions adopted no heat flow at depth and solar insolation at the distance of Callisto based on spatially varying albedo, as described in Section 3.2. The model was evolved over ∼15 Callisto days (1 Callisto day = 16.69 Earth days) in time steps of at least 1/500 of a day to allow for equilibration, and was run up to several thermal skin depths below the surface. The key tunable parameters in the model are the thermal inertia and emissivity. To find the best-fit model for our Band 7 Callisto image, we performed model runs over a wide range of thermal inertias, ranging from = 15 to = 2000, corresponding to values appropriate for very unconsolidated material and water ice, respectively (Ferrari & Lucas 2016). For each model, we integrated the emission at depth along the line of sight to account for subsurface radiation that is probed at ALMA thermal wavelengths. , one using = 2000 and another using = 15. Positive (red) and negative (blue) residuals correspond to surface regions in the data that are warmer and cooler than predicted by the model, respectively. The figure demonstrates that very low thermal inertia models predict too little limb emission relative to disk center, while very high thermal inertia models provide a poor fit for the opposite reason; a mixture of the two models provides a better fit without systematic center-to-limb trends in the residuals.
Callisto's (sub)surface at the viewing geometry of the observations, convolved to the resolution of the ALMA data. The model image is then compared directly to the observations to determine the best-fit thermal inertia and emissivity.

Two-Component Mixtures
Although some models generated using a single thermal inertia provided a better fit for either Callisto's limb or center-of-disk, no single model could fit both parts, regardless of treatment of the surface emission incidence angle dependence (e.g. by using smooth or rough surface implementations). We further considered two-component surface models, which have previously been employed to fit observations of the Galilean satellites, and Callisto specifically. For example, Spencer (1987) reported that Voyager infrared cooling curves of Callisto that displayed an initial, fast cooling followed by a slower cooling rate could be modeled by a surface with a high thermal inertia ∼ 300 overlain by a thin, low thermal inertia layer of ∼ 15. Here, we explored how the systematics in the single model residuals (e.g., rings on the limb in the residuals) could be addressed by a two-component approach.
To create each two-component model, we conducted the following: first, two sets of unconvolved models were generated; the first set adopted thermal inertias ranging from = 15-400 and the second set used ℎ ℎ = 500-2000. Each member of the low thermal inertia set was linearly combined with each member of the high thermal inertia set, forming two-component pairs. For each pair of and ℎ ℎ models, we created 11 linear mixtures that ranged from 0% + 100% ℎ ℎ , to 100% + 0% ℎ ℎ ; the models were mixed in increments of 10% (see Fig. 2 for an example mixture). The spatial alignment of each model pair was checked for proper overlap. After the components were mixed, the two-model was convolved with the synthesized ALMA beam and subtracted from the data. The emissivity for each mixed model was obtained by performing a least-squares analysis, with the off-disk flux masked out for the calculation. Goodness of fit for the thermophysical models was determined using the cost function, which takes into consideration the number of resolution elements (i.e., ALMA beams). Details for applying the cost function to ALMA data are described in Section 3.4 of de Kleer et al. (2021b); we note that the cost function is formally limited to data sets with independent data points.

Albedo Map
Our model does not treat the albedo of Callisto as a free parameter, rather it ingests a bolometric (bond) albedo map derived from spacecraft data, following Trumbo et al. (2017, 2018) and de Kleer et al. (2021a. For this map, we combined information from: (1) a high-resolution (∼1 km) USGS greyscale map of Callisto's surface 3 , with (2) spectral albedos measured at 0.35, 0.41, 0.48, and 0.59 µm obtained at 20 locations by Voyager (Johnson et al. 1983). First, the coordinates of the 20 spectral albedo measurements were used to identify 20 corresponding locations on the USGS map. For each of these locations on the USGS map, an average greyscale value was derived from a 40 × 40 km box matched to the resolution of Voyager. Next, because the 0.35 -0.59 µm range covers only ∼80% of the solar spectrum, we extended the spectral coverage by scaling 0.6 -2.5 µm disk-integrated reflections to our data (Clark & McCord 1980). Then, we computed a weighted, wavelength-integrated albedo for each of the 20 points, derived a linear function relating greyscale value to wavelength-integrated albedo, and then mapped this function to the entire USGS map (Fig. 3). The wavelength-integrated map was then multiplied by a phase integral of 0.51 (Buratti 1991) to convert it to a bolometric albedo map. We replaced poorly sampled regions near the poles in the greyscale map with a global average; the exact values of the albedo map at the poles has negligible impact on the thermal model fits since the viewing geometry of the data is nearly equatorial.

RESULTS & DISCUSSION
We present surface properties derived from an ALMA 0.87 mm (343 GHz) observation of Callisto's leading hemisphere, as well as brightness temperature residual maps created with best-fit thermophysical models. In Section 4.1, we present an integrated measurement, and place our result in context with past, unresolved thermal wavelength observations. In Section 4.2, we describe the single thermal inertia model results. In Section 4.3, we discuss the results from modeling the data with two thermal inertia components. Lastly, in Section 4.4 we detail the connections between regional temperature signatures and the local geology.

Disk-Integrated Brightness Temperature
We determined the disk-integrated flux density of Callisto's leading hemisphere is 8.03 ± 0.40 Jy at 343.5 GHz (0.87 mm), which is equivalent to a brightness temperature of 116 ± 5 K. In Fig. 4, we place our measurement in context with previous diskintegrated quantities obtained over the mm to cm wavelength regime. Overall, our result agrees well with past work, including the SMA measurements by Gurwell & Moullet (ALMA memo #594, Butler 2012) collected in the 0.8-1.3 mm range, as well as an IRAM-PdBI observation at 1.3 mm by Moreno (2007). A separate measurement by Ulich et al. (1984) at 1.3 mm is notably about ∼35 K higher than neighboring results. As seen in Fig. 4, there is an overall decrease in brightness temperature with increasing wavelength, with temperatures centered around ∼150 K at 10 µm that drop to ∼90-100 K longward of 1 cm wavelengths. The trend of brightness temperatures cooling off at longer wavelengths can be attributed to several factors, including the fact that lower frequency observations probe deeper into the subsurface where the diurnal wave is attenuated. Additionally, higher thermal inertia (e.g., denser) materials are often found at depth, and can further bring down daytime brightness temperatures. A similar trend of decreasing brightness temperature with wavelength is observed for the other Galilean satellites (see Fig. 4

Global Properties: Single Thermal Inertia
We initially fit our ALMA image of Callisto's leading hemisphere using a model parameterized by a single and a best-fit emissivity, and found the global properties were poorly constrained using this approach. The single-models that best-fit the observed distribution is given by = 600-1800, with emissivities of 0.91-0.93. If it is the case that Callisto's regolith consists of higher-material, then constraining precise values from the model is difficult because the ALMA data probe materials buried at depths that experience low diurnal wave fluctuations; additionally, increasing produces progressively flatter temperature profiles that are difficult to distinguish between. In addition to varying , we tested other model parameters, including the porosity, treatment of refraction, and ice to dust ratio; none of these changes provided a significantly better fit. Mainly, the systematic residuals and poor fits, along with the previous evidence (Spencer 1987) for a thermally heterogeneous surface, motivated our two-component analysis.

Global Properties: Two Thermal Inertias
Given the difficulties of fitting the observed distribution with a single , we now present the thermal properties acquired by modeling the data with two components. The model parameters that satisfied a cutoff of 2× the minimum chi-squared ( 2 ) allow = 15-50 and ℎ ℎ = 1200-2000, with best-fit global spectral emissivities of 0.89-0.91 (Fig. 5) The acceptable mixing ratios favored the ℎ ℎ component, with ranges varying from 70-50%, compared to 30-50% for . An example of how the residuals change with varying the percentage of / ℎ ℎ is shown in Fig. 2. While the endmember models have deviations from the data up to around ±6 K, the maximum deviations for the best-fitting two-component models are ∼2-3 K. Fig. 2 also demonstrates how the systematic inability of the single-component models to fit the limbs, especially the morning side, is remediated by the two-system.
Context for Callisto's thermal properties derived at millimeter wavelengths is provided by previous work at infrared wavelengths. Using the IRIS instrument onboard Voyager, Spencer (1987) found Callisto's infrared emission was best-fit by a model that adopted two, vertically stratified thermal inertia components: an upper = 15 +2 −2 layer that controlled the initial fast cooling and a lower = 300 +200 −200 that maintained the slower cooling later during the eclipse. As such, our result that a simple, onemodel is not sufficient to match Callisto's thermal millimeter emission is consistent with other work. Similarly, our best-fit overlaps with the upper layer inferred by Spencer (1987), although there is a larger discrepancy with our inferred second component ( ℎ ℎ = 1200-2000). However, comparing the thermal properties derived from IRIS data and 345 GHz ALMA data is not straightforward, given that these wavelengths sense the upper mm of a surface and depths of ∼5 cm into the subsurface, respectively. Importantly, the thermophysical model treatment by Spencer (1987) assumes that the two thermal inertia components are vertically stratified, while ours is implemented with horizontal segregation. Spencer (1987) did perform a simple test of horizontally stratified two-models, but used values that were best-fits for Ganymede and not best-fits for Callisto. Future work using individual ALMA observations at different observing frequencies will provide a more robust evaluation of vertical heterogeneity at millimeter wavelengths.
Additional context for Callisto's thermal properties derived at millimeter wavelengths is provided by similar ALMA observations of the other Galilean satellites. The mapping of Ganymede by de Kleer et al. (2021a) and Europa by Trumbo et al. (2018) Figure 5. Summary of best-fit two-models. Each cell represents the linear mixture results for a combination of (vertical axis) and ℎ ℎ (horizontal axis). The text imposed on each cell gives the range of (shortened to on the plot) that satisfied the 2 cutoff. The actual color of each cell corresponds to the minimum 2 achieved within a given mixture range. E.g., the cell corresponding to the mixture of the = 20 and ℎ ℎ = 1600 models is labeled with " = 30-50%", meaning that tested such two-models to estimate thermal inertia for observations of Ganymede at separate wavelengths, but such models were not preferred to single regimes, which differs from the the proposed explanation for Callisto's distribution. For emissivities, the measured global values for Callisto of 0.89-0.91 are higher than the 0.75 derived for Europa at 223 GHz (Trumbo et al. 2018) and 0.75-0.78 measured for Ganymede at 97.5-343.5 GHz (de Kleer et al. 2021a). This result is reasonable given that the emissivity of ice is lower than rock at millimeter wavelengths and Callisto's surface is not as abundantly ice-rich as its Galilean siblings. Moreover, the spatially-resolved ALMA images of the Galilean moons reveal morphological differences suggestive of their variations in thermal properties (Fig. 6). Among the icy Galileans, Callisto's thermal emission is more homogeneous across varying longitudes than that of either Europa or Ganymede, which is consistent with range of inferred thermal inertias inferred for Callisto reaching higher values.

Local Thermal Residuals
In addition to deriving global surface properties from our thermophysical model, we also highlight terrain-specific thermal anomalies that are robust to model parameter changes (Fig. 7). Note: references to west/east when describing the relative location or orientation of thermal residuals refers to geographic west/east on Callisto, not astronomical west/east.  Figure 7. A side-by-side comparison of the localized thermal residuals on Callisto with an albedo map projected to the same observer sublongitude and sub-latitude. The latitude/longitude grid spacing is 30 • . The ellipse in the lower left corner of each panel represents the ALMA synthesized beam at the time of observing. Left panel: projected albedo map of Callisto with major geologic features identified. Contrast added for ease of identification. Middle panel: Representative residuals from the range of best-fit two thermal inertia models. Right panel: Residual contour lines drawn on the projected albedo map to better illustrate the correlation between certain warm/cold regions and local terrain. Contours are drawn for ± 3, ± 2, and ± 1 K.

Valhalla Impact Basin
A noticeably cool region observed in both the calibrated image (Fig. 1) and the model residuals (Fig. 7) is well centered on the Valhalla impact basin. This thermal residual is colder by ∼8 K than surrounding terrain in the calibrated image, and about by ∼2-3 K compared to model predictions (which incorporate albedo), meaning the SNR of this residual is ∼10 with the 0.2 K measured image rms. Given that the spatial footprint of our ALMA beam at the distance of Callisto is ∼790 km, this thermal residual (which is larger than the beam) covers the entirety of Valhalla's high-albedo ∼ 360 km central zone and part of the ∼ 1900 km inner ridge and trough zone (Moore et al. 2004). The strength of the cold residual tapers off in outer trough terrain of Valhalla, which extends to ∼ 3800 km (Moore et al. 2004). In terms of morphology, this thermal feature may be slightly extended to the southwest, which is a different direction than the ALMA beam elongation and therefore not an artifact of the non-circular resolution element.
A correlation between cold thermal features and impact craters is a phenomenon observed on both Jovian and Saturnian icy satellites. Often, impacts excavate higher thermal inertia materials from the subsurface (e.g., ice-rich, higher density), creating localized terrain that maintains cooler daytime temperatures relative to surface materials that are more responsive to diurnally varying insolation. This is consistent with the thermophysical model residuals at the location of Valhalla improving as higher global thermal inertias were tested, as demonstrated in Fig. 2. Given Valhalla is one of the largest impact features in the Solar System, it is expected that this geologic unit contains relatively more exposed subsurface materials. Impact features on other Galilean moons observed to be thermally-cold via disk-resolved observations include Tros and Osiris on Ganymede (de Kleer et al. 2021a;Brown et al. 2022), and Pwyll on Europa (Trumbo et al. 2017(Trumbo et al. , 2018. Anomalously cold craters have also been observed in the Saturnian system on Titan (Janssen et al. 2016) and Rhea (Bonnefoy et al. 2020). With the results for Valhalla presented here, this general trend is confirmed to extend to the largest class of impacts in the Solar System.

Adlinda, Heimdall, and Lofn
Valhalla is not the only impact feature in our observation co-located with regional thermal anomalies -a collection of three large craters in the southern hemisphere may be the source of the cold residual located near 40 • S and 12 • W. This crater complex includes Adlinda and Heimdall, both of which are multi-ring impact features, and Lofn, which sits in-between/on top of the former two. Adlinda is marked by outer troughs that reach up to ∼ 850 km, however almost 1/3 of its structure is obscured by ejecta from Lofn (Moore et al. 2004). Lofn itself is an exceptionally bright, young anomalous dome crater and is Callisto's largest non multi-ring feature ( ∼ 355 km) (Moore et al. 2004). The morphology of this crater (e.g., unusually shallow relief) and the presence of water-rich deposits suggest the original Lofn impactor struck a relatively liquid/slushy region of Callisto's interior (Schenk 2002;Greeley et al. 2001Greeley et al. , 2000. Heimdall bears outer rings that extend up to ∼ 400 km and center geologic units with albedos high enough to saturate the Galileo images (Greeley et al. 2000). Because these individual craters are below the resolution of this observation, the measured properties are averages between these craters and the surrounding terrains. However, the region can still be seen to possess distinct thermal properties (Fig. 7).
As with the Valhalla cold region, the cold spot near the Adlinda/Lofn/Heimdall complex could be explained by impact-excavated higher-surface products. The residual itself diverges from model predictions by about 3-4 K (so SNR ∼15-20) and is roughly 1.5 beams across and is extended toward the northeast. The high albedos of Lofn and Heimdall, as well as spectroscopic evidence, point toward the presence of increased surface ice relative to the surrounding cratered plains, which could result in locally higher thermal inertia. Ligier et al. (2020) interpret their near-infrared maps of Callisto obtained using SINFONI on the VLT to suggest the small grains (25-100 µm) abundance is uniquely enhanced near Lofn and Heimdall (and at high latitudes in general). For an icy surface such as Callisto's, areas with small-grained ice particles can exhibit ice sintering, which increases grain-grain contacts and thereby increases the effective thermal inertia.
Unlike the Valhalla cold residual, however, this cold residual is not well-centered on the three-crater complex-rather, the temperature minimum is located roughly 12 • (or roughly 0.08 ′′ ) northeast of Lofn's bright central region. Although this region is closer to Callisto's limb, where the resolution tapers off due to geometric foreshortening, the fact that the residual is extended toward the northeast, e.g. not in the direction of beam elongation, suggests the offset may be real. Moreover, in the visible albedo map (Fig. 7), there is similar cratered terrain to the west of this residual that diverges from model predictions on order of the 0.2 K image rms, implying the thermal properties of the region corresponding to the cold residual are unique. A potential relationship to this residual may be found in the Galileo NIMS maps of CO 2 4.2 µm band depths by Hibbitts et al. (2000). The vicinity of Heimdall/Lofn bears locally deep CO 2 band depths, with moderate band depths also extending in a northeastern direction. However, because a link is difficult to confirm we only note the spatial correspondence of these features.

Additional Localized Thermal Anomalies
In addition to the thermal residuals on Callisto that appear connected to specific geologic units (e.g, Valhalla and the Adlinda/Heimdall/Lofn complex), there are other residuals that are not correlated to distinct geologic units. The first region is a warm, L-shaped residual that occupies the area west and south of Valhalla. The vertical component of the L is centered at almost 90 • W and is symmetric across the equator with a divergence from the model of ∼2-3 K, while the horizontal portion of the L is similarly warm, and extends eastward to ∼25 • W. The reason for the L shape may be because the warm region is cropped by the cold Valhalla impact to the northeast. The second region is the colder area that occupies most of Callisto's disk east of Valhalla, starting at ∼25 • W and extending toward the limb. According to the geologic map of Callisto (Greeley et al. 2000), both regions are classified as the same cratered plains terrain, despite exhibiting both cold and warm residuals. This is unlike the other significant thermal residuals shown in Fig. 7, where the units (i.e. large impacts) are linked to only one type of thermal signature (i.e., cold).
Given the morphology and location of these residuals, one speculation is that the ALMA data are sensitive to surface texture alterations induced by the Jovian particle environment. In the Jovian and Saturnian systems, the leading hemispheres of the icy satellites are preferentially bombarded by micrometeorites and neutral particles (such as dust), with the satellites' equatorial to mid-latitude regions most prominently affected, and their trailing hemispheres are preferentially impacted by charged species. Examples of these phenomena include the hemispheric color asymmetry on Europa (Carlson et al. 2009) and "Pacman"-shaped thermally-cold regions on Mimas and Tethys detected by Cassini (Howett et al. 2012(Howett et al. , 2011. Given that the warm L-shape on Callisto's leading hemisphere occupies mid-latitude regions near the satellite's apex of motion, it is possible the local regolith has been texturized to a less-compact material, perhaps due to micrometeorite bombardment since less-compacted regoliths generally have lower thermal inertias and grain contact areas due to lower bulk densities, resulting in faster daytime warming. The micrometeorite bombardment hypothesis is consistent with explanations for similar distributions on Ganymede's leading hemisphere detected with ALMA (de Kleer et al. 2021a). Moreover, as summarized by Moore et al. (2004), solar phase curves produced by Buratti (1991) also suggest Callisto's leading side is less compact than on the trailing side, although Domingue & Verbiscer (1997) using opposition effect data from Thompson & Lockwood (1992) infer the opposite result. Although beyond the scope of the current work, a trailing hemisphere image of Callisto at millimeter wavelengths could confirm if the warm region is actually confined to the leading hemisphere. Even if future observations confirm the equatorial warm region is leading hemisphere-specific, it may originate from several, interacting external modification processes.

CONCLUSION
We present a leading hemisphere thermal image of Callisto obtained at millimeter wavelengths (0.87 mm/343 GHz). Our first approach at modeling the thermal emission of Callisto used a thermal model parameterized by a single thermal inertia, for which we obtained a global best-fit of = 600-1800 (J m −2 K −1 s −1∕2 ) and best-fit emissivities of 0.91-0.93. The single-approach yielded high and systematic residuals that suggested a two-model may describe Callisto's surface more accurately. Our global best-fit two-model adopted = 15-50 and ℎ ℎ = 1200-2000, with acceptable mixing ratios of = 30-50% and ℎ ℎ = 70-50%, and emissivities of 0.89-0.91. Our result that a single-model is not sufficient to accurately model Callisto's surface is consistent with results from Spencer (1987) derived using IRIS infrared data. The ℎ ℎ we derived is higher than properties inferred from the infrared, and the presence of a low-thermal inertia component that may suggest the presence of a very uncompact surface component is consistent with other optical results (e.g., Buratti 1991). Callisto's thermal inertia components as inferred from these ALMA data are higher than those on Ganymede ( ∼ 450-750; de Kleer et al. 2021a) and Europa ( ∼ 95; Trumbo et al. 2017Trumbo et al. , 2018, which agrees with its more thermally-uniform appearance (Fig. 6).
In addition to deriving a global best-fit to the distribution, we also examined localized residuals to determine correlations, if any, with regional geologic features. We found that the Valhalla impact basin was co-located with a cold temperature anomaly ∼3 K, which extends the association between large impact craters/colder temperatures previously observed on other moons up to the multi-ring impact class. In addition, we found that the other primary cold region in our map is in the vicinity of a collection of craters including Adlinda, Heimdall, and Lofn, of which the latter two are among the brightest and youngest large features on the leading hemisphere. However, the cold residual is notably extended away from this suite of craters toward the northeast, which may imply surface textures related to factors other than cratering may be relevant. One residual for which there is no obvious association with a geologic unit is a warm L-shaped region that is nearly centered on the leading hemisphere and is constrained to about ± 30 • latitude. We speculated that this feature, like similar warm residuals observed on Ganymede's leading hemisphere (de Kleer et al. 2021a), might represent a region of lower thermal inertia resulting from micrometeorite bombardment that preferentially texturizes the leading hemisphere. However, such a hypothesis is difficult to confirm in the context of a single thermal observation.
The data presented here provide the first spatially-resolved image of Callisto's surface at ALMA observing frequencies. Future analysis of observations that include imaging of Callisto's trailing hemisphere at 0.87 mm/343 GHz, as well as mapping at other thermally-relevant wavelengths (e.g., 100-343 GHz), will help build the global and vertical-depth coverage needed to better understand the thermal features presented here. Additionally, this work offers ground-based observational context for the future European Space Agency (ESA) JUICE (JUpiter Icy Moons Explorer) mission to the Jovian satellites. The JUICE science program highlights Callisto as a key element in crafting the narrative of Jovian satellite formation and the emergence of potentially habitable environments and ocean-bearing worlds (ESA 2014). ALMA imaging of Callisto in the ∼100-343 GHz range will provide a near-subsurface (∼cm-m depths) complement to the JUICE instrument payload, which includes ground-penetrating radar (RIME, ∼9 MHz, depths of ∼9 km) and the Sub-millimeter Wave Instrument (SWI, ∼530-601 GHz).