Anelastic torsional oscillations in Jupiter's metallic hydrogen region

We consider torsional Alfven waves which may be excited in Jupiter's metallic hydrogen region. The axisymmetric zonal flow fluctuations have previously been examined for incompressible fluids in the context of Earth's liquid iron core. Integrated theoretical models for the deep-seated Jovian dynamo, implementing radial changes of density and electrical conductivity from the molecular hydrogen envelope, have reproduced its strong, dipolar magnetic field. Adopting these models, we find anelastic torsional waves which travel in cylindrical radius on timescales of several years, or longer, in the metallic region. Being excited by the convection that is vigorous in the outer layers of the dynamo region, they can propagate both inwards and outwards and be reflected at a magnetic tangent cylinder, enabling a standing wave, or 'oscillation'. Identifying reflections in data could determine the depth of the tangent cylinder and hence the equatorial radius of the metallic region. Also, this may distinguish the Jovian torsional waves from those in Earth's core, where observational evidence has suggested waves mainly travelling outwards from the rotation axis. These waves can transport angular momentum and possibly give rise to variations in Jupiter's rotation period of magnitude no greater than tens milliseconds. In addition those internal disturbances may be revealed above the metallic region as changes in the zonal flow, with an amplitude of about one tenth of that of the steady zonal jets. The Juno's gravity measurements have indicated that the zonal flows visible at Jupiter's surface penetrate down to a depth of a few thousands of kilometres, so torsional waves in the interior could imprint themselves as surface flow variations.


Introduction
Torsional Alfvén waves (TWs) are a special class of magnetohydrodynamic (MHD) waves, whose transverse motions are confined to cylindrical surfaces aligned with the rotation axis. They are perturbations about the Taylor state [37] expected at leading order when the Coriolis, Lorentz, and pressure gradient forces are in balance in the momentum equation, the so-called magnetostrophic balance. The linear theory for incompressible Boussinesq fluids was introduced by Braginsky [5] and has been applied to Earth's core, in which fluid motions of liquid iron are believed to generate the global, intrinsic magnetic field. The axisymmetric disturbances can propagate in cylindrical radius (denoted by s hereafter), perpendicular to the rotation axis. This enables the waves to transport the angular momentum to other regions, including the rocky mantle and solid inner core, through electromagnetic and/or other couplings (see [30] for an overview). The evidence for such waves within the fluid core has been discussed using core flow models inverted from the observed geomagnetic secular variation (SV): the zonal component was found to exhibit fluctuations with a near six-year period [15,16]. They may also account for a decadal variation of the length-of-day (LOD) of Earth [20]. Such information provides insight on the deep interior by constraining physical quantities that cannot be measured directly, such as the field strength within the core and the electrical conductance of the lowermost mantle.
We here extend the study of TWs to compressible fluids by applying the anelastic approximation, where sound waves are omitted. This is of some interest for geophysical modelling, since there is a density increase of 22 % from the bottom to the top of Earth's fluid outer core [e.g . 22]. Extension to anelasticity is, however, more strongly motivated by a desire to explore the internal dynamics of gas planets and stars, which mostly consist of hydrogen and helium. Jupiter is the largest planet in the solar system and also has the strongest planetary magnetic field, with a surface magnitude of O(10 G), or O(1 mT). Dynamo action is predicted to operate in a metallic hydrogen region situated below a molecular hydrogen envelope. The phase transition is expected to occur continuously between 0.85 and 0.90 R J , with R J being Jupiter's mean radius at the 1 bar level. Adopting an interior model including the transition [13], dynamo simulations for anelastic convection have reproduced Jupiter-like magnetic fields [23,14]. The gas giant is rapidly rotating with a period of 9.925 hours (the System III); changes on the order of tens of milliseconds have been noted [19,32]. The rapid rotation and strong magnetic field in the metallic region yield a force balance with the inertial and viscous forces small compared to the Coriolis and Lorentz forces, the quasi-magnetostrophic balance. Jupiter is therefore a good candidate for hosting TWs.
MHD waves excited within the gas giant may produce decadal variations, as shown below. In-situ observations of Jupiter have the longest history amongst all planets other than Earth, spanning over forty years since the Pioneer epoch in the early 1970s. Coverage is however sparse; although data retrieved from past missions has enabled the construction of global models for the magnetic field, such data was limited to spherical harmonics of degree no higher than seven [8,29]. Ridley & Holme [29] showed time-dependent field models to be preferable to steady models, and attempted to invert the SV to flows at the top of the expected metallic region. The Juno spacecraft is currently orbiting the gas giant and the newly available data displays a magnetic field sampled at the closest point ever to any planetary dynamo region [3,4,24,9]. Over the planned five-year mission it will better define the field in both temporal and spatial resolution. Also, the gravitational sounding has indicated that zonal flows extend thousands of kilometres below the Jovian surface [26,17]. The cloud motion has been tracked for decades by Earth-based telescopes to measure changes to, or periodicities in, the zonal winds [41]. The coloration, brightening, and outbreak events, sometimes leading to global upheavals, have been monitored for more than one hundred years [31,12].
Of more theoretical interest is the nature of excited TWs. Since the Alfvén waves are able to propagate in s inwardly and outwardly, early studies proposed TWs in the form of standing waves and sought wave motions in the form of normal mode solutions [5,46,7], often referred to as torsional oscillations. However, Earth's core flow inversions/assimilation [16] and numerical geodynamo simulations [44,38,35] have illustrated TWs travelling predominantly in an outwards direction with no obvious reflection at the boundaries. This could be explained through preferred excitation [38,39,40] near the tangent cylinder (TC, the imaginary cylinder aligned with the rotation axis that circumscribes the inner core) and dissipation beneath and above the core-mantle boundary (CMB) [33,34]. Studies ignoring dissipation showed that the effect of spherical geometry and variable internal magnetic fields can give rise to asymmetric reflections and hence weaken reflected waves [10]. We shall demonstrate that TWs in the gas giant's metallic region can be reflected from a magnetic TC, which is formed due to the transition to molecular hydrogen, and exhibit its standing nature.

