Hot Spot Offset Variability from Magnetohydrodynamical Thermoresistive Instability in Hot Jupiters

Following H23, we consider the atmospheric layer extending from a pressure of $1.0$ bar at the base to $0.01$ bar at the top, with gravitational acceleration $g_{p}=9.0\rm~{}m~{}s^{-2}$. We assume that the atmosphere is in hydrostatic balance and composed of an ideal gas of pure molecular hydrogen. As the modeled layer represents only 3% of the planetary radius, we adopt the plane parallel approximation where we map the the spherical coordinates $(\phi,\theta,r)$ (longitude, latitude and radius, respectively) onto Cartesian coordinates $(x,y,z)$. We assume that the flow velocity is in the longitudinal direction and depends only on height, i.e. $u_{x}(z)$. incompressibility ($\nabla\cdot{\bf u}=0$) then implies $u_{y}=u_{z}=0$. The magnetic field, satisfying $\nabla\cdot{\bf B}=0$, is

where $B_{0}$ is a background radial field and $B_{x}$ is the toroidal field induced by differential rotation. As discussed by H23, such a radial field component could arise from a misaligned dipole, higher order multipole, or a local dynamo (e.g. Rogers & McElwaine 2017; Dietrich et al. 2022). The simple magnetic field geometry described by Equation (1) can still support torsional Alfvén oscillations impacting zonal flows.

H23 assumed full axisymmetry. We relax this here by allowing the temperature $T$ to depend on $x$ (corresponding to a longitudinal variation in $T$). Note that $u_{x}$ and $B_{x}$ remain axisymmetric: they must be independent of $x$ in order for the velocity and magnetic field to remain divergence free. With these assumptions, the magnetohydrodynamics (MHD) equations (Davidson, 2001) reduce to

These are the same set of equations considered by H23 except for the advective term on the left hand side of Equation (4), which allows for the longitudinal advection of temperature. Since the layer is thin, we assume that other terms involving horizontal gradients, for example in the ohmic dissipation term, are negligible compared to vertical gradients. We use an overbar to indicate quantities that are not time-evolved (see Section 2.3 for further details).

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

models the effects of angular momentum transfer to the equator from higher latitudes, with the pressure $P$ measured in bars, and $\dot{v}$ setting the peak amplitude. Equation (3) is the $x$-component of the induction equation. The temperature-dependent MD is given by

taken from Menou (2012b) 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) considering sodium and potassium only (these elements give the dominant contribution to the ionization). The heat capacity is $c_{p}$, and the thermal conductivity is defined as

where $\sigma$ is the Stefan-Boltzmann constant and $\kappa_{\rm th}$ is the Rosseland mean opacity, operating mostly in the infrared regime, where thermal emissions takes place. We use the thermal conductivity to set the dynamic viscosity

where the Prandtl number $\mathrm{Pr}$ is assumed constant throughout the atmosphere.

In Equation (4) which gives the temperature evolution, we write the irradiation flux as

where $F_{s}(\phi)$ is the incoming flux from the host star at the surface of the planet set by the irradiation temperature, $\kappa_{v}$ is the visible opacity which we keep constant at $4.0~{}\times~{}10^{-4}~{}\rm m^{2}~{}kg^{-1}$ throughout all our simulations as in M12, and the $\sqrt{3}$ factor comes from the exponential in Equation (29) of Guillot (2010).

We now develop a low order expansion in the longitudinal direction, applying the approach¹¹1Tritton (1988) studied convection in a vertical torus heated at the bottom and cooled at the top. Interestingly, for that case an expansion of $T(\phi)$ in $\cos$ and $\sin$ terms leads to a set of equations equivalent to the famous Lorenz equations (Lorenz, 1963). of Tritton (1988). We expand the $\phi$-dependence of temperature and MD as

We also assume for simplicity that the angular variation of the irradiation flux is given by $F_{s}(\phi)=F_{s}(0)(1+\cos\phi)/2$, which has a maximum at the substellar point ($\phi=0$) and drops to zero at the anti-substellar point²²2We note that including a more realistic irradiation profile, $F_{s}(\phi)\propto\cos\phi$ for $|\phi|<\pi/2$ and $0$ for $\pi/2<|\phi|<\pi$ gives a similar low order expansion $F_{s}(\phi)=(F_{s}(0)/\pi)(1+(\pi/4)\cos\phi)$. ($\phi=\pi$). Inserting these expansions into Equation (4), with $\phi=x/R$, and identifying terms that are independent of $\phi$, proportional to $\cos\phi$, or proportional to $\sin\phi$, gives evolution equations for the amplitudes $T_{0}$, $T_{1}$, and $T_{2}$:

The original 2D problem given by Equations (2)–(4) has been reduced to the set of coupled 1D equations (2), (3), and (12)–(14). These can now be solved to find the time-evolution of $u_{x}(z,t)$, $B_{x}(z,t)$, and the temperature components $T_{0}(z,t)$, $T_{1}(z,t)$ and $T_{2}(z,t)$. One complication is that $\eta$ appears in the induction equation (Equation (3)), introducing a $\phi$-dependent term into that equation. To remain consistent with our assumption of axisymmetric $B_{x}$, we use the axisymmetric term $\eta_{0}$ when evaluating those terms in the induction equation.

To obtain the coefficients $\eta_{0}$, $\eta_{1}$, and $\eta_{2}$ at each time step, we first evaluate $\eta(\phi)$ using Equations (6) and (10). A least squares fit of the functional form of Equation (11) to $\eta(\phi)$ then gives

where the sums are over the grid of $n_{\phi}$ $\phi$-values at which $\eta(\phi)$ has been evaluated. Note that since $\eta$ is exponential in temperature, a sinusoidal expansion of $\eta(\phi)$ as given by Equation (11) is not necessarily a good approximation. We find that this approach reproduces $\eta(\phi)$ to within a factor of 6 at worst, but within ten percent during most of the simulation, which is reasonable considering that the amplitude of $\eta$ can change by orders of magnitude during an oscillation (see Appendix A for further details).

The boundary conditions for $u_{x}$ and $B_{x}$ are the same as specified in H23. We impose the inner layer of our domain to be corotating with the core of the planet, setting $u_{x}=0$, and we set the outermost layer to be free of any mechanical stresses, setting $\partial u_{x}/\partial z=0$. At the base of our domain, we impose magnetic stresses to vanish, setting $\partial B_{x}/\partial z=0$ while at the top we consider a magnetic vacuum boundary condition, setting $B_{x}=0$. We assume a constant axisymmetric heat flux $F_{0}$ coming from the hot planetary core, so that the bottom boundary condition for the axisymmetric part of the temperature $T_{0}(t)$ is ${\partial T_{0}}/{\partial z}={-F_{0}}/\bar{\chi}$. We use the flux $\sigma T_{\rm int}^{4}$, for an interior temperature of $T_{\rm int}=150\ \mathrm{K}$, as in M12 and H23 for $F_{0}$. The axisymmetric interior flux thus impose $T_{1}$ and $T_{2}$ components to have constant zero flux, setting $\partial T_{1}/\partial z=\partial T_{2}/\partial z=0$. At the outermost layer, we consider the thermal flux to adjust to the temperature. To do so we consider the thermal flux as

where the $4/3$ factor is introduced to reproduce the results presented by Guillot (2010). With this equation, we can isolate the temperature derivative, and use the longitudinal expansion procedure to get $\partial T_{0}/\partial z$, $\partial T_{1}/\partial z$ and $\partial T_{2}/\partial z$. The free parameters used for the temperature profile are the equilibrium temperature $T_{\rm eq}$ and the thermal opacity $\kappa_{th}$.

We ran a grid of models covering the ranges of parameter values shown in Table 1, using the same initial conditions for all simulations. For our initial conditions, we set the atmosphere to be in solid-body rotation, thus $u_{x}=0$ across all layers, and the background magnetic field is undisturbed, leaving only the radial field, thus $B_{x}=0$. We let the initial steady-state temperature profile obtained by an ODE solver relax in a time-dependent code, while also letting $P,\bar{\rho},\bar{\mu}$ and $\bar{\chi}$ adjust to the changing temperature, but these quantities are kept constant once the relaxation procedure completed. As the fluid is initially under solid rotation, $T_{2}=0$ across the atmosphere. We then use Equation (6) with the initial $\tilde{T}(\phi)$ profile to construct the initial $\tilde{\eta}(\phi)$.

The numerical solver is a straightforward generalization of the scheme presented in H23. The velocity, the magnetic field, and all components of the temperature are advanced in time exactly as in H23. We use 400 grid points equidistant in radius, and 360 grid points in $\phi$ to carry out the fit to $\eta(\phi)$. This numerical grid is fine enough to resolve boundary layers and longitudinal dependencies with good computational speed. We adopt a constant time step, chosen short enough to properly capture the bursts following the TRI. This requires $\Delta t=20\ \mathrm{s}$, with a typical simulation requiring three million time steps. The radial resolution needs to be chosen carefully, because under-resolving the system leads to the appearance of kinks in the profiles that artificially inject energy into the system through enhanced dissipation. These kinks appear once the instability is triggered, therefore they do not fundamentally change the behavior of the system, but they do artificially extend the decay phase of post-burst Alfvén waves. To avoid any numerical issue with our magnetic diffusivity, we have also imposed a maximum of $10^{12}\rm m^{2}~{}s^{-1}$ onto the MD, mostly notable near the beginning of the simulations, where the nightside is at its coldest.

