JIRAM (Adriani et al., 2008, 2016, 2017) combines an imager and a spectrometer into one unit. The spectrometer data is not used in this work, but was discussed, together with the instrument details, in Mura et al. (2020). A dedicated de-spinning mirror compensates for the spacecraft rotation (at approximately 2 rotations per minute). However, due to a misfunctioning of the electronic board, starting from orbit 44 it was decided to stop using the mirror and to operate JIRAM without this feature. Since the exposure time for high-radiance targets, such as the hotspots, is very short (∼few milliseconds), the smearing of the images is small and is corrected via software. During each spacecraft rotation, JIRAM performs an acquisition that includes two images, one in the L (3.3–3.6 µm) and one in the M (4.5–5 µm) spectral ranges. The field of view (FoV) of JIRAM can be tilted within the plane perpendicular to the Juno spin axis by adjusting the timing of the acquisition. JIRAM is not capable of adjusting its FoV in any other direction without requiring a spacecraft reorientation; on some occasions, such reorientation is actually performed to allow JIRAM to observe Io.

The imager channel consists of a single detector with two filters, separated by a 10-pixel inactive strip. The pixel angular resolution (IFOV) is 237.8 µrad, resulting in a FoV of 5.87° by 1.74° for both the L and M bands. The responsivity of both imager channels is 2 × 10⁶ DN/(W m^-2 sr⁻¹ s⁻¹).

The L-band and M-band filters were originally designed to detect, respectively, the auroral and atmospheric emission from Jupiter. When observing Io, the L filter is more sensitive to reflected sunlight, whereas the M filter is more sensitive to thermal emission. However, the L filter is less utilized due to the spectrometer slit receiving no signal if the target is entirely within the L frame. So far (that is, until the last orbit whose data is used for this study, i.e., orbit 49), more than 2000 images of Io have been acquired in the M band, and almost 200 in the L band.

The exposure time of JIRAM for observing Io has been varied over the years, starting with long times and gradually decreasing to shorter times. The reason for this choice was to prioritize the discovery of new hotspots, being aware that this way, the brightest hotspots would lead to signal saturation. JIRAM saturation occurs gradually and, as reported in Adriani et al. (2017), a Digital Number (DN) value below 10,000 ensures linearity of the detector. Mura et al. (2020) suggest not exceeding 12,000 DN to avoid saturation. In most images prior to orbit 37 (Oct. 2021), saturation is very high, and these images should be used only for characterization of hotspot location as in Zambon et al. (2023), while the total radiance of hotspots can be heavily underestimated (in Supplementary Table S1 we give an estimation of the saturation effect for the whole JIRAM dataset until the orbit 49 of Juno). For this reason, we use the best part of the JIRAM database, specifically the latest orbits. In addition to having better spatial resolution, these later orbits exhibit almost negligible saturation signal fractions. Images of Io captured before orbit 5 were solely taken for the purpose of conducting instrumental calibrations and should not be regarded for scientific use. Table 1 provides a summary of the observations utilized in this study. In Figure 1 we show a composition of the L and M images of Io taken during the 43rd orbit of Juno.

To enhance the signal-to-noise ratio and implement a super-resolution algorithm, the images captured both in the M- and L-band were stacked and superimposed. The method is the same used in Mura et al. (2020); here we just report the number of images of each observation (Table 1). The removal of reflected sunlight is performed in a subsequent step (see Section 3.1).

The geometry of JIRAM observations can be reconstructed using routines and data provided by SPICE/NAIF (Acton, 1996). Generally, an accuracy comparable to the angular size of the JIRAM pixel is sufficient (for both the imager and spectrometer, which is approximately 237.8 µrad). However, in some cases, a significant deviation between the images and the reconstructed geometric information is found, suggesting the need for further correction. For example, this can occur in observations of Jupiter’s moons performed after a spacecraft reorientation, resulting in a slight wobbling phenomenon. Regardless of the reason, not investigated here, it is evident that a correction would be beneficial for the images used in this study.