Theory
The theoretical framework of incompressible TWs is well documented [e.g. 5,38,22]. In the light of those studies, we consider anelastic fluids where the Lantz-Braginsky-Roberts formulation [6,27,25] is adopted and explore anelastic TWs within the electrically conducting region of the gas planet. We assume the equilibrium state is close to adiabatic, well-mixed, and hydrostatic with density ρ eq . The velocity perturbations of the waves u are subsonic, so that the continuity equation is ∇ · ρ eq u = 0 .
We assume a basic state dependent only on spherical radius, r, and denote it by subscript 'eq' hereafter. We focus on the rapid dynamics, in which the characteristic timescale is much shorter than the diffusion time. This allows us to begin with the momentum equation where Ω is the rotational angular velocity, j is the current density, B is the magnetic field, S is the entropy, T eq is the temperature in the equilibrium state, e r is the unit vector in the radial direction, andp is a reduced pressure incorporating the density and the gravitational potential. To look at fluctuations corresponding to TWs, we consider the axisymmetric z-independent azimuthal flow by taking averages of the φ-component of the momentum equation over cylindrical surfaces to give where the azimuthal and axial averages are defined as f = (1/2πs) 2π 0 f sdφ and f = (1/h) z+ z− f dz, respectively, with h = z + − z − for any scalar field, f . Outside the TC, the integral is limited by z ± = ± r 2 o − s 2 with r o being the radius of the planet. Hereafter we shall focus on the region outside the TC. With the divergence theorem, the Coriolis force becomes F C = −(Ω/πsh) ∇·ρ eq u dV for geostrophic cylinders. From (1), this term vanishes, implying zero net mass flux across the surfaces. For the magnetostrophic balance where the inertia and F R are negligible, (3) gives F L = 0, i.e. the Taylor state for anelastic fluids. The Lorentz and Reynolds forces may be rewritten as respectively. The second term of each represents the surface term coupling the internal fluid region and the outside, magnetically or dynamically. Since the currents vanish outside the metallic hydrogen zone, the F L surface term will be small, and the average over the cylinder could be taken only over the conducting region. For the stress-free outer boundary used in Jupiter simulations, the F R surface term vanishes also. However, unlike the magnetic term, the molecular non-conducting region can contribute significantly to the F R integral, because the convection-driven velocities are large there, as we shall see in section 4.2.
We now make the ansatz of splitting magnetic field and velocity into their mean and fluctuating parts: The time interval τ is chosen to be significantly longer than the expected wave-period, but not excessively longer to avoid unnecessary computational expense. Here the induction equation for compressible fluids is Recall that we primarily seek the rapid dynamics within the conducting fluid region so we ignore any dissipation; the magnetic diffusion will become substantial as the wave goes up to the poorly-conducting zone, and will damp waves through ohmic dissipation, but the frequencies and the wave-forms in the conducting region should be relatively unaffected by diffusion. We now substitute (6) into the time-derivative of the momentum equation (3). There is a question of whether the TW equation should be expressed in terms of ρ eq u φ /s or u φ /s, because the time-derivative of ρ eq u φ appears in (3), but (6) contains spatial derivatives of u φ , not ρ eq u φ . Here we choose u φ as the dependent variable: we separate u ′ φ into its geostrophic and ageostrophic parts, and ignore the second ageostrophic term, because it is small compared to the first term in our simulations. We then obtain The theory is equivalent to that of the incompressible case [38] but now the effect of compressibility remains in the Lorentz term. The fluctuating components are assumed to be significantly smaller than the mean parts. Moreover, when quadratic quantities comprised of a mean part (i.e. | B 2 s | and | B 2 φ |) are greater than those of cross terms (e.g. | B s B φ |), the equation can be reduced to leave the homogeneous part where U 2 A = B 2 s /µ 0 ρ eq , implying a wave equation for angular velocity in the anelastic case [22]. We note that another possible definition would be B 2 s /µ 0 ρ eq [see Appendix A], but this is less convenient in our formulation. Here the restoring force of the wave is represented by F LR , while the remaining terms of the Lorentz force can be summed up to a term F LD = F ′ L − F LR . This term represents the convection-driven fluctuations, which interact with the magnetic field to drive the TWs through the Lorentz force and to modify their waveforms and/or speeds. Below we see the latter effects but it still minor in our simulation, so we call F LD a forcing term. The waves can also be driven by convective perturbations in the Reynolds force denoted by F ′ R . A perturbation of angular velocity, u ′ φ /s, can propagate in s with Alfvén speed, U A . The speed depends on the magnitude of the background field, B 2 s , and the density, ρ eq , both of which vary with s. This special mode is nondispersive, i.e. the speed is independent of wavenumbers. Since the equation allows both inward and outward propagation, a superposition of those modes could yield standing waves and enable normal mode solutions. However, observational data for Earth and numerical simulations indicate a preference for (outwardly) propagating waves over standing ones (sec. 1).

