# Variability from thermo-resistive instability in the atmospheres of hot
jupiters

Raphaël Hardy Département de Physique, Université de Montréal, Montréal, QC,
H3C 3J7, Canada
Department of Physics and McGill Space Institute, McGill University, Montréal,
QC, H3A 2T8, Canada
Institut de Recherche sur les Exoplanètes (iREx), Université de Montréal,
Montréal, QC H3C 3J7, Canada  Andrew Cumming Department of Physics and McGill
Space Institute, McGill University, Montréal, QC, H3A 2T8, Canada
Institut de Recherche sur les Exoplanètes (iREx), Université de Montréal,
Montréal, QC H3C 3J7, Canada  Paul Charbonneau Département de Physique,
Université de Montréal, Montréal, QC, H3C 3J7, Canada
<EMAIL_ADDRESS>
###### Abstract

The atmosphere of a hot jupiter may be subject to a thermo-resistive
instability, in which the increasing electrical conductivity with temperature
leads to runaway Ohmic heating. We introduce a simplified model of the local
dynamics in the equatorial region of a hot jupiter that incorporates the back
reaction on the atmospheric flow as the increasing electrical conductivity
leads to flux freezing, which in turn quenches the flow and therefore the
Ohmic heating. We demonstrate a new time-dependent solution that emerges for a
temperature-dependent electrical conductivity (whereas a temperature-
independent conductivity always evolves to a steady-state). The periodic cycle
consists of bursts of Alfven oscillations separated by quiescent intervals,
with the magnetic Reynolds number alternating between values smaller than and
larger than unity, maintaining the oscillation. We investigate the regions of
pressure and temperature in which the instability operates. For the typical
equatorial accelerations seen in atmospheric models, we find instability at
pressures $\sim 0.1$–$1\ {\rm bar}$ and temperatures $\approx 1300$–$1800\
{\rm K}$ for magnetic fields $\sim 10\ {\rm G}$. Unlike previous studies based
on a constant wind velocity, we find that the instability is stronger for
weaker magnetic fields. Our results add support to the idea that variability
should be a feature of magnetized hot jupiter atmospheres, particularly at
intermediate temperatures. The temperature-dependence of the electrical
conductivity is an important ingredient that should be included in MHD models
of hot jupiter atmospheric dynamics.

Astrophysical fluid dynamics – Magnetohydrodynamics – Exoplanet atmospheric
variability – Exoplanet atmospheric dynamics

##  1 Introduction

Extreme irradiation from their host stars and tidally locked orbits make for
interesting dynamics in hot jupiters’ (HJs) atmospheres. The temperature
gradient creates a dayside to nightside flow (Komacek & Showman, 2016),
contributing to the driving of equatorial superrotation (Showman & Polvani,
2011; Read & Lebonnois, 2018; Imamura et al., 2020). These strong zonal flows
displace the hottest part of the surface away from the substellar point. The
equatorial superrotation usually leads to eastwards displacement of the hot
spot (Showman & Guillot, 2002; Cooper & Showman, 2005; Showman et al., 2009;
Rauscher & Menou, 2010; Kataria et al., 2016). However, in some HJs, westward
displacements and variability have been observed (Armstrong et al., 2016; Dang
et al., 2018; Jackson et al., 2019; Bell et al., 2019; von Essen et al.,
2020).

One explanation for the westward displacements of the hot spot and variability
involves magnetohydrodynamical (MHD) effects that reverse the direction of the
atmospheric flow (Rogers & Komacek, 2014; Rogers, 2017; Hindle et al., 2019,
2021a, 2021b). Magnetic effects are expected to be important in HJs
atmospheres (Rogers & Komacek, 2014) because the high temperatures lead to
partial ionization of the atmosphere, in particular of alkali metals (Batygin
& Stevenson, 2010; Perna et al., 2010a). Magnetic fields have two effects:
they can lead to dissipation through Ohmic losses, and they can change the
dynamics through $\vec{J}\times\vec{B}$ forces. Simulations have tended to
focus on magnetic drag and the associated dissipation, adding a drag term to
the momentum equation (Perna et al., 2010a) that can significantly slow the
atmospheric flow (e.g. Beltz et al. 2022), as well as heating the atmosphere
(e.g. Batygin & Stevenson 2010). This is appropriate at low magnetic Reynolds
number (lower temperatures), where flux freezing does not hold and Ohmic
dissipation dominates advection. Rogers & Komacek (2014) included the full
dynamics by solving the MHD equations directly, and found that at higher
magnetic Reynolds numbers (temperatures), magnetic forces can lead to time-
variability, including oscillatory flows, and flow reversals (westward flow).
Rogers (2017) and Hindle et al. (2021a) used these results to derive limits on
the magnetic field strength in several HJs.

In this paper, we demonstrate an additional source of variability coming from
the temperature-dependence of the electrical conductivity, which leads to a
temporal variation of the electrical conductivity as the temperature changes
over time. Menou (2012a) pointed out that the strong increase of the
electrical conductivity with temperature allows the development of a thermo-
resistive instability (TRI). A small increase in temperature leads to
increased electrical conductivity, larger induced currents and therefore
enhanced Ohmic heating that further increases the temperature, thus yielding a
positive feedback. Describing the flow dynamics with a parametrized magnetic
drag law, Menou (2012a) found local instability at pressures $\approx 3$–$300\
{\rm mbar}$, magnetic field strengths $B\sim 10\ {\rm G}$, and temperatures
for which the magnetic drag is weak enough to allow $\approx 10\ {\rm
km~{}s^{-1}}$ zonal wind speeds. Global simulations that include magnetic
effects have generally not included a time-dependent electrical conductivity.
Rauscher & Menou (2013) and Beltz et al. (2022) did include a temperature-
dependent magnetic drag force, but none of the studies that solve the full MHD
equations included a time-dependent conductivity.

In this paper, we introduce a simplified model to explore the interplay
between the dynamical effect of the magnetic field and the changes in
electrical conductivity over time as temperature evolves. We find a dynamical
version of the local TRI identified by Menou (2012a) that leads to time-
variability in the flow, either in the form of Alfvén oscillations or short
bursts separated by longer quiescent intervals. The key ingredient is the
temperature-dependent conductivity, since it allows the magnetic Reynolds
number to cycle between small and large values during the oscillation. We then
discuss the application of our results to HJ atmospheric conditions and the
implications for HJ atmospheric dynamics.

##  2 Local model with a time-dependent magnetic diffusivity