A processing method similar to that used for analyzing the images in Mura et al. (2020) and Zambon et al. (2023) has been implemented. However, to avoid arbitrariness in interpreting the limb position, two substantial improvements have been introduced here. On Io’s dayside, especially in equatorial regions, Io can be quite bright, such that even one or two pixels outside the limb can have a significant radiance level, and radiance does not drop to zero sharply between pixels inside and outside the limb. This phenomenon gradually decreases toward the poles. Limb fitting on images without correcting for this effect would systematically shift the reconstructed limb toward the daylight direction. Therefore, an approximate and preliminary reconstruction of geometries was first carried out by using SPICE geometries only, and the image was corrected by applying a Lambertian lighting model; after this correction, the entire limb of the daylight side has an approximately constant radiance level.

The second modification involves introducing a routine that, to avoid the arbitrariness of the human operator’s choice of limb position, identifies it knowing that the dayside limb must be (after the correction mentioned in point 1) a semicircle with a known radius. As long as Io’s dimensions in pixels in JIRAM images do not exceed about 100 pixels, Io can be considered a sphere (for images of observations after those studied here, where Io is seen from closer, such as those in Juno’s orbits 51, 53, 57, and 58, an ellipsoid shape must be used). The routine finds the contour of constant radiance corresponding to a semicircle of the right radius and then determines the position of Io’s center. The improvement over limb fitting done by a human operator can be estimated on the order of one pixel, which is not much, but the method also has the advantage of being deterministic and repeatable.

Finally, as already mentioned, the images of each orbit are grouped and stacked. Limb fitting is performed on the already averaged super-image. Consequently, for orbits 41–49, we have 4 maps in the M band and 4 maps in the L band (as shown in Figure 2, see Section 3.1). The covered regions are very similar, and as a final check, it has been verified that the position of hotspots is the same in all 8 maps (see Supplementary Figure S4 for one example). Note that the coverage of the southern hemisphere is quite poor.

The radiation from the dayside of Io consists of both reflected sunlight and blackbody emission. Far from the hotspots, the equilibrium temperature expected for the surface of Io in its subsolar point within the thermal skin depth (of the order of a few mm) is −130 K (Rathbun et al., 2004), well below the detection threshold of JIRAM, so this component can be neglected. The reflected sunlight, on the other hand, must be considered as it is significantly bright at JIRAM wavelengths. Io is covered with a variety of different materials, including sulfur, sulfur dioxide frost, and various types of silicates, resulting in a diverse and complex surface. We found that an empirical model, similar to the Oren-Nayar one (Oren and Nayar, 1994), with respect to a simple Lambert model, provides a reliable way to describe the way light scatters off the dayside surface of Io. In this first version, spatial variability of albedo is not included in the model, i.e., albedo is assumed uniform (but different at different wavelengths). In the context of this model, a body surface roughness is represented by a statistical distribution of microfacets. These microfacets contribute to the scattered light in different directions, resulting in a diffuse reflection. The model, in fact, incorporates the roughness σ as a parameter (defined below). Furthermore, this model takes into account the incident light direction and the surface normal to calculate the intensity of reflected light. This accounts for the non-uniformity of Io’s surface and allows for more accurate rendering of its appearance under different lighting conditions.

Applying such a model to Io’s surface enables us to simulate the complex interactions of light with the various materials present. This allows generation of realistic values of reflected sunlight, that can be subtracted from the radiances measured by JIRAM to obtain the radiance due to Io’s hotspots only. Eqs 1–6 describe the model: θ _i and θ _r are the incident and reflection zenith angles (i.e., solar and emission angles), α and β are convenient angles defined by equations and 2: α = ⁡ max θ i , θ r (1) β = ⁡ min θ i , θ r (2)

The roughness of the surface is given by the parameter σ, which is the standard deviation of the gradient of the surface elevation and ranges from 0 to ∞. For σ = 0, all facets lie in the same plane and the body is a perfect sphere. The horizontal projections of the incident and reflected beam (azimuthal angles) form an angle φ (φ _i - φ _r ). E ₀ is the irradiance at the sub-solar point, and ε is a parameter to take into account the albedo. Then, three parameters C1, C2 and C3 are defined in the following way: C 1 = 1 − 0.5   σ 2 σ 2 + 0.33 (3) C 2 = 0.45 σ 2 σ 2 + 0.09   sin ⁡ α , cos ⁡ φ > = 0 0.45 σ 2 σ 2 + 0.09   sin ⁡ α − 2 β π 3 , cos ⁡ φ < 0 (4) C 3 = 0.125 σ 2 σ 2 + 0.09   4   α   β π 2 2 (5)