Model description
To explore excitation of TWs in the gas giant, we adopt Jovian dynamo models, which were built by Jones [23] (hereafter referred to as J14) and Dietrich & Jones [11]. Only the essential part of the analysis is shown below: see J14 for the detailed description of the model set-up. The models explore the selfgeneration of magnetic fields by anelastic fluid motions in rotating spherical shells. The equilibrium reference state calculated by [13] was used, and viscous and diffusion terms are taken into account. The reference state density, ρ eq , electrical conductivity, σ eq , and temperature, T eq , arise from a composition comprising of a metallic hydrogen region above a rocky core, r ≥ r c ∼ 6.45 × 10 6 m ∼ 0.09 R J , and its continuous transition to a molecular hydrogen region. The transition begins at r ∼ 0.85-0.90 R J and only the region below a cut-off level, r ≤ r cut ∼ 6.70 × 10 7 m ∼ 0.96 R J , is modelled in our simulations. The density scale height, N ρ = ln [ρ eq (r c )/ρ eq (r cut )], between the core boundary and the cut-off radius is approximately 3.08. Convection is largely driven by a uniform entropy source, which is released as the planet cools; this differs from the present geodynamo, which is primarily driven by buoyancy sources arising from the inner core boundary due to its freezing.
As the electrical conductivity, σ eq , drops by more than five orders from the metallic to the molecular region, a poorly-conducting layer is formed at the outermost part of the shell. Despite compressibility, the Proudman-Taylor constraint still strongly influences fluid motions in the outer layers when electrical conductivity is negligible. This produces a second imaginary cylinder, aligned with the rotation axis, that circumscribes the metallic hydrogen region which we call the magnetic tangent cylinder (MTC) [11], located at s ∼ 0.85-0.90 R J . This is in addition to the traditional 'kinematic' TC found at s = r c = 0.0963 r cut ≡ s tc , circumscribing the solid core. Unlike the kinematic TC, the MTC is not precisely defined, as the conductivity drop occurs over a finite radius range, but this range is thin enough for the MTC concept to be useful here: we hereafter denote the minimal s (∼ 0.885 r cut ) as s mtc . The Jovian core leaves only a small fraction of the domain inside the TC. We shall concentrate on the region outside the TC but inside the MTC, i.e. s tc ≤ s ≤ s mtc .
We select three (models A, E, and I) out of nine models examined by J14, which differ only in model parameters and entropy outer boundary conditions. The chosen models and key quantities are listed in Table 1. The global Rossby number, Ro, quantifying the relative strength of the inertia to the Coriolis force, is shown to be no greater than 5 × 10 −3 . The advective terms therefore only have a minor role, as assumed in the linear theory. The Elsasser number, Λ, measures the ratio of the Lorentz force to the Coriolis force and is found to be approximately 6-10 in our simulations. Model I was reported to reproduce a magnetic field morphology broadly resembling that measured by Juno [24]. Some smaller scale features of Jupiter's magnetic field as revealed by the mission more recently [9] differ from the models, notably in the equatorial asymmetry of the small scale field. However, wave propagation is determined mainly by the large scale magnetic field, and TWs involve averaging over cylinders passing through both hemispheres. So this refinement of the Jovian magnetic field is not likely to affect our results greatly.
The magnetic fields self-generated in our simulations are non-reversing and dipolar during the simulations. They act as the background field for the MHD wave motions discussed below. The propagation speed of TWs is determined by the cylindrically averaged B s field. In figure 1, a solid curve depicts the nondimensional Alfvén speed, U A , as a function of cylindrical radius, s, normalized by the cut-off radius, r cut , for model I. Here the time and length are scaled by the magnetic diffusion time and the shell thickness (D = r cut − r c ), respectively, and the bounds for z-averages are taken at r cut . In the figure, the dashed line represents the Alfvén speed with the density taken to be its constant mid-radius value at r m = (r c + r cut )/2 in the definition of U A . As the density ρ eq decreases with radius r, the z-mean Alfvén speed increases with s (in accordance with the definition of U A ). This gives rise to an increase in U A in the anelastic case with a peak at s/r cut ∼ 0.6. For larger s, the speed gradually weakens as the MTC is approached and crossed. Plots in other cases look similar; models A and E demonstrate a well-defined peak at 0.6 s/r cut 0.7. Table 1 also lists the speeds U A (s mtc ) at the MTC radius and the expected traveltimes τ A from the core boundary s tc to the s mtc . The speeds U A are used for conversion to our dimensional time unit below (sec. 4.1): a Jovian scale U J A is demonstrated on the right-hand side of the axis in fig. 1.

