Magnetohydrodynamical torsional oscillations from thermo-resistive instability in hot jupiters

Working under the MHD approximation (Davidson, 2001), we consider a geometrically-simplified setup within the equatorial plane, depicted schematically in Figure 1. Our solution domain is meant to represent the upper atmosphere of HJs, going from $1$ bar at the base to $0.01$ bar at the top. For a planet with equilibrium temperature $T_{\mathrm{e}q}=1000\ {\rm K}$ and surface gravity $g_{p}=9.0$ m s^-2, this corresponds to a radial extent $r_{0}=0.940$–$0.966~R_{p}$ assuming hydrostatic balance. For simplicity we assume an ideal gas of pure molecular hydrogen. As the layer is geometrically thin, we adopt the plane-parallel approximation in which cartesian coordinates $(x,y,z$) map to the zonal, latitudinal and radial directions in spherical coordinates. We solve the zonal ($x$) components of the MHD equations in the equatorial plane ($y=0$), assuming axisymmetry and that the flow is confined to the equatorial plane. We thus suppress all radial and latitudinal flow displacements, i.e., $u_{y}=u_{z}=0$.

We express the axisymmetric magnetic field as

thus ensuring $\nabla\cdot{\bf B}=0$, and with $B_{0}$ setting the (fixed) strength of the radial magnetic component in the equatorial plane. The assumption of a radial background magnetic field allows us to write down the simplest model possible in the equatorial plane that captures the dynamics of Alfvén oscillations (Figure 1). A global aligned dipole lacks any radial component at the equator, and would require adding the latitudinal dimension to capture Alfvén oscillations. Of course, the true field geometry is likely to be more complicated than an aligned-dipole. A misaligned dipole or higher order multipoles would give radial field at the equator, and the actual magnetic field geometry could also be much more complex if a local dynamo operates in the atmospheric layer (Rogers & McElwaine, 2017; Dietrich et al., 2022). We adopt this simplified model as a starting point to capture the physics of the interaction of the temperature-dependent MD with magnetic torques.

The reduced MHD equations solved are thus:

Equation (2) is the $x$-component of the incompressible Navier-Stokes equation where $\bar{\rho}$ is the density, and $\bar{\mu}$ is the dynamic viscosity. The fluid motions are forced by an acceleration term

with the pressure $P$ measured in bars and $\dot{v}$ setting the peak amplitude. This forcing is meant to capture the effects of angular momentum transfer to the equator from higher latitudes, for example by Rossby and Kelvin waves (Showman & Polvani, 2011).

Equation (3) is the $x$-component of the induction equation with $\eta$ the temperature-dependent MD,

taken from Menou (2012a) and based on the results of Draine et al. (1983). The ionization fraction $\chi_{e}$ is obtained from the Saha equation (Rogers & Komacek, 2014) adopting solar abundances as given in Lodders (2010). Equation (4) is the heat transport equation with specific heat capacity $c_{p}$ and thermal conductivity

where $\sigma$ is the Stefan-Boltzmann constant. The irradiation flux is

where $F_{s}$ is the flux from the host star incident on the atmosphere. In equations (7) and (8), $\kappa_{\rm th}$ is the Rosseland mean opacity (assumed constant), dominated by opacity at infrared-wavelengths, where most of the heat from the planet is emitted, while $\kappa_{v}$ is the visible opacity, where the spectral flux of the host star peaks. The thermal conductivity is used to set the dynamic viscosity via an assumed Prandtl number $\mathrm{Pr}$:

The last two terms on the RHS of Equation (4) are the viscous ($\phi_{\nu}$) and ohmic ($\phi_{\eta}$) dissipation functions.

