Anelastic torsional oscillations in Jupiter's metallic hydrogen region

TL;DR

This paper models anelastic torsional Alfvén waves in Jupiter's metallic hydrogen region, exploring their propagation, reflection, and potential observational signatures to understand Jupiter's internal dynamics and magnetic field.

Contribution

It introduces a model of anelastic torsional waves in Jupiter's metallic hydrogen region, including their propagation, reflection, and observational implications, extending prior incompressible models.

Findings

Waves travel perpendicular to the rotation axis over several years.

Reflections at the magnetic tangent cylinder can form standing waves.

Waves may cause observable variations in Jupiter's rotation and zonal flows.

Abstract

We consider torsional Alfv\'en waves which may be excited in Jupiter's metallic hydrogen region. These axisymmetric zonal flow fluctuations have previously been examined for incompressible fluids in the context of Earth's liquid iron core. Theoretical models of the deep-seated Jovian dynamo, implementing radial changes of density and electrical conductivity in the equilibrium model, have reproduced its strong, dipolar magnetic field. Analysing such models, we find anelastic torsional waves travelling perpendicular to the rotation axis in the metallic region on timescales of at least several years. Being excited by the more vigorous convection in the outer part of the dynamo region, they can propagate both inwards and outwards. When being reflected at a magnetic tangent cylinder at the transition to the molecular region, they can form standing waves. Identifying such reflections in observational data could determine the depth at which the metallic region effectively begins. Also, this may distinguish 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 of milliseconds. In addition these internal disturbances could give rise to a 10% change over time in the zonal flows at a depth of 3000 km below the surface.

Tables

