# Hot Spot Offset Variability from Magnetohydrodynamical Thermoresistive
Instability in Hot Jupiters

Raphaël Hardy Département de Physique, Université de Montréal, Montréal, QC,
H3C 3J7, Canada
Department of Physics and Trottier Space Institute, McGill University,
Montréal, QC, H3A 2T8, Canada
Institut Trottier 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
Andrew Cumming Department of Physics and Trottier Space Institute, McGill
University, Montréal, QC, H3A 2T8, Canada
Institut Trottier de Recherche sur les Exoplanètes (iREx), Université de
Montréal, Montréal, QC H3C 3J7, Canada

###### Abstract

Hot Jupiter atmospheres are possibly subject to a thermoresistive instability.
Such an instability may develop as the ohmic heating increases the electrical
conductivity in a positive feedback loop, which ultimately leads to a runaway
of the atmospheric temperature. We extend our previous axisymmetric one-
dimensional radial model, by representing the temperature and magnetic
diffusivity as a first order Fourier expansion in longitude. This allows us to
predict the hot spot offset during the unfolding of the thermoresistive
instability and following Alfvénic oscillations. We show a representative
simulation undergoing the thermoresistive instability, in which the peak flux
offset varies between approximately $\pm 60^{\circ}$ on timescales of a few
days with potentially observable brightness variations. Therefore, this
thermoresistive instability could be an observable feature of hot Jupiters,
given the right timing of observation and transit and the right planetary
parameters.

magnetohydrodynamics (MHD) – magnetic diffusivity – planets and satellites:
gaseous planets – planets and satellites: atmospheres – planets and
satellites: magnetic fields

##  1 Introduction

Hot Jupiter (HJ) atmospheres are interesting case studies of extreme
atmospheric dynamics. Being tidally locked due to their proximity ($\lesssim
0.1~{}{\rm AU}$) to their host star, the extreme radiative heating on their
dayside sustains large temperature gradients between their locked day and
nightside, which drives equatorial jets (Showman & Polvani, 2011; Komacek &
Showman, 2016; Read & Lebonnois, 2018; Imamura et al., 2020). These jets are
usually prograde in both observations and hydrodynamical simulations (Showman
& Guillot, 2002; Cooper & Showman, 2005; Showman et al., 2009; Rauscher &
Menou, 2010; Kataria et al., 2016). The jets advect the hot spot, the hottest
region of the atmosphere, away from the substellar point either eastward or
westward depending on the direction of the zonal winds. There are, however, a
few exceptions displaying retrograde jets (Armstrong et al., 2016; Dang et
al., 2018; Bell et al., 2019; Jackson et al., 2019; von Essen et al., 2020).

The temperature regime characterizing HJs leads to partial ionization of their
atmospheres, which can then couple to the magnetic field (Rogers & Komacek,
2014). Magnetic effects have been proposed as the cause of retrograde winds
and westward hot spots, as well as atmospheric variability (Rogers & Komacek,
2014; Rogers, 2017; Hindle et al., 2019, 2021a, 2021b; Hardy et al., 2022,
2023). The main source of ionization is alkali metals such as potassium and
sodium, which have low first ionization energies (Perna et al., 2010; Batygin
& Stevenson, 2010). Previous studies have investigated magnetic coupling
between the winds in the upper atmosphere and the magnetic field pervading the
planetary interior, leading to slower prograde equatorial jets, or even
retrograde jets (Rauscher & Menou, 2013; Rogers, 2017; Menou, 2012a; Hardy et
al., 2022, 2023).

As the alkali metals just start being ionized at the temperatures of HJ
atmospheres, the ionization fraction is very temperature sensitive. Thus,
small temperature variations can have very large impact on the magnetic
diffusivity (MD) $\eta$ and magnetic coupling (Perna et al., 2010). Menou
(2012a) showed that with $\eta$ decreasing as temperature increases can drive
a thermoresistive instability (TRI), leading to runaway ohmic heating of the
atmosphere (see also Hubbard et al. 2012 and Price et al. 2012). Rauscher &
Menou (2013) included this effect in 3D atmospheric circulation models under
the assumption that the magnetic Reynolds number Rm remains small and so the
magnetic forces can be treated as a drag term (Perna et al., 2010).

In order to explore the full range of dynamics including the regime of large
Rm where the plasma is strongly-coupled to the field, we developed a
simplified axisymmetric model of the equatorial plane of a HJ (Hardy et al.
2022, 2023, hereafter H22 and H23). The model assumes that angular momentum is
continuously pumped into the equatorial plane from higher latitudes, driving a
zonal flow at the equator (Showman & Polvani, 2011). The radial component of
the magnetic field at the equator supports torsional Alfven waves and
interacts with the flow. This geometry enables a study of the feedback between
ohmic heating, the evolving $\eta$, and the dynamics, while also including the
full radial structure of the atmosphere. We showed that the outcome of the TRI
in HJs with intermediate temperatures (equilibrium temperatures
$T_{\mathrm{eq}}\approx 1000$–$1200\ \mathrm{K}$) is to create bursts of
Alfven oscillations separated by longer periods of quiescence. These time-
dependent bursting solutions are not present when $\eta$ is assumed to be
time-independent; they represent a new class of time-dependent behaviour
driven by the temperature-dependence of $\eta$.

