Morphology of the Io Plasma Torus From Juno Radio Occultations

The jovian moon Io disperses about 1 ton/s of material in the planetary magnetosphere, mainly by sublimation of SO2 from the surface and by its intense volcanic activity. The ejected material supplies the plasma cloud surrounding Jupiter known as Io Plasma Torus (IPT). The radio communication between Juno and the Earth DSN station crosses the IPT near the closest approach. Being a dispersive medium, the IPT introduces a path delay in the signal, which can be analyzed to retrieve the density distribution of electrons. We used radio tracking data from the first 25 orbits to investigate the morphology of the IPT and its variability. We adopted a static and axisymmetric model for the electron density and we updated it including temporal and longitudinal variability. We found that our best fit model must include both variabilities, even though on average the morphology of the IPT agrees with previous analyses. Our results suggest that the density of the outer region of the IPT fluctuates over 50% the average value over a typical time scale of about 420 days.


Density Models
In order to start our investigation, we adopted the empirical model for the electron density given by Phipps and Withers (2017). In their study, they considered the radial distribution of electrons given by Bagenal and Sullivan (1981) from Voyager 1 and modeled each region with a Gaussian-like distribution centered at the distance of local maximum density (they omit the 1/2 factor at the exponent and we will do the same, but we nevertheless use the shorthand Gaussian for this profile). They noticed that the warm torus cannot be well represented by a single Gaussian model, so they split this region into two pieces: the warm and the extended torus (even though their temperature and chemical composition are the same). Subsequently, they extended their model off the centrifugal equator, modeling the vertical profile of the IPT with another Gaussian for each region in order to take into account the different scale heights.
In this paper we considered the coordinates to be expressed in a left-handed centrifugal frame of reference whose symmetry axis is tilted by 6. (System III longitude). We sometimes refer to this frame as "IPT-fixed." Taking cylindrical coordinates in this frame (i.e., z axis aligned with the symmetry axis, r away from Jupiter and  E clockwise), the reference density model reads: where  1,2,3 E i stands for cold torus, ribbon, warm torus and ext stands for the extended torus. The parameters of the model are the peak density of each region N, the radial position and extension R and W and the scale height H. The model can be made more realistic by including the parameter Z which takes into account a possible offset of the torus from the centrifugal equator. An estimate of the values of the parameters and density contour plots can be found in Phipps et al. (2018). However, the above-mentioned model cannot be fully resolved by radio occultations performed by Juno because of the geometry of the observations. Indeed, the shape of the signature in the path delay introduced by the IPT strongly depends on its density gradient in the normal direction from the centrifugal plane and it goes as 2 i E H . Taking the values from Table 3 in Phipps et al. (2018), it is easy to see that the density gradient of the cold torus is about two order of magnitude larger than the one of the ribbon, the warm and the extended torus. These latter three regions on the contrary have similar thickness, thus their signatures are quite mixed up by the occultation geometry. This issue was also pointed out by Phipps et al. (2018), hence we make use of a simplified model with a single region to represent the ribbon, the warm and the extended torus, as done for the analyses of PJ 01, 03, 06, 08 10, 11, 14, and 15 (Phipps et al., 2018(Phipps et al., , 2019(Phipps et al., , 2021. We will call this region the outer torus and the whole model is given by In the Supporting Information S1, we showed that this density profile and the model presented in (Phipps et al., 2018, Appendix A) are both well suited to describe the data within the uncertainty. In turn this also confirms that we are not able to resolve the outer part of the torus in its subregions.
Both Equations 12 and 13 are axially symmetric and static in a centrifugal frame, but in general the torus parameters may be functions of time and longitude in an IPT-fixed frame. This implies that the parameters of those equations can be considered as the "zeroth-order" and that they can be expanded using a Fourier series of time and longitude. One of our main goals is to improve the axisymmetric model including both from the surface , whose neutrals are ionized by collisions with already-existing plasma and photoionization, thus an enhanced volcanic activity may affect the IPT changing its density, temperature and composition (Delamere et al., 2004;Hikida et al., 2020;Koga et al., 2019). Nevertheless, the correlation between Io's volcanism and IPT is very thorny. First of all, the brightness of a volcanic events is not necessarily proportional to its emission of gas and dust. Second, Io hosts different types of volcanic activity such as caldera-like formations, flow-dominated eruptions, explosive events and plumes (Davies, 2001;Lopes & Williams, 2015). Each of them is more or less regular and exhibits different magma, gas and dust supply. Third, the altitude, longitude and latitude of an ejection affect the interaction between the outgassing and the torus, determining if new material is ionized or if it fall back on Io's surface (McDoniel et al., 2019). As a last warning, the response of the IPT to plasma injection of different composition, temperature and density may be complex and it is still debated and under observation Yoshikawa et al., 2017;Yoshioka et al., 2018). Thus any agreement with the periodicity of Loki Patera or typical orbital periods of Io should be taken as a mere indication. With these cautions, we tested a model with T as a free parameter to be determined.
In this study, we used a static axisymmetric profile of the density as reference model and we tested the model obtained by setting  1 . This is the lowest order of the Fourier expansion that includes both temporal and longitudinal variations and we will refer to this model as mixed model.

Analytical TEC Model
In Phipps and Withers (2017), the occultation geometry undergoes two simplifications in order to calculate the TEC. First, the line of sight between Juno and the DSN station was assumed to be parallel to the centrifugal equator, which is reasonable because of the relatively small angle of the centrifugal equator with respect to the ecliptic plane. Nevertheless, the TEC obtained by integration along such lines of sight is insensitive to the radial parameters R and W. Second, they considered the projection of the trajectory of Juno on a meridional plane in order to integrate the density (see, for example, Figure 1 of Phipps et al., 2018). Actu-7 of 25 ally, the path of the radio signal is azimuthally slanted with respect to the surface of the IPT, thus it crosses different longitudes at each instant. Therefore in principle we must take into account the full 3D geometry of the occultations (Phipps et al., 2018(Phipps et al., , 2019(Phipps et al., , 2021. Nonetheless, we show that it is possible to reduce the problem in two dimension exploiting the close proximity of Juno to Jupiter introducing only a small error. Detailed calculation on how such reduction was performed can be found in the Supporting Information S1. First, we performed the projection of the S/C position on a proper meridional plane. The plane we are considering is perpendicular to the centrifugal plane and passes through the G/S position at each instant and the center of mass of Jupiter (see Figure S1). We then evaluated the projection of the S/C position onto this plane and we made use of this approximate trajectory of Juno to carry out the integration. We evaluated the error introduced by such approximation by using the IPT described by Equation 12 and simulating the occultations performed by Juno. For each occultations we found that the relative difference between the TEC computed with the full 3D geometry and the one obtained with the approximation is usually less than 1% and almost never greater than 5%, except far from the occultation, where  0 E TEC , thus making the relative error large even though the absolute TEC difference is small. Given that the uncertainty we estimated for our data set is always greater than 5% of the minimum path delay, this has negligible consequence on the parameter estimation.
Second, we integrated the density profile in Equation 13 using the parametrization where  E is the radial coordinate in the meridional plane away from Jupiter, E m is the angular coefficient of the line of sight between S/C and G/S and E q is the intercept with the z-axis of the IPT-fixed frame. In the end, we were able to obtain an analytic expression for the TEC: b a E f f b f a and sc E R and gs E R are the radial position of the S/C and the G/S respectively. Here  E is the line of sight between the projected S/C position and the G/S. This way we can plug Equation 15 into Equation 4 and obtain a relation between the observed path delay and the TEC along the line of sight between Juno and the DSN station.

Parameter Estimation
Two different technique were used to estimate the typical parameters of the IPT in Equation 13: the single arc fit and the global fit. In the first case, each path delay signature is fitted individually with the delay obtained from Equations 4 and 15, where E m and  E are computed from the position of the S/C and the DSN station as function of time. Using the single arc technique each PJ is characterized by its own set of parameters, which is different from one another. Example of this approach can be seen in Phipps et al. (2019Phipps et al. ( , 2021. The main advantages of this method are its simplicity (there are few parameters to be estimated) and that the resulting residuals are within the uncertainty of the data (see ). Besides, the parameters R and Z in Equations 13 and 15 are geometrically coupled as can be seen schematically in Figure 2. The global fit approach instead tries to fit all the data using a single set of parameters to characterize the IPT. This technique makes the parameters R and MOIRANO ET AL. 10.1029/2021JA029190 8 of 25 Z in Equations 13 and 15 to be decoupled if the line of sight between Juno and the DSN station has different tilt from one PJ to the other. Furthermore, it allowed us to exploit all the timespan and the sector coverage (see Section 2.5) to include temporal and longitudinal variability in the model. Indeed, the model used for a global fit must include all the effects that govern the variability of the IPT, many of which are still unknown. In turns, including all such effects lead to very complex models which may make the fitting algorithm to not converge due to the nonlinear nature of the model. In this study, we show results from both approaches. Since R and Z are coupled, we first need an estimate of 1 R E and 2 R E for the single arc fit using Equation 15. We take these values from Voyager 1 data Bagenal and Sullivan (1981) . These values are obtained from Voyager 1 data (Bagenal & Sullivan, 1981) and the above-mentioned gaussian fit we used to estimate 2 E R for the single arc fits.
During each occultation the radio signal crosses three main plasma sources which affect the path delay: the Earth's ionosphere, the IPM and the IPT. Thus, in order to fit the signature of the IPT, we first need to remove the contribution due to the other two sources. In addition, the DSS-25 of the DSN showed a differential delay between the X and Ka-band during uplink, which strongly affect the path delay estimated by radio tracking data in a dual uplink dual downlink configuration (see Figure 7 of Zannoni, 2020). In this study, we will use an estimate of this delay performed at JPL in November 2018 .
The ionospheric delay is removed by means of GPS data (Thornton & Border, 2000), thus our main concern is removing the plasma contribution due to the IPM, whose main source is the solar wind. Its high variability (Ebert et al., 2014;Matthaeus et al., 1991) makes difficult to determine its impact on the delay. If we consider the density of the solar wind to be a decreasing function only of the radius from the Sun (e.g.,: Köhnlein, 1996), the time derivative of the TEC depends only on the relative motion of the S/C and the G/S. For the time interval of a single PJ (about 6-8 hr), this means that we could, in principle, fit the solar wind contribution with a straight line whose slope depends on the relative velocity of S/C and G/S.

of 25
Actually, the path delay far from the signature of the IPT shows other signatures which cannot be ascribed to a linear drift due to G/S-S/C relative motion in a stationary solar wind. The fluctuations observed in the path delay may be due to specific conditions of the electron content along the line of sight between Juno and the DSN station that are not measurable. We empirically removed this effect by fitting a nine-degree polynomial as already done in previous data analysis of Juno radio occultations (Phipps et al., 2019(Phipps et al., , 2021. We substituted the signature of the IPT with a linear interpolation between the beginning and the end of it and we used this path delay to remove the background. The linear interpolation stabilizes the fit, while leaving missing data can make the fit to exhibit strong fluctuations. In Figure 3, an example of the background removal during PJ17 can be inspected. In this way we calibrated for the solar wind contribution taking into account an unknown variability, maybe due to different periods of solar activity, longitudinal asymmetries in the solar wind or time-specific events. Once we get rid of the background contributions to the path delay, we use a Markov Chain Monte Carlo algorithm coded in MATLAB to fit the data (Haario et al., 2001(Haario et al., , 2006. This code combines the delayed rejection algorithm (DR, which derives from the Metropolis-Hastings algorithm) with the adaptive metropolis algorithm (AM) and it is thus referred as DRAM (Haario et al., 2006).
For the single arc fits we started from a uniform distribution for each parameter. We chose and to keep the model as simple as possible, thus we left a larger space for the parameters of this region to allow them to change more freely. The first and second order coefficients in Equation 14 are normalized to the zero order coefficient and their boundary are set at −0.5 and 0.5. The period T in Equation 14 is left as a free parameter with initial value T = 450 days, which is taken approximately from the result of de Kleer et al. (2019). The boundary are based on the time spanned by our data set: the lower limit is given by twice the maximum time between two subsequent PJs (about 300 days) and half of the time between PJ01 and PJ25 (about 630 days). Summaries of initial values and boundaries used for the global fits can be found in Table 2. To evaluate the weight to use for our measurements, we consider three different contributions: solar wind scintillation, ionosphere and background fluctuations. The uncertainty due to scintillation of the solar wind is given by numerical integration of the uncertainty in the Doppler shift in the X-band and it is a function of the SEP angle (Sun-Earth-Probe) (Iess et al., 2014). The error on the TEC due to the terrestrial ionosphere depends on the elevation of the spacecraft and it ranges from 3 TECU at high elevation to 8 TECU below  10 E . Since we did not included the elevation in computing the line of sight between Juno and DSN station, in this study we adopted the average uncertainty  iono E TEC = 5TECU (Thornton & Border, 2000). The background fluctuations are more evident in some PJs than others. These may be caused by anomalous conditions in the solar wind which may be time-specific and localized. In order to take into account these fluctuations we use the RMS of the residuals obtained from a linear fit to the background plasma. We did not use the 9th degree polynomial used for the background removal because it absorbs also the fluctuations, thus the residuals are very similar among the occultations and so the resulting weights.

Data
Our investigation of the IPT uses PJs dedicated to the gravity experiment. Indeed the closest the flyby is to Jupiter, the more relevant the gravity field becomes on the trajectory of the spacecraft. Besides, given the morphology of the IPT, the gravity experiment and the occultations of the IPT occur in the same time windows. So far, we have made use of 15 PJs (Table 1), which allowed us to span more than 3 years since August 2016 (PJ01) to February 2020 (PJ25). Plots of the signature of the IPT in the path delay data can be seen in Figure 5. The gray area represents the uncertainty estimated as explained in the end of Section 2.4. As can be noticed, the minimum value of the delay after the background removal varies between −300 and −600 mm, which points out strong variations in the TEC of the IPT over the three years of observations. Because Juno is performing occultations every  E 53 days and that each occultation spans a limited longitudinal sector, such variability can be due to longitudinal asymmetry and/or temporal periodicities.
The path delay of PJ01 and 13 in Figure 5 appears shallower compared to the other since they were measured using a X/X-X/Ka link, so they are due only to the plasma crossed in downlink. If the tracking during these PJs had been performed using a X/X-Ka/Ka link, their signature would have been about 2.25 as deep.
The dates (yyyy-mm-dd) refer to the day when the tracking started. b The times are referred to UTC DSN station time ( 3 E t in Figure 1). c The longitude is referred to System III. d Each occultation is     (30 40 ) E from the peak density.
. This factor of has been included in the fitting algorithm, so that we fitted all the data as they were acquired using only the X/X-Ka/Ka link. Nevertheless, in Section 3 we report the actual path delay extracted using the multifrequency technique (Bertotti et al., 1993;Mariotti & Tortora, 2013), therefore the signatures of PJ01 and 13 show only the downlink plasma contribution.
Even though the polar orbits of Juno change gradually in longitude, so far the radio signals have spanned many sectors of the IPT, as can be seen in Figure 4.

Results
In this section, we present the result obtained using the single arc approach and the global fit described in Section 2.4. Establishing the convergence of a Monte Carlo chain is still an open problem, as there is no unique way to determine if the chain has converged to the target distribution. Only a handful of diagnostic Note. The columns 2-4 refer to the mixed model and they represent the initial value (init.), the upper and lower boundaries (UB and LB), the best fit value from the MCMC algorithm with the 1- E uncertainty (best) respectively. The columns 5-7 are the same, but for the axisymmetric model. a The uncertainties were adapted after the stability test.  Cowles & Carlin, 1996). In this study, we considered a small autocorrelation of the chain as diagnostic for the convergence.
The equilibrium parameter distributions we obtained are all univariate and they are gaussian-shaped in nearly all the cases.
In order to quantify the difference between the axisymmetric and other models, we reported the difference between the square root of the mean weighted sum of squared residuals (MWSSR) of the axisymmetric model and the MWSSR of the mixed model. These values can be found alongside the plots in Figure 7. The MWSSR is computed as where i E W , i E O and i E E are the weight, the observed value and the expected value (from the fit) of the i-th data point, while  E L N are the degrees of freedom (i.e., number of data points minus number of parameters). The global solutions undergo a stability test in order to verify the uncertainties (whose details can be found in the Supporting Information S1). We made use of such test to correctly evaluate the uncertainty of our results.

Single Arc Fits
In this section, we present the fit we obtained with the single arc approach described in Section 2.4. The single arc approach is not able to decouple the radial position of the two region of the IPT (   sigmas. Such differences may be due different method for the removal of the background. Indeed, the fiducial range of each occultation we choose may be different from other analyses. Besides, we substituted the signature of the IPT with a straight line before removing the background, which can also affect partially the resulting data. Our results for the offset differs from the of Phipps et al. (2019Phipps et al. ( , 2021 because of the different frame: we adopted a centered dipole tilted by 6.  8 E , while they used an offset dipole tilted by 6.  3 E . Taking into account this difference (which accounts up to ∼0.07 J E R depending on the longitude), the results are compatible within one sigma for the cold torus, while for the outer torus the offsets are in agreement within three sigmas. The only exception is PJ 10, when the offset we obtained is different from the reference analysis (Phipps et al., 2021) by more the three sigmas.

Global Fits
The results reported in the remainder of this section are obtained using the global approach described in Section 2.4. Nevertheless, we retained the offsets 1 E Z and 2 E Z as local parameters, since the longitudinal and temporal variability that we will introduce in Section 3.2.2 involve only the density and the scale height. This choice helped to align the fitted path delay with the signature in the data, which is particularly relevant for the cold torus: indeed, its size (i.e., 1 E H ) can be similar to its offset, therefore neglecting the variability of the offset shown in Figure 6 may strongly affect the determination of the density and scale height of the cold torus. Even though this effect should be less significant for the outer torus, it can still affect the shape of this region, especially for large offset as the ones observed during PJ 01, 21, and 22.

Axisymmetric Model
The axisymmetric model fits the data in Figure 5 using Equation 15, where all the parameters are constant over time and longitude. The results in this section were used as a reference and as input for the mixed models. As can be seen by the plots in Figure 7, such model cannot reproduce well all the observed data. In some occultations, the difference between the observed and fitted path delay is as large as 100 mm (e.g., PJ08,11,18,21,23,and 25). In addition, the width of the signature is very poorly represented by this model in multiple occasions (e.g., PJ14,18,21,22,and 25). Nevertheless, the axisymmetric model can be viewed as an estimate of how the IPT looks like on average during the Juno mission.
As explained in Section 2.3, the global approach allowed us to retrieve the radial position of the IPT by exploiting the little tilt of the S/C-G/S line of sight with respect to the centrifugal plane.
and the same holds for 2 E H . The total density and scale height need to be bounded in order to avoid potential negative values that might artificially result from the fit. In order to include this constraint we adopted a penalty function (Kuri-Morales & Gutiérrez-García, 2002) in our calculations. Therefore any set of parameters explored by the MCMC algorithm that does not satisfy the constraint is graded much less poorly than anyone that provides positive density and scale height. Additional details on this penalty function can be found in the Supporting Information S1. and The mixed model then consists of a longitudinal modulation at fixed longitude (second term in Equation 18), a temporal pulsation (third term) and a longitudinal modulation whose amplitude and phase vary over time (fourth term).
The introduction of both longitudinal and temporal variations in our model led to a remarkable improvement of the fit, as can be noticed by looking at Figure 7. All the PJs except PJ08, 10, and 17 are better fitted using this model compared to the other global fits. This is particularly evident for PJ11, 13, 14, 22, 23, and 25. The overall MWSSR improved by 40% compared with the reference model.
The value of the zero order parameters in Table 2  . The corrections of the first order are quite important, as density fluctuations can be as high as 0.5 times the average, while the scale height can vary by about 0.3 times the average (taking into account both temporal and longitudinal variability). The amplitude we obtained at the first order showed that longitudinal modulation are more important than temporal one at the first order for the density, but they are both relevant for the scale height. Indeed, we obtained that H . This is not in contrast with the fact that we observed PJs with large difference in the minimum path delay at nearly the same longitude. First, we still need to discuss the role of the second order correction. Second, those observations cannot be explained by a purely longitudinal modulation and require a temporal variability of the IPT. At the same time, temporal variability does not rule out longitudinal variation, which can then be retrieved if the model also includes pulsations and/or other time dependant changes in the IPT.
In Figure 8, we plotted the second order amplitude (  Figure 8), we can notice that the peak of the density fluctuation occurs about 50 days later than the peak in scale height. Comparing the maxima of both these curve with the corresponding first order amplitude in Table 3, it is evident that the second order contributions can be as relevant as the first order ones. The phase difference between density and scale height at the second order (bottom panel in Figure 8)  , as can be seen in the bottom panel in Figure 8. The phase Equation 21 in the bottom plot of Figure 8 slowly drifts periodically from ∼  100 , which points out that the density and the scale height are almost anticorrelated.
In the top panels of Figure 9, we showed the color-coded density and scale height computed using Equation 18 with the retrieved parameters (the vertical axis spans a single period of 432 days). As can be seen, the mixed model revealed high temporal and longitudinal fluctuations. In the middle panels of Figure 9 we showed the density and scale height relative fluctuation as functions of time at specific longitudes  and  350 E t days). The anticorrelation seems to occur mostly when the amplitude of the scale height fluctuations is higher.
The period we obtained with the mixed model is   432 2 E T days. This is about 5% less than the periodicity of the volcanic activity of Loki Patera retrieved by de Kleer et al. (2019) and about 6% and 10% less than the evolution timescales for the semimajor axis and eccentricity of Io's orbit, but nevertheless it can be considered as in indication that the IPT undergoes temporal variation on a similar timescale.

Discussion
From the single arc fits we found that each set of parameter greatly varies from one occultation to the other, which point out that the IPT undergoes substantial changes over few weeks. Indeed, the density of the cold torus varies considerably from PJ03 ( . The mean density of the cold torus from the single arc approach was 1,800 3 cm E (weighted with the uncertainty), ranging mostly between 1,500 3 cm E and 25,00 3 cm E . This value is in good agreement with the electron density at about 5.2 R J E from re-analysis of Voyager 1 data (Bagenal, Dougherty, et al., 2017) and about 800 3 cm E higher than the density measured by Galileo (Bagenal et al., 1997). The cold torus exhibit strong enhancements, like during PJ 01, 08, and 21, and depletions, like PJ 03 and 13. The outer torus shows an average density of 2,500 3 cm E ranging mostly between 2,000 3 cm E and 3,200 3 cm E , except during PJ11 and PJ25, when the outer torus appeared much less dense (about 1300 3 cm E ). Such difference during those occultations is evident also by looking at the path MOIRANO ET AL.

10.1029/2021JA029190
19 of 25 delay in Figure 5. The average density of the outer region from the global approach is also consistent with data from Voyager 1 (Bagenal & Sullivan, 1981), but about 600 3 cm E lower than the data from Galileo (Bagenal et al., 1997). This may point out that the density of outer region of the IPT is more stable around 2,500 3 cm E , but occasional enhancement or depletion can occur.
Some of the PJs occurred at similar longitudes (see Figure 4), therefore we were able to observe that PJ03 and 13 exhibited very depleted cold torus during both occultations, which occurred at about  70 E . On the other side, there are evidence of density variation in the cold torus, as can be seeing by comparing PJ11 and 18, which occurred at about  250 E . The outer torus also showed different peak density values at similar longitudes, as can be noticed by comparing PJ11-18 and PJ06-23. This suggest that density variations of the IPT are time-dependant, although this does not rule out the possibility of longitudinal modulations. This is in agreement with the different electron density measured by Voyager 1 with respect to Galileo (Thomas et al., 2004).
In our models, we did not take into account the presence of the ribbon between the cold and the outer torus. This region has limited radial extent (<0.1 R J E ) and high density (3,000-4,000 3 cm E ), although this region showed longitudinal and temporal variability in both density and temperature (Thomas et al., 2004). Besides, this feature does not seems to be always present (Bagenal et al., 1997). The vertical extent of the ribbon is comparable with the outer torus (Bagenal, 1994;Phipps et al., 2018;   Phipps & Withers, 2017), thus its signature is mixed with the signature of the outer regions. In order to approximately estimate the relative contribution of the three regions of the outer torus in Equation 12, we can integrate this equation along lines parallel to the centrifugal equator, which is equivalent to Equation 17 in Phipps and Withers (2017). Taking the values of Voyager 1 from Table 2 of Phipps and Withers (2017) as a reference, the relative contributions to the TEC of the outer torus are approximately 10%, 32%, and 58% for the ribbon, warm torus and extended torus respectively. Although a single gaussian is not well suited to fit the radial profile of the electron density compared to a three-gaussian model, as showed in Figure S3, the best fit values in both cases would lead to a similar TEC for the occulting geometry of Juno. Therefore, the above-mentioned proportions we estimated from Voyager can be used to approximately retrieve the density of the ribbon, the warm and the extended torii as follow: first, we estimate the TEC using   Both results are similar to the values found from previous analysis of the first PJs (Phipps et al., 2019). Nevertheless, the results from the global fit in Figure 7 showed poor agreement with the data. In particular, during PJ 14, 23, and 25 the scale height, which governs the width of the signature, overestimates the length of the occultation, while during PJ18, 21, and 22 it is underestimated.
Using the global fit we were able to include the radial positions where the peak density of the ribbon and the warm torus are expected to occur. In our model, we did not take into account the offset of the magnetic equator with respect to the center of mass, which is given by System III longitude (Connerney, 1993;Dessler, 1983). Including this offset, the torus is displaced by about  0.1 J E R along the offset vector. Nevertheless, the inclusion of such offset doesn't change substantially our results.
The offsets retrieved from the global fit with the axisymmetric model agree with the results from the single arc fit within the uncertainty.
The mixed model, which includes both temporal and longitudinal variability, showed the residuals were improved by about 40%. The predicted amplitude of the density fluctuations is large enough to explain the observed variability of the path delay, whose minima range between −300 and −600 mm. Indeed, we can see that all the values for 2 E N and 2 E H obtained with the single arc approach in Figure 6 fall within the extrema predicted by the mixed model in the top panels of Figure 9. The density of the cold torus we retrieved from both the axisymmetric and the mixed model agrees with previous estimate based on Juno data (Phipps et al., 2018(Phipps et al., , 2019(Phipps et al., , 2021, but they appears higher than the results from Galileo (Bagenal et al., 1997). The average density of the outer torus is also very similar to past in-situ observation of the IPT (Bagenal et al., 1997;Bagenal & Sullivan, 1981), but we retrieved wide fluctuations that can reach up to 70% of the average. Cassini recorded longitudinal variations in the electron density of about 20% (Steffl et al., 2006), while volcanic active periods can lead to a density increase by about 20%-25% (see for example Yoshioka et al., 2018). Our measurements are sensitive to the TEC, therefore it is possible that the fluctuations we retrieved are partially due to variation of the radial extension of the IPT, as the path delay depends on the TEC.
As a rough approximation, we expect the density and scale height to be nearly anticorrelated because they are related to the temperature by The mixed model is the only model we tested which reproduced this feature, at least when the amplitude of the corrections is quite large. Nevertheless, the anticorrelation should not be taken too strictly. The relation between density and temperature  E NT constant holds if a polytropic equation of state is a proper choice for the plasma in the IPT, which relies on the assumption that density and temperature variations take place at constant pressure. Besides, the scale height depends on the temperature anisotropy (i.e., the ratio between perpendicular and parallel temperature of the plasma), on the parallel temperature and on the ambipolar electric field (Thomas, 1992), which in turn depends on the ion composition. While we found no indications that the anisotropy varies substantially over time, comparison between Voyager 1, Cassini and Galileo seems to point out that the ion composition can change at different epochs (Nerney et al., 2017). If such changes occur almost simultaneously with temperature variations, the resulting scale height is caused by both phenomena and not just by temperature variation. In the end, nonstationary injection of gas from Io into the torus may lead to transient states during which the anticorrelation between density and scale height is not mandatory. The anticorrelation should hold more strictly when the supply of plasma is quite steady and the ion composition does not change appreciably.
The best fit period for the mixed model is   432 2 E T days. We re-run the fit with the mixed model changing the initial value of T to check the stability of this result. We chose T = 400, 500, and 550 days and MOIRANO ET AL. days. This period is quite similar to the periodicity of the activity of Loki Patera, which is taken as proxy for the overall volcanic activity on Io (de Kleer et al., 2019;Tsuchiya et al., 2019). In Figure 10, the infrared emission of Loki Patera is plotted as function of time alongside the average density of the outer torus predicted by the mixed model. Far from establishing a definitive connection between the morphology of the IPT and the volcanism of Io, nonetheless the similarity between our result and the periodicity of Loki is quite remarkable. In addition, Hisaki reported increased volcanic activity during the periods March-April 2015 and May-June 2016 . These events were observed alongside the decrease of the rotation period of IPT and increase in thermal electron temperature. Furthermore, the increased activity observed in the early 2015 was concurrent with an increased hot electron fraction (Hikida et al., 2020), which can affect the ionization of the main elements in the IPT (Steffl et al., 2006) and potentially its electron content as a consequence. Yoshioka et al. (2018) retrieved the radial profile of the electron and ion densities that matched the spectrum observed by the EXCEED instrument, onboard Hisaki, in the same period, showing that the electron peak density at ∼6R J E can increase up to  E 20%. The above-mentioned period observed by Hisaki occurred after the peak emissions from Loki Patera, as can be seen in Figure 10. Hikida et al. (2020) and Tsuchiya et al. (2019) suggested that Kurdalagon Patera and Pillian Patera may also be involved in the enhanced volcanic activity detected in 2015. Comparing the epochs observed by Hisaki, the activity of Loki Patera and the present prediction, we noticed that the electron density of the IPT starts to increase after that the emission from Loki increased and that a second minor peak is reached after the periods of enhanced activity recorded by Hisaki. Our results seems to suggest a potential correlation of the electron content of the IPT with the periods of volcanic activity, observed from September 2013 to June 2018. Nevertheless, the dynamics response of the IPT to the time-variable volcanic activity of Io is not straightforward, hence conclusions should be drawn carefully, as we discussed in Section 2.2. Adding future occultations and including direct observations of volcanic activity might strengthen this indication in the future.
The improvement of the residuals, the anticorrelation between density and scale height and the potential agreement between our period and the volcanic activity place the mixed model as the best candidate among our models to describe the variability of the IPT during the Juno mission.

Conclusions
We made use of the path delay measured during 15 occultations of the Io plasma torus to inspect its morphology and variability, which was investigated by a Fourier expansion of a parametric model of the density distribution (Phipps & Withers, 2017). We assumed the IPT to be made by two regions: a cold inner torus and an outer warm torus, both of which were modeled by gaussian profiles in the radial and vertical direction (in cylindrical coordinates in a frame of reference tilted by 6. ).
We integrated the density profile of the IPT along the line of sight between Juno and the DSN station and we obtain an analytical formula for the TEC, which is proportional to the path delay. The integration was carried out in an approximated 2-dimensional geometry, which introduced a small error (compared to the uncertainty) with respect to a full 3D integration. The result can be efficiently employed on iterative algorithms, such the Markov Chain Monte Carlo algorithm we used in this work (Haario et al., 2001(Haario et al., , 2006. The equation we obtained takes into account the tilt of the line of sight with respect to the centrifugal equator, thus it can be used to retrieve the radial position and extent of the IPT, at least in principle. Actually, the radial extension and the peak density are strongly coupled by the geometry of the occultations, thus we needed to fix the radial extension ( 1 E W and 2 E W in Equations 13 and 15).
We analyzed our data using two approaches: the single arc fit and the global fit. From the former we obtained peak density, vertical extension and offset from the centrifugal equator for each occultation. We needed to fix the radial position when using the single arc approach because the radial position and the vertical offset were strongly coupled by the geometry (this coupling is removed in the global fit). This method is the same as used in previous analysis of the IPT morphology during the Juno mission (Phipps et al., 2018(Phipps et al., , 2019(Phipps et al., , 2021. The densities we obtained are in agreement with previous analyses within three sigma for the cold torus and within 1 sigma for the outer torus. Deviations of the average values between previous and this studies may be due to the different fit to the background signal. Indeed, we enforce the background removal to follow a straight line in place of the signature of the IPT, which is different from what used by Phipps et al. (2019Phipps et al. ( , 2021. This in turn may affect the depth of the signature. We also noticed that the offset deviates by about   1 E with respect to the position predicted by a nominal centrifugal equator based on a dipolar magnetic field. As highlighted by Phipps et al. (2020), this suggest that the inclusion of higher-than-dipole moments of the magnetic field and the plasma disk should be taken into account to compute the centrifugal equator. Besides, the latitudinal distribution of plasma between 5 and 10 R J E is a function of the radial distance (Phipps et al., 2020), which should also be considered when evaluating the offset of the IPT.
The residuals obtained from te mixed model are about 40% better than the reference model. Besides, the density and scale height predicted by this model are mostly anticorrelated, even though correlation is still possible when the density and scale height show little fluctuations. This feature is appreciated, as we expect the density to be inversely proportional to the temperature, while the scale height to be proportional to the square root of the temperature. Lastly, the characteristic period we found of about 430 days is somewhat close to the periodicity in volcanic activity and orbital changes of Io recently investigated (de Kleer et al., 2019). Even though we do want to point out any hard evidence between volcanism on Io and the response of the IPT to the mass supply (which requires ad hoc theoretical modeling of the Io-IPT interaction and monitoring of the activity on Io, both of which are beyond the scope of our work), at least we take this as an indication of such potential interaction.
Without delving deeply into the physical interpretation of the results we obtained, we showed that both longitudinal and temporal variations of plasma in the IPT are likely to occur, which should be kept in mind in planning and analysing long-term missions in the Jupiter's inner magnetosphere.

Data Availability Statement
Spacecraft and DSN station data can be retrieved at https://naif.jpl.nasa.gov/naif/data.html. The Doppler data and ancillary information used in this analysis are archived in NASA's Planetary Data System (Buccino, 2016).