To model the local dynamics in the equatorial plane, we consider the evolution
of a zonal flow in the presence of a background radial magnetic field $B_{r}$.
Rather than assume a fixed velocity for the zonal flow, since we are
interested in the dynamics, we assume that it is forced by a net flow of
momentum into the equatorial region, represented by an acceleration $\dot{v}$.
This could correspond to either an eastward or westward acceleration depending
on the overall global dynamics (Showman & Polvani, 2011; Hindle et al.,
2021b).

Working under the MHD approximation (Davidson, 2001), we solve the zonal
components of the axisymmetric MHD equations in the equatorial plane of the
planet. Thus, all radial and latitudinal flow displacements are eliminated
from the system. For the background magnetic field, we consider a purely
radial magnetic field111For the purposes of our simplified model, we consider
a background radial field only. A global aligned dipole field would have only
horizontal components at the equator, but radial field components would be
expected if the dipole is misaligned, or higher order field components are
present. respecting ${\bf\nabla}\cdot{\bf B}=0$. In this geometry, radial
shear drives torsional Alfvén oscillations.

As the thickness of the atmospheres of HJs are very small compared to their
radii, we can use the plane parallel approximation to write our first set of
equations. With these approximations in mind, the MHD equations involving the
local values of the zonal flow of the fluid in the reference frame of the
rotating planet, $U_{\phi}$, toroidal magnetic field, $B_{\phi}$, and
temperature, $T$, are:

$\frac{\partial U_{\phi}}{\partial t}=\frac{B_{r}}{4\pi\rho}\frac{\partial
B_{\phi}}{\partial r}+\dot{v},$ (1) $\frac{\partial B_{\phi}}{\partial
t}=B_{r}\frac{\partial U_{\phi}}{\partial
r}+\eta\frac{\partial^{2}B_{\phi}}{dr^{2}},$ (2) $\rho c_{p}\frac{\partial
T}{\partial t}=-\rho
c_{p}\frac{(T-T_{0})}{\tau}+\frac{\eta}{4\pi}\left(\frac{\partial
B_{\phi}}{\partial r}\right)^{2},$ (3)

where $\eta$ is the magnetic diffusivity, $\rho$ the gas density and $c_{P}$
the heat capacity. In the heat equation (3), the first term on the right-hand
side represents Newtonian cooling to a reference temperature $T_{0}$ on a
timescale $\tau$, and the second term is Ohmic heating. We take the constant
acceleration $\dot{v}$ to be positive (eastward), but our results apply
equally to the westward $\dot{v}<0$ case with the direction of flow reversed.

We now assume that radial variations are on the scale of the thickness of the
layer $H$, writing the radial derivative of a quantity $f$ as $\partial
f/\partial r\sim f/H$. We can then write down a simplified version of the MHD
equations:

$\frac{dU_{\phi}}{dt}=-\frac{B_{r}}{4\pi\rho}\frac{B_{\phi}}{H}+\dot{v}=-\frac{v_{A}^{2}}{H}\frac{B_{\phi}}{B_{r}}+\dot{v},$
(4) $\frac{dB_{\phi}}{dt}=B_{r}\frac{U_{\phi}}{H}-\eta\frac{B_{\phi}}{H^{2}},$
(5) $\rho c_{p}\frac{dT}{dt}=-\rho
c_{p}\frac{(T-T_{0})}{\tau}+\frac{\eta}{4\pi}\left(\frac{B_{\phi}}{H}\right)^{2},$
(6)

where the Alfvén speed $v_{A}$ is given by $v_{A}^{2}=B_{r}^{2}/4\pi\rho$.

To capture the temperature-dependence of the magnetic diffusivity $\eta$, we
write

$\eta\left(T\right)=\eta_{0}\exp\left[E_{\eta}\left(\frac{T_{0}}{T}-1\right)\right],$
(7)

where $\eta_{0}=\eta(T_{0})$ and the parameter $E_{\eta}$ controls the
temperature sensitivity ($d\ln\eta/d\ln T=-E_{\eta}T_{0}/T$). The exponential
form for $\eta(T)$ is similar to the diffusivity from alkali metal ionization,
which involves the thermal ionization fraction (e.g. Menou 2012b).

We write the evolution equations in dimensionless form by defining
$v=U_{\phi}/v_{A}$, $b=B_{\phi}/B_{r}$, $\dot{V}=\dot{v}H^{2}/\eta_{0}v_{A}$,
the magnetic Reynolds number ${\rm Rm}=v_{A}H/\eta_{0}$, and the Eckert number
${\rm Ec}=v_{A}^{2}/c_{p}T_{0}$. We measure time in units of the Ohmic
diffusion time $H^{2}/\eta_{0}$, and temperature with respect to the reference
temperature $T_{0}$, setting $\eta_{0}t/H^{2}\rightarrow t$,
$\eta_{0}\tau/H^{2}\rightarrow\tau$, and $T/T_{0}\rightarrow T$. In these
units, the Alfvén timescale $H/v_{A}$ is $1/{\rm Rm}$. If we estimate the
cooling timescale as $H^{2}/\kappa_{T}$, with $\kappa_{T}$ the thermal
diffusivity, then the parameter $\tau\approx\eta_{0}/\kappa_{T}$. The
dimensionless form of equations (4)–(6) are then:

$\frac{dv}{dt}=-{\rm Rm}~{}b+\dot{V},$ (8) $\frac{db}{dt}={\rm
Rm}~{}v-\tilde{\eta}{b},$ (9) $\frac{dT}{dt}=-\frac{(T-1)}{\tau}+{\rm
Ec}~{}\tilde{\eta}b^{2},$ (10)

where $\tilde{\eta}(T)~{}=~{}\eta(T)/\eta_{0}$ is given by equation (7). These
equations allow us to solve for the evolution of $b$, $v$ and $T$ given the 5
parameters ${\rm Rm}$, ${\rm Ec}$, $E_{\eta}$, $\tau$, and $\dot{V}$.

We solve the nonlinear coupled ordinary differential equation system defined
by equations (7) to (10) as an initial value problem. Towards this end we use
a solver from SciPy222The 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.

In the absence of forcing ($\dot{V}=0$), equations (8) and (9) describe damped
Alfvén waves with (dimensionless) frequency given by $\omega^{2}={\rm Rm}^{2}$
and damping timescale $1/\tilde{\eta}$. Note also that with $\dot{V}=0$, the
limit of magnetic drag is recovered by setting $db/dt=0$, in which case
$b={\rm Rm}v/\tilde{\eta}$ and the momentum equation reduces to $dv/dt=-({\rm
Rm}/\tilde{\eta})^{2}v$, representing a drag force with associated timescale
$(\tilde{\eta}/{\rm Rm})^{2}$, in agreement with equation (1) of Rauscher &
Menou (2013).