In this paper, we extend the model of H23 by relaxing the assumption of
axisymmetry. We do this by using a low order Fourier expansion in the
azimuthal angle $\phi$, following the approach of Tritton (1988) who studied
the problem of convection in a torus. This expansion enables us to include the
variation in $\eta$ with $\phi$ in an approximate way, allowing a study of the
displacement angle of the hot spot (temperature maximum) during the different
phases of the bursts created by the TRI.

The paper is organized as follows. We describe the model in Section 2, and
discuss in detail the properties of a single a representative solution,
followed by Section 3 where a briefer discussion of behavior variations across
the model’s parameter space. We close the paper (Section 4) by discussing
possible observational signatures of the TRI, including variations in the hot
spot offset and the atmospheric temperature.

##  2 Extended One Dimensional Model in the Equatorial Plane

###  2.1 Model Setup

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

${\bf B}(z,t)=B_{x}(z,t){\hat{x}}+B_{0}{\hat{z}},$ (1)

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

$\frac{\partial u_{x}}{\partial
t}=\frac{B_{0}}{{\mu}_{0}{\bar{\rho}}}\frac{\partial B_{x}}{\partial
z}+\frac{{\bar{\mu}}}{{\bar{\rho}}}\frac{\partial^{2}u_{x}}{\partial
z^{2}}+a_{x},$ (2) $\frac{\partial B_{x}}{\partial t}=B_{0}\frac{\partial
u_{x}}{\partial z}+\frac{\partial\eta}{\partial z}\frac{\partial
B_{x}}{\partial z}+\eta\frac{\partial^{2}B_{x}}{\partial z^{2}},$ (3)
$\displaystyle\frac{\partial T}{\partial t}+u_{x}\frac{\partial T}{\partial
x}$
$\displaystyle=\frac{1}{{\bar{\rho}}c_{p}}\frac{\partial\bar{\chi}}{\partial
z}\frac{\partial T}{\partial
z}+\frac{\bar{\chi}}{{\bar{\rho}}c_{p}}\frac{\partial^{2}T}{\partial
z^{2}}+\frac{1}{{\bar{\rho}}c_{p}}\frac{\partial F_{\rm irr}}{\partial z}$ (4)
$\displaystyle+\frac{{\bar{\mu}}}{{\bar{\rho}}c_{p}}\left(\frac{\partial
u_{x}}{\partial
z}\right)^{2}+\frac{\eta}{{\mu}_{0}{\bar{\rho}}c_{p}}\left(\frac{\partial
B_{x}}{\partial z}\right)^{2}.$

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

$a_{x}=\dot{v}\exp({-P/\mathrm{bar}})~{},$ (5)

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

$\eta(T)=0.023\frac{\sqrt{T}}{\chi_{e}}{\rm m^{2}~{}s^{-1}},$ (6)

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

${\bar{\chi}}=\frac{16\sigma T^{3}}{3\kappa_{\rm th}{\bar{\rho}}},$ (7)

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

${\bar{\mu}}=\frac{\Pr\bar{\chi}}{c_{p}}.$ (8)

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

$F_{\rm irr}=F_{s}(\phi)\exp(-\sqrt{3}\kappa_{v}P/g),$ (9)

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).

###  2.2 Longitudinal Expansion

We now develop a low order expansion in the longitudinal direction, applying
the approach111Tritton (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

$\tilde{T}(\phi)=T_{0}+T_{1}\cos\phi+T_{2}\sin\phi,$ (10)
$\tilde{\eta}(\phi)={\eta_{0}}+{\eta_{1}}\cos\phi+{\eta_{2}}\sin\phi.$ (11)

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
point222We 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}$:

$\displaystyle\frac{\partial T_{0}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{1}{{\bar{\rho}}c_{p}}\frac{\partial\bar{\chi}}{\partial
z}\frac{\partial T_{0}}{\partial
z}+\frac{\bar{\chi}}{{\bar{\rho}}c_{p}}\frac{\partial^{2}T_{0}}{\partial
z^{2}}+\frac{1}{2}\frac{1}{{\bar{\rho}}c_{p}}\frac{\partial F_{\rm
irr}}{\partial z}$ (12)
$\displaystyle+\frac{{\bar{\mu}}}{{\bar{\rho}}c_{p}}\left(\frac{\partial
u_{x}}{\partial
z}\right)^{2}+\frac{\eta_{0}}{{\mu}_{0}{\bar{\rho}}c_{p}}\left(\frac{\partial
B_{x}}{\partial z}\right)^{2},$ $\displaystyle\frac{\partial T_{1}}{\partial
t}$ $\displaystyle=$
$\displaystyle-\frac{u_{x}T_{2}}{R}+\frac{1}{{\bar{\rho}}c_{p}}\frac{\partial\bar{\chi}}{\partial
z}\frac{\partial T_{1}}{\partial
z}+\frac{\bar{\chi}}{{\bar{\rho}}c_{p}}\frac{\partial^{2}T_{1}}{\partial
z^{2}}$ (13)
$\displaystyle+\frac{1}{2}\frac{1}{{\bar{\rho}}c_{p}}\frac{\partial F_{\rm
irr}}{\partial
z}+\frac{\eta_{1}}{{\mu}_{0}{\bar{\rho}}c_{p}}\left(\frac{\partial
B_{x}}{\partial z}\right)^{2},$ $\displaystyle\frac{\partial T_{2}}{\partial
t}$ $\displaystyle=$
$\displaystyle\frac{u_{x}T_{1}}{R}+\frac{1}{{\bar{\rho}}c_{p}}\frac{\partial\bar{\chi}}{\partial
z}\frac{\partial T_{2}}{\partial
z}+\frac{\bar{\chi}}{{\bar{\rho}}c_{p}}\frac{\partial^{2}T_{2}}{\partial
z^{2}}$ (14)
$\displaystyle+\frac{\eta_{2}}{{\mu}_{0}{\bar{\rho}}c_{p}}\left(\frac{\partial
B_{x}}{\partial z}\right)^{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

$\eta_{0}=\frac{1}{n_{\phi}}\sum_{\phi}\eta(\phi),$ (15)
$\eta_{1}=\frac{\sum_{\phi}\eta(\phi)\cos(\phi)}{\sum_{\phi}\cos^{2}(\phi)},$
(16)
$\eta_{2}=\frac{\sum_{\phi}\eta(\phi)\sin(\phi)}{\sum_{\phi}\sin^{2}(\phi)},$
(17)

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).

###  2.3 Boundary Conditions, Initial Conditions, and Choice of Parameter
Values

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

$-\chi\frac{\partial T}{\partial z}=\frac{4}{3}\sigma T^{4},$ (18)

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}$.