Internal dynamics: zonal flow fluctuations and their excitation
The time averaged components of azimuthal velocity, u φ , show one very strong prograde jet outside the MTC and rather incoherent mean alternating flows within it (Figure 6 of J14). In spite of the presence of generated magnetic fields and anelasticity, axisymmetric zonal flows inside the MTC still retain a significant fraction of the z-independent part of the flow. By removing the mean part, we identify fluctuations of azimuthal flows, u ′ φ , which are of interest here. Figure 2 displays contours of u ′ φ in s-t space for the three runs. In each diagram, white curves indicate the calculated Alfvén speed, U A , to compare with the computed fluctuations. A dimensional time is shown on the top of each image (details in the following subsection). Run A ( fig. 2a) shows that some disturbances emerge near s/r cut 0.6 and move outward to the poorly conducting layer; they can also be found to travel inwards towards the core boundary. Their propagation speeds fit well with the predicted U A , suggesting that they are anelastic torsional Alfvén waves. They become more evident when filtered (see Appendix B). Travelling TWs are found in Earth-like Boussinesq models [44,38,39,35]; they mostly originate in the vicinity of the TC, where vigorous convection occurs near the solid inner core. No obvious standing TWs have been found in geodynamo simulations to date. Figure 2b displays contours of u ′ φ for model E, where the relative strength of the viscosity to the Coriolis force is decreased. We see significant fluctuations repeatedly occurring at an outer radius, s/r cut 0.6. Interestingly, beneath the MTC, waves appear to form a node at around s/r cut ∼ 0.65 and ∼ 0.8, so figure 2b shows evidence of standing waves being excited in the Jovian models. We also find propagating features at a later time, t 0.003. There are signatures of reflection, highlighted by the white lines, around the MTC at, for instance, t ∼ 0.0038.
A simple one-dimensional model of Alfvén waves propagating into a region where the diffusivity increases over a transition region was considered (not shown; see Appendix C for the uniform diffusivity case). It shows that incident waves whose wavelength is shorter than, or comparable to, the thickness of the transition region are absorbed by diffusion, whereas waves with a wavelength longer than the transition thickness are mostly reflected. The theory also shows there is no phase change in u ′ φ , so a red patch in figure 2b should reflect into a red patch, as seen in the figure. When a wave packet reaches the MTC, the shorter wavelength components comparable to the (rather thin) transition region thickness are absorbed, while the longer wavelength components are reflected. This contrasts with the circumstances around Earth's CMB, which is a hard boundary of the fluid; there a combination of the viscous dissipation and magnetic dissipation across the CMB controls the behaviour [34].
In model I, where the entropy flux at r cut is a given constant, the nature of reflections from the MTC has been studied. Figure 2c shows the interaction of the blue feature with the MTC at 0.0008 t 0.0010. Note that poor resolutions and/or improper filters may make the reflecting nature of the waves less clear (see Appendix B). To describe the time evolution, we also present profiles of u ′ φ in figure 3. A trough is initiated at t ∼ 0.0006 and s/r cut ∼ 0.75. As time evolves, it eventually grows, while the waveform becomes sharper (fig. 3a). This suggests a nonlinear influence on the TWs, arising from the terms F LD and/or F R [e.g. 43]. At t ∼ 0.00095, the trough reflects at around s mtc but also passes through the transition zone. The patterns of the incident and reflected waves are compared in Fig. 3b, which illustrates a positive reflection. However, there is a superposition of continuously excited waves, so the amplitude and shape vary in our nonlinear simulation, and it is hard to determine the reflection coefficient or the phase shift accurately. An abrupt change in the U A profile may also yield reflections [e.g. 1]. Forward simulations of linear, nondissipative TWs in spherical geometry reported internal reflections where the gradient of a background magnetic field was steep [10]. The role of the background velocity on the reflections in our simulations has not yet been elucidated.
The excitation mechanism of the waves is investigated in figures 4a and b which display the forcing terms F ′ R and F LD , respectively, for run I. The Reynolds force is found to be important in the outer regions, 0.6 r cut s s mtc , whereas the Lorentz force is more evenly spread throughout the region and so more dominant in the interior. This is understandable as the cylinders defining the wave motion which have larger s have a greater proportion of area in the vigorously convecting outer layers. The term F ′ R better matches the location and time at which u ′ φ disturbances begin to travel than F LD . This is in spite of the small Rossby number; such an initiation was pointed out in Boussinesq cases [38]. The convective motions are most vigorous in the outer layers of our Jupiter models, similar to Fig. 6d of J14; the density stratification enhances convective velocities to get the heat flux out. This produces a convergence of the Reynolds stress, particularly through the ρu s u φ term, and continuously forces fluctuations, ∂u φ /∂t, by almost-hydrodynamic Rossby waves (not shown), which are much faster than the Alfvén waves. TWs in model A are also predominantly excited by F ′ R ; it is rather mixed with Lorentz terms F LD in model E. Models for Earth, by comparison, have shown that the driving of TWs is possible by either the Reynolds or Lorentz force. Geodynamo simulations have often shown the Reynolds force to be the largest contributor [38]. However, as parameters are moved towards their Earth-like values, geodynamo and magnetoconvection simulations [39] display a growing influence of the Lorentz force. This is to be expected as the role of the magnetic field increases as a balance closer to magnetostrophy is achieved. The models for the Jovian dynamo discussed in this work use moderate values of the Ekman number, E, representing the ratio of the viscous force to the Coriolis force, and the resulting Elsasser number is smaller than was possible in [39]. It is not yet clear whether TWs in Jupiter will be primarily driven by Reynolds or Lorentz force, and possibly both will be significant. In Earth, Lorentz force will dominate, but in Jupiter convective velocities increase with radius. It is possible strong convection in the outer regions could provide a significant contribution to the driving over the whole MTC, even though the density in these upper regions is small. Further simulations at lower E and greater field strengths are necessary to decide this driving mechanism issue.