We impose corotation with the core of the planet at the base, $u_{x}=0$, and set vanishing mechanical stresses at the top boundary, $\partial u_{x}/\partial z=0$. For the magnetic field, we set $\partial B_{x}/\partial z=0$ at the lower boundary, so as to ensure vanishing magnetic stress. At the outer boundary, we opt for a vacuum boundary condition, implying $B_{x}=0$. For the temperature, we emulate a constant flux of energy coming from the core, $F_{0}$, and a flux coming from the host star, $F_{s}$, which is absorbed at depths controlled by the visible opacity, $\kappa_{v}$, as prescribed by Equation (8). We therefore set the inner temperature boundary condition at a constant flux corresponding to the sum of the two flux sources. The upper boundary is kept at a constant temperature $T_{0}$.

We set the heat flux from the core as $\sigma T_{\mathrm{int}}^{4}$ where $T_{\rm int}=150$ K for all models, following M12, and then parametrize the thermal profile with $T_{\rm eq}$ and $\kappa_{\rm th}$. For each choice of parameters, we set the values of $F_{s}$, $T_{0}$, and $\kappa_{v}$ by comparing with the temperature profiles presented in M12 which are given as a function of $T_{\rm eq}$. We solve the steady state form of Equation (4) with only the first three terms, i.e. we ignore any viscous and ohmic heating, jointly with the hydrostatic balance equation, using a ODE solver from SciPy¹¹1The function is solve_ivp from the integrate module. The documentation can be found at https://docs.scipy.org/doc/scipy/reference/generated/ scipy.integrate.solve_ivp.html, a Python package, using an explicit Runge-Kutta method of order 5. We then determine the values of $T_{0}$ and $F_{s}$ as a function of $T_{\mathrm{eq}}$ by fitting the temperature profiles shown in M12, and adopt a value $\kappa_{v}=6.93\times 10^{-4}$ m² kg^-1 for all models that reproduces the pressures at which the heat is deposited in M12.

The same initial conditions are used throughout: solid body rotation at the lower boundary rate $u_{x}=0$, and purely radial magnetic field, i.e. $B_{z}=B_{0}$ and $B_{x}=0$. The initial temperature profile is found jointly with the boundary conditions, as described in the previous paragraph, but we also let the solution from the ODE solver relax on our simulation grid. During this relaxation process, the barred variables in Equations (2)–(4) are allowed to adjust to the dissipationless steady-state temperature profile, but remain fixed thereafter. The initial MD is then simply given by Equation (6).

The numerical solution of Equations (2)–(4) requires the MD $\eta$ and dissipation functions $\phi_{\nu}$, $\phi_{\eta}$ to be advanced in time. Because of the strongly nonlinear dependency of $\eta$ on $T$, and of the dissipation functions $\phi_{\nu}$ and $\phi_{\eta}$ on flow and field gradients, simply evaluating these quantities at the current (known) time step $n$ when advancing Equations (2)–(4) to $n+1$ can lead to runaway numerical instabilities. Consequently, we opted to subdivide each time step into four substeps, as depicted schematically in Figure 2. In substep 1, we first advance the magnetic field and azimuthal velocity from time step $n$ to $n+1$ using a mid-point value for the MD, so as to ensure second-order accuracy in the spirit of centered finite differences. Equations (2) and (3) are advanced using a Crank-Nicolson scheme for the linear terms and Euler explicit for the nonlinear term. In substep 2, we use these newly calculated values of $B$ and $u$ at $n+1$ to do a mid-point ($n+1/2$) evaluation of the viscous and ohmic dissipation functions, namely the last two terms on the RHS of Equation (4). Substep 3 then advances $T$ from $n$ to $n+1$ via Crank-Nicolson, with mid-point evaluation of all source terms. The final substep 4 extrapolates $T$ at mid-step, from $n+1$ to $n+3/2$, via a truncated Taylor series, i.e., at every grid point $k$ we write

the temperature derivative at $n+1$ being computed by evaluating the RHS of Equation (4) at $n+1$ for the linear terms, and $n+1/2$ for the dissipation functions. The MD at $n+3/2$ is then computed via Equation (6). This completes the $n\to n+1$ time step, and the process can now repeat anew to advance from $n+1$ to $n+2$.