We first show a representative simulation exhibiting periodic bursts of Alfvén oscillations triggered by the TRI. To best illustrate the features of a burst and subsequent decaying Alfvén waves, we have chosen the simulation with parameters $B_{0}=30$ G, $\dot{v}=0.006\rm~{}m~{}s^{-1}$, $\Pr=0.01$, $\kappa_{\rm th}=0.0008\rm~{}m^{2}~{}kg^{-1}$ and $T_{\rm eq}=1000$ K.

Figure 1 shows the time series of the longitudinal velocity and magnetic field component, as well as the temperature and magnetic Reynolds number (${\rm Rm}=u_{x}H/\eta_{0}$, with $H$ the local pressure scale height) at the center of the domain. In addition, we also show the longitudinal position of the maximum of $T(\phi)$ in the center of the domain (bottom panel). On the right-hand side of the figure, we zoom in to show the Alfvénic behavior of the velocity and magnetic field. The general behavior is similar to that found in H23. In panel (c) of Figure 1, we see that the hottest point at mid-depth reaches longitudinal displacement of $76^{\circ}$ during the build-up phase. However, once the TRI is triggered, the position of the hottest point varies greatly, even reaching westward positions for over a day at a time.

Figure 2 shows the ratios of the main timescales in the system during the simulation. Panel (a) indicates that thermal diffusion overwhelms advection almost throughout the simulation, with the exception of the Alfvénic oscillations at depth, where they share similar values. As thermal diffusion is almost always the leading heat transport mechanism to a varying degree throughout the simulation, finding the displacement of the temperature maximum is not simply a matter of time-integrating the velocity. As diffusion overwhelms advection, the hottest point does not reach large longitudes as heat diffuses quickly before it can be advected. In panel (b), comparing the advection and Alfvén timescales, we do not have a dominant timescale throughout. During build-up, advection dominates, but during and shortly after the TRI is triggered, the Alfvén timescale decreases to become slightly smaller than the advection time. The ratio of the Alfvén timescale to thermal diffusion timescale is not shown, as it is always $\gg 1$.

Figure 3 shows the depth-dependence of key variables as a function of time. Panels (a), (b) and (c) show the evolution of the velocity, magnetic field and the temperature. These are similar to the results from H23. We see strong heating of the atmosphere at depth during the oscillatory phase, with the TRI initially starting at a pressure of $\approx 0.2$ bars and propagating downwards towards the 1 bar level. Panel (f) shows the evolution of the magnetic Reynolds number. As the magnetic Reynolds number crosses unity, the TRI is triggered and as it dissipates through Alfvénic oscillations, the system returns to a quiescent build-up phase evolving towards a steady state, until the next TRI is triggered. The magnetic field strength during the oscillation phase can become large enough that magnetic pressure is significant in some regions; we discuss this further in Appendix B.

The longitudinal variations of temperature are summarized in panel (d), which shows the temperature difference between coldest and hottest spots at each depth, and panel (e), which shows the angular position $\phi_{T_{\mathrm{max}}}$ of the hottest point at each depth. The rapid vertical thermal diffusion prevents strong variations of $\phi_{T_{\mathrm{max}}}$ with depth, except at the highest pressures where the angular displacements lag those at higher altitudes. The largest offsets occur when the advection timescale is closest in size to the thermal diffusion timescale (Figure 2). Were advection to be the dominant heat transport mechanism, we would expect an isothermal ring to form around the planet; in the opposite limit of negligible advection, we would expect the hottest point to remain at the substellar point. As the simulation moves between these limits, the offset can go beyond the terminator when advection is strongest, but returns on the eastern dayside when diffusion dominates again. The largest temperature difference across the surface is seen in the decaying phase, where $\Delta T$ reaches $\approx 310\ \mathrm{K}$ at $P\approx 0.2$ bars, lasting for about 10 days after the burst of oscillations. This is much smaller than the change in the axisymmetric component of the temperature $T_{0}$ during the burst, so that the heating from the TRI is close to being axisymmetric.

