In Fig. 1, we superimpose a steady, axisymmetric, zonal flow profile on a background map of the magnetic field^{2} derived from Juno magnetic field observations^{11} from the spacecraft’s first 33 orbits. The flow is dominated by an equatorial jet, which induces intense secular variation in the vicinity of the Great Blue Spot (the region of concentrated field at the equator) as the magnetic field associated with this spot is swept eastwards. Owing to its dominant role in generating the secular variation^{1,2,3,12}, a recent set of orbits by the Juno spacecraft^{13} was targeted at this region.

To begin, we produce a new model including the magnetic field observations from these targeted passes (and other subsequent orbits over other regions of the planet). The model is produced using the same method as the model in Fig. 1 (Methods and ref. ^{2}). One pass, PJ02 (in which PJ stands for perijove), did not acquire any data, so the number of data-yielding orbits is 41 compared with 32 orbits for the earlier model; note we refer to these models in terms of the last orbit used, that is, the 33-orbit model (Fig. 1) and the 42-orbit model. The resulting 42-orbit model has a global misfit of 492 nT compared with 411 nT for the 33-orbit model (for comparison, the root-mean-square (r.m.s.) field strength of the observations is 282,000 nT); within a box around the spot (Fig. 2) the misfit is 934 nT compared with 675 nT (where the r.m.s. field strength is 393,000 nT) and the maximum speed of the equatorial jet is 0.64 cm s^{−1} compared with 0.86 cm s^{−1}. Thus, the fits we obtain to the 42-orbit dataset are poorer than those to the earlier 33-orbit dataset, especially near the spot, indicating that a steady flow performs worse as the time interval spanned by the passes increases. We may have expected, instead, that with the addition of these later passes that the misfit would decrease because these passes are at higher altitude over the spot and hence sample weaker field. Except for the southern hemisphere south of 30 °S, where the flow resolution is poor^{2}, the flow profiles are broadly similar; however, the equatorial jet speed is reduced by 26% in the 42-orbit solution, suggesting that the flow may be changing with time. The pattern of residuals in Fig. 2a lends additional support to this possibility: we can identify pairs of passes that are spatially adjacent but separated in time that have oppositely signed residuals over the spot, notably PJ19 and PJ36, and PJ24 and PJ38. Oppositely signed residuals will result for adjacent passes if the actual flow speed at the time of the passes is greater than the steady flow solution for one pass and smaller for the other.

To examine the possibility that the flow speed is varying, we allow the flow to vary in amplitude on a pass-by-pass basis. We do this by applying a velocity scale factor to the flow for each pass (Methods). The velocity scale factor does not change the flow profile, instead it simply scales its amplitude. By doing so, we find the adjusted flow speed that gives the best fit for a particular pass, but for a different pass that flow speed will probably be different. These adjusted flow speeds represent the average flow speed from the baseline epoch of 2016.5 for each particular pass. In Fig. 2b we show the residuals after applying velocity scale factors to each pass. The residuals are reduced, especially for passes to the west of the spot. The misfit within the box is 721 nT, a variance reduction of 40% from the steady flow solution. This variance reduction can be considered as the maximum that can be achieved simply by varying the flow speed. However, this variation is only physically reasonable if we can find a time-varying flow consistent with the pass-by-pass velocity scale factors, in other words a time-varying flow that yields the corresponding average flow speed for each pass. It is possible, instead, that the different velocity scale factors (or average flow speeds) are mutually inconsistent.