With forcing included ($\dot{V}>0$), equations (8)–(10) have a steady state
solution

$b=\frac{\dot{V}}{{\rm Rm}};\hskip
28.45274ptv=\tilde{\eta}~{}\frac{\dot{V}}{{\rm Rm}^{2}},$ (11)

where $\tilde{\eta}$ is evaluated at the steady-state temperature for which
$dT/dt=0$ (i.e. the value of $T$ for which the right hand side of equation
(10) vanishes, which yields $\bm{T-1=}{\rm\bf
Ec}\tilde{\bm{\eta}}\bm{\tau}\bm{b}^{2}$). The physical balance represented by
equation (11) results from the driven zonal flow winding the radial field
until the Ohmic dissipation associated with the induced toroidal magnetic
component balances the rate of kinetic energy input to the layer.

For a temperature-independent diffusivity, the evolution of the driven system
is always to the steady-state given by equation (11), although the nature of
the evolution depends on the magnetic Reynolds number. At high ${\rm Rm}$ the
system moves towards the steady-state via a series of damped Alfvén
oscillations, whereas at low ${\rm Rm}$ the system evolves smoothly to the
steady state. The situation is very different when the magnetic diffusivity
depends strongly on temperature, as we discuss in the next section.

##  3 Dynamical Thermo-Resistive Instability

We now show that equations (8)–(10) have unstable time-dependent solutions
when the magnetic diffusivity depends strongly on temperature. We first
present some representative evolutions and then discuss the regions of the
parameter space where the system is unstable.

###  3.1 Examples of time-dependent solutions

Figure 1 shows three numerical integrations of equations (8)–(10) for
different choices of parameters. In each case, we start the system from rest
($v=0$), with no toroidal magnetic field ($b=0$), and at the reference
temperature ($T=1$). The top panel of Figure 1 shows an example in which the
magnetic diffusivity does not depend on temperature ($E_{\eta}=0$). In this
case, with low magnetic Reynolds number (${\rm Rm}=0.1$), the system evolves
smoothly to the steady state given by equation (11).

The middle panel of Figure 1 shows the dramatically different behaviour that
occurs when the magnetic diffusivity is temperature-sensitive. This
integration has the same parameters as the top panel except that we now set
$E_{\eta}=10$. Instead of a steady-state solution, we find sustained Alfvénic
oscillations. While not a pure Alfvén wave due to the temperature-dependent
Ohmic dissipation term, the solution shows the signature $\pi/2$ phase lag
between $b$ and $v$ expected for Alfvén waves.

Figure 1: Time evolution of velocity $v$, toroidal field $b$, temperature $T$
and effective magnetic Reynolds number ${\rm Rm}_{\rm eff}={\rm
Rm}|v|/\tilde{\eta}$, showing the three different outcomes: (a) the system
reaches steady-state (${\rm Rm}$=0.1, $\tau$=1, ${\rm Ec}$=1, $E_{\eta}$=0,
$\dot{V}$=0.1); (b) oscillations (${\rm Rm}$=0.1, $\tau$=1, ${\rm Ec}$=1,
$E_{\eta}$=10, $\dot{V}$=0.1) and (c) bursts (${\rm Rm}$=0.01, $\tau$=10,
${\rm Ec}$=1, $E_{\eta}$=10, $\dot{V}$=0.01).

It is important to note that during the oscillation, the effective Reynolds
number ${\rm Rm}_{\rm eff}={\rm Rm}|v|/\tilde{\eta}$ reaches values above $1$,
even though the input parameter ${\rm Rm}<1$ (see red curve in Figure 1). In
other words, the temperature increase is enough to reduce the magnetic
diffusivity to a level where ${\rm Rm}_{\rm eff}$ becomes larger than unity.
This temporary crossing of the effective magnetic Reynolds number from ${\rm
Rm}_{\rm eff}<1$ to ${\rm Rm}_{\rm eff}>1$ seems to be a requirement for
sustained oscillations to occur. We find that oscillations do not occur for
other choices of parameters where the system either never reaches ${\rm
Rm}_{\rm eff}$ values above unity, or has ${\rm Rm}_{\rm eff}$ values that
stay above unity. In these cases, the system either evolves smoothly to steady
state (${\rm Rm}_{\rm eff}<1$) or shows damped oscillations to a steady-state
(${\rm Rm}_{\rm eff}>1$), similar to the situation with $E_{\eta}=0$.

Figure 2: Angular frequency $\omega/\omega_{A}$ of the oscillation cycle as a
function of $\dot{V}$ and $\tau$ for (a) ${\rm Rm}=0.1$, ${\rm Ec}=1.0$ and
(b) ${\rm Rm}=0.01$, ${\rm Ec}=1.0$ and as a function of ${\rm Ec}$ and $\tau$
for (c) ${\rm Rm}=0.01$, $\dot{V}=0.1$ and (d) ${\rm Rm}=0.01$,
$\dot{V}=0.01$. The magnetic diffusivity sensitivity is kept constant for all
panels at $E_{\eta}=10.0$. The black lines represent the isocontours of the
angular frequency $\omega/\omega_{A}$. The solid red lines represent the
instability criterion of equation (15); the system should be oscillating above
it. The dashed red lines represent the analytic estimate of equation (14) for
the boundary between bursts and Alfvénic oscillations.

The bottom panel of Figure 1 shows a different behaviour in which the time-
dependence takes the form of regular short-duration bursts. This case is
intermediate between the solutions with Alfvénic oscillations and those that
evolve smoothly to a steady-state. Compared to the previous example, this case
has a longer cooling timescale, slower forcing, and lower ${\rm Rm}$. The long
phase between bursts, lasting for more than 10 Alfvén times, can be understood
as a quasi-steady-state in which the system slowly evolves along a sequence of
steady-state solutions satisfying equation (11) with increasing temperature as
the system slowly heats up. Eventually, the temperature becomes large enough
that the TRI occurs, driving the system to ${\rm Rm}_{\rm eff}>1$. The
magnetic field then undergoes a rapid growth as it is now coupled to the flow,
the resulting large magnetic torque then leads to a rapid drop in the
velocity, an Alfvén oscillation follows, and eventually the cycle restarts. In
the example shown, the system oscillates for only half a period after the
onset of the burst, but in other cases multiple damped oscillations can occur.
Depending on the parameters, we find that the interval between bursts varies
smoothly from large to smaller values comparable to the Alfvén period, at
which point the bursts smoothly transition to Alfvénic oscillations.