Parameters | Values
---|---
$B_{0}$ (G) | [1, 100]
$\dot{v}$ (m s-2) | [0.0001, 0.01]
$\kappa_{\rm th}$ (m2 kg-1) | [0.0001, 0.001]
$T_{\rm eq}$ (K) | [1000, 1200]
Pr | [0.01]
Table 1:  Ranges of the parameters used in the simulations. See H23 for
further justifications for the parameter ranges.

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)$.

Figure 1: An example of bursting behavior driven by TRI. We show the time
series of (a) velocity and magnetic field, (b) temperature and magnetic
Reynolds number, and (c) angle of the maximum temperature. All quantities are
measured in the center of the domain. The model shown has 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. The right panels
zoom in on the gray shaded area indicated in the left panels. Figure 2: Ratios
of timescales in the system during oscillations; (a) the advection timescale
($\tau_{\rm adv}=R/\absolutevalue{u_{x}}$) versus the thermal diffusion
timescale ($\tau_{\rm diff}=H^{2}/\chi$), (b) the advection timescale versus
the Alfvén timescale ($\tau_{\rm Alf}=H/u_{A}$). The length scale used to
measure these is the pressure scale height, $H$, for $\tau_{\rm Alf}$ and
$\tau_{\rm diff}$ which we update with the axisymmetric temperature value at
every time step and the radius is used for $\tau_{\rm adv}$.

###  2.4 Numerical Scheme

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.

Figure 3: Spatiotemporal evolution of (a) the velocity, (b) the magnetic
field, (c) temperature, (d) the longitudinal temperature difference between
the hottest and coldest points at each depth, (e) angular offset of the
hottest point at each depth, and (f) magnetic Reynolds number for the same
simulation as in Figure 1.

##  3 Results

###  3.1 A Representative Simulation

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 4: Depth-longitude temperature maps during three different instants;
(a) initial condition, (b) first build-up phase and (c) during the first burst
of the simulation. The thin black lines each represent an isocontour separated
by 50 K from each others, and the thicker black line track the hottest point
at each layers.

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.

Figure 5: Day and nightside angle-averaged thermal flux at the surface (full
and dotted lines respectively), normalized by the flux value on the dayside at
the start of the simulation (top panel), and the angular offset of the maximum
in the thermal flux (bottom panel). We show results for simulations with a
range of different $B_{0}$ values. The representative solution from §3.1
($B_{0}=30\ \mathrm{G}$) is plotted in orange.

###  3.2 Dependence on input parameters

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}$.

##  4 Discussion

We have presented here a simple model which encapsulates the longitudinal
dependency of the temperature in a HJ atmosphere susceptible to the TRI.
Building on the works of H23, this model has allowed us to characterize the
reaction of the hot spot offset to our imposed acceleration and the TRI. We
were able to simulate the hot spot displacement during a burst triggered by
TRI from 0.01 bar to 1 bar. These simulations suggest that it may be possible
to detect the occurrence of this instability in the atmospheres of these gas
giants. Our simulations indicate that the TRI could produce significant hot
spot offset variations, although somewhat larger than the offset range found
from Spitzer’s phase curves by Bell et al. (2021). Indeed, the simulation
presented in this work showed peak flux offset ranging from $66^{\circ}$
eastward to $56^{\circ}$ westward. The recurrence time and Alfvénic
oscillation periods found in H23 are still valid for this longitudinally-
extended model, with our representative solution being on the shorter end of
the period spectrum. The reanalysis of phase curves of 16 planets done by Bell
et al. (2021) shows offsets between $-38.7^{\circ+3.2}_{~{}-3.2}$ (CoRoT-2b)
and $43.4^{\circ+5.4}_{~{}-6.1}$ (HD209458b) in the colder regime ($T_{\rm
irr}<2500$ K). Comparing our results to the model presented in Komacek et al.
(2017), our offset results for oscillating systems are comparable to the
fastest superrotating equatorial jets presented in that study
($8\rm~{}km~{}s^{-1}$), and long drag timescales ($\tau_{\rm
drag}=10^{6}-10^{7}$ s). Indeed, Komacek et al. (2017) are expecting offsets
of around $60-70^{\circ}$ with these physical characteristics. Figure 6 shows
drag timescale throughout the representative solution. Outside of the TRI, the
viscous drag is dominating, with associated timescale in the range
$10^{5}-10^{8}$ s, depending on depth. During TRI, magnetic drag is the
biggest contributor, dropping the drag timescale down to $10^{3}-10^{4}$ s.
The repercussions of this drop is seen in the hot spot offset. Indeed, after
the decaying oscillations and before build-up starts, the offset is near
$0^{\circ}$, as predicted by models with small drag timescales in Komacek et
al. (2017). Thus, our offset values seem reasonable considering our simplified
model.