Additionally, due to the boundary conditions imposed on the magnetic field, the ohmic heating is at its maximum at the top boundary. This raises issues, as we are keeping the temperature fixed at this end of the domain. To avoid buildup of a very thin thermal boundary layer, we have introduced a progressive ramp-down to zero of the viscous and ohmic heating terms in the outer 20% of the grid, following a $\propto-z^{3}$ dependency so as to maintain continuity of the second derivative in temperature. We have verified that this procedure does not affect significantly the TRI behavior or character of the resulting Alfvénic oscillations.

We use 100 grid points equally spaced in $z$ for all simulations. This was found to provide an acceptable compromise so as to adequately resolve boundary layers without becoming computationally prohibitive. Our time-stepping algorithm uses a constant time step which allows for suitable time resolution during the burst phase. To achieve this, we have used $\Delta t=10$ s for our simulations. We therefore need a few million time steps to complete a simulation. This represents around 2 core-days on a 4.90 GHz CPU.

The free parameters for the simulations are the radial magnetic field strength, $B_{0}$, the velocity forcing amplitude, $\dot{v}$, the thermal opacity, $\kappa_{\rm th}$ and the equilibrium temperature, $T_{\rm eq}$. We opted to keep the Prandtl number small and constant throughout all simulations as ohmic heating is the focus of this study and HJ atmospheres are highly inviscid in our pressure range (Rogers & Komacek, 2014; Beltz et al., 2022). Table 1 lists the ranges in different parameter values that were used in our parameter space exploration.

The choice of parameter values is guided by the results of H22. H22 presented a local model with a similar magnetic field and flow geometry to that considered here, that gave estimates for timescales and instability criterion as a function of pressure. Rewritten in dimensional units, the instability criterion that H22 derived is

where $E_{\eta}$ is the temperature-sensitivity of the MD, as defined in H22. We compute the instability criterion of Equation (11) across the domain, under the initial conditions, to ascertain in a first approximation which free parameter combinations satisfy Equation (11) and may result in an unstable system. We find that the H22 instability criterion is a good predictor of instability in our 1D model, and the TRI is robustly present. The fact that radial variations are included slightly changes the region of parameter space conducive to the TRI as compared to the local model results of H22. In addition, the behaviour is more complex since the local model parameters vary significantly across the radial direction, and velocity and magnetic field perturbations interact across a large pressure range.

We first examine in detail a representative simulation exhibiting recurrent oscillatory bursts triggered by the TRI. Figure 3 shows the time series of velocity, magnetic field perturbation, temperature and magnetic Reynolds number at the center of the domain through a series of seven recurrence cycles of oscillatory bursts (reminiscent of Figure 1 (c) of H22). The local magnetic Reynolds number is defined as ${\rm Rm}=u_{x}H/\eta$, where $H$ is the local pressure scale height. The close-ups on the right hand side of Figure 3 focus on a single burst, and highlight the Alfvénic oscillations and their signature $\pi/2$ phase lag between velocity and magnetic field, as well as the exponential relationship between Rm and temperature coming from the ionization fraction, $\chi_{e}$ (c.f. Equation (6)). As pointed out by H22, the bursting behaviour is associated with a transition from $\mathrm{Rm}<1$ during the build-up phase between bursts to $\mathrm{Rm}>1$ during the oscillatory phase.

Figure 4 shows the phase space $(u_{x},B_{x},T)$ trajectory of the center point of the domain. The system starts at low velocity and low temperature and moves away from the $u_{x}=B_{x}=0$ plane as the velocity increases. The trajectory veers abruptly upwards once the critical temperature for TRI is reached. The initial runaway in temperature is as discussed by M12, but is halted, ultimately, by the dynamical back-reaction of the magnetic field on the flow, similar to the behavior seen in the local model of H22. After the oscillations have damped and the system cooled towards its equilibrium temperature, slow acceleration resumes, building up shear and inducing an azimuthal magnetic field, eventually triggering TRI and another oscillatory burst. A marked separation of timescale exists between the slow buildup phase, typically lasting weeks to months and in which the system spends most of its time, and the very rapid oscillatory burst and decay, spanning a few days.