###  3.2 Domains of instability and analytic estimates

Panels (a)-(b) of Figure 2 show the frequency of the oscillation cycle as a
function of the forcing strength $\dot{V}$ and cooling time $\tau$, for two
different values of ${\rm Rm}$. The other parameters are fixed at
$E_{\eta}=10$ and ${\rm Ec}=1$. The instability occurs in the region with
$\omega>0$ (green shading), at small $\tau$ and large $\dot{V}$. In the region
where $\omega=0$ (pale yellow), the solution goes to steady-state and
sustained oscillations do not occur.

Since the sustained oscillations require the effective magnetic Reynolds
number to cross above and below one, we require ${\rm Rm}<1$ for instability.
Indeed, in Figure 2, we see that the instability region shrinks in size as
${\rm Rm}$ increases. It is noteworthy that we see an extended region of
bursts (with $\omega<\omega_{A}$) only for the lower ${\rm Rm}$ case; at
higher ${\rm Rm}$, the frequencies are close to $\omega_{A}$ throughout the
unstable region. As a different slice of parameter space, panels (c)-(d) of
Figure 2 show the frequency in the ${\rm Ec}$–$\tau$ plane, this time with
fixed ${\rm Rm}=0.01$ and $\dot{V}=0.1$. This shows that there exists a
considerable region parameter space within which we find sustained
oscillations, even when ${\rm Ec}$ takes values $\ll 1$.

We can understand the location of the unstable region with simple analytic
arguments. First, in the long build-up phase of the oscillation, ${\rm Rm}\ll
1$ and magnetic drag is weak so that the velocity increases linearly with
time, $v\approx\dot{V}t$. Assuming steady-state in equation (9), the
corresponding magnetic field is $b\approx{\rm Rm}\dot{V}t/\tilde{\eta}$. We
can thus rewrite equation (10) as

$\frac{dT}{dt}=-\frac{(T-1)}{\tau}+\frac{{\rm Ec}{\rm
Rm}^{2}(\dot{V}t)^{2}}{\tilde{\eta}}$ (12)

during the build-up phase. Direct numerical integration of equation (12) shows
good agreement with our integrations of the full equations until the moment
when the TRI begins.

The recurrence timescale is set by the time taken to reach the TRI. Making a
linear perturbation $T\rightarrow T+\delta T$ in equation (12) gives an
expression for the recurrence time as the time at which $d\delta T/dt$
transitions from being negative to positive (i.e. when $d\delta T/dt=0$),

$\tau_{\mathrm{rec}}^{2}=\frac{\tilde{\eta}T^{2}}{\dot{V}^{2}\tau{\rm
Rm}^{2}{\rm Ec}E_{\eta}}.$ (13)

Comparing this timescale to the Alfvén timescale, $\tau_{A}=1/{\rm Rm}$, the
transition from bursts to Alfvénic oscillations occurs when

$\dot{V}^{2}\tau=4\frac{\tilde{\eta}T^{2}}{{\rm Ec}E_{\eta}},$ (14)

plotted as the dashed line in Figure 2 (we have assumed small temperature
changes, $T\approx 1$ and $\tilde{\eta}\approx 1$, when plotting the dashed
line and have inserted a prefactor of 4 to better agree with the numerical
results in this limit).

Figure 2 shows that the system becomes stable at low values of $\dot{V}$. This
instability boundary corresponds well with the location where the recurrence
time $\tau_{\mathrm{rec}}$ starts to become longer than the time the system
takes to reach the steady-state described by equation (11), $\tau_{ss}$. Since
$v\approx\dot{V}t$ in the build-up phase, we find that
$\tau_{ss}\approx\tilde{\eta}/{\rm Rm}^{2}$. Requiring
$\tau_{ss}>\tau_{\mathrm{rec}}$ gives the criterion for instability as

$\dot{V}^{2}\tau>\frac{{\rm Rm}^{2}T^{2}}{\tilde{\eta}E_{\eta}{\rm Ec}}.$ (15)

We find that in our numerical solutions, the factor $T^{2}/\tilde{\eta}\approx
4$. The solid line in Figure 2 shows this relation assuming that
$T^{2}/\tilde{\eta}=4$, and gives a good match to the instability boundary.

As well as stabilizing at low $\dot{V}$, we also find that the oscillations
stabilize at large $\tau$. Very approximately, values of the cooling time
$\tau\gtrsim\tau_{A}=1/{\rm Rm}$ stabilize the oscillations, since the
temperature then does not change sufficiently during the oscillation period to
enable the effective magnetic Reynolds number to alternate between values
larger and smaller than unity.

The value of $E_{\eta}$ is kept constant throughout the panels of Figure 2.
The criteria presented in equations (14) and (15) show that the effects of
$E_{\eta}$ and ${\rm Ec}$ on the system are degenerate as they always appear
in pairs in the equations, and so our results can be rescaled to a different
value of $E_{\eta}$ as needed.

Finally, for comparison we also derive the instability criterion from Menou
(2012a) in our notation. The equivalent to our equation (12) for the
instability of Menou (2012a) is

$\frac{dT}{dt}=-\frac{(T-1)}{\tau}+\frac{{\rm Ec}{\rm
Rm}^{2}V^{2}}{\tilde{\eta}}$ (16)

where we have taken the weak drag limit with the wind velocity $V$ specified
as a parameter and $b={\rm Rm}\,V/\tilde{\eta}$. In this case, the instability
criterion ($d\delta T/dt>0$) is

$V^{2}>\frac{\tilde{\eta}T^{2}}{{\rm Rm}^{2}E_{\eta}{\rm Ec}\,\tau}.$ (17)

This is equivalent to equation (15) with $V=\dot{V}\tilde{\eta}/{\rm Rm}^{2}$,
i.e. the steady-state velocity under forcing $\dot{V}$.

Figure 3: Values of ${\rm Rm}$, $E_{\eta}$, ${\rm Ec}$, and $\tau$ shown as a
function of temperatures and pressures relevant for HJ atmospheres. The
magnetic Reynolds number is taken from equation (19). The temperature
sensitivity of the magnetic diffusivity is defined as $E_{\eta}\equiv-
T_{3}d\ln\eta/d\ln T$. The Eckert number is taken from equation (18). The
radiative cooling time is taken from Table 2 of Showman et al. (2008) with a
linear interpolation at higher temperatures and was turned into dimensionless
form with the Ohmic time $t_{\rm Ohm}$ from equation (20). We assume a surface
gravity equal to that of Jupiter and a magnetic field of $B=10\ {\rm G}$ (note
that ${\rm Rm}\propto B$, ${\rm Ec}\propto B^{2}$; $E_{\eta}$ and $\tau$ are
independent of $B$). The black lines represent the isocontours of the
presented parameters in each panel.