Figure 6: Total drag timescale defined as $\tau_{\rm
drag}=[\tau_{\nu}^{-1}+\tau_{\rm mag}^{-1}]^{-1}$, where
$\tau_{\nu}=H^{2}/\nu$ is the viscous drag timescale and $\tau_{\rm
mag}=\mu_{0}\rho\eta_{0}/B_{0}^{2}$ is the magnetic drag timescale. Viscous
drag timescale ($\tau_{\nu}$) is shorter during build-up, but as TRI occurs,
$\tau_{\rm mag}$ becomes much smaller and the magnetic drag timescale
dominates the flow dynamics. The black line denotes where
$\tau_{\nu}=\tau_{\rm mag}$.

Keating et al. (2019) have inferred the temperatures on the day and nightside
for 12 exoplanets. The temperature contrast between the dayside and the
nightside varies between 150 K for HD189733b and 1040 K for WASP-18b. We note
however that only HD189733b and HD209458b have dayside temperatures below 1500
K. The latter has a temperature contrast of 189 K. The rest of the sample
presented in that paper is too hot to be conducive to the TRI. We note that
the analysis presented by Bell et al. (2021) find a temperature contrast of
287 K for HD189733b. The representative solution of Section 3.1 has a maximum
temperature contrast of 310 K after the TRI burst. However, during build-up as
the system nears steady-state before the TRI triggers, the temperature
contrast is around 50 K only. Thus, assuming the temperature contrast of
Keating et al. (2019) are from steady-state planets, our model seems to
underestimate this characteristic of the atmosphere. Our first order Fourier
expansion restricts the way we can heat the atmosphere in our modeling. As our
heating is modulated by $(1+\sin(\phi))/2$, we are also heating parts of the
nightside, reducing the temperature contrast between day and nightside.

The flux variation caused by the thermal runaway of the TRI, as shown in
Figure 5, would certainly be the most conspicuous observable telltale sign of
the TRI. Indeed, for our representative solution, the outward thermal flux
momentarily grows fivefold. Therefore, any planets with a much larger
luminosity than expected may be undergoing TRI. However, the short-comings of
our model must again be emphasized. The dynamical variables, $u_{x}$ and
$B_{x}$, are assumed axisymmetric. The longitudinal velocity is expected to
have many small scale structures interacting with higher and lower latitudes,
breaking symmetry. Large scale flows coming from higher and lower latitudes
would also break symmetry. As for the magnetic field, only a perfectly aligned
dipole would be axisymmetric. With such asymmetry, we would assume to be
harder for the TRI to trigger globally as turbulence, local instabilities,
etc. may slow down winds, thus reducing ohmic heating. On the other hand,
turbulence may enhance ohmic heating through shearing. Thus, the exact effects
of a more complex model on the TRI is unclear. TRI may rather operate locally
and significantly reduce the predicted values of Figure 5. In addition, while
the model presented in this work requires ohmic heating to reach the critical
magnetic Reynolds number, such heating could be secondary if considering the
increase in temperature at the morning terminator. Indeed, as the cold gas of
the nightside arrives on the dayside, stellar flux could be the heat source
required to trigger TRI at the morning terminator.

Simulations with larger magnetic field strengths tend to show shorter
recurrence periods and smaller oscillation amplitudes, and vice versa, for a
given set of parameters. From H23, we recall that beyond a certain magnetic
field strength, for a given set of $T_{\rm eq}$, $\dot{v}$ and $\kappa_{\rm
th}$, the system is no longer susceptible to the TRI, and instead reaches a
steady state. Therefore, as the field strength increases, the offset range
diminishes. Moreover, as argued in Section 3.2, the other observable
characteristic of TRI, the outgoing flux, should also favor weaker field, as
stronger fields show smaller variations in their outgoing thermal flux. The
estimated dipole field strength by Yadav & Thorngren (2017) for HJs averages
around 100 G, but as our input field strength represent the radial component,
and as argued in H23, a 100 G dipole field yields a radial component up to
17.4 G with an inclination of $10^{\circ}$ in respect to the rotation axis. At
the maximum value predicted by Yadav & Thorngren (2017), $\approx 250$ G, and
the same inclination, we would get a radial field strength of 43.4 G. Thus,
given the right set of parameters, all field strength proposed by Yadav &
Thorngren (2017) are potential candidates for TRI.