Rescaling to the dimensional unit
To examine whether signals due to TWs may be detectable in observational data, we first convert the nondimensional time in our simulations to a dimensional unit. This issue arises because the current numerical models are still limited to accessible parameters [35,2], and cannot solve the equations with variables spanning many orders of magnitude, which is the case within the gas giant. Amongst several choices we here choose the magnetic field scale by equating the field strength at the magnetic outer boundary in simulations to the observed outer boundary value [38]. Since the density is quite well-known, this gives the Alfvén speed and hence a conversion between dimensionless and dimensional time.
Jovimagnetic models show the magnitude of the radial component to be no greater than 60 G, or 6 mT, on a surface of r ∼ 0.85 R J : in the equatorial region it is seemingly no greater than 30 G and weaker than 1 G for large regions [9]. Taking 30 G as a reasonable maximal field magnitude at our MTC radius, s mtc , and the density, ρ eq , of 8.53 × 10 2 kg m −3 , then an Alfvén speed, U J A , at this radius is approximately 9.16 × 10 −2 m s −1 . By matching this value with those of our simulations, our dimensional time unit τ unit is calculated through DU A /U J A , where D is the shell thickness of 6.06×10 7 m. From this we calculate dimensional versions, τ J and τ J A , of the analyzed interval, τ , and the TW traveltime, τ A , respectively. Values for each run are listed in Table 2. While time units vary from 6.1 to 8.8 thousand years (as do analyzed time windows τ J from 31 to 44 years), the traveltimes τ J A all fall within a 9-13 year window. A difficulty for our scaling arises when converting the averaged azimuthal velocities into dimensional units. The typical convective velocity at 0.85 R J is believed to be around 10 −2 m s −1 [24], but the surface equatorial zonal flow is nearly 100 m s −1 . Simulations do get zonal flows that are larger than the convective flow, but they cannot yet reach the 10 4 ratio due to the enhanced viscosity in the models, so it is uncertain how the axisymmetric azimuthal flow should be scaled. Taking the unit of velocity as D/τ unit gives 3.15 × 10 −4 m s −1 for run A, and this is the unit used for the averaged azimuthal flow in the figures. This gives a rather large convective velocity estimate of 0.5 m s −1 but a reasonable estimate of the zonal flow at 0.85 R J of about 2 m s −1 . If we use the dimensionless time unit which puts the convective flow at 10 −2 m s −1 , the amplitude of the azimuthal flow is reduced by a factor of around 50.

Length-of-day variation
Fluctuations in axisymmetric zonal flows produce variations in the angular momentum of the metallic hydrogen region which can be transferred to other parts of the planet. This may produce fluctuations of the rotation period of the gas giant, namely LOD: this is often defined with the magnetic field (System III) that is generated in the metallic region. In Earth, in contrast, the LOD is fixed to the reference frame of the mantle. Earth's LOD variation with a period of nearly six years with amplitude O(10 −4 s) has been identified; its origin could be angular momentum exchange between the fluid core and the rocky mantle through MHD waves ( [15]; sec. 1). One may envisage an analogous coupling in Jupiter between the deeper conducting metallic region and the overlying transition-molecular envelopes, as well as a Jovian LOD fluctuation.
We evaluate the influence by calculating the axial angular momentum change that is deduced from the axisymmetric disturbances in our metallic hydrogen region, and those outside the region, In figure 5 the solid and dotted curves display the time evolutions of δσ and δσ omtc in model E, respectively. The δσ of the conducting region shows a quasiperiodic variation, corresponding to the flow oscillations ( fig. 2b). The evolution is almost perfectly anti-correlated with the change δσ omtc of the outermost transition zone, as it should be since total angular momentum is conserved. Of interest is the coupling mechanism across the MTC: our simulations indicate this is provided by Reynolds stress and Maxwell stress. Using our standard time unit and the density (ρ eq (r m ) = 2.56 × 10 3 kg m 3 ), we convert the dimensionless δσ of maximum amplitude 38.7 to its dimensional version, δσ J with maximum amplitude ∼ 1.58 × 10 32 N m s. Assuming a value of I = 2.56 × 10 42 kg m 2 for the moment of the inertia [13] and a daily period of P = 3.57 × 10 4 s for the planet, the change δσ J is equivalent to a period δP of approximately 12.5 ms. Here δσ J = IδΩ = −2πIδP /P 2 , where Ω is the angular velocity. If we use the longer dimensionless time unit which puts the convective flow at 10 −2 m s −1 , the amplitude of the period change is reduced by a factor of around 50 (sec. 4.1). In table 2 we give this alternative scaling, which for model E gives about 0.25 ms, in brackets.
Jupiter's mean LOD is determined to a precision more accurate than seconds using the System III. These estimates largely rely on measurements of the decametric radio emission from the magnetosphere since the 1950s. Decadal averages of the observed radio rotation period shows its changes on the order of tens of milliseconds; this remains the subject of some debate [19,32,29]. The measurements may reflect a time-varying SV due to unsteady convective flow, rather than changes being solely due to TWs. Our estimates indicate that TWs could be a part of LOD fluctuations, but separating the convective flow-induced changes from the TW changes will not be easy.