Now it is possible to define the two radiances, L1 (accounting for direct reflection) and L2 (accounting for multiple reflection of light between facets): L 1 = E 0 cos ⁡ θ i e   ε π   C 1 + C 2   cos ⁡ φ   tan ⁡ β + C 3   1 − cos ⁡ φ   tan α + β 2 (6) L 2 = E 0   cos ⁡ θ i e   0.17   ε 2 π   σ 2 σ 2 + 0.13   1 − ⁡ cos ⁡ φ   2   β π 2 (7)

Finally, the radiance L of the light reflected by the surface is just: L = L 1 + L 2 (8)

Note that we include a free parameter e, not present in the original version of the model, to allow a better fit of the data to the model. The model parameters are constrained using data taken from orbits 41, 43, 47 and 49 in both L and M bands. All hotspots have been removed from images to fit only the contribution from sunlight. In Supplementary Figures S1, S2 we show the data before the fit is performed; in Supplementary Figure S3 we show the model results, and the images after subtraction in one of the four cases. The parameters obtained from the fits are: ε = 0.19, σ = 0.88 and e = 0.90 for the L band; and ε = 0.34, σ = 1.5, e = 0.80 for the M band. Note that because we introduced an extra parameter e, the parameters ε and σ have values that do not reflect exactly the original meaning. In particular, the albedo or the Bond albedo should be calculated by using the formula. The result of this removal procedure, for orbits 41, 43, 47 and 49, is shown in Figure 2.

Some hotspots are localized within depressions (paterae), a situation in which the crater edges can effectively shield a portion of the emitted radiation as the emission angle becomes large. It is thus conceivable that Lambert’s cosine law may be an inaccurate approximation for the radiance variability of hotspots. In our dataset, we have a sample of four images taken at closely spaced time intervals (from orbit 37), during which the Juno probe has shifted its position slightly, albeit still significantly, relative to Io, presenting each hotspot under a different emission angle. Given the short duration of the overall observation (15 min), it is fully reasonable to assume that the observed radiance variability for each hotspot is solely attributable to changing observation conditions. Although no hotspot changes its emission angle over a very large range, the ensemble of hotspots covers an overall emission angle range practically from zero to 90°. Once the reflected solar component is removed, it is therefore quite straightforward to perform a best-fit of the radiance variability law as a function of the emission angle θ _r , which is assumed to be cos ^α (θ _r ), where the exponent α is a parameter that is presumably greater than one, due to the aforementioned concept (radiation emitted from depressed areas, and radiance that decreases faster than cos ^α as α increases). Indeed, as shown in Figure 3, where the measured radiances are plotted on the x-axis and the estimated radiances on the y-axis, the law reproduces the variability of hotspots quite well, and the value of α obtained from the best fit (1.25), is indicative, as stated earlier, of a prevalence of emission from slightly depressed surfaces compared to the mean ground level. This concept has been proposed by de Kleer and de Pater (2016), who noted that correcting the radiance with cos ¹ (θ _r ) may be appropriate for some hot spots but not for others (see their Figure 12). It is important to have an accurate law for the emission angle, as in the continuation of this study, we will seek to assess the temporal variability of hotspots, but these are observed in situations of slightly different emission angles (orbits 41, 43, 47, and 49 have similar vantage points, but not exactly the same), and therefore precise correction is necessary. Because this effect is attributed to the shadowing of the hotter surface by colder patera walls, we can safely assume that the same function is valid for the L band as well.