As the true physical boundary conditions on the different temperature
component is unclear, we have tested many different combinations, before
settling on the one presented in this work. We tried having $T_{1}$ and
$T_{2}$ fixed to zero at both boundaries. At the bottom, this would represent
a core with an axisymmetric temperature. At the top, it would mimic an
efficient thermal conductor. For the axisymmetric part, we kept it constant at
the top and kept the internal planetary flux, as in H23. However, with these
assumptions, the temperature contrast had a maximum around half that
characterizing our representative simulation, i.e., much smaller than
anticipated from other numerical models and observations. Keeping the
temperature constant at the outermost layer also has the effect of radically
increasing the outgoing thermal flux. The angle averaged thermal flux coming
out of the atmosphere would increase tenfold during TRI. Thus, the flux-
temperature relation was chosen for the top boundary condition on the
temperature components, and an insulating bottom boundary condition for
$T_{1}$ and $T_{2}$ was also adopted.

While we only presented results from simulations with Prandtl number of 0.01,
we tested larger values, i.e. giving more importance to viscosity. A larger
Prandtl number has the effect of reducing the drag timescale, thus reducing
the flow velocity, but also enhancing viscous heating. This shifts the region
of parameter space susceptible to TRI to lower temperature and stronger
forcing, as to counteract the effect of stronger viscous drag. However,
velocities of around $10~{}\rm km~{}s^{-1}$ seemed to still be required to
trigger TRI. For example, Pr=0.1 yields periodic TRI at equilibrium
temperatures of 700-800 K, putting these atmospheres closer to the warm
Jupiters regime rather than the hot Jupiters one. Moreover, enhancing the
Prandtl number beyond the expected value from micro physics could also be a
way to take into account omitted physical phenomenon, such as flow out of the
equatorial plane, shocks, and turbulence.

Finally, we want to point out that the magnetic Prandtl number (${\rm
Pm}=\nu/\eta_{0}$) always remain below unity. For the representative solution
of Section 3.1, it reaches a maximum value of $1.4\times 10^{-1}$, and goes
down to $\approx 10^{-7}$ at depth during the quiescent phase, when $\eta_{0}$
is at its largest. Thus, ohmic dissipation is always greater than viscous
dissipation.

In summary, together with H22 and H23, our results show that a temperature
dependent MD can impact the atmospheric dynamics drastically. Thus, in the
regimes where the magnetic Reynolds number is around unity, it would be
imperative to further extend models to implement such dependence. A complete
3D model or simulation would be necessary to incorporate all physical
mechanisms susceptible to impact, positively or negatively, the TRI. Such a
model would also lead to better predictions of the observable offset
variability caused by the TRI, but also the changes in luminosity associated
with the thermal runaway of the instability.

This work was supported by the Natural Sciences and Engineering Research
Council of Canada (NSERC) Discovery grants RGPIN-2018-05278, RGPIN-2023-03620,
and RGPIN-2024-04050. R.H., P.C., and A.C. are members of the Centre de
Recherche en Astrophysique du Québec (CRAQ). A.C. is grateful to the Isaac
Newton Institute for Mathematical Sciences, Cambridge, for support and
hospitality during the programme “Anti-diffusive Dynamics: from sub-cellular
to astrophysical scales” supported by EPSRC grant no EP/R014604/1.

##  Appendix A Accuracy of the sinusoidal assumption for the angular
dependence of the magnetic diffusivity

As discussed in the main text, because $\eta(T)$ is strongly non-linear in $T$
(Equation (6)), a sinusoidal $T(\phi)$ as assumed in Equation (10) does not
give $\eta(\phi)$ that is also sinusoidal. Therefore, when we expand
$\eta(\phi)$ in the form of Equation (11), with $\sin\phi$ and $\cos\phi$
terms only, we introduce an inaccuracy.

To assess the magnitude of this error, Panel (a) of Figure 7 compares the true
$\eta(\phi)$ corresponding to $T(\phi)$ obtained using Equation (6) (solid
lines in color) and the corresponding Fourier fit $\tilde{\eta}(\phi)$
(Equation (11)) (black dotted lines). We show curves at various times during
the oscillation cycle of the TRI, during which $\eta$ varies by more than
three orders of magnitude as the atmospheric temperature changes. We note how
at some epochs, the fitted diffusivity may reach negative values at some
longitudes. These negative values happen when the range of the diffusivity is
largest, making it difficult to fit a sinusoidal curve onto it. The relative
deviations between the real profiles and their fit are plotted in panel (b).
Considering the large range over which $\eta$ varies over the oscillation, we
find that the $\eta(\phi)$ profiles are reproduced reasonably by the
sinusoidal approximation, with largest deviations reaching 250% locally, but
usually remaining around 20%.

Panels (c) and (d) of Figure 7 are similar to panels (a) and (b) except now we
compare the longitudinally-averaged component $\eta_{0}$ with the true
$\eta(\phi)$ profile. This is important because, as discussed in section 2, we
use $\eta_{0}$ to replace $\eta$ in the induction equation. We underestimate
the true $\eta$ by more than a factor of six at the substellar point at some
phases of the Alfvénic oscillations, down to 8% right before TRI. This leads
to a small change in the instability boundary in parameter space compared to
H23, as discussed in section 3.2.

##  Appendix B Ratio of magnetic pressure to gas pressure