##  4 Application to hot jupiters

###  4.1 Parameter values

In the previous section we showed that instability occurs only for particular
ranges of the parameters ${\rm Ec}$, $E_{\eta}$, ${\rm Rm}$, $\tau$, and
$\dot{V}$. We now estimate the expected values of these parameters in the
atmosphere of a HJ. Values of ${\rm Ec}$, $E_{\eta}$, ${\rm Rm}$, and $\tau$
are shown in Figure 3 as a function of pressure and temperature.

Figure 3 shows that the temperature sensitivity of the magnetic diffusivity is
indeed large, and lies within a narrow range of values. We define
$E_{\eta}\equiv T_{3}d\ln\eta/d\ln T$, where $T_{3}=T/1000\ {\rm K}$, and
calculate the magnetic diffusivity $\eta=230\sqrt{T/{\rm K}}/x_{e}\ {\rm
cm^{2}\ s^{-1}}$ including contributions to the ionization fraction from
species up to Ni (following Menou 2012a). We find $E_{\eta}$ takes values
between $\approx 18$ and $29$ for $10^{-3}<P<10\ {\rm bar}$ and $1000<T<2500\
{\rm K}$. Indeed for $T\lesssim 1800\ {\rm K}$, where the ionization fraction
is dominated by potassium, $x_{e}\propto T^{3/4}\exp(-25.2/T_{3})$ (e.g. Perna
et al. 2010b equation (1)), giving $E_{\eta}\approx 25$ almost independent of
pressure.

The Eckert number ${\rm Ec}=v_{A}^{2}/c_{P}T$ is an important parameter for
the instability because it determines whether Ohmic dissipation is able to
heat the gas sufficiently quickly. The magnetic field strength of HJs is
uncertain. We scale our results to $B=10\ {\rm G}$ in this section, although
note that field strengths exceeding $\sim 100\ {\rm G}$ may arise in HJ
atmospheres (Reiners & Christensen, 2010; Cauley et al., 2019; Yadav &
Thorngren, 2017). The Alfvén speed is $v_{A}=B/(4\pi\rho)^{1/2}\approx 570\
{\rm cm\ s^{-1}}\ (B/10\ {\rm G})\,T_{3}^{1/2}(P/1\,{\rm bar})^{-1/2}$. For
the heat capacity, we assume pure molecular hydrogen with an ideal gas
equation of state, giving $\rho=2m_{p}P/k_{B}T\approx 2.42\times 10^{-5}\ {\rm
erg\ cm^{-3}}\ T_{3}^{-1}\,(P/1\,{\rm bar})$ and
$c_{P}\approx(7/2)k_{B}/2m_{p}\approx 1.45\times 10^{8}\ {\rm erg\ K^{-1}\
g^{-1}}$. Therefore ${\rm Ec}\approx B^{2}/14\pi P$, or

${\rm Ec}\approx 2.3\times 10^{-6}\ \left(\frac{B}{10\rm
G}\right)^{2}\left(\frac{P}{1\rm bar}\right)^{-1}.$ (18)

We see in Figure 3 that ${\rm Ec}$ takes larger values at low pressure,
suggesting that instability is more likely there, as also emphasized by Menou
(2012b) for the thermo-resistive instability.

The magnetic Reynolds number ${\rm Rm}=v_{A}H/\eta$ can be written in terms of
the ionization fraction as

${\rm Rm}\approx 0.1\,\left(\frac{x_{e}}{10^{-7}}\right)\
T_{3}\left(\frac{B}{10\rm G}\right)\,\left(\frac{P}{1\rm bar}\right)^{-1/2},$
(19)

where we take the length scale $H$ to be the pressure scale height $H\approx
k_{B}T/2m_{p}g\approx 1.6\times 10^{7}\ {\rm cm}\ T_{3}\,(g_{J}/g)$,
normalized to the surface gravity of Jupiter $g_{J}=GM_{J}/R_{J}^{2}$. Figure
3 shows that ${\rm Rm}$ can take on a large range of values since $x_{e}$ is
so sensitive to the local temperature. For example at $T\approx 1200\ {\rm
K}$, $x_{e}\sim 10^{-12}$–$10^{-10}$ for pressures in the range $10\ {\rm
bar}$–$1\ {\rm mbar}$, while for $T\approx 2000\ {\rm K}$, $x_{e}\sim
10^{-8}$–$10^{-7}$ in the same pressure range.

For the radiative cooling timescale, we follow Menou (2012b) and refer to the
tabulated values in Showman et al. (2008) (Table 2 in that paper). For
$T\approx 1200\ {\rm K}$, the approximate cooling timescales given there are
$\sim 3000\ {\rm s}$ ($\sim 10^{6}\ {\rm s}$) for a pressure of $1\ {\rm
mbar}$ ($1\ {\rm bar}$). For $T\approx 2000\ {\rm K}$, these timescales become
$\sim 300\ {\rm s}$ ($\sim 10^{5}\ {\rm s}$). Our dimensionless parameter
$\tau$ is the ratio of the cooling timescale to the Ohmic diffusion time
$t_{\rm Ohm}\approx H^{2}/\eta$, which we define using the pressure scale
height as the lengthscale,

$t_{\rm Ohm}\approx 3.5\times 10^{3}\ {\rm s}\
\left(\frac{x_{e}}{10^{-7}}\right)\
T_{3}^{3/2}\left(\frac{g}{g_{J}}\right)^{-2}.$ (20)

Proportional to $x_{e}$, the Ohmic time spans a large range of values, ranging
from $\sim 10^{-2}\ {\rm s}$ to hours depending on the local physical
conditions. Normalizing the radiative cooling time by the Ohmic time, we
obtain $\tau$, shown in the last panel of Figure 3.

Perhaps the most uncertain parameter is the forcing $\dot{V}$ since this
depends on the overall flow dynamics (and indeed the back reaction of the
dynamics in the equatorial plane on the global flow). The units of $\dot{V}$
are $v_{A}/t_{\rm Ohm}$ which is

$\displaystyle\frac{v_{A}}{t_{\rm Ohm}}$ $\displaystyle\approx$ $\displaystyle
0.16\ {\rm cm\ s^{-2}}\ \left(\frac{B}{10\rm
G}\right)\,T_{3}^{-1}\left(\frac{P}{1\rm bar}\right)^{-1/2}$ (21)
$\displaystyle\times\left(\frac{x_{e}}{10^{-7}}\right)^{-1}\left(\frac{g}{g_{J}}\right)^{2},$