Flow change above the metallic region
Unlike terrestrial planets, the gas giant's composition may allow for a direct inference of deep-origin perturbations via observation of surface flows above the metallic region. Figure 6a and b portray contours of the fluctuating zonal flow u ′ φ on the cut-off surface, r cut , in the northern and southern hemisphere for model E, respectively. The latitude-dependent data is displayed in s-t space to enable comparison with the figures and wave speeds shown earlier. Note that the plots include the region outside the MTC. The amplitude is scaled by the maximum of the mean zonal flow u φ on the same surface, which is indeed the maximal speed of the prograde, equatorial jet reproduced in the simulation (sec. 3.2). Figs. 6c-d demonstrate the same data filtered to remove all periods outside the range from t = 0.00063 to 0.0025, i.e. from 5.6 years to 22 years in the dimensional unit for 30 G.
In both hemispheres we find corresponding fluctuations on the surface, although they look much noisier than the internal wave motions. The filter used in figures c-d helps to visualize the wave signals more clearly. The variations are found to be almost symmetric with respect to the equator, which is a consequence of the predominantly z-independent flow. Standing oscillations and both equatorward and poleward propagation are seen at mid and high latitudes where s/r cut s mtc , whereas the equatorial region s/r cut s mtc features only equatorward migration. We interpret this as partial transmission through the MTC and absorption within the resistive, transition layer (sec. 3.2; Appendix C). Such an abrupt change in zonal flow fluctuations on spherical surfaces signifies the location of the MTC. Thus it can act as an indicator of the location where the transition begins, i.e. the magnetic dissipation becomes significant. This is however hard to identify when searching in a few snapshots only, as examined for zonal wind profiles; exploration in θ-t space -sometimes called a Hovmöller diagram -is essential for the identification.
On the surface the maximum of the fluctuating velocity is 11 % of the mean velocity. Other models show analogous fractions of 12-15 %; the values are listed in table 1. Converting to dimensional units as before, the fluctuation amplitudes at r cut ∼ 0.96 R J are found to be about 0.1-0.2 m s −1 and 3-5 ×10 −3 m s −1 for the equatorial field of 30 G and 0.6 G, respectively, whereas the mean velocities are an order greater (table 2).
Deep seated zonal flows, and possibly small-scale convective plumes, may penetrate into the stably-stratified weather layer unless the stratification is too strong [47,36,18]. The gravitational harmonics recently obtained by the Juno mission seemingly support this. Provided flows are governed by the thermal wind balance and some radial structures, the cloud-top zonal winds would stretch to ∼ 0.95-0.97 R J , possibly reaching part of the metallic hydrogen zone [26,17]. The steady equatorially-symmetric part would slow to less than 6 m s −1 at about 0.96 R J [17]. This is compatible with our conversion range of 1.5-0.03 m s −1 (see u φ in table 2). At lower depths flow models inverted from the jovimagnetic SV suggest a velocity of order 10 −2 m s −1 at the top of the expected conducting region (given 0.85 R J ), not far from those estimated with scaling properties based on available heat fluxes [29,45]. Those deep dynamics may arise to set the thermodynamical conditions at the cloud deck and to trigger visible photochemical changes.
Earth-based campaigns have monitored the long term variability of the cloud and/or atmospheric appearance of the gas planet. Global upheavals are known to comprise of recurrent activities spreading over several latitudinal bands and to occur at intervals of decades, irregularly in most cases [31,12]. At some epochs, 5-or 10-years periodicities in jetstream outbreaks or fades/revivals were recognized at the North Temperate Belts (NTB), 23-35 • N. Recent datasetsprimarily collected with the Hubble Space Telescope between 2009-2016 -have been updated to identify the most relevant change in zonal winds near 24 • N of about 10 m s −1 and 5-7 year periods at a few lower-latitudes [41]. These seem interesting in comparison to the internal flows simulated earlier. First, the amplitude of such disturbances is ∼ 10 %, or less, of the 150 m s −1 stable jet [41]. Second, the NTB latitudes correspond to cylindrical radii s ∼ 0.82-0.92 R J , which likely lie in the outermost part of the metallic region and the transition zone. At these radii/latitudes, TWs in our models exhibited oscillations with periods of several years or longer and sometimes at irregular intervals.