Figure 4 shows the depth-longitude temperature maps during three distinct instants. Panel (a) shows the initial condition of the atmosphere. We see a temperature contrast of about 1000 K as the fluid is at rest. In panel (b), the advection has become important and the temperature is almost homogeneous in longitude. The hot spot is also advected eastward, with the smallest deviation at lower pressure, where the thermal diffusion is fastest, as shown in Figure 2. Panel (c) shows the large rise in temperature caused by the TRI and the shift in hot spot cause by the reversal of the wind. While at lower pressure, the hot spot is near substellar, while at depth the hottest region is well into the western side of the dayside of the planet.

Figure 5 (a) shows the time-evolution of the thermal flux coming out of the atmosphere for different values of $B_{0}$. We first evaluate the outwards flux at the surface $F_{\mathrm{out}}(\phi)=-\bar{\chi}\partial\tilde{T}(\phi)/\partial z|_{z=z_{\mathrm{top}}}$ as a function of $\phi$, and then integrate $\int_{-\pi/2}^{\pi/2}F_{\mathrm{out}}(\phi)\cos\phi\ d\phi$ to approximate observing the dayside of the planet (e.g. during secondary eclipse). The same procedure is used for the nightside. The kinetic energy from the forcing, which is ultimately transformed into heat through magnetic field induction followed by ohmic heating, enhances the thermal flux by factors of $2$–$3$ compared to the initial advection free value. The dotted lines representing the angle-averaged thermal flux of the nightside is indistinguishable from the dayside at the start of TRI, but as irradiation is smaller, this region can cool down faster.

Figure 5 (b) shows the longitudinal position of the thermal flux peak. The faster vertical thermal diffusion at low pressure means that the offset of the hot spot at the surface is reduced compared to the offset of the hottest point at depth (Figures 3 and 4). Nevertheless, the behavior in time remains similar: it grows toward an equilibrium displacement of $60$–$70^{\circ}$ eastwards during the quiescent phase, oscillates around $0^{\circ}$ after the onset of TRI, reaching in some cases more than $50^{\circ}$ westward. Once the flow decouples from the field and the oscillations damp down, the hot spot displacement then relaxes back toward equilibrium over about 10 days.

The set of parameters chosen for the simulation just discussed are only one of many choices that lead to TRI. With four different control parameters that can be varied, the behavior of the different solutions changes. Increasing the radial magnetic field strength reduces the recurrence period of the TRI, by reducing the velocity needed to trigger the instability, heating being faster with a stronger magnetic field. As such, the kinetic energy reservoir is smaller when the TRI triggers for larger $B_{0}$ values, thus the thermal flux outgoing from the atmosphere is reduced, as seen in Figure 5 (a). The hot spot offset gets smaller with stronger $B_{0}$ values, as seen in Figure 5 (b), reflecting the quenching action of the Lorentz force on the zonal flow. However, if the field is too strong, then the system simply reaches steady-state, as shown in H23. Increasing $B_{0}$ also reduces the number of Alfvén oscillations as the field is stiffer. The acceleration $\dot{v}$, when increased, also reduces the recurrence time as more energy is injected into the system per unit time. This energy is converted into heat and the critical temperature where ${\rm Rm}=1$ is reached faster. However, if the acceleration is too large, then the system will remain too hot and recurrent TRI will not be possible; whereas if it is too small, then the critical temperature will never be reached. The equilibrium temperature dictates the initial magnetic Reynolds number value. Considering only equilibrium values compatible with TRI, a hotter atmosphere will reach TRI faster than a colder one with the same parameters, therefore reducing the recurrence period. As for the thermal opacity $\kappa_{\rm th}$, a larger value will also help the atmosphere to reach TRI faster as heat is dissipated more slowly, therefore also reducing the recurrence period. The quantitative effect of $B_{0}$ and $T_{\mathrm{eq}}$ can be seen in Figure 10 of H23.

Although we find very similar behavior to H23, for a given set of parameters the precise evolution is different. As discussed in section 2, we use the longitudinally-averaged MD $\eta_{0}$ in the induction equation, which is larger than in the fully axisymmetric model of H23, where only substellar conditions were considered. As we are also considering the colder nightside, $\eta_{0}$ is larger than $\eta$ at the substellar point. Thus, the system needs to reach larger velocities to generate large enough values of $B_{x}$ to compensate for the larger diffusivity, which reduces the value of the magnetic Reynolds number, making it harder for the system to reach the critical value of ${\rm Rm}=1$, triggering the instability. Only if the system is able to reduce its temperature after the TRI through an outgoing thermal flux at the upper boundary, thus going back to Rm smaller than unity, can the system generate recurrent bursts. Overall, the larger MD shifts the instability region to larger values of $T_{\rm eq}$, $\kappa_{\rm th}$ and $\dot{v}$.