We examine whether such a flow exists by fitting the pass-by-pass velocity scale factors with a simple sinusoidally varying flow model with a single period and no damping (Methods). We omit PJ01 from this analysis as that orbit passes over the spot less than 2 months after the baseline epoch and thus is insensitive to variations in the flow (the flow would advect the spot by less than 0.05° during those 2 months). The best-fit solution is shown in Figs. 2c and 3: it has a period of 3.8 years and results in a variance reduction within the box of 24.8%. As expected, the variance reduction on a pass-by-pass basis varies substantially (Fig. 3b), as those passes with velocity scale factors that differ substantially from unity will have their fit enhanced more than a pass with a factor close to unity. Note that Fig. 3 shows the residuals to the radial component of the field, rather to the three components of the magnetic field, as the radial component is more readily interpreted in terms of changes in the flow speed. In a few cases, though, other components of the field show a much larger reduction in misfit than the radial component, most particularly *B*_{ϕ} (the east component of the magnetic field) for PJ24. In other words, there is not necessarily a one-to-one correspondence between the residuals in Fig. 2 and the variance reductions in Fig. 3. Comparing Fig. 2a with 2c, we can see that the residuals of the pairs of passes discussed earlier (PJ19 and 36, and 24 and 38) are much reduced. For most passes, the red bars in Fig. 3b (the normalized misfits to the sinusoidal model) are below the grey line corresponding to 1 (the normalized misfit of the 42-orbit steady flow model), but two passes (PJ26 and PJ37) stand well-above the grey line indicating that they are fit worse by the sinusoidal model than by the 42-orbit steady flow model. These two passes are the most easterly passes within the box. PJ37 requires a flow speed almost 15% more rapid than that of PJ36 and PJ38, which though temporally adjacent to PJ37 are not spatially adjacent to PJ37, indicating that additional spatial complexity in the flow may be required. PJ26 is, instead, fit by a slower flow than the sinusoidal model arguing instead for additional temporal complexity. Additional complexity could take the form of more than one wave being present or wave damping. In case our results are skewed by these two passes, we repeat the sinusoidal fit omitting them, as shown by the light red curve in Fig. 3a. The fit to most of the remaining passes, in particular PJ24 and the targeted passes (PJ36, PJ38, PJ39, PJ41 and PJ42) is improved. The period of the sinusoidal fit is changed by only a small amount from 3.8 to 4.1 years.

The period of roughly 4 years suggests that this is a torsional oscillation or Alfvén wave rather than, for example, a MAC (magnetic-Archimedean-Coriolis) wave^{14}, which would have a much longer period. Torsional oscillations have also been proposed as the origin of cloud level variability in Jupiter on subdecadal timescales^{15}: the zonal shear associated with a torsional oscillation may modulate the heat flux from the deep interior, which may in turn result in variability of observed infrared emissions at cloud level. The wave speed of torsional oscillations is determined by the r.m.s. value of the component of the magnetic field, \({\bar{B}}_{{\rm{s}}}\), perpendicular to the rotation axis^{10} (where the average is taken over longitude and the latitude band of interest). For an equatorial belt of ±10° about the equator (the latitudinal extent of the deep equatorial jet), we find \({\bar{B}}_{{\rm{s}}}=0.6\) mT at 0.9 *R*_{J}. This corresponds to an Alfvén wave speed of 10^{−2} ms^{−1}.

The period of the oscillation depends, of course, on its wavenumber *k*, for which we have no direct observation. If the cloud level variability is due to torsional oscillations, then the wavenumber can be estimated from the length scale of those variations, yielding dimensionless wavenumbers *k**R*_{J}/2π in the range 10 to 15 (ref. ^{15}). Here, however, we are examining a single equatorial fluctuation rather than a set of torsional oscillations spanning a wide range of latitudes. For the equatorial jet, a dimensionless wavenumber of 10 could be considered (although this would be based on the azimuthal extent of the jet rather than its wavenumber in the *s* direction) yielding a period of roughly 15 years: that is, four times longer than that found here.

However, our estimate of \({\bar{B}}_{{\rm{s}}}\) may be too small: the field is most probably stronger at depths below 0.9 *R*_{J}, but the field below that depth cannot be reliably estimated from the externally observed potential field owing to the rapid increase of electrical conductivity with depth^{16}; and second, intense, small scale magnetic fields (which will be geometrically attenuated in the observations at satellite altitude) may serve to increase \({\bar{B}}_{{\rm{s}}}\) further.

A period of 4 years corresponds to a field strength \({\bar{B}}_{{\rm{s}}}\approx 3\) mT, similar to the field strength associated with the spot itself, so the wave may instead be a localized Alfvén wave propagating along the field lines associated with the spot (which are largely in the *s* direction), rather than an axisymmetric torsional oscillation, in which case a superimposed longer period torsional oscillation may then also be excited.