Io’s volcanic centers erupt lavas that are mafic in composition, based on their color and flow morphologies, and they have temperatures consistent with basaltic, or possibly ultramafic, compositions (e.g., Stansberry et al., 1997; Williams et al., 2001). Early analysis of Galileo data (McEwen et al., 1998; Davies et al., 2001) was suggestive of ultramafic lavas, which erupt at higher temperatures than basaltic lavas, up to −1800 K for the Pillan hotspot. However, subsequent analysis of the same data by Keszthelyi et al. (2007) indicates a lower limit of −1600 K, consistent with basaltic to lunar-like compositions. Mura et al. (2020), using JIRAM data, report a large range of temperature values that can fit the observed spectra (from 250 to more than 1,000 K). Some of these temperature values can be the effect of non-uniform filled pixel, that is they are somewhat an average of the hot and cold regions inside a pixel (where the hot regions are usually smaller but emitting more radiant power in comparison with colder regions). Note that the size of the pixel of the spectrometer of JIRAM is exactly equal to the pixel of the imager. Hence, we safely assume, or at least not exclude, that the same variability of detected temperature will be valid for the imager data. This is crucial when estimating the total power for a hotspot. A method to infer the total power from the M-band spectral radiance has been proposed and used in a recent paper (Davies et al., 2024). We note that, because of the variability of the detected temperature, this method may have large uncertainty. In Supplementary Figure S5 we show the ratio between the total power emitted from a black body, and the radiance at the M band center (for a grey body, with emissivity smaller than one, the ratio does not change). In the next Section, from the L/M ratio we will show that this range is completely covered by the data of our observation. A similar variability of temperatures can be found in Mura et al. (2020). In conclusion, estimating the total power only from the M band may lead to uncertainties up to one order of magnitude (and such uncertainty does not have a linear dependence with temperature); hence, we prefer to implement a 2-band total power estimation.

The estimation of the hotspot temperature values on Io through the analysis of JIRAM spectra was first obtained by Mura et al. (2020). However, the data selected for our current study consists solely of dual-band images, because JIRAM spectra are not ideal for this study. In fact, as Juno gets closer to Io, Io gets larger (and faster) inside JIRAM’s field of view, and this implies that the spatial coverage of the spectrometer is reduced (after orbit 43, a malfunctioning of the de-spinning mirror prevented the use of the spectrometer anyway). However, it is possible to determine the temperature of a grey body source using only the radiance values in the L and M bands. As depicted in Fig. 6, the ratio between L and M band radiance has an almost linear relationship with temperature. This simulation employs the actual band ranges for L and M filters to integrate the Planck function at various temperatures. Notably, some uncertainties may exist in this estimation, such as potential variations in the emissivity at 3.45 and 4.8 µm. The provided value serves as an indicative average temperature within the image pixel, resembling an “effective temperature” as outlined in Eq. 1 of Lopes and Spencer (2007). It is crucial to acknowledge that this assumes uniform emissivity across the ∼3–∼5 µm range (i.e., a grey body in this spectral range), which might not be the case. Consequently, the temperature derived in this manner can be less precise than that derived from the full spectrum. Nevertheless, this estimation surpasses the accuracy of using a single band, which is influenced by factors such as the fill factor, indicating the partial occupancy of a pixel by the source, and uncertainties in emissivity magnitude.

This method requires that each hotspot has a radiance spectrum that is indicative of one temperature only: as a matter of fact, most spectra in Figure 3A of Mura et al. (2020) do not show bimodal distributions. However, to evaluate the significance of the temperature that is estimated in this way, one can consider the background temperature and the imagers’ minimum detectable signal, which at 1 s of exposure time is 5 10⁻⁶ W sr⁻¹m⁻² and increases as the exposure time decreases. Such value, in the M band and for usual exposure times, corresponds to about 150 K; in practice, JIRAM would be able to measure the surface temperature in regions away from hotspots only if it could observe areas completely free of hotspots with very long exposure times and with high resolution. However, this did not occur during the observational campaign. In practice, any temperature above 500 K and any filling factor greater than 1 in ten thousand cause the radiance in the 2–5 micron range, and therefore especially in the individual L and M bands, not to be influenced by anything other than the hottest source. Consequently, our temperature retrieval is possible. Once the temperature retrieval is performed, it is possible to calculate the filling factor that provides the correct values for the two radiances, L and M, and verify that the found value is within an acceptable range (i.e., not fitting the temperature of a too small portion of the pixel).

After the temperature of a hotspot is estimated (i.e., the temperature of a grey body that best fit the radiance in the two bands) it is trivial to compute the total power by integrating the Planck function (over the wavelength and solid angle) and multiplying by the hotspot area (see next Section for details).