Concluding remarks and discussions
We have demonstrated, through our anelastic models, that torsional Alfvén waves could be excited in Jupiter's metallic hydrogen region, which is expected to exist below a molecular hydrogen layer. The axisymmetric MHD disturbances can propagate in cylindrical radius on timescales of Alfvén speeds in a medium with a variable equilibrium density. In the Jovian dynamo models we adopted, waves were excited at an outer radius of the conducting region, at which nonaxisymmetric convective motions were vigorous and the resulting stresses drove the axisymmetric fluctuations. TWs were found to travel both towards the molecular region and towards the core boundary. Modes propagating outwards were found to be partially transmitted to the poorly conducting layer but could also be reflected around the MTC. Consequently, this results in waves travelling inwards, back into the deeper interior of the conducting region. Since convection perturbs the fluid at all times with timescales of Rossby waves, it is able to continuously supply a source for TWs travelling in both directions. Provided that their amplitude and timing match, a superposition of the opposed propagation enables the formation of standing waves, as observed in our model E. It is however not straightforward; our simulations illustrated a mixture of travelling, reflecting, and standing waves. Indeed this could be the case in reality; this differs from conventional arguments where propagating waves and normal mode oscillations are treated as being mutually exclusive. The variability in the form and location of the waves here is in contrast to recent work displaying TWs in the geodynamo [44,39,35]. In such models the waves appear to be preferably excited at a location with vigorous convection near the inner core and they do not reflect upon impact with the CMB.
A key requirement for reflection here is the existence of the MTC, which is created by the drastic decrease of the electrical conductivity in the gas giant. The MTC may act as an interface for the waves approaching the magneticallydissipative fluid layer, which enables reflection as well as transmission. This significantly differs from the corresponding region in Earth's core, in the vicinity of the CMB, which is typically modelled as a hard bound for the fluid motion. The interface created by the varying conductivity may allow reflection of waves to be a feature within Jupiter. Whilst the size of the dynamo region is currently hard to define, detecting reflections from data may enable us to infer the radius where the transition from metallic to molecular hydrogen indeed begins. This is analogous to how seismology has constrained the structure of the deep Earth.
TW traveltimes across the metallic region were estimated at several years, provided that the time units in our simulations were chosen so that the Alfvén speed at the equator at our MTC level matches that suggested by jovimagnetic models and adiabat density models. An equatorial radial field of maximal strength 30 G yields traveltimes of 9-13 years; longer timescales are feasible when a weaker field is implemented.
The fluctuations of zonal flows yielded the exchange of angular momentum between the metallic hydrogen region and the overlying transition-molecular regions. With the same time units being adopted, the magnitudes could be translated to variations in the LOD no greater than 10 −2 s. Alterations in Jupiter's radio rotation period [19] might partly embrace the true LOD changes besides the magnetic SV.
Our simulations also demonstrated the wave motions identified through zonal flows on a spherical surface above the metallic region. The surface fluctuations were sizable, up to 15 % of the maximal amplitude of the steady zonal component. The reflecting and/or transmitting nature across the MTC could be projected upwards from the metallic region. Juno's gravity measurements have constrained interior models by identifying that visible surface zonal flows penetrate downwards significantly [26]. Assuming this deep origin, variations at the cloud deck could display partial evidence of TWs.
Another source for detection is the jovimagnetic SV which is inferred at the top of the metallic region. The projection from internal wave motions to the SV is rather complicated, as discussed in Earth's fluid core. TWs may contribute to the occurrence of geomagnetic jerks; they cannot however account for all phenomena alone (see [28] for a review). To link wave motions to the magnetic SV, zonal flow fluctuations seen at the surface are required to linearly induce rapid variations in the surface magnetic field. Here caution is necessary since such magnetic changes can be initiated by other terms in the induction equation: effects of the other flow component [42] and of nonlinear couplings with the field may remain important [21]. Nevertheless, an increase of spatial and temporal coverage in magnetic data is expected to better resolve the SV and inverted flow models on the top of the metallic region. The ongoing Juno magnetic measurements, coupled with theoretical studies, will offer a promising route to develop our knowledge on the dynamics in the dynamo region. Table 1: Dimensionless input and output parameters of dynamo simulations selected from J14 [23]. The columns E, Ra, and H list the Ekman number, the Rayleigh number, and the volumetric entropy source, respectively. The kinetic and magnetic Prandtl numbers are fixed at P r = 0.1 and P m = 3, respectively. Constant entropy on both boundaries, rc and r cut , for models A and E; for model I, constant entropy on the inner boundary and fixed entropy-flux outer boundary. In all three models, a stress-free outer boundary, a no-slip inner boundary, and electrically insulating conditions both at both boundaries are used. The columns τ , Ro, and Λ list the analyzed time interval, the Rossby number (equal to urmsE/P m in our scaling with urms being the rms flow vigor), and the Elsasser number (equal to the rms field strength, Brms, in our scaling), respectively, where averages over the whole volume are taken. These are used to yield the Lehnert number Le = Brms/( √ ρeqµ 0 DΩ) = ΛE/P m of 6.6-7.3×10 −3 and the Alfvén number A = urms √ ρeqµ 0 /Brms = Ro/Le of 0.45-0.62. The columns U A (s mtc ), τ A , δσ, u ′ φ (r cut ), and u ′ φ / u φ (r cut ) represent, respectively, the Alfvén speed at the MTC, the traveltime across the conducting region (from the core boundary to the MTC), the maximal amplitude of the axial angular momentum fluctuation in the metallic region, the maximal fluctuating zonal velocity, and the fraction of the maximal fluctuation part to the maximal time-averaged part at the cut-off boundary.  Table 2: Dimensional parameters of the simulations. The first column τ unit lists the time unit which is represented in years and is used to convert our dimensionless time to its dimensional version. The assumption is made of an Alfvén speed of 9.16 × 10 −2 m s −1 (1.83 × 10 −3 m s −1 ), or, equivalently, a radial field strength of 30 G (0.6 G), at the equator at the top of the metallic region ∼ 0.85 R J ; calculated as described in the main text. The columns τ J , τ J A , δσ J , δP , u ′ φ (0.96 R J ), and u φ (0.96 R J ) present the dimensional version of the analyzed interval τ , traveltime τ A , maximal angular momentum δσ, maximal LOD variation, maximal zonal flow fluctuation u ′ φ at the surface r cut ∼ 0.96 R J , and maximal mean zonal flow u φ at the same surface, respectively.         (d) Figure 6: (a-b) Zonal flow fluctuations, u ′ φ , on the surface of the cut-off boundary, r cut ∼ 0.96 R J , for run E. In northern (a) and southern (b) hemispheres. The ordinates represent the cylindrical radius s/r cut so that white curves indicate the same phase paths as those shown in fig. 2b. The horizontal dashed lines denote the range of the MTC radii, s/R J ∼ 0.85 and 0.90, corresponding to latitudes ∼ 32 • and 26 • on the surface, respectively. The amplitude is scaled by the maximum of the steady part, u φ , corresponding to about 1.4 m s −1 in the dimensional unit; see table 2. (c-d) Same as figures a-b but with periods outside 0.00063≤ t ≤0.0025, or 5.6 yrs ≤ t J ≤ 22 yrs, filtered from the data.  1). The columnsÛ A (s mtc ),τ A , and δσ represent, respectively, the Alfvén speed at the MTC, the traveltime across the conducting region (from the core boundary to the MTC), and the maximal amplitude of the axial angular momentum fluctuation in the metallic region.  fig. 2b]. White curves indicate phase paths of the Alfvén speed,Û A = B 2 s /µ 0 ρeq , in the present formulation. The horizontal dashed lines indicate a range of the MTC radius, s/r cut ∼ 0.89 and 0.94. Note that the amplitude ρequ ′ φ is small outside the MTC than in the metallic interior, whereas u ′ φ is substantial outside the MTC (see fig.2b). Appendix C. Alfvén waves approaching a resistive zone We consider a Cartesian, one-dimensional model for Alfvén waves approaching a resistive layer. Let x = 0 the interface between a perfectly conducting fluid (for negative x) and a weakly conducting one (for positive x). They are permeated by a uniform background magnetic field B 0 in the x direction. For simplicity we assume an incompressible fluid with ρ 0 being constant density. We then suppose the variables B = B 0êx + b y (x)ê y , u = u y (x)ê y (C.1) to rewrite the equations of induction and momentum as ∂b y ∂t = B 0 ∂u y ∂x + ∂ ∂x η(x) ∂b y ∂x and ∂u y ∂t respectively. Here the magnetic diffusivity, η(x) = 1/µ 0 σ(x), varies in x: it is set zero for x < 0 and to a constant nonzero value η 0 for x > 0. At the interface the field and velocity are continuous, i.e. the continuity condition across x = 0 is required for b y and ∂b y /∂x. Eq. (C.2) may be reduced to where the Alfvén speed V A = B 0 / √ ρ 0 µ 0 . Now we seek solutions of the form b y = e iωt e −ikx + Re ikx for x < 0 (C.5) b y = T e iωt e λx for x > 0 (C. 6) where λ, R and T are complex and k 2 = ω 2 /V 2 A . For x > 0, substituting (C.6) into the respective wave equation (C.4) gives (C.7) When the waves travel quickly so that ω ≫ V 2 A /η 0 , the valid solution is So the continuity condition on b y and ∂b y /∂x implies, respectively, We hence obtain the reflection coefficient, For ω ≫ k 2 η 0 , this yields R → −1 and T → 0, i.e. nearly perfect reflection. From (C.2), this is equivalent to positive reflection of u y across the interface. When the approximations are inappropriate, it gives rise to partial reflection and partial transmission.