Because our model precludes vertical motions, there is no consequence to the
magnetic pressure $P_{\mathrm{mag}}$ exceeding the gas pressure
$P_{\mathrm{gas}}$. In reality this could possibly lead to a buoyancy-driven
instability (e.g. Newcomb 1961; Parker 1975; Fan 2021; see H23 for further
discussion). We do find that many models reach
$P_{\mathrm{mag}}/P_{\mathrm{gas}}$ ratios larger than unity in some small
regions of the atmosphere. As an example, Figure 8 shows
$P_{\mathrm{mag}}/P_{\mathrm{gas}}$ (the inverse plasma-$\beta$) for the
representative model discussed in Section 3. It can be seen that
$P_{\mathrm{mag}}/P_{\mathrm{gas}}$ reaches $\approx 3$ near the base of the
layer during the Alfven oscillations. Investigation of the effects of this,
for example whether instabilities have time to grow and the subsequent
vertical transport, would be an interesting avenue for future investigation.

Figure 7: (a) $\eta(\phi)$ corresponding to $T(\phi)$ obtained using Equation
(6) (solid lines in color) and the corresponding Fourier fit
$\tilde{\eta}(\phi)$ (Eq. 11) (black dotted lines). (b) The fractional
difference between $\eta(\phi)$ and $\tilde{\eta}(\phi)$. The values are taken
equidistant in time during a burst and the following decaying Alfvén waves
(right-hand side of Figure 1). (c) Similar to panel (a) but now comparing
$\eta(\phi)$ with the longitudinal average $\eta_{0}$. (d) The fractional
difference between $\eta(\phi)$ and $\eta_{0}$. Figure 8: Ratio of magnetic to
gas pressure in $\log_{10}$ as a function of time and pressure, for the
solution presented in Section 3.1.