i.e., the ratio of the acceleration in cgs units to equation (21) gives our
parameter $\dot{V}$. Values of acceleration can be estimated from Figure 9 of
Rogers & Komacek (2014), where the mean accelerations as a function of
latitude for some of their models are shown. Taking a value of $5\times
10^{-8}\ {\rm g\ cm^{-2}\ s^{-2}}$ for the volumetric force from their Figure
9, and the corresponding pressure and temperature of $\approx 10\ {\rm mbar}$
and $\approx 1300\ {\rm K}$ to calculate the density, we obtain an
acceleration $\approx 0.3\ {\rm cm\ s^{-2}}$ ($\approx 0.3\ {\rm km\ s^{-1}\
day^{-1}}$). For this value of acceleration, and taking the ionization
fractions above, we find $\dot{V}\sim 10^{-6}$–$10^{-5}$ at $T\approx 1200\
{\rm K}$, and $\dot{V}\sim 0.01$ at $T\approx 2000\ {\rm K}$. Figure 2 of
Rogers (2017) gives about a factor of ten smaller acceleration, showing
velocity variations of $\sim 10^{4}\ {\rm cm\ s^{-2}}$ on timescales of $\sim
10^{6}\ {\rm s}$, suggesting accelerations of $0.01\ {\rm cm\ s^{-2}}$ at the
mbar pressure level.

###  4.2 Unstable regions as a function of pressure and temperature

We present in Figure 4 the oscillation frequency and instability regions for
different choices of $B$ and $\dot{v}$. The red region has a magnetic field
strength of $B=40\rm~{}G$ and an acceleration of $\dot{v}=1\rm~{}cm~{}s^{-2}$,
the green region has $B=10\rm~{}G$ and also has $\dot{v}=1\rm~{}cm~{}s^{-2}$.
The blue region has $B=3\rm~{}G$ and $\dot{v}=1.4\rm~{}cm~{}s^{-2}$. The
values of the magnetic fields are the same as in Figure 1 of Menou (2012a),
and the values of $\dot{v}$ have the same relative differences. The values of
the necessary parameters are calculated as in Figure 3 for each choice of $B$.
For $\dot{v}$, we applied an exponential decay with respect to pressure,
multiplying the values of $\dot{v}$ above by $\exp(-P/{\rm bar})$, and used
equation (21) to calculate the dimensionless parameter $\dot{V}$.

Figure 4: Frequency $\bm{f}$ in day-1 of the instability domain for the
atmospheres of HJ with magnetic field strength of $B=10~{}\rm G$ and
$B=40~{}\rm G$ and an acceleration of $\dot{v}=1~{}\rm cm~{}s^{-2}$ in green
and red respectively. The blue region represents the results for $B=3~{}\rm G$
and $\dot{v}=1.4~{}\rm cm~{}s^{-2}$. The parameters used to solve equations
(8)–(10) were obtained the same way as in Figure 3. For $\dot{V}$, we applied
the exponential decay, $\exp(-P/{\rm bar})$ and used the scaling presented in
equation (21) to transform $\dot{v}$ to its dimensionless counterpart,
$\dot{V}$. Frequencies are recast in dimensional units using the Alfvén time,
as given by equation (22). The black lines represent the edges of the unstable
regions. Figure 5: Thermo resistive instability domain using Equation (17)
with $\tilde{\eta}T^{2}=1$, combined with Equation (7) of Menou (2012a). The
default magnetic drag law ($n_{s}=n_{e}=1$) was used, as well as the
exponential decay at higher pressures; $\exp(-P/{\rm bar})$. The green and red
regions correspond to input zonal velocities of $V_{\phi}=7\rm~{}km~{}s^{-1}$
with magnetic fields of $B=10$ and $40~{}\rm G$, respectively. The blue region
represent the instability region for $V_{\phi}=10\rm~{}km~{}s^{-1}$ and
$B=3~{}\rm G$. These parameters were chosen to replicate Figure 1 of Menou
(2012a) with equation (17). The black lines represent the edges of the
unstable regions.

The instability regions are well defined within the mbar to a few bars
pressure levels, and within 1250 to 2300 K. The unstable regions are shaded
according to the frequency. We note that none of the simulations for these
parameters reach sustained Alfvénic oscillations, and are all within the burst
regime. The bursts usually take the form of a series of decaying Alfvén waves
after reaching the critical point. Only the hottest part of the blue region
comes close to the Alfvénic regime, where oscillations reach values slightly
above $0.75~{}\omega/\omega_{A}$.

The stripes that can be seen in Figure 4 are related to the number of Alfvén
oscillations that occur following a burst. We find that the sequence of damped
Alfvén oscillations terminates, and the growth of velocity due to $\dot{V}$
resumes, either at a maximum or minimum of velocity, which corresponds to the
time in the oscillation when $b=0$ and ohmic heating vanishes. If the system
reaches the growth phase when $v$ is at a minimum, then it needs to reach the
critical velocity from negative values, like in panel (c) of Figure 1.
However, if the growth phase starts when $v$ is at a maximum and therefore
positive, the system will have a head start in reaching the critical velocity,
reducing the overall period and increasing the frequency of the TRI. As the
parameters are varied, there are discrete transitions in the number of half-
oscillations following the burst, leading to the jumps in frequency that can
be seen in Figure 4.

We have also plotted the instability regions corresponding to the analysis of
Menou (2012a) (given by equation [17]) for different input zonal velocities
and magnetic field strengths in Figure 5. We note that we used different
values for the radiative cooling timescales, shifting our results to higher
temperatures compared to Menou (2012a). However, our results still show the
same behaviour; a higher magnetic field results in an instability region at
lower temperatures (red), and a weaker field yields an instability at higher
temperatures (blue). The green and blue instability domains are divided into
two distinct regions. These separations would disappear if we were to relax
our criterion by only $10\%$.

Both Figure 4 and 5 show unstable regions. We have represented the analogous
regions with the same colours. There are many resemblances to note. First, a
stronger magnetic field will shift the TRI to lower temperatures as the
interaction between flow and the field gets stronger with the latter. A larger
acceleration $\dot{v}$, or zonal velocity $V_{\phi}$, makes the instability
region larger. The overall shape the the regions are also similar. Indeed,
they all resemble a titled droplet, meaning that our results and those of
Menou (2012a) agree to some extent about the areas of instability surrounding
a point in the $P-T$ space.