Table 1. Table 1: Dimensionless input and output parameters of dynamo simulations selected from J14 [ 23 ] . The columns E 𝐸 E , R ​ a 𝑅 𝑎 Ra , and H 𝐻 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 𝑃 𝑟 0.1 Pr=0.1 and P ​ m = 3 𝑃 𝑚 3 Pm=3 , respectively. Constant entropy on both boundaries, r c subscript 𝑟 c r_{\textrm{\scriptsize c}} and r cut subscript 𝑟 cut r_{\textrm{\scriptsize 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 τ 𝜏 \tau , R ​ o 𝑅 𝑜 Ro , and Λ Λ \Lambda list the analysed time interval, the Rossby number (equal to u rms ​ E / P ​ m subscript 𝑢 rms 𝐸 𝑃 𝑚 {u_{\textrm{\scriptsize rms}}E}/{Pm} in our scaling with u rms subscript 𝑢 rms u_{\textrm{\scriptsize rms}} being the rms flow vigour), and the Elsasser number (equal to the rms field strength, B rms subscript 𝐵 rms B_{\textrm{\scriptsize rms}} , in our scaling), respectively, where averages over the whole volume are taken. These are used to yield the Lehnert number L ​ e = B rms / ( ρ eq ​ μ 0 ​ D ​ Ω ) = Λ ​ E / P ​ m 𝐿 𝑒 subscript 𝐵 rms subscript 𝜌 eq subscript 𝜇 0 𝐷 Ω Λ 𝐸 𝑃 𝑚 Le={B_{\textrm{\scriptsize rms}}}/{(\sqrt{\rho_{\textrm{\scriptsize eq}}\mu_{0}}D\Omega)}={\sqrt{{\Lambda E}/{Pm}}} of 6.6-7.3 × 10 − 3 absent superscript 10 3 \times 10^{-3} and the Alfvén number A = u rms ​ ρ eq ​ μ 0 / B rms = R ​ o / L ​ e 𝐴 subscript 𝑢 rms subscript 𝜌 eq subscript 𝜇 0 subscript 𝐵 rms 𝑅 𝑜 𝐿 𝑒 A={u_{\textrm{\scriptsize rms}}\sqrt{\rho_{\textrm{\scriptsize eq}}\mu_{0}}}/{B_{\textrm{\scriptsize rms}}}={{Ro}/{Le}} of 0.45-0.62. The columns U A ​ ( s mtc ) subscript 𝑈 𝐴 subscript 𝑠 mtc {U_{A}(s_{\textrm{\scriptsize mtc}})} , τ A subscript 𝜏 A \tau_{\textrm{\scriptsize A}} , δ ​ σ 𝛿 𝜎 \delta\sigma , u ϕ ′ ¯ ​ ( r cut ) ¯ subscript superscript 𝑢 ′ italic-ϕ subscript 𝑟 cut {\overline{u^{\prime}_{\phi}}(r_{\textrm{\scriptsize cut}})} , and u ϕ ′ ¯ / u ϕ ~ ¯ ​ ( r cut ) ¯ subscript superscript 𝑢 ′ italic-ϕ ¯ ~ subscript 𝑢 italic-ϕ subscript 𝑟 cut {{\overline{u^{\prime}_{\phi}}}/{\overline{\widetilde{u_{\phi}}}}(r_{\textrm{\scriptsize 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.

Run	$E$	$Ra$	$H$	$\tau$	$Ro$	$\Lambda$	$U_{\text{A}}\left(s_{\text{mtc}}\right)$	${\tau }_{\text{A}}$	$\delta \sigma$	$\overline{u_{\varphi }^{\prime }}\left(r_{\text{cut}}\right)$	$\overline{u_{\varphi }^{\prime }}/\overline{\widetilde{u_{\varphi }}}\left(r_{\text{cut}}\right)$
A	$2.5\times {10}^{-5}$	$1.1\times {10}^7$	1.5	0.005	0.0037	5.5	291	0.0016	37.4	772	0.145
E	$1.5\times {10}^{-5}$	$2.0\times {10}^7$	1.4	0.005	0.0041	8.8	421	0.00121	38.7	662	0.106
I	$1.5\times {10}^{-5}$	$2.0\times {10}^7$	1.4	0.005	0.0035	10.3	347	0.00120	42.2	686	0.115

Table 2. Table 2: Dimensional parameters of the simulations. The first column τ unit subscript 𝜏 unit {\tau_{\textrm{\scriptsize 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 9.16 superscript 10 2 9.16\times 10^{-2}  m s -1 , or, equivalently, a radial field strength of 30 G, at the equator at the top of the metallic region ∼ 0.85 ​ R J similar-to absent 0.85 subscript 𝑅 J \sim 0.85\,R_{\textrm{\scriptsize J}} ; calculated as described in the main text. The columns τ J superscript 𝜏 J {\tau^{\textrm{\scriptsize J}}} , τ A J subscript superscript 𝜏 J A {\tau^{\textrm{\scriptsize J}}_{\textrm{\scriptsize A}}} , δ ​ σ J 𝛿 superscript 𝜎 J {\delta\sigma^{\textrm{\scriptsize J}}} , δ ​ P 𝛿 𝑃 \delta P , u ϕ ′ ¯ ​ ( 0.96 ​ R J ) ¯ subscript superscript 𝑢 ′ italic-ϕ 0.96 subscript 𝑅 J {\overline{u^{\prime}_{\phi}}(0.96\,R_{\textrm{\scriptsize J}})} , and u ϕ ~ ¯ ​ ( 0.96 ​ R J ) ¯ ~ subscript 𝑢 italic-ϕ 0.96 subscript 𝑅 J {\overline{\widetilde{u_{\phi}}}(0.96\,R_{\textrm{\scriptsize J}})} present the dimensional version of the analysed interval τ 𝜏 \tau , traveltime τ A subscript 𝜏 A \tau_{\textrm{\scriptsize A}} , maximal angular momentum δ ​ σ 𝛿 𝜎 \delta\sigma , maximal LOD variation, maximal zonal flow fluctuation u ϕ ′ ¯ ¯ subscript superscript 𝑢 ′ italic-ϕ \overline{u^{\prime}_{\phi}} at the surface r cut ∼ 0.96 ​ R J similar-to subscript 𝑟 cut 0.96 subscript 𝑅 J r_{\textrm{\scriptsize cut}}\sim 0.96\,R_{\textrm{\scriptsize J}} , and maximal mean zonal flow u ϕ ~ ¯ ¯ ~ subscript 𝑢 italic-ϕ {\overline{\widetilde{u_{\phi}}}} at the same surface, respectively. Figures in brackets ( ) are based on the longer time unit mentioned in the text: an Alfvén speed of 1.83 × 10 − 3 1.83 superscript 10 3 1.83\times 10^{-3}  m s -1 , or a field strength of 0.6 G, at r ∼ 0.85 ​ R J similar-to 𝑟 0.85 subscript 𝑅 J r\sim 0.85\,R_{\textrm{\scriptsize J}} is assumed.

Run	${\tau }_{\text{unit}}$ [${10}^3$ yrs]	${\tau }^{\text{J}}$ [yrs]	${\tau }_{\text{A}}^{\text{J}}$ [yrs]	$\delta {\sigma }^{\text{J}}$ [${10}^{32}$ N m s]	$\delta P$ [${10}^{-3}$ s]	$\overline{u_{\varphi }^{\prime }}\left(0.96R_{\text{J}}\right)$ [m s-1]	$\widetilde{\overline{u_{\varphi }}}\left(0.96R_{\text{J}}\right)$ [m s-1]
A	6.10	30.5	9.7	2.20 (0.044)	18 (0.35)	0.24 (0.0049)	1.7 (0.033)
E	8.81	44.1	10.7	1.58 (0.032)	13 (0.25)	0.14 (0.0029)	1.4 (0.027)
I	7.26	36.3	13.3	2.09 (0.042)	17 (0.33)	0.18 (0.0036)	1.6 (0.031)

Table 3. Table 1: Dimensionless output parameters of the simulations for the current formulation (cf. table  1 ). The columns U ^ A ​ ( s mtc ) subscript ^ 𝑈 𝐴 subscript 𝑠 mtc {\hat{U}_{A}(s_{\textrm{\scriptsize mtc}})} , τ ^ A subscript ^ 𝜏 A \hat{\tau}_{\textrm{\scriptsize A}} , and δ ​ σ ^ 𝛿 ^ 𝜎 \delta\hat{\sigma} 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.

Run	${\hat{U}}_{\text{A}}\left(s_{\text{mtc}}\right)$	${\hat{\tau }}_{\text{A}}$	$\delta \hat{\sigma }$
A	316	0.0015	38.6
E	452	0.00116	41.45
I	373	0.00111	41.43

Equations

$$\nabla\cdot\rho_{\textrm{\scriptsize eq}}\mbox{\boldmath$u$}=0\;.$$

$$\frac{\partial\mbox{\boldmath$u$}}{\partial t}+(\mbox{\boldmath$u$}\cdot\nabla)\mbox{\boldmath$u$}+2\mbox{\boldmath$\Omega$}\times\mbox{\boldmath$u$}=-\nabla\hat{p}+\frac{1}{\rho_{\textrm{\scriptsize eq}}}\mbox{\boldmath$j$}\times\mbox{\boldmath$B$}-\hat{\mbox{\boldmath$e$}}_{r}\frac{dT_{\textrm{\scriptsize eq}}}{dr}S\;,$$

$$\frac{\partial }{\partial t}⟨{\rho }_{\text{eq}}\overline{u_{\varphi }}⟩$$

$$F_{\text{L}}$$

$$\text{and}{F}_{\text{R}}$$

$$\displaystyle\mbox{\boldmath$u$}(s,\phi,z,t)$$

$$\displaystyle\mbox{\boldmath$B$}(s,\phi,z,t)$$

$$\frac{\partial\mbox{\boldmath$B$}}{\partial t}=\mbox{\boldmath$B$}\cdot\nabla\mbox{\boldmath$u$}-\mbox{\boldmath$u$}\cdot\nabla\mbox{\boldmath$B$}-(\nabla\cdot\mbox{\boldmath$u$})\mbox{\boldmath$B$}\;.$$

$$\frac{\partial^{2}}{\partial t^{2}}\langle\rho_{\textrm{\scriptsize eq}}\overline{u^{\prime}_{\phi}}\rangle={\color[rgb]{0,0,0}\frac{\partial^{2}}{\partial t^{2}}\langle\rho_{\textrm{\scriptsize eq}}(\langle\overline{u^{\prime}_{\phi}}\rangle+\overline{{u^{\prime}}^{\textrm{\scriptsize a}}_{\phi}})\rangle}=\langle\rho_{\textrm{\scriptsize eq}}\rangle\frac{\partial^{2}\langle\overline{u^{\prime}_{\phi}}\rangle}{\partial t^{2}}\;+\;\left\langle\rho_{\textrm{\scriptsize eq}}\frac{\partial^{2}\overline{{u^{\prime}}^{\textrm{\scriptsize a}}_{\phi}}}{\partial t^{2}}\right\rangle\;,$$

$$\frac{{\partial }^2⟨\overline{u_{\varphi }^{\prime }}⟩}{\partial t^2}$$

$$\frac{{\partial }^2}{\partial t^2}\frac{⟨\overline{u_{\varphi }^{\prime }}⟩}{s}=\frac{1}{s^{3}h⟨{\rho }_{\text{eq}}⟩}\frac{\partial }{\partial s}\left(s^{3}h⟨{\rho }_{\text{eq}}⟩U_{\text{A}}^2\frac{\partial }{\partial s}\frac{⟨\overline{u_{\varphi }^{\prime }}⟩}{s}\right)\equiv \frac{\partial }{\partial t}\frac{F_{\text{LR}}}{s⟨{\rho }_{\text{eq}}⟩},$$

$$\delta \sigma =2\pi {\int }_{s_{\text{tc}}}^{s_{\text{mtc}}}{\int }_{z_{-}}^{z_{+}}h⟨{\rho }_{\text{eq}}⟩s^2⟨\overline{u_{\varphi }^{\prime }}⟩dzds,$$

$$\delta {\sigma }_{\text{omtc}}=2\pi {\int }_{s_{\text{mtc}}}^{r_{\text{cut}}}{\int }_{z_{-}}^{z_{+}}h⟨{\rho }_{\text{eq}}⟩s^2⟨\overline{u_{\varphi }^{\prime }}⟩dzds.$$

$$\frac{{\partial }^2}{\partial t^2}\frac{⟨{\rho }_{\text{eq}}\overline{u_{\varphi }^{\prime }}⟩}{s}=\frac{1}{s^{3}h}\frac{\partial }{\partial s}\left(s^{3}h{{\hat{U}}_{\text{A}}}^2\frac{\partial }{\partial s}\frac{⟨{\rho }_{\text{eq}}\overline{u_{\varphi }^{\prime }}⟩}{s}\right),$$

$$\delta \hat{\sigma }=2\pi {\int }_{s_{\text{tc}}}^{s_{\text{mtc}}}{\int }_{z_{-}}^{z_{+}}{s}^{2}h⟨{\rho }_{\text{eq}}\overline{u_{\varphi }^{\prime }}⟩dzds$$

$$\mbox{\boldmath$B$}=B_{0}\hat{\mbox{\boldmath$e$}}_{x}+b_{y}(x)\hat{\mbox{\boldmath$e$}}_{y}\;,\quad\mbox{\boldmath$u$}=u_{y}(x)\hat{\mbox{\boldmath$e$}}_{y}$$

$$\frac{\partial b_y}{\partial t}=B_0\frac{\partial u_y}{\partial x}+\frac{\partial }{\partial x}\eta \frac{\partial b_y}{\partial x}\text{and}\frac{\partial u_y}{\partial t}=\frac{B_0}{{\mu }_0{\rho }_0}\frac{\partial b_y}{\partial x},$$

$$\frac{{\partial }^2{b}_y}{\partial t^2}$$

$$b_y$$

$${\lambda }^2=-\frac{{\omega }^2(V_{\text{A}}^2-\mathrm{i}\omega {\eta }_0)}{V_{\text{A}}^4+{\omega }^2{\eta }_0^2}.$$

$$\lambda =-(1+\mathrm{i})\sqrt{\frac{\omega }{2{\eta }_0}}.$$

$$1+\mathcal{R}=\mathcal{T}\text{and}\mathrm{i}k(-1+\mathcal{R})=-(1+\mathrm{i})\sqrt{\frac{\omega }{2{\eta }_0}}\mathcal{T}.$$

$$\mathcal{R}=\frac{\mathrm{i}k-(1+\mathrm{i})\sqrt{\omega /2{\eta }_0}}{\mathrm{i}k+(1+\mathrm{i})\sqrt{\omega /2{\eta }_0}}.$$

Full text

††journal: Earth and Planetary Science Letters

1 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 [38] 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, gravitational, topographic, and viscous couplings (see [32] for an overview). The evidence for such waves within the Earth’s 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 [16, 17]. They may also account for a decadal variation of the length-of-day (LOD) of Earth [21]. 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.

Here we extend the study of TWs to compressible fluids by applying the anelastic approximation, in which sound waves are ignored. 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]. The 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 $\sim 10\,\textrm{G}$, or $1\,\textrm{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_{\textrm{\scriptsize J}}$, with $R_{\textrm{\scriptsize J}}$ being Jupiter’s mean radius at the 1 bar level. Adopting an interior model including the transition [14], dynamo simulations for anelastic convection have reproduced Jupiter-like magnetic fields [23, 15]. 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 [20, 34]. The rapid rotation and strong magnetic field in the metallic region give rise to a force balance in which the viscous forces are small compared to the Coriolis and Lorentz forces: the quasi-magnetostrophic balance. Jupiter may therefore be 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 have enabled the construction of global models for the magnetic field, such data was limited to spherical harmonics of degree no higher than seven [9, 31]. Ridley & Holme [31] 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 sample the field closer to its source than for any other planetary dynamo so far [3, 4, 24, 10]. 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 [27, 18]. The cloud motion has been tracked for decades by Earth-based telescopes to measure changes to, or periodicities in, the zonal winds [42]. The colouration, brightening, and outbreak events, sometimes leading to global upheavals, have been monitored for more than one hundred years [33, 13].

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, 45, 7], often referred to as torsional oscillations. However, Earth’s core flow inversions/assimilation [17] and numerical geodynamo simulations [43, 39, 37] have found TWs travelling predominantly in an outwards direction with no obvious reflection at the boundaries. This could be explained through preferred excitation [39, 40, 41] 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) [36, 35]. 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 [11]. 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. This leads to the formation of standing torsional waves.

2 Theory

The theoretical framework of incompressible TWs is well documented [e.g. 5, 39, 22]. In the light of those studies, we consider anelastic fluids where the Lantz-Braginsky-Roberts formulation [6, 28, 26] 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 $\rho_{\textrm{\scriptsize eq}}$. The velocity perturbations of the waves $u$ are subsonic, so that the continuity equation is

[TABLE]

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

[TABLE]

where $\Omega$ is the rotational angular velocity, $j$ is the current density, $B$ is the magnetic field, $S$ is the entropy, $T_{\textrm{\scriptsize eq}}$ is the temperature in the equilibrium state, \hat{\mbox{\boldmathe}}_{r} is the unit vector in the radial direction, and $\hat{p}$ is a reduced pressure incorporating the density and the gravitational potential. Hereafter we suppose \mbox{\boldmath\Omega}=\Omega\hat{\mbox{\boldmathe}}_{z} with \hat{\mbox{\boldmathe}}_{z} being the unit vector in the direction of rotation axis. To look at fluctuations corresponding to TWs, we consider the axisymmetric $z$-independent azimuthal flow by taking averages of the $\phi$-component of the momentum equation over cylindrical surfaces to give

[TABLE]

where the azimuthal and axial averages are defined as $\overline{f}=({1}/{2\pi})\int_{0}^{2\pi}f\,d\phi$ and $\langle f\rangle=({1}/{h})\int_{z_{-}}^{z_{+}}f\,dz$, respectively, with $h=z_{+}-z_{-}$ for any scalar field, $f$. Outside the TC, the integral is limited by $z_{\pm}=\pm\sqrt{r_{\textrm{\scriptsize o}}^{2}-s^{2}}$ with $r_{\textrm{\scriptsize 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_{\textrm{\scriptsize C}}=-({\Omega}/{\pi sh})\int\nabla\cdot{\rho}_{\textrm{\scriptsize eq}}\mbox{\boldmathu}\,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_{\textrm{\scriptsize R}}$ are negligible, (3) gives $F_{\textrm{\scriptsize L}}=0$, i.e. the Taylor state for anelastic fluids. The Lorentz and Reynolds forces may be rewritten as

[TABLE]

respectively. The second term of each represents the surface term across an interface between the internal fluid region and the outside, magnetically or dynamically. Since the currents vanish outside the metallic hydrogen zone, the $F_{\textrm{\scriptsize 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_{\textrm{\scriptsize R}}$ surface term vanishes also. However, unlike the magnetic term, the molecular non-conducting region can contribute significantly to the $F_{\textrm{\scriptsize R}}$ integral, because the convection-driven velocities are large there, as we shall see below.

We now make the ansatz of splitting magnetic field and velocity into their mean and fluctuating parts:

[TABLE]

where $\widetilde{f}=({1}/{\tau})\int f\,dt$ with $\tau$ being a time interval, $f^{\prime}=f-\widetilde{f}$, $f^{\textrm{\scriptsize a}}=f-\langle\overline{f}\rangle$, and $\widetilde{f^{\prime}}=0$, $\langle\overline{{f}^{\textrm{\scriptsize a}}}\rangle=0$. The time interval $\tau$ 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

[TABLE]

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 waveforms 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 $\langle\rho_{\textrm{\scriptsize eq}}\overline{u_{\phi}}\rangle/s$ or $\langle\overline{u_{\phi}}\rangle/s$, because the time-derivative of $\langle\rho_{\textrm{\scriptsize eq}}\overline{u_{\phi}}\rangle$ appears in (3), but (6) contains spatial derivatives of $u_{\phi}$, not $\rho_{\textrm{\scriptsize eq}}u_{\phi}$. Here we choose $\langle\overline{u_{\phi}}\rangle$ as the dependent variable: we separate $u^{\prime}_{\phi}$ into its geostrophic and ageostrophic parts,

[TABLE]

and ignore the second ageostrophic term, because it is small compared to the first term in our simulations. We then obtain

[TABLE]

where the non-fluctuating part of the first term of the right hand side sums up to zero, i.e. the Taylor constraint. The theory is equivalent to that of the incompressible case [39] but now the effect of compressibility remains in the Lorentz term. The fluctuating components are assumed to be significantly smaller than the mean parts. Following section 3.1 of [39], we separate the Lorentz terms in (8) into a restoring force part $F_{\textrm{\scriptsize LR}}$ and a driving part $F_{\textrm{\scriptsize LD}}$, where the restoring force part is the term coming from the axisymmetric geostrophic part of $u^{\prime}_{\phi}$ and $F_{\textrm{\scriptsize LD}}$ contains the remaining Lorentz terms. The $F^{\prime}_{\textrm{\scriptsize R}}$ part then corresponds to the driving of the TW by the Reynolds forces. If the Lorentz and Reynolds driving terms are omitted, we obtain the homogeneous free oscillation TW equation

[TABLE]

where ${U_{\textrm{\scriptsize A}}^{2}}={\langle\overline{\widetilde{B_{s}^{2}}}\rangle}/{\mu_{0}\langle{\rho}_{\textrm{\scriptsize eq}}\rangle}$, implying a wave equation for angular velocity in the anelastic case [22]. We note that another possible definition would be $\langle\overline{\widetilde{B_{s}^{2}}}/\mu_{0}\rho_{\textrm{\scriptsize eq}}\rangle$ [see A], but this is less convenient in our formulation. As the restoring force of the wave is represented by $F_{\textrm{\scriptsize LR}}$, the remaining terms of the Lorentz force can be summed up to a term $F_{\textrm{\scriptsize LD}}=F^{\prime}_{\textrm{\scriptsize L}}-F_{\textrm{\scriptsize 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 they are minor in our simulation, so we call $F_{\textrm{\scriptsize LD}}$ a driving term. As we will see below, the waves can also be driven by convective perturbations in the Reynolds force denoted by $F^{\prime}_{\textrm{\scriptsize R}}$.

A perturbation of angular velocity, ${\langle\overline{u^{\prime}_{\phi}}\rangle}/{s}$, can propagate in $s$ with Alfvén speed, $U_{\textrm{\scriptsize A}}$. The speed depends on the magnitude of the background field, ${\widetilde{B_{s}^{2}}}$, and the density, ${\rho}_{\textrm{\scriptsize 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).

3 Numerical simulations

3.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 developed by Dietrich & Jones [12]:see J14 for the detailed description of the model set-up. The models explore the self-generation of magnetic fields by anelastic fluid motions in rotating spherical shells. The equilibrium reference state calculated by [14] was used, and viscous and diffusion terms are taken into account. The reference state density, ${\rho}_{\textrm{\scriptsize eq}}$, electrical conductivity, ${\sigma}_{\textrm{\scriptsize eq}}$, and temperature, ${T}_{\textrm{\scriptsize eq}}$, arise from a composition comprising of a metallic hydrogen region above a rocky core, taken in this model as $r\geq r_{\textrm{\scriptsize c}}\,\sim 6.45\,\times\,10^{6}$ m $\sim 0.09\,R_{\textrm{\scriptsize J}}$, and its continuous transition to a molecular hydrogen region. The transition begins at $r\sim 0.85$-$0.90\,R_{\textrm{\scriptsize J}}$ and only the region below a cut-off level, $r\leq r_{\textrm{\scriptsize cut}}\sim 6.70\times 10^{7}$ m $\sim 0.96\,R_{\textrm{\scriptsize J}}$, is modelled in our simulations, the cut-off being required for numerical reasons. The density scale height, $N_{\rho}=\ln{[{\rho}_{\textrm{\scriptsize eq}}(r_{\textrm{\scriptsize c}})/{\rho}_{\textrm{\scriptsize eq}}(r_{\textrm{\scriptsize 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, ${\sigma}_{\textrm{\scriptsize 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. The constraint is relaxed in the conducting region and 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) [12], located at $s\sim 0.85$-$0.90\,R_{\textrm{\scriptsize J}}$. This is in addition to the traditional ’kinematic’ TC found at s=r_{\textrm{\scriptsize c}}{\color[rgb]{0,0,0}=0.0963\,r_{\textrm{\scriptsize cut}}}\equiv s_{\textrm{\scriptsize 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: for our purposes, we denote $s_{\textrm{\scriptsize mtc}}$ as the minimal $s$ at which the magnetic diffusion term becomes comparable to the other terms in (9), which in our model is at $\sim 0.89\,r_{\textrm{\scriptsize cut}}$. 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_{\textrm{\scriptsize tc}}\leq s\leq s_{\textrm{\scriptsize 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\times 10^{-3}$. The Elsasser number, $\Lambda$, is a dimensionless measure of the magnetic field strength and is found to be approximately 6-10 in our simulations. These yield Alfvén numbers ranging from 0.45 to 0.62 (see the caption of table 1) and the Alfvén speed is faster than the rms velocity overall. 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 [10] 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_{\textrm{\scriptsize A}}$, as a function of cylindrical radius, $s$, normalised by the cut-off radius, $r_{\textrm{\scriptsize cut}}$, for model I. Here the time and length are scaled by the magnetic diffusion time and the shell thickness ($D=r_{\textrm{\scriptsize cut}}-r_{\textrm{\scriptsize c}}$), respectively, and the bounds for $z$-averages are taken at $r_{\textrm{\scriptsize cut}}$. In the figure, the dashed line represents the Alfvén speed with the density taken to be its constant mid-radius value $\rho_{\textrm{\scriptsize eq}}(r_{\textrm{\scriptsize c}}/2+r_{\textrm{\scriptsize cut}}/2)$ in the definition of $U_{A}$. The anelastic Alfvén speed has a peak at ${s}/{r_{\textrm{\scriptsize cut}}}\sim 0.6$. At ${s}/{r_{\textrm{\scriptsize cut}}}\lesssim 0.6$, the density ${\rho}_{\textrm{\scriptsize eq}}$ decreases with radius $r$ and the $z$-mean Alfvén speed increases with $s$ because ${U_{\textrm{\scriptsize A}}^{2}}={\langle\overline{\widetilde{B_{s}^{2}}}\rangle}/{\mu_{0}\langle{\rho}_{\textrm{\scriptsize eq}}\rangle}$. At ${s}/{r_{\textrm{\scriptsize cut}}}\gtrsim 0.6$, the density decrease effect is countered by the drop in $\langle\overline{\widetilde{B_{s}^{2}}}\rangle$ due to the field morphology, so for larger $s$, the speed gradually decreases as the MTC is approached and crossed. Profiles of $U_{A}$ are similar for the other simulations explored here, with peaks at $0.6\lesssim{s}/{r_{\textrm{\scriptsize cut}}}\lesssim 0.7$. Table 1 also lists the speeds ${U_{\textrm{\scriptsize A}}(s_{\textrm{\scriptsize mtc}})}$ at the MTC radius and the expected traveltimes $\tau_{\textrm{\scriptsize A}}$ from the core boundary $s_{\textrm{\scriptsize tc}}$ to the $s_{\textrm{\scriptsize mtc}}$. The speeds $U_{\textrm{\scriptsize A}}$ are used for conversion to our dimensional time unit below (see details in sec. 4.1): a Jovian scale ${U_{\textrm{\scriptsize A}}^{\textrm{\scriptsize J}}}$ is shown on the right-hand side of the axis in fig. 1.

3.2 Internal dynamics: zonal flow fluctuations and their excitation

The time averaged components of azimuthal velocity, $\langle\overline{\widetilde{u_{\phi}}}\rangle$, 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, $\langle\overline{u^{\prime}_{\phi}}\rangle$, which are of interest here. Figure 2 displays contours of $\langle\overline{u^{\prime}_{\phi}}\rangle$ in $s$-$t$ space for the three runs. In each diagram, white curves indicate the calculated Alfvén speed, $U_{\textrm{\scriptsize A}}$, to compare with the computed fluctuations. A dimensional time is shown on the top of each image (details in the following section). Run A (fig. 2a) shows that some disturbances emerge near ${s}/{r_{\textrm{\scriptsize cut}}}\gtrsim 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_{\textrm{\scriptsize A}}$, suggesting that they are anelastic torsional Alfvén waves. They become more evident when filtered (see B). Travelling TWs are found in Earth-like Boussinesq models [43, 39, 40, 37, 41]; 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 $\langle\overline{u^{\prime}_{\phi}}\rangle$ 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_{\textrm{\scriptsize cut}}\gtrsim 0.6$. Interestingly, beneath the MTC, waves appear to form a node at around $s/r_{\textrm{\scriptsize cut}}\sim 0.65$ and $\sim 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\gtrsim 0.003$. There are signatures of reflection, highlighted by the white lines, around the MTC at, for instance, $t\sim 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 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 $\langle\overline{u^{\prime}_{\phi}}\rangle$, 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 [35].

In model I, where the entropy flux at $r_{\textrm{\scriptsize 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\lesssim t\lesssim 0.0010$. Note that poor resolution of observational data and/or improper filters over them may make the reflecting nature of the waves less clear (see B). To describe the time evolution, we also present profiles of $\langle\overline{u^{\prime}_{\phi}}\rangle$ in figure 3. A trough came into existence at $t\sim 0.0006$ and ${s}/{r_{\textrm{\scriptsize cut}}}\sim 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_{\textrm{\scriptsize LD}}$ and/or $F_{\textrm{\scriptsize R}}$. At $t\sim 0.00095$, the trough reflects at around $s_{\textrm{\scriptsize 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 the shape vary in our nonlinear simulation, and it is hard to determine the reflection coefficient or the phase shift accurately. An abrupt change in a $U_{\textrm{\scriptsize 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 [11]. 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^{\prime}_{\textrm{\scriptsize R}}$ and $F_{\textrm{\scriptsize LD}}$, respectively, for run I. The Reynolds force is found to be important in the outer regions, $0.6\;r_{\textrm{\scriptsize cut}}\lesssim s\lesssim s_{\textrm{\scriptsize 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^{\prime}_{\textrm{\scriptsize R}}$ better matches the location and time at which $\langle\overline{{u^{\prime}}_{\phi}}\rangle$ disturbances begin to travel than $F_{\textrm{\scriptsize LD}}$. This is in spite of the small global Rossby number; such an initiation was pointed out in Boussinesq cases [39]. 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 $\rho u_{s}u_{\phi}$ term, and continuously forces fluctuations, $\partial u_{\phi}/\partial t$, by almost-hydrodynamic Rossby waves, which are found to be faster than the Alfvén waves in the simulations by at least a factor 10 (not shown). TWs in model A are also predominantly excited by $F^{\prime}_{\textrm{\scriptsize R}}$; it is rather mixed with Lorentz terms $F_{\textrm{\scriptsize 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 [39]. However, as parameters are moved towards their Earth-like values, geodynamo and magnetoconvection simulations [40] 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 [40]. 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.

4 Application to Jupiter

4.1 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. Current numerical models are limited to numerically accessible parameters [37, 2], so parameters relating to diffusive processes have artificially increased values. No choice of dimensional units correctly scales all the physical processes involved, so we choose scalings which aim to get the most important aspects of TWs right. We choose the magnetic field scale by equating the field strength at the magnetic outer boundary in simulations to the observed outer boundary value [39]. 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\sim 0.85\,R_{\textrm{\scriptsize J}}$: in the equatorial region it is seemingly no greater than 30 G and weaker than 1 G for large regions [10]. Taking 30 G as a reasonable maximal field magnitude at our MTC radius, $s_{\textrm{\scriptsize mtc}}$, and the density, ${\rho}_{\textrm{\scriptsize eq}}$, of $8.53\times 10^{2}$ kg m*-3*, then an Alfvén speed, $U_{\textrm{\scriptsize A}}^{\textrm{\scriptsize J}}$, at this radius is approximately $9.16\times 10^{-2}$ m s*-1*. By matching this value with those of our simulations, our dimensional time unit $\tau_{\textrm{\scriptsize unit}}$ is calculated through ${DU_{\textrm{\scriptsize A}}}/{U_{\textrm{\scriptsize A}}^{\textrm{\scriptsize J}}}$, where $D$ is the shell thickness of $6.06\times 10^{7}$ m. From this we calculate dimensional versions, ${\tau^{\textrm{\scriptsize J}}}$ and ${\tau_{\textrm{\scriptsize A}}^{\textrm{\scriptsize J}}}$, of the analysed interval, $\tau$, and the TW traveltime, ${\tau_{\textrm{\scriptsize 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 analysed time windows ${\tau^{\textrm{\scriptsize J}}}$ from 31 to 44 years), the traveltimes ${\tau^{\textrm{\scriptsize J}}_{\textrm{\scriptsize 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 at depth should be scaled. Taking the unit of velocity as ${D}/{\tau_{\textrm{\scriptsize unit}}}$ gives $3.15\times 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 mean zonal flow at $0.96\,R_{J}$ of about $2$ m s*-1* (table 2). If we use a longer 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. We prefer the shorter time unit, as we believe that future less diffusive models will have a higher ratio of zonal flow to convective flow, allowing a convective velocity of $10^{-2}$ m s*-1* with a zonal flow of $\sim 2$ m s*-1* at $0.96\,R_{J}$.

4.2 Length-of-day variation (LOD)

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 $\mathcal{O}(10^{-4}\,\textrm{s})$ has been identified; its origin could be angular momentum exchange between the fluid core and the rocky mantle through MHD waves ([16]; 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,

[TABLE]

and those outside the region,

[TABLE]

In figure 5 the solid and dotted curves display the time evolutions of $\delta\sigma$ and $\delta\sigma_{\textrm{\scriptsize omtc}}$ in model E, respectively. The $\delta\sigma$ of the conducting region shows a quasi-periodic variation, corresponding to the flow oscillations (fig. 2b). The evolution is almost perfectly anti-correlated with the change $\delta\sigma_{\textrm{\scriptsize 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 both the magnetic and dynamic terms play a role (not shown); it is however uncertain how the coupling arises. TWs with a short wavelength in the $s$-direction will be damped out by magnetic diffusion as soon as they leave the metallic hydrogen region, but longer wavelength TWs are damped less rapidly as a single wavelength could extend right across the transition region (e.g. C), allowing the TWs to be seen at the surface of the planet.

Using our standard time unit and the density ({\color[rgb]{0,0,0}{{\rho}_{\textrm{\scriptsize eq}}(r_{\textrm{\scriptsize c}}/2+r_{\textrm{\scriptsize cut}}/2)}}=2.56\times 10^{3}\,\textrm{kg}\,\textrm{m}^{3}), we convert the dimensionless $\delta\sigma$ of maximum amplitude 38.7 to its dimensional version, $\delta\sigma^{\textrm{\scriptsize J}}$ with maximum amplitude $\sim 1.58\times 10^{32}$ N m s. Assuming a value of $I=2.56\times 10^{42}\,\textrm{kg}\,\textrm{m}^{2}$ for the moment of the inertia [30] and a daily period of $P=3.57\times 10^{4}$ s for the planet, the change $\delta\sigma^{\textrm{\scriptsize J}}$ is equivalent to a period $\delta P$ of approximately 13 ms. Here $\delta\sigma^{\textrm{\scriptsize J}}=I\delta\Omega={-2\pi I\delta P}/{P^{2}}$, where $\Omega$ 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 [20, 34, 31]. 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.

4.3 Flow change above the metallic region

Unlike terrestrial planets, gas giants may allow deep-origin perturbations to be observed at the surface. Figures 6a and b show contours of the fluctuating zonal flow $\overline{u^{\prime}_{\phi}}$ on the cut-off surface, $r_{\textrm{\scriptsize 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. The amplitude is scaled by the maximum of the mean zonal flow $\overline{\widetilde{u_{\phi}}}$ on the same surface, which is the maximal speed of the prograde, equatorial jet reproduced in the simulation (sec. 3.2). Figs. 6c-d show 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 visualise 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. Oscillations and both equatorward and poleward propagation are seen at mid and high latitudes where $s\lesssim s_{\textrm{\scriptsize mtc}}$, whereas the equatorial region $s\gtrsim s_{\textrm{\scriptsize mtc}}$ features only equatorward migration. We interpret this as partial transmission through the MTC and absorption within the resistive, transition layer (sec. 3.2; C). The 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 $\theta$-$t$ space - sometimes called a Hovmöller diagram - is essential for the identification.

On the surface of our cutoff level the maximum of the fluctuating velocity in model E 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_{\textrm{\scriptsize cut}}\sim 0.96\,R_{\textrm{\scriptsize J}}$ are found to be about 0.1-0.2 m s*-1* and 3-5 $\times 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).

The surface zonal flows may extend deep into the interior [e.g. 8, 25, 19]. Jupiter’s gravitational harmonics have recently been obtained by the Juno mission, Kaspi et al. [27] and Guillot et al. [18], providing evidence that the zonal flows do go down to $\sim$ 0.95-0.97 $R_{\textrm{\scriptsize J}}$. To date, observational constraints on the speed of deep zonal flows are weak: [18] show that the zonal flow in the deep interior must be less than 10% of the cloud-level value, while [27] suggest that the zonal flows could fall off exponentially with an e-folding depth of 1000-3000 km to be consistent with the Juno gravity data. This is compatible with our higher velocity conversion range of $\sim 1.5$ m s*-1* (see ${\overline{\widetilde{u_{\phi}}}}$ 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_{\textrm{\scriptsize J}}$), not far from those estimated with scaling properties based on available heat fluxes [31, 44]. Those deep dynamics may set the thermodynamical conditions at the cloud deck and 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 recurrent activities spreading over several latitudinal bands and occur at intervals of decades, irregularly in most cases [33, 13]. At some epochs, 5- or 10- years periodicities in jetstream outbreaks or fades/revivals were recognised at the North Temperate Belts (NTB), 23-35∘ N. Recent datasets - primarily 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 [42]. These seem interesting in comparison to the internal flows simulated earlier. First, the amplitude of such disturbances is $\sim$ 10 %, or less, of the 150 m s*-1* stable jet [42]. Second, the NTB latitudes correspond to cylindrical radii $s\sim$ 0.82-0.92 $R_{\textrm{\scriptsize 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.

5 Concluding remarks and discussion

We have demonstrated, through our anelastic models, that torsional Alfvén waves could be excited in Jupiter’s metallic hydrogen region. 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 the outermost part of the conducting region, where nonaxisymmetric convective motions were vigorous and the resulting stresses drove the axisymmetric fluctuations. TWs were found to travel both outwards and inwards. Modes propagating outwards were found to be partially transmitted to the poorly conducting layer but could also be reflected around the MTC. This results in waves travelling inwards, back into the deeper interior of the conducting region. Since convection perturbs the fluid at all times it is able to continuously supply a source for TWs travelling in both directions. If their amplitude and timing matches, a superposition of the opposed propagation enables the formation of standing waves, as observed in our model E. Our simulations suggest there may be a mixture of travelling, reflecting, and standing waves in giant planets. Our results suggest TWs in giant planets behave rather differently from those in the Earth’s core [43, 39, 40, 37, 41]. In geodynamo 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 magnetically-dissipative fluid layer, which enables reflection as well as transmission. 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 an exchange of angular momentum between the metallic hydrogen region and the overlying molecular regions. With the time units adopted, the waves could give rise to variations in the LOD no greater than $10^{-2}$ s. Alterations in Jupiter’s radio rotation period might partly be due to true LOD changes as well as the magnetic SV. We also note that uncertainties in the chosen scaling, which arose from the limitation of the current numerical dynamo models, could affect our LOD variation estimates.

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 at $r\sim 0.96R_{J}$. The reflecting and/or transmitting nature across the MTC could be projected upwards from the metallic region to the surface. Juno’s gravity measurements have constrained interior models by identifying that visible surface zonal flows penetrate downwards significantly [27, 18]. Assuming this deep origin, variations at the cloud deck could display some evidence of TWs.

Another possible way to detect TWs is from 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 the context of Earth’s fluid core. TWs may contribute to the occurrence of geomagnetic jerks; they cannot however account for all phenomena alone (see [29] for a review). 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.

Acknowledgments

We acknowledge support from the Japan Society for the Promotion of Science (JSPS) under Research Activity Start-up No. 17H06859, as well as from the Science and Technology Facilities Council of the UK, STFC grant ST/N000765/1. This work was undertaken on ARC2, part of the High Performance Computing facilities at the University of Leeds, UK. Also this work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This work was partly made during the visit at the Lorentz Center Leiden. We also thank Amy Simon for helpful discussions. Comments by Thomas Gastine and Johannes Wicht helped to improve the manuscript.