## 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
  * Bell et al. (2021) Bell, T. J., Dang, L., Cowan, N. B., et al. 2021, MNRAS, 504, 3316, doi: 10.1093/mnras/stab1027
  * 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)
  * Dietrich et al. (2022) Dietrich, W., Kumar, S., Poser, A. J., et al. 2022, MNRAS, 517, 3113, doi: 10.1093/mnras/stac2849
  * Draine et al. (1983) Draine, B. T., Roberge, W. G., & Dalgarno, A. 1983, ApJ, 264, 485, doi: 10.1086/160617
  * Fan (2021) Fan, Y. 2021, Living Reviews in Solar Physics, 18, 5, doi: 10.1007/s41116-021-00031-2
  * Guillot (2010) Guillot, T. 2010, A&A, 520, A27, doi: 10.1051/0004-6361/200913396
  * Hardy et al. (2023) Hardy, R., Charbonneau, P., & Cumming, A. 2023, ApJ, 959, 41, doi: 10.3847/1538-4357/ad0968
  * Hardy et al. (2022) Hardy, R., Cumming, A., & Charbonneau, P. 2022, ApJ, 940, 123, doi: 10.3847/1538-4357/ac9bfc
  * 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, 922, 176, doi: 10.3847/1538-4357/ac0e2e
  * Hindle et al. (2021b) —. 2021b, ApJ, 916, L8, doi: 10.3847/2041-8213/ac0fec
  * Hubbard et al. (2012) Hubbard, A., McNally, C. P., & Mac Low, M.-M. 2012, ApJ, 761, 58, doi: 10.1088/0004-637X/761/1/58
  * 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
  * Keating et al. (2019) Keating, D., Cowan, N. B., & Dang, L. 2019, Nature Astronomy, 3, 1092, doi: 10.1038/s41550-019-0859-z
  * Komacek & Showman (2016) Komacek, T. D., & Showman, A. P. 2016, ApJ, 821, 16, doi: 10.3847/0004-637X/821/1/16
  * Komacek et al. (2017) Komacek, T. D., Showman, A. P., & Tan, X. 2017, ApJ, 835, 198, doi: 10.3847/1538-4357/835/2/198
  * Lodders (2010) Lodders, K. 2010, Astrophysics and Space Science Proceedings, 16, 379, doi: 10.1007/978-3-642-10352-0_8
  * Lorenz (1963) Lorenz, E. N. 1963, Journal of the Atmospheric Sciences, 20, 130, doi: 10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2
  * 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
  * Newcomb (1961) Newcomb, W. A. 1961, Physics of Fluids, 4, 391, doi: 10.1063/1.1706342
  * Parker (1975) Parker, E. N. 1975, ApJ, 198, 205, doi: 10.1086/153593
  * Perna et al. (2010) Perna, R., Menou, K., & Rauscher, E. 2010, ApJ, 719, 1421, doi: 10.1088/0004-637X/719/2/1421
  * Price et al. (2012) Price, S., Link, B., Epstein, R. I., & Li, H. 2012, MNRAS, 420, 949, doi: 10.1111/j.1365-2966.2011.19807.x
  * 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
  * 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. (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
  * Tritton (1988) Tritton, D. J. 1988, Physical fluid dynamics /2nd revised and enlarged edition/ (Osford Science), 246–254
  * 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


# Progress in the development of a KITWPA for the DARTWARS project

M. Borghesi  C. Barone  M. Borghesi  S. Capelli  G. Carapella  A. P. Caricato
I. Carusotto  A. Cian  D. Di Gioacchino  E. Enrico  P. Falferi  L. Fasolo  M.
Faverzani  E. Ferri  G. Filatrella  C. Gatti  A. Giachero  D. Giubertoni  V.
Granata  A. Greco11footnotemark: 1 C. Guarcello  D. Labranca  A. Leo  C. Ligi
G. Maccarrone  F. Mantegazzini  B. Margesin  G. Maruccio  C. Mauro  R. Mezzena
A. G. Monteduro  A. Nucciotti  L. Oberto  L. Origo  S. Pagano  V. Pierro  L.
Piersanti  M. Rajteri  A. Rettaroli  S. Rizzato  A. Vinante  M. Zannoni  INFN
- Laboratori Nazionali di Frascati, Via Enrico Fermi, I-00044, Frascati, Italy
INFN – Sezione di Lecce, Via per Arnesano, I-73100 Lecce, Italy  University of
Salento, Department of Physics, Via per Arnesano, I-73100 Lecce, Italy
University of Milano Bicocca, Department of Physics, Piazza della Scienza ,
I-20126 Milano, Italy  INFN - Milano Bicocca, Piazza della Scienza, I-20126
Milano, Italy  University of Salerno, Department of Physics, Via Giovanni
Paolo II, I-84084 Fisciano, Salerno, Italy  INFN - Napoli, Salerno group, Via
Giovanni Paolo II, I-84084 Fisciano, Salerno, Italy  University of Sannio,
Department of Science and Technology, via Francesco de Sanctis, I-82100,
Benevento, Italy  University of Sannio, Department of Engineering, Corso
Garibaldi, I-82100 Benevento, Italy  INO-CNR BEC Center, Via Sommarive,
I-38123 Povo, Italy  University of Trento, Department of Physics, Via
Sommarive, I-38123, Povo, Trento, Italy  Fondazione Bruno Kessler, Via
Sommarive, I-38123, Povo, Trento, Italy  INFN - TIFPA, Via Sommarive, I-38123,
Povo, Trento, Italy  IFN-CNR , Via Sommarive, I-38123 Povo, Trento, Italy
INRiM - Istituto Nazionale di Ricerca Metrologica, Strada delle Cacce, I-10135
Turin, Italy  Politecnico di Torino, Corso Duca degli Abruzzi, I-10129 Turin,
Italy  INFN - Torino, Via Pietro Giuria, I-10125 Turin, Italy

###### Abstract

DARTWARS (Detector Array Readout with Traveling Wave AmpliﬁeRS) is a three
years project that aims to develop high-performing innovative Traveling Wave
Parametric Ampliﬁers (TWPAs) for low temperature detectors and qubit readout
(C-band). The practical development follows two diﬀerent promising approaches,
one based on the Josephson junctions (TWJPA) and the other one based on the
kinetic inductance of a high-resistivity superconductor (KITWPA). This paper
presents the advancements made by the DARTWARS collaboration to produce a
first working prototype of a KITWPA.

##  1 Introduction

Nowadays, many applications rely on the faithful detection of low power
microwave signals at cryogenic temperatures. This is especially true in the
field of the superconducting quantum computation, where quantum-limited noise
microwave amplification is paramount to infer the qubit state with high
fidelity.

For these applications, the goals are to maximize the signal to noise ratio of
extremely feeble microwave signals while allowing a broad readout bandwidth.
The latter is extremely relevant in all the applications where the devices are
required to be multiplexed over large bandwidths. To solve this problem,
parametric amplification with superconducting circuits, a well known technique
used for low noise amplifiers, will be exploited and developed to its
technical limits.

DARTWARS [1] (Detector Array Readout with Traveling Wave AmplifieRS) aims to
develop two types of traveling wave parametric amplifiers: TWJPA and KITWPA.
The technical goal is to achieve a gain value around 20 dB, such that the
added noise of a commercial semiconductors low temperature amplifiers (HEMT)
could be neglected, with a high saturation power (around -50 dBm), and a
quantum limited or nearly quantum limited noise ($T_{N}<600$ mK) in C-band.
These features will lead to the readout of large arrays of detectors or qubits
with virtually no noise degradation.

##  2 Parametric amplification and KITWPA

A traveling wave parametric amplifier is the physical implementation of a
parametric oscillator. By carefully modulating the properties of the
oscillating system, the oscillator will absorb energy from the surroundings,
resulting in an exponential growth of the oscillation amplitudes. A parametric
amplification of a current signal $I$ can thus be achieved by making the
signal travel through a transmission line (TL) made of a non-linear inductors
$L(I)$, resulting in a wave equation which is equal to the one of the
parametric oscillator. In particular, KITWPA exploits the nonlinear kinetic
inductance of a superconducting coplanar waveguide (CPW) to generate
parametric interaction between a strong RF input pump at a frequency
$\omega_{p}$ and signal at a frequency $\omega_{s}$ ($\omega_{s}<\omega_{p}$)
, resulting in a signal amplification and generation of idler(s) product(s)
$\omega_{i}$ [2].

The inductance of the superconducting line comprises a geometric inductance
$L_{g}$ and a kinetic inductance $L_{k}$. The former depends on the geometry
of the circuit, while the latter depends on the supercunducting material and
it is related to the inertia of the Cooper pairs. The kinetic inductance makes
the overall inductance of the line non linear, and along the direction $x$ of
the CPW can be modeled as

$L(x,t)=L_{0}(x)\Bigg{[}1+\Big{(}\frac{I(x,t)}{I^{*}}\Big{)}^{2}\Bigg{]}$ (1)

where $L_{0}(x)$ is the linear kinetic inductance and $I^{*}$ is a constant
scaling factor.

The fact that a superconducting line has zero $dc$ resistance at $T<T_{c}$
suggests the motivation for the use of superconducting parametric amplifiers:
a very low power dissipation and an added noise that approaches the minimum
set by quantum mechanics.

One of the goal of DARTWARS is to develop a KITWPA amplifiers that operate in
a three-wave mixing (3WM) fashion when biased with a dc current $I_{d}$. In
this scheme, the pump exchange its energy with a signal and idler products
when the following phase-matching conditions are satisfied

$\omega_{p}=\omega_{s}+\omega_{i}\\\ ;\\\ k_{p}-k_{s}-k_{i}=\frac{\chi
I_{p0}^{2}}{8}(k_{p}–2k_{s}-2k_{i})$ (2)

where ${k_{p},k_{s},k_{i}}$ are the wavenumbers of the corresponding signals
and $\chi=1/(I^{*2}+I_{d}^{2})$ and $I_{p0}$ is the pump amplitude at the
amplifier input.

##  3 Preliminary results from the KITWPA fabrication

To achieve the phase matching condition of eq (2), DARTWARS will exploit a
lumped-element transmission line for the KITWPA design. In this configuration,
the characteristic impedance of the TL is periodically loaded to implement
both a wide stopband at $3\omega_{p}$ and a narrow stopband near $\omega_{p}$.
Placing the pump tone below this narrow stopband allows to easily match the
conditions of eq (2) while suppressing higher pump harmonics, which would
result in a gain and noise degradation. In addition, the lumped elements TL
allows to achieve a high gain with a shorter TL (tens of cm compared to few
meters of a periodic loading TL).

The material of choice to fabricate the KITWPAs is NbTiN and the films are
being produced at the FBK facility. The critical temperature of the film
produced can be seen in Figure 1.

Figure 1: Normalized transition curves of the NbTiN films fabricated at FBK.
The different colors correspond to different film depositions.

Even if the measured critical temperatures suggest that new sputter targets
will be required to achieve the best literature values ($\sim 15$ K), the
results of the sputtering fabrication processes are encouraging, indicating
that a wide range of $L_{k}$ can be easily achieved with the current setup.

Figure 2: Comparison between the two different methods adopted to measure the
kinetic inductance of the supeconducting films.

To characterize the produced films, the NbTiN has been patterned into lumped
element resonators [3]. The kinetic inductance was assessed on three different
wafers with two techniques. The first technique, labeled as $L_{k}$-$Sim$,
estimates $L_{k}$ as the value required to match the measured resonant
frequency of the KID with the one extracted from simulation, which takes into
account only $L_{g}$. The second approach, indicated as $L_{k}$-$T_{c}$,
computes $L_{k}$ from

$L_{k}=\frac{\hbar}{1.76\pi k_{b}T_{c}}R_{n}$ (3)

The two techniques are comparable, as shown in Figure 2.

Figure 3: Example of different resonant profiles due to the different powers
of the RF signal.

Finally, the maximum achievable non-linearity of the film
$\frac{I_{c}}{I^{*}}$ was assessed from eq (1). The resonance frequency
$f_{0}$ was measured by changing the power $P$ of the RF probe tones (Fig 3),
while the conversion between $P$ and $I^{2}$ was assessed through Sonnet and
Qucs simulations. The following relation was used to evaluate $L_{k}$: $(2\pi
f_{0})^{-2}=(L_{k}+L_{g})C$.

Figure 4: $I^{*}$ measurement of the resonator with the highest kinetic
inductance

Figure 4 shows the results achieved for one of the produced resonators.
Overall, the maximum achieved non-linearity from the different wafers produced
was roughly 0.25, close to the literature value of 0.34 [4]. Upon these
promising results, the first KI-TWPA prototypes are currently beeing produced.

## Acknowledgments

This work is supported by the Italian Institute of Nuclear Physics (INFN),
within the Technological and Interdisciplinary research commission (CSN5), by
the European Union’s H2020-MSCA Grant Agreement No. 101027746, and by the
Joint Research Project PARAWAVE of the European Metrology Programme for
Innovation and Research (EMPIR). PARAWAVE received funding from the EMPIR
programme co-financed by the Participating States and from the European
Union’s Horizon 2020 research and innovation programme.

## References

  * [1] A. Giachero, C. Barone, M. Borghesi, G. Carapella, A. Caricato, I. Carusotto, W. Chang, A. Cian, D. Di Gioacchino, E. Enrico, et al., Detector array readout with traveling wave amplifiers, arXiv preprint arXiv:2111.01512 (2021).
  * [2] R. P. Erickson, D. P. Pappas, Theory of multiwave mixing within the superconducting kinetic-inductance traveling-wave amplifier, Physical Review B 95 (10) (2017) 104506.
  * [3] P. K. Day, H. G. LeDuc, B. A. Mazin, A. Vayonakis, J. Zmuidzinas, A broadband superconducting detector suitable for use in large arrays, Nature 425 (6960) (2003) 817–821.
  * [4] M. Malnou, M. Vissers, J. Wheeler, J. Aumentado, J. Hubmayr, J. Ullom, J. Gao, Three-wave mixing kinetic inductance traveling-wave amplifier with near-quantum-limited noise performance, PRX Quantum 2 (1) (2021) 010302.