There are also striking differences between Figure 4 and 5. Arguably the most
important disagreement is that with a larger field (red) our instability
region is smaller than for a smaller field (green). This boils down to the
difference in nature of the instabilities. The instability presented by Menou
(2012a) and in Figure 5 is considered activated when a thermal runaway occurs.
The model presented in this paper is considered unstable when an oscillatory
steady state is reached. In our model a strong magnetic field is less
conducive to producing oscillations, but is more inclined to produce a thermal
runaway in the thermo-resistive framework of Menou (2012a). Our dynamical
model’s unstable regions are also shifted towards higher pressures and lower
temperatures. There is also much more overlap between the different coloured
regions in Figure 4 than in Figure 5. While $\dot{v}$ and $V_{\phi}$ play
analogous roles in the two models, pushing the system to the critical point,
their fundamental natures are very contrasting.

##  5 Discussion

It is no surprise that a sheared fluid threaded by a magnetic field has an
oscillatory response to steady forcing. However, for a temperature-independent
magnetic diffusivity, these oscillations are transient: Ohmic losses damp the
oscillations and the system evolves to a steady state in which the forced
winding of the field lines is balanced by magnetic diffusion (e.g. for a
similar situation in stellar spin down see Charbonneau & MacGregor 1993). We
find that a temperature-dependent magnetic diffusivity introduces new
solutions. The local model introduced in §2 shows that the interplay of Ohmic
diffusion and the temperature-dependent diffusivity leads to long-lived
sustained oscillations in which the magnetic Reynolds number alternates
between large and small values during an oscillation cycle (Figure 1).

This dynamical thermo-resistive instability may be relevant for some regions
of HJ atmospheres (Figure 4). We find that sustained oscillations occur when
(1) the magnetic Reynolds number ${\rm Rm}<1$, so that there is opportunity to
generate an induced field and Ohmic heating; (2) the cooling time is short
compared to the Alfvén timescale, so that the gas temperature can respond to
Ohmic heating on the timescale of an Alfvén oscillation, and (3) for values of
${\rm Ec}=v_{A}^{2}/c_{P}T$ approaching unity, so that Ohmic heating is able
to raise the temperature efficiently. These requirements lead to unstable
regions at intermediate temperatures $T\approx 1500\ {\rm K}$ and lower
pressures $P\lesssim 0.1\ {\rm bars}$, for lower magnetic field strengths
$B\lesssim 10\ {\rm G}$.

The timescale for the oscillations is set by the Alfvén time across a pressure
scale height $t_{A}=H/v_{A}$, which is

$t_{A}\approx 10.5\ {\rm h}\ T_{3}^{1/2}\left(\frac{B}{10\rm
G}\right)^{-1}\,\left(\frac{P}{1\rm
bar}\right)^{1/2}\left(\frac{g}{g_{J}}\right)^{-1}.$ (22)

The sustained oscillations can take the form of modified Alfvén oscillations
with period $\approx t_{A}$ (e.g. middle panel of Figure 1), or bursts of
oscillations separated by long intervals of duration $\gg t_{A}$ (e.g. lower
panel of Figure 1), which could lead to interesting observable effects. For
the parameters expected in HJ atmospheres, we find that the solutions are in
the burst regime. There are two competing effects that determine the
oscillation frequency. One is that the Alfvén time is shorter at lower
pressures, as described by equation (22). The other is that the oscillations
become closer to the Alfvenic regime at higher pressure, which shortens the
oscillation period at higher pressures. The resulting oscillation frequencies
are shown in Figure 4 for different cases. Were such oscillations observed, it
would be possible to use our results to put constraints on the strength of the
magnetic field of the planet observed, knowing the depth and temperature of
the observed region, as well as the period of the oscillation. Possible
degeneracies may arise from such an analysis, as Figure 4 suggests, therefore
our model may only be able to impose an upper limit on the field’s strength.

Our work generalizes the analysis of Menou (2012a) to include the dynamical
effect of the magnetic field. Whereas Menou (2012a) assumed a given background
wind velocity, we have instead assumed a given background acceleration that
drives the flow. This is motivated by the behaviour seen in global dynamical
models, which show net eastward or westward accelerations in the equatorial
region (eastward in hydrodynamic models, westward in magnetized cases; Hindle
et al. 2021b). The resulting instability criterion (equation [15]) is
equivalent to the instability criterion of Menou (2012a) with the wind
velocity taken to be the steady-state velocity under the applied forcing
(equation [11]). Putting in the units, this is a velocity
$V=\dot{V}\eta/v_{A}^{2}=\dot{V}t_{\rm drag}$, i.e. the velocity acquired in a
time

$t_{\mathrm{drag}}\approx 2.5\ {\rm d}\
\left(\frac{x_{e}}{10^{-7}}\right)^{-1}T_{3}^{-1/2}\left(\frac{P}{1\rm
bar}\right)\left(\frac{B}{10\rm G}\right)^{-2},$ (23)

is the magnetic drag timescale (see equation [1] of Rauscher & Menou 2013).
This helps to understand why we prefer lower values of $B$: to achieve the
$\sim{\rm km/s}$ velocities that Menou (2012a) requires for instability with
an acceleration $\sim 1\ {\rm cm\ s^{-2}}$ requires a drag timescale of $\sim
10^{6}\ {\rm s}$. With a constant velocity rather than acceleration, Menou
(2012a) found that increasing $B$ favoured instability; we find the opposite
here (compare equations [15] and [17] and Figures 4 and 5).

Our results add support to the idea that variability should be a feature of
magnetized HJ atmospheres, in particular at intermediate temperatures
$T\approx 1500\ {\rm K}$. Even the highest frequency simulations represent a
period of oscillation of over 10 days. Such timescales are longer than typical
integration times for observations and should therefore be resolved in
observations of HJ atmospheres. However, the timescale of variability could be
much shorter at lower pressures or higher magnetic field strengths (equation
(22)). In particular, the decaying Alfvén waves following bursts could be
challenging to resolve as the Alfvén time is smaller than the period between
bursts.

Spatially-extended MHD simulations of atmospheric dynamics including a
temperature-dependent $\eta$ are needed to go beyond our local model and
address whether the TRI plays a role in generating observable variability. In
particular, the forcing acceleration $\dot{V}$ is put in by hand in our local
model, but arises as part of the global atmospheric dynamics. Indeed, as
explained in detail by Hindle et al. (2021b), global Lorentz forces can lead
to a reversal of the mean flow at the equator in magnetized atmospheres, i.e.
changing the sign of $\dot{V}$ compared to hydrodynamic models. Similarly, our
local analysis shows that the TRI is unstable only within a range of
pressures; understanding the interaction of the low pressure unstable region
with the deeper stable regions requires going beyond a local analysis.