Figure 5 displays the temporal evolution of the $T$–$P$ profile over one recurrence period. The profile changes rapidly once the TRI is triggered. The temperature at the higher pressures can exceed 3000 K at burst peak, as compared to $\sim$1300 K between bursts²²2In this example and other cases that reach peak temperatures close to $3000\ \mathrm{K}$, the temperature gradient is steep enough that convection would be expected to occur (not included in our calculation). However, the radial temperature gradient is likely enhanced by the fixed temperature outer boundary condition, and would be reduced with a more realistic outer boundary.. Such temperatures are expected at these pressures for equilibrium atmospheres with $T_{\rm eq}\approx 2500$ K.

We find that the (time-evolving) gradient of MD strongly influences the character of bursts and oscillations developing in the simulations. Figure 6 shows time series for three simulations differing only in their treatment of the diffusivity. In blue, we show a simulation based on the equations as presented in §2.1. In green, we show a simulation without the spatial gradient of the MD, i.e. the second term on the right hand side of Equation (3) has been artificially set to zero. In orange, we show a simulation in which the MD is time-independent, fixed to its initial profile.

While all three simulations show the same initial growth in velocity, the fixed diffusivity simulation (orange) simply grows to a steady-state, as in the local model of H22 and akin to most previous works treating magnetic effects as a drag force. The other two simulations, with time-varying MD, both develop oscillations, but bursting behavior only occurs if the (time-evolving) spatial MD gradient is retained. This can be traced to large spatial gradients of magnetic fields building up in the domain, enhancing dissipation and leading to very rapid rise of the temperature (the “burst”), subsequently transitioning to decaying Alfvén waves. These large magnetic gradients, building up in the lower part of the domain, can be seen in Figure 7 (top row), where time-evolving profiles of the velocity and magnetic field are shown during a single post-burst Alfvénic oscillation. Such strong gradients are absent in the simulation with the $\nabla\eta$ term artificially zeroed out (bottom row). The velocity and magnetic field also reaches much smaller values when the $\nabla\eta$ term is removed. For example, at the center of the domain the peak magnetic field is $\approx 5\ \mathrm{kG}$ for the full simulation (Figure 3), but drops to $50\ {\rm G}$ when the $\nabla\eta$ term is removed, and $4\ {\rm G}$ for the fixed $\eta$ case.

With our geometrical setup and simplifications (i.e., ${\bf B}\equiv B_{x}(z){\hat{x}}+B_{0}{\hat{z}}$ and $\nabla\eta\equiv\,$d$\eta/$d$z$), the gradient term $-\nabla\eta\times(\nabla\times{\bf B})$ associated with the ohmic dissipation term in the $x$-component of the induction equation becomes equivalent to $(\nabla\eta\cdot\nabla)B_{x}$, which is structurally akin to advection of $B_{x}$ by a “flow” $\nabla\eta$. This is sometimes referred to as diamagnetic pumping. In our simulations the field is advected inwards, as $\nabla B_{x}$ is negative during build-up, while $\nabla\eta$ is positive. When the Alfvén oscillation starts, large gradients build up near the inner boundary, which lead to the very sharp increase in temperature due to a magnetic field varying very rapidly with depth, producing large amount of ohmic heating.

We find that certain choices of parameters are needed for sustained oscillations/bursts. For example, they occur only for equilibrium temperatures $T_{\mathrm{eq}}\lesssim 1200\ {\rm K}$. This can be understood in terms of the magnetic Reynolds number, since a key characteristic of sustained oscillations is that the system needs to move between weak and strong flux freezing, i.e. between ${\rm Rm}<1$ and ${\rm Rm}>1$ (H22). Figure 8 shows an estimate of the magnetic Reynolds number in the $T$–$P$ plane (for this estimate we assume a characteristic velocity of 5 km s^-1, a rough average of the velocities in our oscillating solutions). Figure 8 suggests that the maximum equilibrium temperature allowing TRI-driven oscillatory bursts is around 1200 K, since hotter models have $\mathrm{Rm}>1$ throughout the atmosphere (in agreement with Dietrich et al. (2022), who also find Rm$=1$ at around $1300$ K). Any system above this equilibrium temperature will be characterized by strong flux freezing from the beginning and will never be able to oscillate around ${\rm Rm}=1$. Conversely, colder atmospheres below around 1000 K may never reach ${\rm Rm}>1$ as they do not experience enough heating, at least under reasonable velocity forcing.

