NASA's Juno spacecraft, launched in 2011, has been monitoring Jupiter since its orbital insertion on 4 July 2016. Among its nine instruments is a six‐frequency microwave radiometer (MWR) designed to measure Jupiter's microwave emission emanating from depths ranging from a few tenths of a bar down to several kilobars. As a passive remote sensor that can avoid much of the interference of synchrotron radiation that stems from Jupiter's radiation belts and confounds earth‐based observations at longer wavelengths, MWR can infer details of the Jovian atmosphere, namely, the distribution of atmospheric constituents that either have significant opacity in the range of frequencies spanned by MWR (ammonia, primarily) or affect the lapse rate (water). The ultimate objective of this instrument is to determine a global water (and therefore oxygen) abundance to sufficient precision to discriminate among various competing theories of Jupiter's origins and, by extension, those of our solar system (Atreya et al., 2019; Gautier et al., 2008; Helled & Lunine, 2014; Li et al., 2020; Mousis et al., 2012; Owen et al., 1999; Wong et al., 2008).
A key step toward achieving this goal is the determination of the angular and positional dependence of Jupiter's thermal emission. Retrieval of global water abundance, in particular, relies heavily on an accurate description of limb darkening (Janssen et al., 2005). The determination of these quantities is achieved through a deconvolution of the measured antenna temperatures. This work describes (1) the formulation of the deconvolution, (2) validation of the resulting nadir brightness temperatures and limb darkening which are subsequently used in atmospheric retrievals, and (3) selected results from MWR‐favorable orbits. Section 2 begins with a formulation of the deconvolution as a linear problem mapping a set of coefficients describing Jupiter's thermal emission to a set of measured antenna temperatures and validates the choice of basis set used to represent the linear mapping with an atmospheric model. This section also presents a heuristic algorithm to mitigate the effects of nonthermal sources through the use of a synchrotron model and subsequent downselection of MWR measurements to a subset useful for the deconvolution. Section 3 describes the approach used to solve the linear problem and validate its solution, and section 4 shows results from all MWR‐favorable orbits between July 2016 and April 2018, covering latitudes between 60°S and 60°N.
The Juno MWR instrument (Janssen et al., 2017) comprises six receivers whose center frequencies are approximately equally spaced in log‐frequency over a range of 0.6 GHz (50 cm, channel 1) to 21.9 GHz (1.37 cm, channel 6). Each receiver measures the power received through the antenna in a narrow (∼3.5%) bandwidth centered on its nominal frequency and averaged over contiguous time intervals of 100 ms. In the microwave region it is customary to scale the directionally dependent radiant intensity to be the brightness temperature in units of Kelvin through the use of the Rayleigh‐Jeans approximation (Janssen, 1993). In this case the power collected by each receiver antenna can be expressed as an antenna temperature: [Image Omitted. See PDF]where G(ϑ, φ) is the antenna gain function, normalized to unity over an integration over the two‐dimensional unit sphere, and is the sky brightness distribution. (The gain function G(ϑ, φ) defined in this way is a factor of 4π smaller than the traditional definition of antenna gain and is proportional to the beam pattern, whose maximum value is normalized to unity.) The antenna gain patterns are approximately Gaussian with beamwidths of approximately 21° (channels 1 and 2) and 11 to 12° (channels 3–6). Measurements of G(ϑ, φ), represented on a 1° by 1° grid, were made pre‐launch and are described in more detail in Janssen et al. (2017). The noise for each measurement varies from about 1.0 K (channel 1) to 0.2 K (channel 6), with an independent systematic error in absolute calibration of approximately 2% for each channel (Janssen et al., 2017). The spacecraft spins at a fixed rate of approximately 2 rpm, and Jupiter and its environs are continuously scanned as the Juno spacecraft travels from north to south through periapsis. The mapping obtained on Jupiter depends on orbit and spacecraft orientation as discussed below.
Three sources contribute to the antenna temperature: observed Jovian thermal emission, T(J ) (the term of interest), the galactic background, T(g), and a term, T(s), which encapsulates other effects such as synchrotron emission, aurorae, or lightning that are not accurately modeled. [Image Omitted. See PDF]The principal objective of the deconvolution process is to determine brightness temperatures T(B)(r, μ) due to Jovian black body emission as functions of positions, r, on Jupiter and the cosine of the emission angle, μ. The first step, then, is to separate as well as possible the three contributors to the measured antenna temperatures. The convolved galactic background was measured during the mission's cruise phase, prior to Jupiter orbital insertion (JOI) and is assumed to be accurately modeled. The contribution stemming from synchrotron emission is more problematic. Discussion of how to handle this term is deferred to section 2.5. For the time being, we shall assume that both synchrotron emission and the galactic background can be accounted for and subtracted appropriately from the measured data, so that the thermal contribution to the antenna temperature, T (J), can be isolated. Analogous to Equation 1, the quantity T (J) may be expressed as a convolution of the gain function and the spatially dependent brightness temperature T (B) arising only from Jupiter's atmospheric thermal emission: [Image Omitted. See PDF]where r is the intercept point of the direction vector (ϑ, φ) with Jupiter's one‐bar equipotential surface and μ is the cosine of the angle formed by the direction vector and the surface normal at this intersection. The intercept points are determined using the SPICE library (Acton, 1996; Acton et al., 2018) along with ancillary data appropriate for the Juno mission. For a set of measurements indexed by i, this convolution is discretized into a summation over a set of direction vectors, indexed by j, in a rotating reference frame in which the antenna boresight is fixed: [Image Omitted. See PDF]where rij denotes the intercept point on Jupiter in direction j, μij denotes the cosine of the emission angle, and δAj denotes the area of the 1° by 1° patches comprising the angular integration over the unit sphere. In this work deconvolutions are performed on a per‐orbit basis, so that each measurement set consists of all the measurements within a single orbit that survive a screening process (described in section 2.5). The brightness temperatures are then expanded over sets of positional basis functions {hp(r)} indexed by planetary grid point p and angular basis functions {fk(μ)} indexed by k: [Image Omitted. See PDF]where ckp are the coefficients we would like to solve for. The planetary grid points are generated using a Healpix (Górski et al., 2005) spherical tessellation comprising 49,152 equal‐area pixels, where is the Healpix resolution parameter, which must be a power of two. The pixel centers lie on 255 unique latitudes, and the spacing between latitudes coarsens from 0.6° at the equator to 0.8° toward the poles. The positional basis consists of bilinear interpolation functions, each localized to a unique pixel such that it attains a maximum value of unity at the grid point associated with the pixel and vanishes outside a region bounded by its neighboring grid points. The choice of angular basis functions is discussed in section 2.3.
The antenna temperature in Equation 4 can thus be expanded in terms of the coefficients ckp: [Image Omitted. See PDF]
The matrix element Mi, kp represents the contribution to the ith measured antenna temperature from the kth angular coefficients at planetary pixel p. Equation 6 can be more compactly expressed in matrix form as a linear relation: [Image Omitted. See PDF]where M is a linear mapping from a spatial and angular description of Jupiter's thermal emission to a set of measured intensities and depends only on the particulars of the instrument (e.g., gain patterns) and on its position and orientation for those data included in the measurement set, and we have explicitly added a term to reflect the fact that actual measurements contain instrument noise. An optimal estimation solution to Equation 7 can expressed as a minimization of the sum of two terms with respect to the coefficient vector c: [Image Omitted. See PDF]where cp is an a priori estimate for c. The quantity Sm represents the covariance of the measurement noise, and Sc, effectively a tuning parameter, is an estimate of the parameter covariance. The limiting case where implies no prior information and hence no regularization, while setting Sc to a tunable diagonal matrix and cp to zero is equivalent to Tikhonov regularization in its simplest form (Rodgers, 2000). The solution to Equation 8 is given by (Rodgers, 2000; Tarantola, 2005) [Image Omitted. See PDF]
The measurement covariance Sm is determined by both the characteristics of the antennae and the integration algorithm used to determine Juno MWR's antenna temperatures from raw counts generated from each receiver's voltage‐to‐frequency converter. A radiometer noise model has been developed (Janssen et al., 2017) to simulate the instrument's response to the antennae's received power, and a Monte Carlo simulation to determine the measurement covariance matrix was performed. However, it was found that the resource requirements to establish off‐diagonal components of Sm such that a matrix inversion of Sm was sufficiently accurate were prohibitive. In this work, we approximate Sm by ignoring its off‐diagonal components. For most channels, this should be a good approximation; the Monte Carlo simulations found the largest off‐diagonal elements to be at least a factor of 5 smaller than the diagonal for the first five channels, though off‐diagonals for the shortest wavelength antenna, channel 6, could be as large as half of the diagonal values. This diagonal approximation is tantamount to the assumption that the measurement error is dominated by white noise, so that errors associated with measurements made at different times are uncorrelated. In this case, the diagonal components of the measurement covariance (i.e., the variances) were found to be well‐approximated within the set of data used in the deconvolution by a quadratic dependence on measured intensity: [Image Omitted. See PDF]The coefficients αi in Equation 10 used in this analysis are provided in Table 1. Details of the fitting are provided in the supporting information.
TableLeast‐Squares Fit of Parameters Describing Measurement Variance in Equation 10 From Monte Carlo Simulation| Channel | Frequency (GHz) | α0 (K2) | α1 (K) | α2 |
| 1 | 0.5998 | 2.618 × 10−1 | −9.557 × 10−5 | 1.301 × 10−6 |
| 2 | 1.2476 | 1.249 × 10−1 | 2.603 × 10−4 | 5.392 × 10−7 |
| 3 | 2.5966 | 8.016 × 10−2 | 1.458 × 10−4 | 1.468 × 10−7 |
| 4 | 5.2153 | 5.773 × 10−2 | 1.207 × 10−4 | 2.652 × 10−8 |
| 5 | 10.004 | 1.933 × 10−2 | 5.146 × 10−5 | 2.566 × 10−10 |
| 6 | 21.898 | 2.725 × 10−2 | 8.775 × 10−5 | −1.563 × 10−7 |
Equation 6 describes the mapping of a parameterized form of the brightness temperature T(B)(r, μ) to a set of observations. This parameterization of the brightness temperature is represented at each relevant region on Jupiter by a linear combination of angular basis functions. Different basis sets have been used in the literature in the context of modeling limb darkening of stellar atmospheres. Typically, some combination of integer and half‐integer powers of the cosine of the emission angle is used (e.g., Claret, 2000; Heyrovský, 2007). For an ideal instrument in the limit of vanishing measurement error, increasing the number of basis functions can potentially yield a more accurate representation of the angular dependence of the thermal emission. In practice, however, instrument noise constrains the resolvability of this functional dependence and limits the improvement in fit one might expect from using larger bases. Furthermore, increasing the number of degrees of freedom describing the angular variation of the thermal emission beyond a minimal value increases the effective degeneracy of the problem and further ill‐conditions M as its columns become no longer linearly independent to within measurement noise. Therefore, a trade‐off exists: One should select the smallest basis set that is still capable of fitting observational data to within measurement error.
For guidance, we evaluate the suitability of angular basis sets by analyzing their performance on a simple, adiabatic description of Jupiter's atmosphere based on an equilibrium cloud condensation model (Weidenschilling & Lewis, 1973). Details of the model, the Jupiter Atmospheric Radiative Transfer Model (JAMRT) are described in Janssen et al. (2017). The solid curves in Figure 1 illustrate the angular dependence of equatorial brightness temperatures for Juno's six channels as computed by JAMRT assuming a moist adiabat. Lower boundary conditions include deep atmosphere volume mixing ratios of NH3 = 351 ppm and H2O = 2,500 ppm, the maximum likelihood solutions given in Li et al. (2020) for Jupiter's equatorial deep constituent abundances. These values correspond to 2.76 times protosolar NH3 and 2.7 times protosolar H2O (Asplund et al., 2009). Abundances of other constituents (He, CH4, H2S) are as cited in Atreya et al. (2019). The temperature is pinned to a reference value of 132.79 K at the 0.5‐bar level (Seiff et al., 1998). Although the radiative‐transfer component of JAMRT does account for refraction and non‐plane parallel transmission, relative uncertainties in these quantities magnify uncertainties in path length at near‐limb viewing geometries, so that relative errors in brightness temperatures as μ → 0 may be quite large. The computed brightness temperatures are fit over the range μ > 0.6 with a quadratic polynomial using a constant weighting function evenly spaced in angular space. These quadratic fits, denoted by dashed lines in Figure 1, by eye demonstrate good agreement with the modeled values over the range of emission angles over which they were fit but show significant deviations at larger angles ( ).
Enlarge this image.1 Figure. Equatorial brightness temperatures as functions of cosine of emission angle (μ) for Juno's six channels computed by JAMRT (solid). The model assumes a moist adiabat, deep atmosphere mixing ratios of NH3=351 ppm, H2O = 2,500 ppm, and an upper boundary condition of T(P=0.5bar)=132.79 K. Dashed curves are least‐squares quadratic fits using a weighting function that is evenly spaced in angular space and constant and nonzero over the interval 0.6 < μ < 1. A quadratic fit is appropriate for small emission angles but insufficient to cover the entire range of possible emission angles.
At an absolute minimum, at least two fitting parameters are needed, one to describe an intensity at a particular emission angle and another to describe variation with emission angle. It was found through experimentation with actual Juno MWR data that a third independent parameter describing some sort of curvature with respect to emission angle could also be consistently retrieved. However, various attempts to admit additional parameters not substantially correlated with the other three were unsuccessful owing to the limitations that instrument noise places on the resolution of the retrieval. It is clear from Figure 1, however, that for none of the six channels can a simple quadratic polynomial fit the entire span of admissible emission angles. Since the MWR data are described by a convolution with a beam pattern of significant width (up to 21° for the two lowest channels; Janssen et al., 2017), it is possible for emission from larger angles to make non‐negligible contributions to the measured brightness temperatures even with fairly restrictive sampling criteria.
To reduce the errors incurred from quadratic fits to TB(μ) over the full range of emission angles we introduce shape functions ξ(μ) for each channel given by the ratio of a modeled brightness temperature to its extrapolated fit: [Image Omitted. See PDF]where are brightness temperatures determined from a model atmosphere (ultimately chosen to be the one described in Figure 1) and is the least‐squares quadratic estimate of the model where the fit is performed only over the interval 0.6 < μ ≤ 1.0. For small emission angles (μ > 0.6), ξ(μ) is equal to unity (up to small fitting errors less than 1.5 × 10−4) but at larger emission angles accounts for deviations from a quadratic fit. For each channel and for each physical location that is modeled, the angular dependence of the thermal emission is fit to [Image Omitted. See PDF]where (corresponding to an emission angle of 37°) is a conveniently chosen reference value for the cosine of the emission angle. The specific form of this quadratic expression is chosen to provide a more intuitive meaning for the coefficients: c0 is the nadir brightness temperature; c1 is the absolute limb darkening when ; and c2 represents the additional reduction of the brightness at (or 53°) over that of a linear extrapolation. It follows, then, from the expression in Equation 12 that the angular basis functions are given by [Image Omitted. See PDF]
Since ξ(μ) depends on the details of the atmospheric temperature and mixing ratio profiles, which are unknown, it would not appear at first glance to be a useful quantity. However, we shall demonstrate that for the subset of measurements that are sampled for the deconvolution, the dependence of ξ(μ) on atmospheric properties is sufficiently weak that any reasonable profile may be selected to reduce bias. Figure 2 shows the range of values for each channel that ξ(μ) assumes over a reasonable spectrum of model atmospheric boundary conditions: deep ammonia and water volume mixing ratios span values of 305 to 395 ppm and 0 to 7,000 ppm, respectively; the reference temperature at 0.5 bar ranges from 130 K to 135.6 K; the lapse rate follows both dry and moist adiabats; and the local gravitational acceleration, which linearly affects the lapse rate, spans all latitudes from equator to pole. The bounds of deep ammonia and water abundances were informed by an estimation of the probability distribution of equatorial values of these quantities presented in Li et al. (2020). Despite the large variations in model atmospheres, ξ(μ) remains constrained to fairly narrow bands of values at each emission angle, especially for channels 2 through 6. The thickness of these bands serves as rough metric of the sensitivity of the thermal emission at high emission angles to atmospheric variation once emission at smaller angles is determined. Channel 6 shows almost no variation with atmospheric model, while channel 1 shows the greatest sensitivity among all the channels. Nevertheless, channel 1's variation is still less than 10% for all emission angles less than 81° (μ ≈ 0.15). For the remainder of this work, we employ the shape function corresponding to the model atmosphere described in Figure 1.
Enlarge this image.2 Figure. Extent of shape functions for all six channels using a wide range of model atmospheres: deep ammonia and water mixing ratios span values of 305 to 395 ppm and 0.0 to 7,000 ppm, respectively; the reference temperature at 0.5 bar ranges from 130 K to 135.6 K; the lapse rate follows both dry and moist adiabats; and the local gravitational acceleration spans all latitudes from equator to pole. Despite the atmospheric variation, shape profiles follow fairly narrow envelopes. Profiles selected for this work are represented by the dashed curves and correspond to the specific conditions indicated in Figure 1. (The range of shape profiles for channel 6 is much narrower than the others and nearly coincides with the selected profile for channel 5.)
To evaluate how well this specific choice of shape function can fit actual data in our synthetic tests we introduce a measurement covariance matrix which incorporates information about the emission angle distribution of the set of measurements included in the sampled MWR data. We define an angular contribution function, dE/dθ, whose integral over all emission angles equals the sum of all energy received by MWR over some set of measurements: [Image Omitted. See PDF]This quantity characterizes the angular emission distribution that contributes to the total received power over some measurement set. A more detailed discussion of how the MWR measurements are sampled is deferred to Appendix A, but the primary criterion that determines the sampling is that at least 99% of the gain pattern, G(ϑ, φ), intersects the planet. For the purposes of the analysis in this section, we also restrict the measurement set to those measurements whose footprint maxima lie within 25 degrees latitude of perijove in order to focus on the region where MWR's resolution is highest—the angular distribution varies with latitude, and measurements at high latitudes do not strongly affect the deconvolution of brightness temperatures near perijove. There have been two general categories of spacecraft orientation during Juno orbits, (1) MWR orbits in which the spin vector of the spacecraft is normal to the orbital plane so that smaller emission angles are favored and (2) gravity
orbits in which the high gain antenna is pointed toward the earth in order to perform gravity science. The spacecraft spin axes for early gravity orbits (before perijove 10) were also coincidentally favorable for MWR so that the two orbits share many characteristics of their viewing geometries. Though not a primary objective, data from portions of gravity orbits can be shown to be useful for MWR analysis. Figure 3 illustrates the angular contribution function for the first 12 perijoves. The contributions to antenna temperature from MWR orbits (solid) and also early gravity orbits (dashed) are weighted toward smaller emission angles (especially at low latitudes), while later gravity orbits (perijove 10 and later) show significant contributions from larger emission angles.
Enlarge this image.3 Figure. Distribution of emission angles (dE/dθ in Equation 14) contributing to each of the first 12 perijoves for both MWR (solid) and gravity (dashed) type orbits. Only measurements in which the footprint latitude (i.e., latitude of the intersection of the boresight vector with Jupiter's 1‐bar equipotential surface) lies within 25 degrees of perijove are included.
We begin with Equation 7, in which the rows of M correspond to emission angle “observations” evenly spaced in θ and modeled by JAMRT, and the columns correspond to basis function. The unconstrained weighted least‐squares solution of Equation 7 is given by [Image Omitted. See PDF]where the measurement covariance matrix, Sm, is taken to be a diagonal matrix that scales inversely with the density of observations at a given emission angle. That is, is proportional to the angular contribution function. Note that the solution in Equation 15 is invariant with respect to a constant scaling of the covariance matrix and is therefore independent of the absolute magnitude of the measurement noise. From this solution, we determine the limb darkening [Image Omitted. See PDF]which is the fractional reduction (in percent) of the intensity of the thermal emission as a function of emission angle. One subsidiary goal of the Juno MWR's objective of quantifying Jupiter's global water abundance is a determination of R(45) in channels sensitive to water to within 0.1% (Janssen et al., 2017).
Figure 4 compares the errors in fits to synthetic data generated using geometries for both a typical MWR orbit (perijove 5) and a typical gravity orbit (perijove 11). Each axis plots |R(fit)(45) − R(model)(45)|, the absolute difference of the limb darkening at 45° (in percent) between “true” (i.e. modeled by JAMRT) and fit values for an ensemble of atmospheres within the bounds specified in Figure 2. On the ordinate axis, no shape function is employed, while on the abscissa the shape function based on the atmospheric boundary conditions specified in Figure 1 is used in the fitting. Each individual data point represents a unique model atmosphere from the ensemble, and each color represents a different channel. The cluster of points labeled PJ5 corresponds to the MWR orbit, while data outside this cluster refer to results for the gravity orbit. The tight clustering of perijove 5 data in Figure 4 shows that the fits for the MWR orbit are much better than for the gravity orbit. The contrast stems from the difference in angular contribution function shown in Figure 3. The maximum biases in limb darkening are 0.009% for the case where a shape function is used and 0.027% for the case where no shape function is used. These values are safely within the bounds of MWR's target of 0.1%. For this type of orbit, the use of a shape function is salutary but not critical. In contrast, for the orbits whose angular distribution function is similar to that of the gravity orbit, not using a shape function leads to unacceptably large biases. Figure 4 shows that for all model atmospheres in a gravity orbit the fits for channel 2, for example, are all biased between −0.64% and −0.46%, while using the shape function constrains the bias for all channels except channel 1 to less than 0.11%.
Enlarge this image.4 Figure. Comparsion of limb darkening fit in percent with (x‐axis) and without (y‐axis) the use of a shape profile for a nearly ideal perijove (PJ5) and a gravity orbit (PJ11). Channels are differentiated by color and each data point represents a model atmosphere within the bounds specified in Figure 2. The cluster of points near the origin and labeled by “PJ5” represents the use of a synthetic measurement set whose emission angle distribution is that of perijove 5, a typical MWR orbit; the remaining points correspond to the emission angle distribution of perijove 11, a typical gravity orbit. As a guide to the eye, faint gray dashed lines indicating unbiased fits for the cases with and without shape profile are superimposed. For the gravity orbit, employing a shape function reduces bias in fitted limb darkening.
Each cluster of points roughly follows a line of unity slope. This means that the spread of the error in the fits over the ensemble of modeled atmospheres is not altered by the introduction of the shape function. Rather, the key benefit of the shape function is to reduce the overall bias of the solution (effectively, the distance of a cluster of points to the dotted lines in Figure 4). Interestingly, channel 1 is not strongly biased even without the adoption of a shape function, so its use for this channel does not impart a significant improvement. However, channel 1 does exhibit a larger spread of errors across the range of atmospheres in the model ensemble than the other channels do. The larger spread was found to be largely due to JAMRT's sensitivity to water abundance: the largest values of |R(fit)(45) − R(model)(45)| in Figure 4 correspond to extreme values of the model ensemble's deep water abundance (0 or 7,000 ppm). Our water opacity model (Bellotti et al., 2016) suggests that the water opacity may become comparable to that of ammonia at depths to which channel 1 is still sensitive (although uncertainties in the water absorption are quite large at those depths). Our model's increasing water opacity at depth and the wide range of water abundances spanned in the model ensemble result in a correspondingly larger spread in errors in channel 1.
An analysis of the fit of the nadir brightness temperatures follows a narrative similar to that of limb darkening. Relative errors in nadir brightness temperatures for the MWR orbit are small for all model atmospheres (less than 0.066%) even without the use of a shape function. Using a shape function for the gravity orbit reduces the fitting bias from 0.61% to less than 0.10% for channels 2 through 6 while slightly reducing maximum fitting errors for channel 1 from 0.28% to 0.24%.
In sum, we have demonstrated a simple three‐term representation of the angular dependence of Jupiter's thermal emission that is able to characterize a fairly wide range of model atmospheres under conditions similar to those encountered in MWR‐type orbits to a precision that satisfies mission requirements. Additionally, the introduction of a shape function partially mitigates the bias stemming from the unfavorable geometries associated with later gravity orbits. Caution is advised, however, for channel 1, whose greater sensitivity to atmospheric variation negates much of the improvement in bias correction that the shape function provides in the other channels.
The positional basis functions presented in section 2.2 are bilinear interpolation functions, each localized to one of 49,152 equal‐area Healpix pixels. This bijective correspondence between basis function and pixel leads to a rather large basis set as roughly 30,000 pixels make some nonzero contribution to the measured antenna temperatures in a typical orbit, roughly three to eight times the number of observations included in the measurement set (depending on channel and orbital geometry). Direct use of such a large set of basis functions leads to a underdetermined linear system of equations that must be strongly regularized. However, the size of the basis set can be drastically reduced by grouping basis functions together and assuming the coefficients do not vary over the pixels associated with the grouped basis functions. A new set of positional basis functions {HP(r)} may be derived from the original ones through simple summation: [Image Omitted. See PDF]where P is a superpixel comprising a set of pixels p over which the thermal emission is treated as uniform. Thus, the measured thermal emission in Equation 6 may expressed as [Image Omitted. See PDF]where for all pixels p contained in a superpixel P. This means that we can effectively coarsen the operator M by simply summing over selected columns of M.
Optimizing the agglomeration of pixels into larger superpixels is nontrivial and depends on two key observations. First, there may be extended regions encompassing many pixels where Jupiter's thermal emission is not expected to vary greatly, so that over these regions the coefficients can be treated as constant, and a single physical basis function may be used. For example, at lower latitudes, Jupiter's visually bright and dark axisymmetric features, known as zones and belts, are suggestive of far greater variation latitudinally than longitudinally, so that it may be sensible to group pixels into superpixels encompassing many longitudes. Second, the extent of the superpixels should ideally be roughly commensurate with the effective resolving power of the instrument, which may vary considerably owing to orbit geometry, spacecraft spin orientation, and the degree to which data contaminated by synchrotron radiation are screened out. The large eccentricity of Juno's orbit alone changes MWR's effective resolution by nearly two orders of magnitude within a single orbit: for the orbits considered in this study, Juno at perijove is only a few thousand kilometers from Jupiter's one‐bar level, while at high latitudes the spacecraft is much further away (e.g. up to 100,000 kilometers at 75°S depending on details of the orbit). Furthermore, the spatial resolution of the coefficients varies by angular basis function, coarsening with the degree of the associated fitted term. In this work we do not attempt a thorough optimization of pixel grouping. Instead, in section 3, we explore the performance of a couple of variant coarsenings.
The Juno spacecraft takes approximately two hours to travel from pole to pole, over which time roughly 60,000 MWR science data measurements are made per channel. However, only a small fraction of these measurements is useful for mapping Jupiter's thermal emission. There are two general categories of effects that must be screened out so as to limit bias. One is geometric; the other is pollution of the measured intensities from unmodeled sources such as synchrotron radiation, plasma‐induced waveguiding from aurorae, and lightning events. In this section we discuss how these effects are mitigated.
Because Juno spins and is therefore pointed away from Jupiter at least half of the time (and even more so as Jupiter's angular diameter shrinks the further the spacecraft is from the planet), a majority of the MWR science data measurements do not contribute to the deconvolution. It was demonstrated in section 2.3 that the fit of the thermal emission to a quadratic dependence on cosine of the emission angle is improved when the distribution of emission angles that contribute to a given measurement is weighted toward smaller emission angles. Also, geometric uncertainties associated with planet shape and refraction‐induced bending of the optical path are magnified at high emission angles. A trade‐off, therefore, exists between incorporating more data to reduce statistical variance and limiting the data set to geometries that support a good fit to the chosen bases. To address this trade‐off we require that the fraction of each channel's antenna pattern that does not intersect Jupiter not exceed some threshold value, which we take in this study to be 1%. That is, we select only for those measurements that satisfy [Image Omitted. See PDF]This particular threshold value concomitantly enforces a reasonable cutoff for emission angle. During the first perijove, for channel 1 the maximum emission angle passing this criterion is 54° (though this value varies by channel owing to different channels' beamwidths). Typical values for the fraction of measurements exceeding 60° emission angle range from 0% for channels 1 and 2 (whose beamwidths are relatively large) to 2.5% for channel 6, which has the narrowest beamwidth. (An example of the relationship between this criterion and emission angle is demonstrated in the Supplemental Materials.)
Juno's microwave radiometer is sensitive, especially at longer wavelengths, to sources other than black body emission. One such source is lightning (Brown et al., 2018), which manifests itself, primarily in channel 1 but also to a lesser extent in channel 2, as anomalous spikes in measured antenna temperatures. These events are identified through an iterative scheme which first partitions the timeseries into a set of independently analyzed spins. For each measurement a smoothing spline is fit over the entire spin (with the measurement removed) using weights at each observation commensurate with the expected instrument noise level for that observation. The difference between the antenna temperature and the value of the smoothing spline at the time of the measurement yields a residual for each measurement in the spin. If the largest positive residual for a given measurement is greater than the 4σ noise level, the measurement is identified as a lightning event and removed from the observation set. The process iterates until no more lightning events are identified. Failure to account for lightning in the deconvolution was found to produce differences in R(45) of up to 0.3% in channel 1 at northern midlatitudes where lightning is prevalent.
The Juno mission was designed to avoid much of the synchrotron radiation that muddles earth‐based measurements at longer microwave wavelengths. Nevertheless, for particular spacecraft orientations and latitudes it remains the largest non‐atmospheric component of MWR data. Arising from the centripetal acceleration of high‐energy charged particles in Jupiter's magnetic field, the observed synchrotron radiation decays (very roughly) by a factor of 5 (Bolton et al., 2002; de Pater, 2006) with each channel number increment (and corresponding doubling of frequency). Though most problematic for channel 1, its signature is evident in Jupiter‐looking data in the first four channels. We primarily mitigate the effects of synchrotron emission by removing affected observations from the measurement set used to perform the deconvolution. A heuristic algorithm has been developed to screen synchrotron emission, which is described in detail in Appendix A.
Despite our best efforts to exclude synchrotron‐contaminated measurements from the dataset used to perform the deconvolution, synchrotron emission nevertheless can make a non‐negligible contribution to the measured intensities via an antenna's back and side lobes even when the antenna is pointed directly at Jupiter and away from the synchrotron source. Indeed, the greatest source measured by channel 1 is due to synchrotron emission in the equatorial region with maximal values ranging between 4000 K and 6300 K depending on orbit, far dwarfing Jupiter's nadir thermal emission of 800–900 K. In such cases (generally when Juno is within about 30° of the equator), the effect of synchrotron emission cannot be merely screened but must be accounted for explicitly by subtracting its contribution from the measured antenna temperature. To account for synchrotron radiation, we use a semi‐empirical model developed by (Adumitroaie et al., 2016; Connerney et al., 2018; Levin et al., 2001; Santos‐Costa et al., 2017) which computes the synchrotron emission produced by the motion of electrons trapped in the Jovian magnetic field, as it is experienced following a particular line of sight. Originally applied to model Earth‐based observations, the model has been modified for the Juno mission to simulate the synchrotron emission as seen from an arbitrary spacecraft position in the vicinity of Jupiter and has been further enhanced to allow ingestion of externally produced electron distributions, such as those obtained through higher fidelity physics‐based transport simulations (Santos‐Costa & Bolton, 2008). The model is currently not accurate enough to simply subtract the calculated synchrotron contribution from the measured intensities for the purposes of the deconvolution. However, informed by off‐planet measurements, it can be adjusted to yield a presumably better on‐planet estimate. Details of how the back and side lobe contribution is determined are presented in Appendix B.
Using the fine resolution of the Healpix grid described earlier leads to an ill‐conditioned mapping with strongly correlated coefficient covariances. This problem can be mitigated by using some form of regularization that correlates closely spaced points. One method, described earlier, is to group pixels into superpixels over which the coefficients are assumed to be constant. This grouping should be based on the estimated spatial resolution of the mapping M in Equation 7. We note that the variation in thermal emission measured by the radiometer is much greater meridionally than zonally, especially at longer wavelengths and at lower latitudes. The greater meridional variation compared to zonal variation, a result of Jupiter's strong Coriolis force, mirrors many other observed properties, including colorful belts and zones, east‐west jets, high altitude hazes, the upper‐tropospheric concentrations of gases such as ammonia and phosphine, turbulence power spectra, and the distribution of lightning flashes. Furthermore, for nearly all orbits thus far, the spin‐plane has been oriented in such a way that longitudinal coverage has been limited, except in the polar regions. Therefore, we start by assuming the thermal emission is cylindrically symmetric, so that pixels sharing the same latitude are grouped together. We then estimate the latitudinal resolution with a Backus‐Gilbert approach (Backus & Gilbert, 1968; Backus et al., 1970).
Like the regularized minimum residual solution of Equation 8, Backus‐Gilbert involves a trade‐off between two quantities. Whereas Equation 8 entails a compromise between residual and prior information, Backus‐Gilbert attempts to strike a balance between variance and resolution. The motivation behind Backus‐Gilbert is that the variance of the solution due to measurement noise is inversely related to the resolution of the solution. For any given latitude and specified variance, Backus‐Gilbert finds an optimal linear combination of coefficients for which the “spread”, a measure of the latitudinal extent of the linear combination, is minimized. Details of the implementation in the context of this work are described in Appendix C. This linear combination, normalized to unity over an integration in latitude, may be interpreted as a resolution function (or averaging kernel). Figure 5a illustrates resolution functions for all three angular basis functions at three different latitudes (±45° and 0°) for the longest wavelength channel during the first orbit under the constraint that the standard deviation of the limb darkening at 45° due to measurement noise be fixed at 0.1%. Additional assumptions are made that the uncertainties in the coefficients are uncorrelated, and that they contribute equally to the overall limb darkening uncertainty. These assumptions are not strictly valid, but they enable a useful estimate of the instrument resolution. Figure 5a shows that the coefficients are well‐resolved at the equator but broaden significantly away from perijove due to the increasing distance of the spacecraft from the planet. We note that the resolution function is not symmetric about the equator; it is narrower at +45° than it is at −45° because for this orbit, perijove occurs at a sub‐spacecraft latitude that is just north of the equator (∼ 3.8° N). The differing line types show that the resolution decreases in order of increasing angular basis function: the nadir brightness temperature (c0) exhibits better resolution than the absolute limb darkening (c1), which in turn is more finely resolved than the measure of non‐linearity (c2). This ordering of resolution holds for all MWR‐like orbits but is not necessarily true for other types of orbits in which the distribution of emission angles may be peaked far away from nadir, thereby reducing the resolution of the nadir brightness.
Enlarge this image.5 Figure. (a) Examples of the resolution function for each angular basis function at ± 45° and 0° for channel 1, perijove 1 assuming a limb‐darkening uncertainty of 0.1%. The resolution is reduced for higher order coefficients and away from perijove. (b) Resolution function width for channels 1 and 3 for perijove 1 for c0 (solid), c1 (dashed), and c2 (dotted). (c) channel 1 resolution function width for c0 and all MWR favorable orbits. Large variations are due to reduced data set size at latitudes where synchrotron contamination is large.
We define a width of the resolution function to be the minimum range of latitudes over which the integrated resolution function equals 0.68, equivalent to the area within one standard deviation of the maximum of a normal distribution. This width represents an effective resolution (in degrees) of the coefficients achievable under the constraint of a 0.1% uncertainty in limb darkening at 45° due to measurement noise. Figure 5b plots this latitudinal resolution for channels 1 and 3 as a function of latitude for the first perijove. These two channels are representative of the others not shown here: Channel 2 is very similar to channel 1, and channels 3 through 6 have nearly identical widths. There are two reasons for this clustering. First, the beam widths of channels 1 and 2 (∼21°) are much larger than those of channels 3 through 6 (∼11° to 12°) (Janssen et al., 2017). Their respective maximal resolutions are consequently poorer. Second, the lower channels are more strongly affected by the screening of measurements contaminated by synchrotron radiation. For the MWR orbits, the number of measurements used in the deconvolution of channel 1 data is typically a factor of two smaller than those of the higher channels. The reduced data set size of the lower channels further limits the resolution of the instrument. Additionally, because synchrotron contamination is more prevalent at off‐nadir angles, the resolution of angular coefficients associated with higher powers of μ decreases more rapidly than the nadir brightness c0, as evidenced by the local peak of c2 around −20° in Figure 5b.
Figure 5c shows the variation by perijove of the latitudinal resolution for the nadir brightness c0 of channel 1. For most orbits, a resolution finer than 5° is found over a range of −30° to 50°, but certain orbits (perijoves 3, 4, 6, and 9) exhibit regions where synchrotron pollution has substantially limited the instrument's resolving power.
The Backus‐Gilbert approach described in the previous section and in the appendix yields a solution to Equation 7 in terms of a transformed set of coordinates , where the rows of A are the aforementioned resolution functions. It is tempting to try to invert A to determine both the fine‐grid coefficients c and the residual T (J) − Mc as a measure of the goodness of the fit. However, A is very ill‐conditioned; coefficients derived is this manner exhibit large, unphysical fluctuations in regions where we try to resolve on a length scale finer than measurements allow. We consider two approaches to determining the brightness coefficients, a straightforward weighted least‐squares solution on a coarse grid and a regularized one.
In the first approach, the angular coefficients are grouped into latitude bins whose widths are commensurate with the those of the resolution functions as described in the previous section. That is, the grid spacing at any given latitude is approximately equal to the width of the resolution function determined by specifying a 0.1% uncertainty in the limb darkening at 45°. In this case no regularization is applied, and the weighted least‐squares solution is given by Equation 15. The difficulty with an unregularized approach is in finding precisely the correct resolution for each coefficient that simultaneously avoids overfitting to instrument noise while eliminating latitude‐dependent features in the residual, all while constrained to a Healpix grid that allows only for discrete changes in resolution.
The second approach employs a prior estimate and regularizes with an appropriate (diagonal) choice of Sc in the optimal estimation solution (Equation 9). Different priors were tested, including a linear interpolation of the coarse‐grid approach as well as using perijove‐averaged values (to be discussed in section 4), but the derived coefficients were not found to be strongly sensitive to choice of prior, particularly within about 20 degrees of perijove. For this work, we use the Backus‐Gilbert‐derived coefficients as a priori estimates. These coefficients represent a weighted average value over some distribution associated with each latitude, but in this case we treat them as prior values localized at each fine latitude grid point. The degree of regularization is governed by the magnitude of Sc. The maximum curvature point of the so‐called “L‐curve” is an often used heuristic to determine Sc (Hansen & O'Leary, 1993), but we found that this choice resulted in an overly regularized solution. Instead, we set to be a diagonal matrix, whose nonzero values are proportional to the mean nadir brightness temperature used in the data set for each channel. The proportionality factors, which represent the degree to which the coefficients are allowed to vary before the regularization term becomes comparable to residual term, are taken to be 5 × 10−3 for the nadir brightness c0 and 1.25 × 10−3 for the other two angular basis functions.
It is also important to note that the reduction of the deconvolution to one dimension (i.e., independent of longitude), while appropriate at low to midlatitudes, becomes less useful in the polar region for two reasons. First, the spacecraft is further from the planet at high latitudes so that the range of longitudes that contribute to a measurement increases. Second, at high latitudes the regions over which the model assumes the physical properties to be constant become annuli in the limit of full polar coverage rather than compact areas. For this reason, we restrict figures in this work to latitudes less than 60° and defer analysis of polar regions to future work.
Figure 6 comprises panel plots of the latitudinal dependence of various quantities derived from the deconvolution of channel 3 measurements during the first perijove. The upper four panels (a–d) show the output of the deconvolution. The first three are the coefficients themselves, and the fourth, derived directly from the first three, is R(45), the limb darkening at 45°. For this channel, the limb darkening is qualitatively similar to c1 because the quadratic term c2 is relatively small, and the variation in nadir brightness c0 across latitudes is less than 15%. The coarse grid solution is represented by red error bars whose horizontal extents indicate the latitudinal coarseness of the solution and whose vertical extents represent the one‐sigma uncertainties due to instrument noise alone. It is important to realize that these uncertainties account only for variance due to instrument noise and therefore may underestimate the total error in regions where the fit is poor. For this data set, a priori Backus‐Gilbert coefficients (indicated by the dashed lines) are nearly indistinguishable from the regularized solution (shown as solid blue lines). The good agreement is due to the relatively small spread of the resolution function for this perijove as shown in Figure 5.
Enlarge this image.6 Figure. Example of deconvolution from channel 3 measurements during the first perijove for both regularized (blue) and unregularized, coarse grid (red) solutions. Horizontal/vertical extents of the error bars indicate grid coarseness and 1σ standard deviation arising from instrument noise, respectively. Backus‐Gilbert solutions (dashed line), used as a prior value, are nearly indistinguishable from the regularized solution for this case. Panels (a), (b), and (c) plot the three angular coefficients, and panel (d) shows the limb darkening at 45°. Panels (e) and (f) plot individual residuals relative to the instrument noise level (≈0.38 K for this measurement set) as well as moving averages of the residuals (multiplied by two to better visualize spatial variation). The regularized solution tracks the unregularized solution in the top four panels and also eliminates spatial features in the residuals of panel (e). A goodness of fit metric, the reduced χ2, as defined in Equation 20 is displayed in panel (g) and shows a few regions (∼−39°, ∼15°, and ∼37°) where the fit is not optimal due to unmodeled longitudinal variation.
The remaining panels illustrate the goodness of the fits. Panels (e) and (f) plot residuals of each measurement for the two cases, normalized to the instrument noise level, which from Table 1 is approximately 0.38 K for this measurement set. Abscissa values denote planetocentric latitude intersect points of the antenna boresight with the planet. We loosely term these footprint latitudes, though one should not lose sight of the fact that the measurements encompass contributions from footprints of variably sized and shaped neighborhoods containing the intersect points. The solid curves are a moving average of the residuals over a 1.4° latitude interval (multiplied by a factor of two to place them on the same scale as the residuals). For the coarse‐grid case, spatial structure is clearly evident in the moving averages, indicating that this solution is underresolved in at least one of the coefficients. It is important to keep in mind that the Backus‐Gilbert widths discussed in section 3.1 serve only as an estimate of the spatial resolution needed to reduce the instrument noise‐induced spread of R(45) to the 0.1% level but do not necessarily represent the ultimate resolving power of the instrument. We did not find it generally possible to reduce residuals to the level of instrument noise within the constraints imposed by Healpix discretization without either overfitting or some sort of regularization. In contrast, the regularized solution in Figure 6f exhibits locally averaged residuals much closer to zero.
Both residual panels also evince latitudinal variability in the spread of the residuals. A measure of this spread in shown in the bottom panel (g) which plots an effective “local” reduced χ 2, which we define as [Image Omitted. See PDF]wherein the squares of the normalized residuals are summed over the same 1.4° latitude interval as in panels (e) and (f) and divided by an effective number of degrees of freedom per bin, N(ϕ) − ν, where N(ϕ) is the number of measurements within the bin and ν is the effective number of independent model parameters associated with the bin. We estimate ν to be total number of coefficients solved in the coarse grid solution multiplied by the fractional size of the bin (1.4/180). For this perijove and channel, typical values of N range from 60 to 90, depending on latitude, and ν ≈ 1.45. Values of χ 2 within a range of 1 ± 2N −1/2 are indicative of a good fit. Values below this range would indicate overfitting, but there is no evidence for overfitting in these data. It is clear, however, that neither model fully describes Jupiter's thermal emission to the extent the measurements allow, as evidenced by the elevated values of χ 2 at certain latitudes (e.g., at ∼−39°, ∼15°, and ∼37°). These latitudes correspond to regions where the MWR instrument both is able to resolve and actually detects longitudinal structure. We note that while larger values of χ 2(ϕ) indicate a fit that is not optimal, they do not necessarily indicate the retrieved values are biased: Enhanced values of χ 2(ϕ) are visible around 15° in panel (g) that do not correspond to significant departures from zero in the averaged mean values in panel (f). The converse, however, is true: Biased solutions (e.g., near −39° ) are reflected in elevated values of χ 2(ϕ).
It is possible for the deconvolution algorithm to reduce χ 2 so that it is close to unity everywhere, but additional constraints must be imposed to reduce the space of solutions that are still consistent with the instrument uncertainty model. In the regularized case described above, Healpix pixels at the same latitude are grouped into a single superpixel over which the angular coefficients are constant. To allow for longitudinal variation, we tighten the grouping criterion as follows. We first define a spatial contribution function ηp for all pixels p, whose sum over p equals the total energy measured by the instrument over the measurement set used in the deconvolution: [Image Omitted. See PDF]Comparison with Equation 6 shows that [Image Omitted. See PDF]where the summation is over the three angular basis indices k and the measurement index i. The coefficients are taken from the longitudinally independent regularized solution so that for fixed k the coefficients ckp have the same value for all pixels p at the same latitude. This spatial contribution function is effectively a measure of the total energy measured by MWR that is emitted by a given pixel over the data set used in the deconvolution. At each of the 255 unique latitudes of the Healpix grid, the pixel with maximum contribution function is identified. If the contribution function exceeds some threshold value η∗, the pixel is added as a singleton to the set of spatial basis functions. Otherwise, its nearest neighbors (at the same latitude) are added one at a time to a grouped superpixel until the sum of the contributions exceeds η∗. If, at a given latitude, pixels with nonzero contribution functions remain, the process iterates until they are exhausted. Upon its conclusion, the algorithm yields a set of single latitude pixel groupings, each of which contributes at least η∗ to the total energy received by MWR. Setting η∗ to 30% of the maximum value for ηp typically yields ∼1,000 groupings, or an average of roughly four longitudinal points per latitude.
Solving Equation 9 using the pixel grouping described above yields a two‐dimensional map of the coefficients. However, because the total number of coefficients becomes comparable to the number of measurements, it becomes much more difficult to disentangle correlations among the coefficients. Figure 7 provides examples of two extreme cases in the solution of Equation 9. The top panel plots the spread in nadir brightness c0 (in orange), which is allowed to vary in both latitude and longitude. In this case, the coefficients c1 and c2 are held fixed, and a prior determined in Figure 6 (shown in blue) is assumed. The coefficient covariance for c0 is floated so that the total χ 2 is approximately 1.1. The bottom panel shows that the local reduced χ 2 is everywhere near unity, which means that allowing for longitudinal variation enables the existence of a solution consistent with the characteristics of the instrument. The middle panel plots R(45) (in purple), using the same prior as in the top plot. However, in this case the nadir brightness is fixed, and only c1 and c2 are allowed to vary. The resulting reduced χ 2, illustrated in the bottom panel, is also diminished compared to the longitudinally invariant case and matches very closely with the case shown in the top panel. Both solutions represent valid fits to the measurements; additional prior information is needed to favor one over the other or, as seems more likely, some intermediate combination of the two. Because of the ambiguity inherent in disentangling longitudinal limb darkening variations from those of nadir brightness, we favor the 1D latitudinal solutions at low and mid latitudes for the remainder of this analysis, where solutions in regions where the local χ 2 is large are given less weight.
Enlarge this image.7 Figure. (a) Nadir brightnesses using the 1‐D solution (blue) as a prior for a longitudinally dependent solution (orange) in which c1 and c2 are held fixed. The vertical spread in orange points at a given latitude is a measure of the possible longitudinal variation at that latitude. (b) 45° limb darkening using the 1‐D solution (blue) as a prior for a longitudinally dependent solution (purple) in which the nadir brightness c0 is held fixed. (c) Admitting additional degrees of freedom to allow for longitudinal variation reduces the spread in residuals for both cases (a) and (b) so that χ2 ∼ 1, but how much to attribute to variation in nadir brightness rather than limb darkening is not resolvable in this case.
Figure 8 shows the local reduced χ 2 as defined by Equation 20 for all MWR‐favorable orbits and for each of the channels. Under most circumstances, χ 2 ∼ 1, indicating a good fit to the measurements to within instrument noise. However, for particular orbits, channels, and latitudes, χ 2 can be quite large. In these cases, the measurements do not offer a consistent story (to within the expected noise level). This may be due to imperfect synchrotron screening (e.g., channel 1, PJ4, near 35°N), auroral effects (e.g., channels 1, 2, and 3, PJ5, north of 50°) (Hodges et al., 2020), or spatial (i.e., longitudinal) inhomogeneities. Large values of χ 2 do not necessarily indicate that the retrieved coefficients are biased, only that there is the potential for them to be so. Figure 8 also demonstrates that higher channel number generally correlates with more instances where χ 2(ϕ) is large, especially in the north and south equatorial belts. This reduction in quality of fit at shorter wavelengths is due to a combination of greater actual longitudinal variation (see, e.g., de Pater et al., 2016; Fletcher, Kaspi, et al., 2020) and the smaller beam sizes which increase the resolvability of smaller features. In contrast, at the longest wavelengths, the quality of fit generally improves because there is less resolvable longitudinal variation in the thermal emission. One example of large longitudinal variation evident in Figure 8 that can be correlated with other observations is illustrated in channels 5 and 6 in perijove 4 in the planetocentric latitude range of −14° to −9°. Four weeks before Juno passed over this region, a mid‐Southern Equatorial Belt convective outbreak was observed in the vicinity of this location in contemporaneous measurements by the Hubble Space Telescope, Atacama Large Millimeter/Submillimeter Array, Very Large Array, the Very Large Telescope, and the Gemini telescope (de Pater et al., 2019; Wong et al., 2020).
Enlarge this image.8 Figure. Values of local reduced χ2 as defined by Equation 20 for each channel. Regions where χ2 ≫ 1 show where residual distributions are much wider than would be expected from instrument noise alone and are indicative of either longitudinal variation or contamination from nonthermal sources.
Figure 9 shows deconvolved nadir brightnesses and 45° limb darkening for all MWR‐type orbits and all channels. Black lines denote weighted averages over all the perijoves, in which the coefficients are weighted by both 1/χ 2(ϕ) and the square root of the spatial contribution function Equation 22 to account for the effective number of measurements contributing to the deconvolution at a given latitude. Latitude ranges where the goodness of fit is deemed suspect (χ 2(ϕ) > 3) are depicted by dotted lines.
Enlarge this image.9 Figure. Deconvolved nadir brightness temperatures and 45° limb darkening for all MWR‐favorable orbits through PJ12. Solid/dashed curve segments correspond to regions where χ2 < 3 and χ2 > 3, respectively. Black curves denote averaged values weighted by 1/χ2 and the relative contribution to the measurement set. For reference, shaded regions denote the range of values spanned by the JAMRT model atmospheres as detailed in Figure 2 in section 2.3. The darker shading additionally constrains the deep ammonia abundance of the model to 351 ppm.
For a sense of comparison to a simple adiabatic model, ranges of JAMRT‐derived values are superimposed as shaded regions. Two different shades are shown. The lighter shading corresponds to the full range of atmospheric parameters used to validate the choice of angular basis functions in section 2.3, while the darker shading additionally constrains the NH3 deep atmosphere abundance to 351 ppm. The relative sizes of the light and dark shaded regions serve as a crude visual indicator of a coefficient's sensitivity to deep atmosphere ammonia abundance relative to its sensitivity to other atmospheric parameters. A similar model, but which also allowed for a more complex characterization of Jupiter's atmosphere, including perturbations to JAMRT's simpler adiabatic vertical profiles, formed the basis of a detailed analysis of the implications for Jupiter's water abundance in the equatorial region in work by Li et al. (2020).
It is worth noting that for channel 1, at all latitudes the 45° limb darkening is substantially lower than the entire range of JAMRT's results, except for a region just north of the equator, where it is only a few tenths of a percent smaller. Outside of a ∼15°‐wide region about the equator, R(45) is at least a full percent smaller than JAMRT's lowest value. This large discrepancy suggests inaccuracies in our current opacity models. There are two principal suspects: free‐free absorption and water opacity. At sufficiently high temperatures, thermal collisions ionize atoms, resulting in a plasma. The degree of ionization depends on the atoms' abundances and ionization energies. Owing to their relatively low ionization energies, alkali metals such as sodium and potassium are expected to be the primary contributors in the pressure and temperature regime that channel 1 is sensitive to. The mixture of conducting electrons and neutrals in the deep atmosphere can create an effective “wall” of opacity at some critical temperature, over which black‐body radiation is not able to propagate upward (Bellotti, 2018). This effect is not currently taken into account in JAMRT. The second source of the limb darkening discrepancy arises from JAMRT's water opacity model. The model (Bellotti et al., 2016) is based on laboratory measurements that were taken under Jovian‐like conditions up to pressures of 100 bar (Karpowicz & Steffes, 2011b), a level far beyond which JAMRT must extrapolate in order to model channel 1 measurements. Furthermore, this model and others (Karpowicz & Steffes, 2011a, 2011b) indicate that the relative opacity of water to ammonia increases substantially at the high pressures and temperatures that contribute to the emission observed by channel 1. Further advances to align the model with observation may necessitate additional laboratory measurements of high pressure water absorption.
The narrow width of the channel 6 gray bar, compared to the actual range of observed variation in nadir brightness temperature and R(45), reveals the effects of simplifying assumptions within the model, especially in the upper troposphere, the only region to which channel 6 is sensitive. In this region, the model's thermal emission depends almost entirely on the ammonia concentration, fixed to the saturation vapor pressure, and an adiabatic temperature profile, which is constrained to a small range (±2.8 K) of reference temperatures at 0.5 bar centered at a nominal value of 132.79 K (Seiff et al., 1998). In reality, kinetic temperatures at these altitude levels vary spatially and temporally (Fletcher et al., 2016; Simon‐Miller et al., 2006), and substantial ammonia subsaturation may be widespread (de Pater et al., 2016). An attempt to account for the observed variation in thermal emission in the upper troposphere using upper‐tropospheric temperature and ammonia variations derived from TEXES observations is discussed in Fletcher, Orton, et al. (2020).
The greatest outlier among the orbits shown is PJ7 (orange) during which Juno passed over the Great Red Spot (GRS), effects of which are visible over a range of latitudes from 30°S to 10°S. Channels 1 and 2 exhibit significantly higher nadir brightness temperatures relative to the averaged value. This elevated emission largely disappears in channel 3. At higher channels, a bifurcation becomes increasingly prominent, in which a relatively warm region south of the GRS coexists with a relatively cold region to the north. Apart from an enhancement in channels 1 and 2, the 45° limb darkening is consistent with variations seen in the other orbits.
Another interesting feature is a latitudinally confined cold region between 10° and 15° observed in channel 1 that appears in a small number of perijove passes (PJ4, PJ5, and possibly PJ7). Interestingly, there is no similar signature seen at that latitude in the other channels for those perijoves. Because this feature is absent in channel 2, it is unlikely that it is due to synchrotron emission or some aurora‐like phenomenon. Rather, it suggests structure that only exists very deep ( 100 bar) in Jupiter's atmosphere.
A periodic, latitudinal structure in the nadir brightness temperature is apparent from −60° to −25° in channels 5 and 6 and mimics what is already known about the upper tropospheric temperature and ammonia distribution (Fletcher, Kaspi, et al., 2020). This structure largely vanishes in channel 4 but re‐emerges in channels 1, 2, and 3 in a fashion that is anti‐correlated with its higher altitude behavior. Figure 10 expands on the nadir brightness temperatures in this latitudinal range. The top six panels display gravity‐adjusted nadir brightness temperatures in K for MWR's six channels relative to their mean values over the range −60° to −25°. The gravity adjustment accounts for the change in thermal emission stemming from the increase in lapse rate attendant with the poleward increase in effective value of gravity. It is implemented as a multiplicative scaling factor of , where ϕ is latitude, and TB is the nadir brightness temperature derived from the same model atmosphere used in section 2.3. For a compositionally uniform deep atmosphere with an adiabatic profile and boundary conditions as in the model atmosphere described in section 2.3, the gravity adjustment scales the brightness temperatures at a given latitude as though the local value of gravity were equal to its value at the equator. This adjustment successfully detrends the poleward enhancement of emission seen in Figure 9; variations over this span of latitudes is less than 4 K for all channels. Additionally, the amplitudes of the oscillations ( 2 K) do not differ strongly by channel, although it is possible that the deconvolution could smooth some channels more than others. The bottom panel shows zonal wind speeds from Hubble OPAL observations (dashed line) (Simon et al., 2015; Wong et al., 2020). Cyclonic belts, indicated by gray bands, are defined by prograde jets on the equatorward side and retrograde jets on the poleward side. They generally show greater brightnesses at high altitude (channels 5 and 6) but lower brightnesses at depth (channels 1, 2, and 3). The fact that the deep‐sensing MWR channels for these latitudes appear to be anticorrelated with those channels sensing the upper troposphere is something that will be investigated in future work, but they imply that either (i) ammonia (or another opacity source) is enhanced in belts at depth, which is opposite to what is seen in the belts in the upper troposphere; or that (ii) kinetic temperatures are cooler in the belts at depth, opposite to what we see in the belts in the upper troposphere; or (most likely) (iii) a combination of those two factors. In either case, this result shows that there's a complex vertical structure to the belts and zones of Jupiter (e.g., see review by Fletcher, Kaspi, et al., 2020), with zonal organization persisting to the deepest levels sensed by MWR. Modeling to disentangle the thermal and compositionally effects is ongoing.
Enlarge this image.10 Figure. The top six panels display each channel's nadir brightness (in K) relative to their means over the plotted southern temperate latitudes. The brightnesses are averaged over all the perijoves and adjusted to account for gravity (as explained in the text). Channels 5 and 6 are anti‐correlated in latitude with channels 3, 2, and 1. Regions of cyclonic wind shear (gray bands) correlate with higher brightness temperature at low pressures (channels 5 and 6) and with lower brightness temperature at depth (channels 1, 2, and 3). The northern‐most (right‐most) cyclonic shear region is co‐located with Jupiter's South Temperate Belt (Fletcher et al., 2016). The bottom panel shows zonal wind speeds from Hubble OPAL observations (dashed line) (Simon et al., 2015; Wong et al., 2020).
This work has detailed a procedure for extracting the angular dependence of Jupiter's thermal emission from measurements made by the Juno Microwave Radiometer for angles less than ∼45°. A key component of this effort has been the development of a model representation that is both sufficiently descriptive to characterize the emission to within instrumental errors yet minimal enough not to suffer from excessive degeneracies. A heuristic screening algorithm has also been developed to identify and cope with the subset of MWR data contaminated by undesired synchrotron emission and auroral effects. We expect our procedure to benefit from the improved models that are being developed to quantify these effects. Nadir brightness temperatures and limb darkening determined from the application of our procedure to data from the orbits favorable to MWR thus far were demonstrated, and prominent salient features were identified. The physical implications of some of these features have been discussed in a number of papers (Fletcher, Orton, et al., 2020; Hodges et al., 2020; Li et al., 2017, 2020; Zhang et al., 2020), and we expect more to follow.
A four‐step heuristic has been developed to screen synchrotron emission. In the first step the antenna temperature time series is partitioned into independently treated spins. For each spin the largest interval for which the second time derivative of the antenna temperature (computed using a smoothing spline to account for instrument noise) is less than a small threshold value of 2 K/s 2 is retained, while other measurements are discarded. (It was found empirically that using smaller threshold values had the potential to remove uncontaminated data that represented real spatial variations in the Jupiter's thermal emission.) A timeseries corresponding to a portion of one spin is illustrated in Figure A1a, during which time Juno's sub‐spacecraft latitude moves from 47.5°S to 48.5S°. As it spins, the spacecraft traces a path across Jupiter scanning from north to south. The area shaded in blue corresponds to the time interval during which 99% of the beam falls on the planet. On either side of this window, the spacecraft is pointed further away from Jupiter and sees increasing synchrotron emission, peaking at 1750 K and 380 K on the northern and southern horizons, respectively (outside the scope of the figure). Measurements that satisfy the curvature criterion are denoted by either blue or green points. The blue points are measurements that are ultimately discarded because they fail at least one other criterion of the screening procedure, while green points successfully survive further screening. A lone instance of lightning is also observed during this spin (marked in red) and removed from the measurement set.
Enlarge this image.A1 Figure. Examples of various screening criteria applied to perijove 1, channel 1. (a) Points correspond to individual measurements of the antenna temperature time series for an 8 s portion of one spin of the Juno spacecraft. The solid curve is a smoothing spline fit to these points with lightning events (red points) removed. The blue rectangle corresponds to the time window wherein 99% of the beam falls on the planet. Measurements that satisfy the curvature criterion d2T/dt2 < 2 K/s2 but fail other criteria are in blue. Measurements that pass all filters are in green. (b) Fore (orange) and aft (purple) quadratic fits of the antenna temperatures. The fit to aft data exhibits an unphysical concave‐up curvature. Points marked by an “x” are discarded because they are contaminated by synchrotron. Final fit of all surviving measurements is indicated by the black line. (c) Residuals with respect to the final quadratic fit are plotted as a function of a signed cosine of emission angle, where fore/aft views are considered positive/negative.
The second stage in the synchrotron emission screening inspects for systematic outliers among measurements viewing the same neighborhood but at different angles. The measurements that have survived prior winnowing are binned into subsets of equal size and nearly equal boresight footprint latitude, so that bin members are expected to have approximately the same thermal emission properties. A trade‐off exists in the choice of bin size: inclusion of too many measurements leads to greater inhomogeneity in the physical properties among its members and therefore hampers the search for outliers, while too few measurements may cause random instrument noise to find too many spurious outliers. It was found that roughly 40 measurements per bin offered a good compromise. These measurements are sorted, starting with the largest emission angle on the side of the time interval for which the synchrotron emission is smallest. (In Figure A1a, this would correspond to the rightmost green dot.) The first four measurements are assumed to be uncontaminated by synchrotron and included in the observation set. A least‐squares quadratic fit is made. If all remaining measured antenna temperatures exceed the values extrapolated from the quadratic fit, those remaining points are discarded on the suspicion that they might be contaminated by synchrotron radiation. Otherwise, the size of the observation set is incremented by one and the process iterates until every measurement in this latitude bin has either been discarded or included in the observation set.
Figure A1b plots the measured antenna temperature as a function of the cosine of the emission angle for the measurements contained in a latitude bin of width 1.1° and centered at 48°S. Measurements eliminated in the second stage of the synchrotron screening are marked by “×” s. Quadratic fits to aft‐only and fore‐only measurements are indicated by purple and orange dotted lines, respectively, The fit to the screened data is shown in black and is nearly identical to the fore fit. The fore and aft fits, by contrast, differ significantly – the fit to aft‐measurements shows an unphysical concave‐up curvature, indicative of synchrotron contamination. Figure A1c shows residuals with respect to the screened fit as a function of the signed cosine of the emission angle, where positive/negative indicate southward‐looking (aft) / northward‐looking (fore) observations, respectively. Most of the points removed from the final screening process, again indicated by “×” s, show large residuals, while the root mean square of the retained measurements relative to instrument noise equals 1.0, signifying a consistency in the measurements in this bin.
The third stage reverts to a spin‐based perspective. For each spin, measurements are partitioned into fore and aft views, and for each view any measurement with emission angle greater than the largest that has already been discarded is also removed from the measurement set. This filtering criterion stems from the observation that absent synchrotron emission between the spacecraft and planet, the magnitude of synchrotron radiation should be an increasing function of emission angle. This step flags additional synchrotron contamination that may have eluded detection in the previous stage.
The final stage also employs a spin‐based perspective. It is shown in B1 that the synchrotron emission can make a non‐negligible contribution to the measured antenna temperatures via an antenna's back and side lobes, especially for the lower channels. That section also demonstrates how this contribution is estimated. Because this back lobe contribution is expected to have a large relative uncertainty, we apply a filter to remove large variations within a spin. Specifically, we discard measurements where the difference between the estimated back lobe synchrotron contribution and the minimum estimated synchrotron contribution within a spin is greater than twice the instrument noise.
The presence of very large synchrotron emission in channels 1 and 2 can make a non‐negligible contribution to the measured intensity through an antenna's back and side lobes even when the antenna is pointed directly at Jupiter and away from the synchrotron source. To account for this, we use the synchrotron model described in section 2.5 to generate skymaps of radiation intensity for each measurement as seen by each channel. These skymaps are then convolved with the antenna gain functions as in Equation 3 to compute an estimate of the synchrotron contribution to each measurement. The resultant time series yields a qualitatively good description of the magnitude of synchrotron emission but is not quantitatively accurate enough to use unmodified in the deconvolution. We demonstrate by example how the model is used to estimate the synchrotron radiation's contribution to the measured antenna temperature via the back and side lobes.
Figure B1 illustrates a segment of a timeseries within a single Juno spin during the first perijove for channel 1. This window is partitioned into off‐planet (orange), on‐planet (blue), and limb (gray) regions. The gray‐blue and gray‐orange interfaces mark where 90% and 10% of the antenna gain pattern falls on the planet, respectively. Raw MWR measurements of the brightness for channel 1 are indicated by black points, and the synchrotron model is denoted by the dashed curve. In the (orange) off‐planet region, the black‐body contribution to the antenna's backlobe is negligible, so that the measured intensities are due almost entirely to synchrotron radiation. In this region, the model qualitatively matches the measured data; peak values are typically within a factor of two of each other, and the ephemeris times of extremal points are approximately aligned. In the (blue) on‐planet region, which contains the measurements used in the deconvolution, the synchrotron contribution comes from the antenna patterns' back lobes and is therefore quite small, more than two orders of magnitude less than the thermal contribution. However, because a precision of ∼10−3 is desired for limb darkening, this back lobe contribution is important to take into account. We compute estimated values of synchrotron intensity for each of the three aforementioned regions. In the on‐planet region, we assume the same temporal profile as the model but scaled by the ratio of the measured and modeled intensities integrated over the entire off‐planet region during a single spin. That is, the model is scaled by [Image Omitted. See PDF]
Enlarge this image.B1 Figure. Example of synchrotron accounting for a portion of a time series corresponding to a single Juno spin, where orange/gray/blue regions denote off‐planet/limb/on‐planet pointing. Black points represent actual MWR measurements; the dashed/solid gray curves denote the original/adjusted synchrotron models. Green dots represent the synchrotron contribution to the measurements used in the deconvolution. Hatched regions indicate the extent over which limb and on‐planet solutions are significantly mixed.
In the limb region, the decay of the antenna temperature is largely driven by Jupiter's occultation of the synchrotron. Therefore, we expect the time decay of this region to be approximately proportional to that of the model and use a proportionality factor determined by a best fit least‐squares scaling for a short one second interval just outside the limb. For the limb to planet transition, this proportionality factor is given by [Image Omitted. See PDF]where t is in seconds. For the planet‐to‐limb transition, the integration interval is [tlimb, tlimb + 1]. These integration intervals are illustrated in Figure B1 by orange points that represent modeled synchrotron emission which must be scaled up to best fit the corresponding measurements (black points). Finally, in the off‐planet region the thermal contributions to the back lobes are ignored, and the measurements are assumed to consist entirely of synchrotron emission. The scaled model, of course, does not generally produce values that are consistent at region interfaces. To preserve continuity, the scaled models are mixed across a region interface at time t0 with a smooth mixing function: [Image Omitted. See PDF]where Δt is the time required to advance from 75% of planet beam coverage to 90%. The hatched regions in Figure B1 illustrate where the limb and on‐planet regions are mixed to within one decay time constant (76%). The solid curve indicates the final scaled model with superimposed green dots highlighting the estimated synchrotron contribution to the measurements used in the deconvolution.
The Backus‐Gilbert method of solving inverse problems is a well‐known tool in earth science. One limitation is its relatively large computational cost compared to other methods. Restricting the solution space to one‐dimension (i.e., neglecting longitudinal dependence), however, limits its computational expense and makes it a more attractive option. The formulation and derivation of the approach is treated in many sources (Backus & Gilbert, 1968; Backus et al., 1970; Rodgers, 2000; Tarantola, 2005). Here, we provide only enough detail relevant to the problem at hand.
Historically, the dependent variable (in this case latitude) in Backus‐Gilbert has been treated as a continuous quantity, but in this work we represent it as a discretized quantity defined by the 255 unique latitudes of our Healpix grid. We define to be an “unresolved” version of the coefficients c defined in Equation 7 and represent its components as linear combinations of the measurements so that [Image Omitted. See PDF]where T(J) is a column vector consisting of the observed antenna temperatures in the measurement set and Q is some linear transformation determined by properties we discuss shortly. The quantity is intended to be an averaging kernel (or resolution function), so that its ℓth row essentially specifies a smoothing of cℓ, as illustrated in Figure 5a. In order for A to serve as a resolution function, we specify that for all rows ℓ the elements of A sum to unity: [Image Omitted. See PDF]where R, formed from summing over the columns of M, is a vector of dimension equal to the size of the measurement set. Note that in Equation C2, the index of summation j indicates a sum over both superpixels, each identified with a unique latitude, and the three angular basis functions. Backus and Gilbert defined the notion of a “spread” at ℓ given by [Image Omitted. See PDF]
[Image Omitted. See PDF]where qℓ is the ℓth row of Q and the matrix Wℓ is given by [Image Omitted. See PDF]
The quantity ρℓj represents a “distance” between coefficient indices ℓ and j, which represent not only latitude but also angular basis function. If ℓ and j correspond to the same basis function, we define ρℓj as the difference between their corresponding latitudes. However, in keeping with the philosophy of the spread as representing a measure of physical resolution only, we apply a penalty to the “distance” between indices belonging to different angular basis functions by adding a large constant (e.g., 1,000) in such cases: [Image Omitted. See PDF]
Since can be viewed as a linear functional qℓ operating on T (J), its covariance transforms as [Image Omitted. See PDF]
The fundamental trade‐off in the Backus‐Gilbert approach is between minimizing the variance of the coefficients (Equation C7) and minimizing their spread (Equation C4). The trade‐off can be quantified by introducing a parameter λ, such that we seek to determine [Image Omitted. See PDF]
The solution to Equation C8 can be determined using the method of Lagrange multipliers and is given by [Image Omitted. See PDF]
The matrix Q is constructed by solving for all ℓ, so that the “unresolved” coefficients can be computed from Equation C1.
Using the definitions of the limb darkening and the coefficients defined by Equation 16 and Equation 12, the uncertainty of the limb darkening is easily shown to be [Image Omitted. See PDF]
Assuming that the coefficients are uncorrelated and that their associated uncertainties (δci/c0, ) are equal yields an estimate of the target covariance. The MWR target value of at 45° yields an uncertainty of δci/c0 ≈ 5.4 × 10−4. The trade‐off parameter λ is adjusted iteratively to align the uncertainty with the covariance specified in Equation C7.
We would like to thank Edwin Sarkissian who helped develop and maintain the pipeline software that generates the Juno MWR antenna temperatures from raw instrument counts. This research was supported by funds from NASA distributed to the Jet Propulsion Laboratory, California Institute of Technology, the University of Michigan at Ann Arbor, the SETI institute, the Georgia Institute of Technology, the University of California at Berkeley, and the Southwest Research Institute. Fletcher was supported by a Royal Society Research Fellowship and European Research Council Consolidator Grant (under the European Union's Horizon 2020 research and innovation programme, grant agreement no. 723890) at the University of Leicester. Wind speeds in Figure 10 are based on NASA/ESA Hubble Space Telescope observations associated with program GO‐15262, with support provided by NASA through grants from the Space Telescope Science Institute (operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5‐26555).
Juno MWR data are available through the Planetary Atmospheres Node of the Planetary Data System (
© 2020. This work is published under http://creativecommons.org/licenses/by/4.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Details
; Levin, Steven 1
; Orton, Glenn 1
; Brown, Shannon T 1
; Adumitroaie, Virgil 1
; Janssen, Michael 1
; Wong, Michael H 2
; Fletcher, Leigh N 3
; Steffes, Paul 4
; Cheng, Li 5
; Gulkis, Samuel 1
; Atreya, Sushil 6
; Misra, Sidharth 1
; Bolton, Scott 7
1 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA, USA
2 SETI Institute, Mountain View, CA, USA; Department of Astronomy, University of California Berkeley, Berkeley, CA, USA
3 School of Physics and Astronomy, University of Leicester, Leicester, UK
4 Georgia Institute of Technology, Atlanta, GA, USA
5 Department of Astronomy, University of California Berkeley, Berkeley, CA, USA
6 Climate and Space Sciences and Engineering, University of Michigan, Ann Arbor, MI, USA
7 Southwest Research Institute, San Antonio, TX, USA