It is important to note that previous work on MHD models of HJ atmospheres has
shown that variable flows arise even with a time-independent $\eta$ profile.
Rogers & Komacek (2014) found that as the influence of the magnetic field
increased (either because of a hotter atmosphere or stronger magnetic field),
the MHD effects changed from a drag force, to oscillatory mean flows, to a
reversed (westward) mean flow. Hindle et al. (2021b) explain the reversal of
the mean flow in terms of Lorentz forces associated with the toroidal magnetic
fields generated by zonal flows. In addition, Rogers & McElwaine (2017) showed
that the spatial gradient of $\eta$ can lead to magnetic field generation in
an atmospheric dynamo. The temperature-dependence of $\eta$ is a further
important ingredient to include in future MHD models of HJ atmospheric
dynamics.

We thank the anonymous reviewer for helpful comments that improved this
manuscript. This work was supported by NSERC Discovery grants. RH, AC, and PC
are members of the Centre de Recherche en Astrophysique du Québec (CRAQ).

## References

  * Armstrong et al. (2016) Armstrong, D. J., de Mooij, E., Barstow, J., et al. 2016, Nature Astronomy, 1, 0004, doi: 10.1038/s41550-016-0004
  * Batygin & Stevenson (2010) Batygin, K., & Stevenson, D. J. 2010, ApJ, 714, L238, doi: 10.1088/2041-8205/714/2/L238
  * Bell et al. (2019) Bell, T. J., Zhang, M., Cubillos, P. E., et al. 2019, MNRAS, 489, 1995, doi: 10.1093/mnras/stz2018
  * Beltz et al. (2022) Beltz, H., Rauscher, E., Roman, M. T., & Guilliat, A. 2022, AJ, 163, 35, doi: 10.3847/1538-3881/ac3746
  * Cauley et al. (2019) Cauley, P. W., Shkolnik, E. L., Llama, J., & Lanza, A. F. 2019, Nature Astronomy, 3, 1128, doi: 10.1038/s41550-019-0840-x
  * Charbonneau & MacGregor (1993) Charbonneau, P., & MacGregor, K. B. 1993, ApJ, 417, 762, doi: 10.1086/173357
  * Cooper & Showman (2005) Cooper, C. S., & Showman, A. P. 2005, ApJ, 629, L45, doi: 10.1086/444354
  * Dang et al. (2018) Dang, L., Cowan, N. B., Schwartz, J. C., et al. 2018, Nature Astronomy, 2, 220, doi: 10.1038/s41550-017-0351-6
  * Davidson (2001) Davidson, P. A. 2001, An Introduction to Magnetohydrodynamics (Cambridge University Press)
  * Hindle et al. (2019) Hindle, A. W., Bushby, P. J., & Rogers, T. M. 2019, ApJ, 872, L27, doi: 10.3847/2041-8213/ab05dd
  * Hindle et al. (2021a) —. 2021a, ApJ, 916, L8, doi: 10.3847/2041-8213/ac0fec
  * Hindle et al. (2021b) —. 2021b, ApJ, 922, 176, doi: 10.3847/1538-4357/ac0e2e
  * Imamura et al. (2020) Imamura, T., Mitchell, J., Lebonnois, S., et al. 2020, Space Sci. Rev., 216, 87, doi: 10.1007/s11214-020-00703-9
  * Jackson et al. (2019) Jackson, B., Adams, E., Sandidge, W., Kreyche, S., & Briggs, J. 2019, AJ, 157, 239, doi: 10.3847/1538-3881/ab1b30
  * Kataria et al. (2016) Kataria, T., Sing, D. K., Lewis, N. K., et al. 2016, ApJ, 821, 9, doi: 10.3847/0004-637X/821/1/9
  * Komacek & Showman (2016) Komacek, T. D., & Showman, A. P. 2016, ApJ, 821, 16, doi: 10.3847/0004-637X/821/1/16
  * Menou (2012a) Menou, K. 2012a, ApJ, 754, L9, doi: 10.1088/2041-8205/754/1/L9
  * Menou (2012b) —. 2012b, ApJ, 745, 138, doi: 10.1088/0004-637X/745/2/138
  * Perna et al. (2010a) Perna, R., Menou, K., & Rauscher, E. 2010a, ApJ, 719, 1421, doi: 10.1088/0004-637X/719/2/1421
  * Perna et al. (2010b) —. 2010b, ApJ, 724, 313, doi: 10.1088/0004-637X/724/1/313
  * Rauscher & Menou (2010) Rauscher, E., & Menou, K. 2010, ApJ, 714, 1334, doi: 10.1088/0004-637X/714/2/1334
  * Rauscher & Menou (2013) —. 2013, ApJ, 764, 103, doi: 10.1088/0004-637X/764/1/103
  * Read & Lebonnois (2018) Read, P. L., & Lebonnois, S. 2018, Annual Review of Earth and Planetary Sciences, 46, 175, doi: 10.1146/annurev-earth-082517-010137
  * Reiners & Christensen (2010) Reiners, A., & Christensen, U. R. 2010, A&A, 522, A13, doi: 10.1051/0004-6361/201014251
  * Rogers (2017) Rogers, T. M. 2017, Nature Astronomy, 1, 0131, doi: 10.1038/s41550-017-0131
  * Rogers & Komacek (2014) Rogers, T. M., & Komacek, T. D. 2014, ApJ, 794, 132, doi: 10.1088/0004-637X/794/2/132
  * Rogers & McElwaine (2017) Rogers, T. M., & McElwaine, J. N. 2017, ApJ, 841, L26, doi: 10.3847/2041-8213/aa72da
  * Showman et al. (2008) Showman, A. P., Cooper, C. S., Fortney, J. J., & Marley, M. S. 2008, ApJ, 682, 559, doi: 10.1086/589325
  * Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564, doi: 10.1088/0004-637X/699/1/564
  * Showman & Guillot (2002) Showman, A. P., & Guillot, T. 2002, A&A, 385, 166, doi: 10.1051/0004-6361:20020101
  * Showman & Polvani (2011) Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71, doi: 10.1088/0004-637X/738/1/71
  * von Essen et al. (2020) von Essen, C., Mallonn, M., Borre, C. C., et al. 2020, A&A, 639, A34, doi: 10.1051/0004-6361/202037905
  * Yadav & Thorngren (2017) Yadav, R. K., & Thorngren, D. P. 2017, ApJ, 849, L12, doi: 10.3847/2041-8213/aa93fd