Figure 9 shows the different regimes in the $T_{\mathrm{eq}}$–$B$ plane. We plot the recurrence period of the TRI, and take $\dot{v}=0.01$ m s^-1 and $\kappa_{\rm th}=0.0005$ m² kg^-1. Instability and resulting oscillations are favored by weaker radial magnetic fields and lower equilibrium temperatures. The cooler temperatures are expected from Figure 8. Weaker radial fields require less forcing to be stretched azimuthally against the opposing Lorentz force, which is necessary to generate electrical currents and ohmic heating, the main heating contributor. As expected from H22, we have identified two different regimes in this figure: sustained oscillations, and systems that reach steady-state either by decaying Alfvénic oscillations or directly reach steady-state. We can see the recurrence period of the oscillations steadily shrinking as the systems approach the larger field and temperature values, even across the black line delimiting the two regimes.

In the simulations presented so far, we assumed for simplicity that the opacity $\kappa_{\mathrm{th}}$ has a constant value. In fact, the opacity varies strongly throughout the atmosphere, as can be seen in Figure 10 where we plot the opacity in the $T$–$P$ plane using the Freedman et al. (2008) opacity tables for dust free gas at solar metallicity. Comparing with the typical atmospheric temperature profiles shown in the Figure, we see that the opacity varies by an order of magnitude within our pressure range, and would also be changing locally during a burst as well as temperature changes.

To check the effect of opacity variations, we ran a simulation with a depth-dependent opacity profile, but still time independent, using the initial temperature profile to calculate $\kappa_{\mathrm{th}}$ as a function of depth. We used the following parameters: $B_{0}=30$ G, $T_{\rm eq}=1000$ K, $\dot{v}=0.01$, $\mathrm{Pr}=0.01$. As a comparison baseline, we also ran a simulation with a constant thermal opacity of $\kappa_{\rm th}=0.0005$ m² kg^-1 ($\log_{10}\kappa_{\rm th}=-3.3$), which corresponds to the opacity at around 0.5 bar for this equilibrium temperature. We find that the overall effect of the depth-dependent opacity is that the outer part of the envelope is not as insulating as with a constant opacity. Therefore, higher velocities are needed to reach the critical temperature, thus increasing the recurrence time of the TRI. Otherwise, the two simulations are qualitatively very similar.

HJs may have larger metallicities than solar (Welbanks et al., 2019). Freedman et al. (2008) also presents the values of the mean Rosseland opacity for a metallicity of $[M/H]=0.3$. On average, for temperatures above 1000 K, the opacity with enhanced metallicity is larger by a factor of 1.7. Furthermore, the opacity presented by Freedman et al. (2008) are for dust and cloud free gas. However, the temperature range at which the TRI operates is in a regime where both dust and clouds are expected, which would enhance the opacity of the atmosphere. The opacity presented in Figure 10 is to be considered as a lower limit of the expected opacity in the atmosphere of a HJ.

Moreover, the opacity enters the calculation of the thermal conductivity $\bar{\chi}$ (see Equation 7), which already varies with depth by a few orders of magnitude via its dependence on density and temperature. Introducing a depth dependence on opacity only adds another order of magnitude in the range of thermal conductivity ($\bar{\chi}\sim 1/\kappa_{\rm th}\bar{\rho}$; Equation 7), and both $\kappa_{\rm th}$ and $\bar{\rho}$ go down as pressure goes down). While not an insignificant variation, considering the uncertainties regarding this quantity in HJs atmospheres, the use of a constant thermal opacity is an acceptable approximation.