© The Authors 2026

1 Introduction

Ultrahot Jupiters (UHJs) are tidally locked gas giants with day-side gas temperatures that exceed 2000–3000 K, resulting in partially ionized atmospheres under intense stellar irradiation (e.g., Helling et al. 2019b). At sufficiently high gas temperatures, alkali metals and other species thermally ionize, leading to a weakly ionized plasma in the atmosphere (Lodders 2003; Rodríguez-Barrera et al. 2015; Helling et al. 2019b; Dietrich et al. 2022). Since the atmospheric gas behaves as a plasma on the dayside, it interacts electromagnetically with the planetary magnetic field. The resulting magnetic drag force has been suggested as a key mechanism for damping zonal jets and contributing to localized heating in these atmospheres (e.g., Perna et al. 2010a; Batygin et al. 2013). The role of magnetic fields in shaping the atmospheric circulation of UHJs remains a challenging problem and has significant implications for the (U)HJ’s climate dynamics and wind patterns.

A major uncertainty is the knowledge about the actual magnetic field strengths of (U)HJs. While scaling laws based on self-sustained dynamo simulations can predict fields of a few to hundreds of Gauss, these are very rough estimates due to limited knowledge of planetary interiors and dynamo processes (see, e.g., Christensen & Aubert 2006; Christensen et al. 2009; Yadav et al. 2013; Reiners & Christensen 2010; Kilmetis et al. 2024; Elias-López et al. 2025). Indirect observational estimates of hot Jupiter magnetic fields from star-planet interactions suggested magnetic field strengths of 20 G to 120 G (Scharf 2010; Vidotto et al. 2010, 2011; Cauley et al. 2019; Savel et al. 2024).

Several studies addressed magnetic drag in general circulation models (GCMs) using a universal drag timescale applied throughout the planetary atmosphere (e.g., Komacek & Showman 2016; Tan & Komacek 2019; Carone et al. 2020) or simplified diffusivity-based descriptions: Perna et al. (2010a,b) parametrized the Lorentz force acting on the plasma with a drag timescale and estimated a uniform Rayleigh friction at each pressure level, neglecting horizontal variations of the diffusivity in the atmosphere. The scalar drag timescale was applied as a Rayleigh friction term on the horizontal winds. This constant Rayleigh friction term is referred to as “uniform drag” throughout this work. These models do not capture the different conditions of the day- and nightside, with the nightside possibly being much less ionized and decoupled from the magnetic field and the dayside being partially coupled to the field.

Rauscher & Menou (2013) and Beltz et al. (2022b) improved this model by introducing a spatially varying but directionally isotropic drag coefficient, calculated from local thermodynamic conditions (e.g., gas temperature and density) and updated at each time step. Beltz et al. (2022b) referred to this as “active magnetic drag,” a term that is used throughout this work. In these models, the drag is implemented as a scalar damping term in the momentum equation and is applied only to the zonal (east–west) wind component. In several follow-up studies, Beltz and collaborators explored the dynamical and observational consequences of magnetic drag in (U)HJ atmospheres using the spatially varying drag description in their GCM (Beltz et al. 2022a, 2023; Beltz & Rauscher 2024; Beltz et al. 2025). These studies extensively investigate the effect of planetary magnetic fields on atmospheric circulation, phase curves, high-resolution emission, and transmission spectra for different (U)HJs, as well as the effects on HJs with eccentric orbits. While these models have provided valuable insight into the importance of magnetic field effects, they neglect the anisotropic response of the partially ionized gas to magnetic fields. Based on the work of Rauscher & Menou (2013) and Beltz et al. (2022b), Christie et al. (2025) derived a drag parametrization from an explicit calculation of the Lorentz force in a dipolar field geometry and showed that including anisotropic drag in the zonal, meridional, and vertical direction (due to Hall and Pedersen currents) modifies the flow in HJ atmospheres differently compared with models that only apply zonal drag. Furthermore, Rogers (2017) demonstrated the effects of magnetic fields on the winds of HAT-P-7 b using spherical 3D anelastic magnetohydrodynamic (MHD) code with the same spatially varying resistivity that is linked to temperature-dependent ionization as Rauscher & Menou (2013). Their simulations showed that wind variability in the atmosphere is linked to the magnetic field strength and coincides with the position of the hotspot in the atmosphere. Assuming that the observed variable winds on HAT-P-7 b are due to magnetism, a minimum magnetic field strength of 6 G can be constrained.

In the study presented here, the parametrized magnetic drag from the differences between the charged particle velocity and the neutral wind velocity is derived. The magnetic drag on the atmospheric gas bulk arises from collisions between the bulk neutral flow and the ionized component in the atmospheric gas under the influence of a magnetic field. While ions and electrons experience electromagnetic forces (Lorentz forces), neutrals do not, leading to momentum transfer that can act as a drag force on the flow. Both Pedersen and Hall drag components are included and applied as momentum sink terms in the horizontal wind equations. The corresponding frictional heating is computed self-consistently from the work done by the drag force on the neutral flow. This approach enables a more straightforward and physically motivated way to implement the magnetic drag within the GCM code ExoRad (Carone et al. 2020; Schneider et al. 2022b; Baeyens et al. 2024). This work adapts ionospheric MHD formulations to the exoplanet hot Jupiter regime, where weak ionization, rapid rotation, and strong irradiation combine to produce a unique plasma environment. The approach of Christie et al. (2025) and the one presented in this work are based on the same underlying non-ideal MHD framework for weakly ionized atmospheres and rely on Pedersen and Hall conductivities or drags to describe magnetic coupling in weakly ionized atmospheres. However, they differ in their formulation and implementation way within the GCM. Christie et al. (2025) describe the magnetic coupling as a drag on the neutral flow based on the Lorentz force, whereas here it is implemented via Pedersen and Hall drag terms derived from drift velocities between the charged particles and the neutrals. Furthermore, in both approaches, the energy equation for the neutral gas takes into account the corresponding feedback, where the kinetic energy dissipated by the drag terms is introduced as a local heating source. The two approaches lead to similar parametrizations of the anisotropic drag for the parameter space considered. The ultrahot regime of WASP-18 b explored in this work allows the Hall term to become significant, leading to flow asymmetries that were not seen in the parameter space considered by Christie et al. (2025).

The physical basis for separating magnetic drag into Pedersen and Hall components is well established in terrestrial and planetary ionospheric electrodynamics. The modeling of MHD and electrodynamic coupling between neutral atmospheres and ionospheric currents has been a central topic in the Solar System atmospheric modeling. The thermo-sphere–ionosphere–electrodynamics general circulation model (TIE-GCM; Richmond et al. 1992; Richmond & Thayer 2000) and related thermosphere-ionosphere GCMs (e.g., Rees 1989; Wang et al. 2004) solve the coupled neutral-ion dynamics and electrodynamics self-consistently and include ion drag and frictional heating in the momentum and energy equations. Additionally, Zhu et al. (2005) derived expressions for the ion drag and Joule heating for the neutral atmosphere in the thermosphere dependent on the magnetic field, the Pedersen and Hall conductivities, the ion cyclotron frequency, and the ion-neutral collision frequency. These approaches show how ionospheric currents create frictional forcing on neutral winds and localized heating in planetary upper atmospheres. Similar studies with coupled magnetosphere–ionosphere–thermosphere systems for Jupiter and Saturn demonstrate that field-aligned currents and ionospheric torques can reshape thermospheric winds and energetics, and that the spatial variations of conductivity and currents are important for the resulting drag distribution (Bougher et al. 2005; Yates et al. 2020; Müller-Wodarg et al. 2012). Neither the model presented here nor the model of Christie et al. (2025) explicitly solve the charged particle dynamics or couples the atmosphere to the magnetosphere, as was done for Earth and for Solar System giant planet GCMs. Nevertheless, these parameterizations are crucial, since they provide a computationally feasible way to capture the main effects of a planetary magnetic field on the wind structures and heat transport in the (U)HJ atmospheres, and allow for a direct comparison to observations. Such models are hence more advanced than simplified Rayleigh drag approaches and represent an important step toward coupled neutral–plasma GCMs for exoplanets.

The study and analysis of the magnetic drag and frictional heating is presented for the atmosphere of UHJ WASP-18 b. WASP-18 b is used as an example to present and discuss a new modeling approach that enables one to consider anisotropic magnetic drag and its resulting effect on the exoplanet climate state. The ultimate aim is to provide a model that allows one to derive the global magnetic field strength, which is otherwise inaccessible for extrasolar planets. WASP-18 b is an UHJ with an equilibrium temperature of 2400 K and a mass of ∼10 Jupiter masses. It orbits an F6-type star at a close distance of 0.02 AU (Hellier et al. 2009; Southworth 2010). Recent Hubble Space Telescope (HST) and Transiting Exoplanet Survey Satellite (TESS) phase-curve observations revealed an inefficient day–night heat circulation (≳96% of the heat remains on the dayside) and a small hotspot offset 3^◦–5^◦ eastward (Arcangeli et al. 2019; Shporer et al. 2019). Dayside thermal emission spectrum obtained with the Near-Infrared Imager and Slitless Spectrograph (NIRISS) on the James Webb Space Telescope (JWST) further exposed water emission features (Coulombe et al. 2023). These observed features suggest the presence of an effective magnetic drag and super-solar metallicity (Deline et al. 2025; Arcangeli et al. 2019) and makes it a suitable laboratory for studying magnetic effects in exoplanet atmospheres. The approach presented here offers a more accurate parametrization of the magnetic drag and heating, and opens the way for more accurate predictions of circulation patterns, temperature distributions, and potentially observable signatures such as phase curves and wind speeds. Magnetic drag on UHJs not only affects winds and temperatures but also where and how clouds may form, and can be observed in these extreme atmospheres (see, e.g., Helling et al. 2021; Kennedy et al. 2025). This study investigates the following research questions:

• To what extent disrupts magnetic drag the equatorial super-rotation of gas giants? Can it change the overall flow pattern?

• How does magnetic drag modify the heat redistribution between the day- and nightsides?

• At which gas pressures does magnetic drag become dynamically important for the atmosphere?

The results show that the equatorial jet is weakened by a magnetically coupled atmosphere (independent of the details of magnetic drag treatments), anisotropic drag that traces the locally changing ionization in particular may lead to a partial decoupling of the day-night flow on the morning terminator in the upper atmosphere.

The paper is organized as follows. Section 2 summarizes the approach that addresses the research questions, and gives an overview of the GCM code ExoRad, which physics is implemented, and the physical and numerical parameters (including boundary layers and magnetic drag implementation) used for the WASP-18 b simulations. Section 4 presents the analysis of the plasma parameters (including the degree of ionization, magnetic Reynolds number) and the detailed derivation of the anisotropic drag parametrization, which is implemented in ExoRad. In Sects. 5 and 6, the results of the simulations with different drag treatments are shown and their influence on the wind pattern and gas temperature distribution are compared, and the limitations of the model are discussed. The effect on the climate parameters such as the day-night gas temperature difference, terminator temperature difference, jet spread, and width are explored. In Sect. 7, a summary of this study is presented. In the appendix, an overview is given of the active drag treatment based on the models of Perna et al. (2010a) and Rauscher & Menou (2013) and the drag treatment of Christie et al. (2025).

2 Approach

This study investigates how magnetic drag affects the atmospheric dynamics of UHJs that show extreme day and night differences, exploring the example of WASP-18 b. WASP-18 b was chosen for this study because it is one of the hottest close-in extrasolar gas giants (Parmentier et al. 2018) with an extreme irradiation environment. The high dayside temperatures reaching a maximum brightness temperature of about 3000 K between 1 bar and 10⁻² bar (Coulombe et al. 2023) lead to high thermal ionization of alkali metals and other species (including Al, Fe, and Ti; see Helling et al. 2019a), resulting in noticeable electrical conductivities, and thus strong magnetic coupling in the atmosphere. Furthermore, its high mass and short orbital period hint at the presence of a strong planetary dipole magnetic field (Elias-López et al. 2025), making it a useful test case for studying magnetic drag effects on atmospheric circulation. Our approach to answer our research question was divided into four steps: (1) The atmospheric thermal structure of WASP-18 b obtained from ExoRad simulations was analyzed to evaluate where magnetic effects are expected to be significant. Therefore, the ionization fraction (Eq. (23)) was calculated from equilibrium chemistry using ExoRad gas temperature and gas pressure profiles, and it was applied to evaluate magnetic coupling through plasma parameters such as the magnetic Reynolds number (Eq. (22)) and the electron plasma frequency (Eq. (19)). Furthermore, the analysis allowed us to identify where nonideal MHD effects such as the Hall effect might become important (Sect. 6). (2) From the three-fluid description of the atmospheric gas (Eqs. (24)–(26)), including electrons, ions, and neutral particles coupled through neutral-plasma collisions (Eqs. (41) and (42)), a parameterization of the magnetic drag in a partially ionized atmosphere was derived that can be used in a GCM. The parameterization separates between Pedersen (dissipative) and Hall (non-dissipative) components, leading to an anisotropic drag (Eq. (36)). The derivation allows one to couple local plasma properties (ionization fraction, collision frequencies, and magnetic field strength) to momentum exchange with the neutral atmosphere. (3) The anisotropic drag parameterization was implemented in ExoRad. For comparison, other known magnetic drag models are summarized (Appendices A and B) and runs with no drag and with simplified drag parameterizations (“uniform” and “active drag”) were performed to analyze the consequences of different magnetic drag treatments on the atmospheric dynamics, using WASP-18 b as example. (4) Finally, the atmospheric dynamics, wind speeds, and gas temperatures across the different drag treatments were compared. These results were analyzed further through parameters that characterize the planet’s climate state as introduced by Plaschzug et al. (2026): the day–night and evening-morning terminator temperature difference, wind jet speed, and width. Additionally, a scale analysis of Ohm’s law was performed to evaluate the relative importance of Pedersen and Hall currents and ambipolar diffusion at different gas pressure levels. This is important to determine whether the drag is isotropic or directionally dependent, and at what altitudes magnetic forces can modify the large-scale circulation. In weakly ionized atmospheres, ions and electrons drift relative to the neutrals because only the charged species are directly affected by the magnetic field. This mechanism is usually referred to as ambipolar diffusion in astrophysics, and it allows the magnetic field to slip through neutral gas via collisions between charged and neutral particles, which transfer momentum between the plasma and the neutrals. The collisional momentum exchange can result in frictional heating and may be important for energy dissipation (e.g., Khomenko & Collados 2012; Hillier 2024). Identifying which magnetic process dominates is important for understanding the physical regime of the atmosphere and the effect of magnetic forces on the atmospheric flow.

3 3D atmosphere model: ExoRad

To model the 3D atmospheric hydro- and thermodynamics of WASP-18 b under different magnetic drag treatments, the 3D GCM ExoRad (Carone et al. 2020; Schneider et al. 2022b) was employed. It uses the hydrodynamical core MITgcm (Adcroft et al. 2004) to solve the Navier-Stokes equations for a hydrostatic atmosphere on a rotating sphere, complementing it with physical parametrizations suitable for a highly irradiated, tidally locked gas giant.

Model setup. An outline of the MITgcm used here is in place, given that various papers are applying the same core with their individual developments (see, e.g., Steinrueck et al. 2025). The MITgcm model applied here solves the hydrostatic primitive equations (HPEs; Showman et al. 2009) using the finite-volume method to discretize the momentum equations in space and the second-order Adams-Bashforth time integration scheme (Lilly 1965) for explicit time stepping on a rotating sphere in a C32 cubed-sphere grid, corresponding to a horizontal resolution of 128 × 64 in longitude (ϕ) and latitude (θ; 2.8^◦ × 2.8^◦). The dynamical core uses a staggered Arakawa C grid (Arakawa & Lamb 1977; Showman et al. 2009). The equation of state is the ideal gas law. The vertical coordinate is defined as the gas pressure, p_gas (computed in Pa). However, for clarity, all pressure values discussed in the text and shown in figures are expressed in bar. The modeled atmosphere extends from 10⁻⁵ bar to 700 bar. Pressure (p) coordinates are preferred in atmospheric modeling because they align better with hydrostatic balance. The full set of equations (horizontal momentum, vertical momentum under hydrostatic equilibrium, mass continuity, and thermodynamic energy equations) in p coordinates, which are solved in ExoRad, are

$$\frac{d{\underset{\_}{v}}_h}{dt}=-f\widehat{\underset{\_}{k}}\times {\underset{\_}{v}}_h-{\nabla }_{p_{\text{gas}}}\Phi +{\underset{\_}{F}}_v,$$(1)

$$\frac{\partial \Phi }{\partial p_{\text{gas}}}=-\frac{1}{{\rho }_n},$$(2)

$${\nabla }_{p_{\text{gas}}}·{\underset{\_}{v}}_h+\frac{\partial \omega }{\partial p_{\text{gas}}}=0,$$(3)

$$\frac{d\Theta }{dt}=\frac{\Theta }{T_{\text{gas}}}\frac{q}{c_p},$$(4)

with the time, t [s], the horizontal velocity component of the neutral gas on pressure surfaces, v_{h =} (u, v, 0), the zonal velocity, u [m/s], the meridional velocity, v [m/s], the gas temperature, T_gas [K], the mass density of the neutral gas, ρ_n [kg m⁻³], the horizontal gradient, ${\nabla }_{p_{\text{gas }}}$ evaluated on constant pressure surfaces, the vertical velocity in pressure coordinates, $\omega =\frac{d{p}_{\text{gas}}}{dt}[\text{Pa}/\text{s}]$, the geopotential on constant-pressure surfaces, Φ [m²s⁻²], the Coriolis parameter, f = 2Ω sin ϕ [s⁻¹], the planetary rotation rate, Ω [rad s⁻¹], the local vertical unit, $\underset{\_}{\widehat{k}}$, and the external sink term in the horizontal momentum, F_v, which comprises Rayleigh damping (see Eqs. (7) and (9)) and frictional drag (Eq. (36)). In the MITgcm, the gas temperature is represented by the potential temperature, Θ [K], which is a measure of entropy¹ (Showman et al. 2009):

$$\Theta =T_{\text{gas}}{\left(\frac{p_0}{p_{\text{gas}}}\right)}^{R_d/c_p},$$(5)

where R_d [J kg⁻¹K⁻¹] is the specific gas constant, c_p [J kg⁻¹K⁻¹] is the specific heat at constant pressure, and p₀ [Pa] the reference pressure (p₀ = 7 × 10⁷ Pa). The thermodynamic energy equation includes the total heating rate per unit mass, q = q_rad + q_fric,deep + q_fric [W kg⁻¹]. $q_{\text{rad}}=g\frac{\partial F^{\text{net}}}{\partial p_{\text{gas}}}\left[{\text{Wkg}}^{-1}\right]$, where g [ms⁻²] is gravity and F^net [Wm⁻²] the total bolometric flux (Eq. (18)), represents radiative heating and cooling, calculated from the stellar irradiation by solving the radiative transfer equation. q_fric,deep [W kg⁻¹] is the frictional heating due to deep friction (Eq. (11)). q_fric [W kg⁻¹] accounts for frictional heating resulting from the dissipation of kinetic energy into thermal energy (Eq. (43) for anisotropic drag, Eq. (A.16) for active drag, and Eq. (6) for the uniform drag approach). The total (or Lagrangian) derivative is here defined by $\frac{d}{dt}=\frac{\partial }{\partial t}+{\underset{\_}{v}}_h\cdot {\nabla }_{p_{\text{gas}}}+\omega \frac{\partial }{\partial p_{\text{gas}}}$.

The conservation of radiative energy – the balance between incoming stellar flux and outgoing planetary flux – was tested for ExoRad in Schneider et al. (2022a) for the UHJ WASP-76 b assuming no interior heat flux (T_int = 0 K). The model was found to preserve radiative energy conservation to 99.9% even for very long simulation times (86 000 days). Kinetic energy and angular momentum was found to be conserved within 0.5% accuracy for a simulation time of 2000 days (Carone et al. 2020).

Different drag treatment configurations. The effect of magnetic drag was explored and different parameterizations were compared: “no drag,” uniform drag, active drag, and anisotropic drag. The no drag model represents a purely hydrodynamic atmosphere, where magnetic effects are neglected and the flow evolves freely under the balance between radiation, pressure gradients, and planetary rotation. In the uniform drag model, a constant drag timescale of τ_drag,uniform = 10⁴ s is applied to both the horizontal momentum equation via F_v=−v_h/τ_drag,uniform and energy equation with the heating rate per unit mass given by

$$q_{\text{fric,uniform }}=\left(v_{\varphi }^2+v_{\theta }^2\right)/{\tau }_{\text{drag,uniform}}$$(6)

(see e.g., Tan & Komacek 2019). This simplification represents a spatially uniform magnetic coupling between the atmospheric flow and the planetary interior. Such uniform drag approaches are commonly used in HJ atmospheric circulation studies as a phenomenological representation of magnetic effects, often without specifying the underlying physical damping mechanism. Using the inverted form of Eq. (A.15), we can estimate the effective local field strength that is required to achieve the prescribed τ_drag,uniform = 10⁴ s, given the local gas density, ρ_n, and the magnetic diffusivity, η: $B_{\text{effective}}\sim \sqrt{\frac{{\mu }_0{\rho }_n\eta }{{\tau }_{\text{drag, uniform}}}}$. Since ρ_n and η vary with p_gas, ϕ, and θ, B_effective is likewise spatially variable. For example, B_effective = 0.1 G is obtained at p_gas = 10⁻³ bar and B_effective = 19 G at p_gas = 1 bar (for ϕ = 0^◦ and θ = 0^◦). Physically, the fixed uniform drag acts as a steady resistivity that slows down the horizontal flow everywhere, on both the hot, partially ionized dayside and the cooler, weakly ionized nightside. In reality, however, magnetic effects should vary strongly with temperature and ionization. In contrast, the active and anisotropic drag models take into account a more realistic, spatially dependent treatment of magnetic coupling. The drag strength depends on the local ionization fraction, x_e (Fig. 4, right), and magnetic field geometry, representing how magnetic forces interact more efficiently in hotter, more ionized regions (Sect. 4 and Appendix A). The active drag was derived from the assumption that the weakly ionized atmospheric winds experience a bulk Lorentz force drag (with the detailed formulation given in Appendix A) and this effect is parameterized with a drag timescale that is described here by ${\tau }_{\text{drag,active}}=\frac{{\mu }_0{\rho }_n\eta }{B_{\text{dip,}r}^2}$ with the radial component of the dipolar magnetic field, B_dip,r (Eq. (29)). The external sink term in the zonal momentum is implemented in ExoRad as F_v,ϕ = −u/τ_drag,active and the frictional heating as q_fric,active = u²/τ_drag,active. In the active drag case, the magnetic coupling acts only on the zonal winds, under the assumption that these dominate the large-scale circulation. However, because meridional winds can be comparable in magnitude, this simplification may miss important effects. The potential limitation of the active drag was discussed in Rauscher & Menou (2013) and was investigated further in Christie et al. (2025). The anisotropic drag model, which is described in Sect. 4.2, includes drag in zonal and meridional directions, creating a more complete representation of how magnetic drag modifies the flow across field lines and determines the global wind pattern. In the anisotropic drag approach, the drag is decomposed into Pedersen and Hall components, resulting in the drag force that both damps and redirects the flow relative to the magnetic field geometry

Numerical parameters. The vertical space covers 47 vertical cells. A logarithmic spacing with 41 grid cells was applied between 10⁻⁵ bar and 100 bar to resolve the radiative and dynamically active regions, while linear spacing with six grid cells was used between 100 bar and 700 bar with a step size of 100 bar. The vertical resolution of the model resolves the atmosphere with about two vertical grid layers per scale height, which is important to accurately capture vertical dynamics. The dynamical time step ∆t =25 s was chosen in the model to satisfy both the Courant–Friedrichs–Lewy (CFL) condition (set by horizontal advection and grid spacing) and the shortest physical timescales introduced by parameterized source terms. Furthermore, a fourth-order Shapiro filter (Shapiro 1970) with a dampening timescale of 25 s was applied that horizontally smooths grid-scale noise (Showman et al. 2009; Carone et al. 2020). Although the application of smoothing methods is routinely used in GCMs, they still might have a non-negligible effect on the atmospheric dynamics (e.g., Heng et al. 2011; Polichtchouk et al. 2014; Skinner & Cho 2021). The simulations ran for a total of 1000 planetary days to allow the temperature structure to stabilize and the system to approach a steady-state solution.

Boundary conditions. The top boundary of the computational domain acts as a solid boundary (ω = 0 Pa/s), while the bottom boundary is an impermeable surface placed at p_bottom = p₀ = 700 bar located below the region of interest. At both the top and bottom boundaries, a free slip boundary condition was applied to the horizontal velocity. To dampen nonphysical gravity wave reflection at the top of the computational volume, a sponge layer (ghost cells) was applied between p_gas = 10⁻⁴ bar and p_gas = 10⁻⁵ bar (Carone et al. 2020; Beltz et al. 2022b; Schneider et al. 2022b). Within this layer, Rayleigh friction was applied to the zonal velocity, u, which acted to bring u closer to its longitudinal mean ū [m/s] via

$$\frac{du}{dt}=-\tilde{k}(u-\overline{u}),$$(7)

where the strength parameter $\tilde{k}\left[{\text{s}}^{-1}\right]$ is a function of p_gas:

$$\tilde{k}\left(p_{\text{gas }}\right)=k_{\text{top }}\max {\left[0,1-\frac{p_{\text{gas }}^2}{p_{\text{sponge }}^2}\right]}^2.$$(8)

The parameters k_top [s⁻¹] and p_sponge [Pa] determine the intensity and the location of the Rayleigh friction in the sponge layer. The default values in this model are k_top = 1728 × 10³ s⁻¹ and p_sponge = 10 Pa (Deline et al. 2025). Since the domain uses a cubed sphere grid, Eq. (7) was computed by first converting the cubed sphere grid values to geographic coordinates, averaging the deprojected u within 20 latitude bands, and finally mapping the resulting ū back onto the cubed sphere grid (Schneider et al. 2022b). This upper sponge layer isolates potentially nonphysical solutions caused by unresolved boundary effects from the interior test volume (see, e.g., also Sect. 2.3. in Helling et al. 2004). Consequently, the numerical solution in the region p_gas < 10⁻⁴ bar was excluded from the analysis of the atmospheric dynamics presented in Sect. 5. Similarly, a deep-layer sponge drag was applied at the bottom boundary to ensure numerical stability, following Carone et al. (2020). With this method shear flow instabilities and nonphysical changes in the flow pattern of the simulation domain are avoided without directly affecting the observable atmosphere because the drag only modifies the winds in the deepest layers (p_gas > 490 bar). This deep-layer sponge drag dissipates the horizontal wind, v_h, via

$$\frac{d{\underset{\_}{v}}_h}{dt}=-k_{deep}{\underset{\_}{v}}_h.$$(9)

The parameter k_deep [s⁻¹] is defined by

$$k_{\text{deep}}=k_{\text{bottom}}\max \left[0,\frac{p_{\text{gas}}-4.9\times {10}^7 \text{Pa}}{p_0-4.9\times {10}^7 \text{Pa}}\right],$$(10)

with the control parameter k_bottom = 86 400 s⁻¹, which dissipates fast wind jets at depth (Carone et al. 2020). The kinetic energy dissipated at the bottom via friction is converted into thermal energy, approximating Ohmic dissipation (Rauscher & Menou 2013; Carone et al. 2020). The corresponding heating rate per unit mass of the neutral gas is given by

$$q_{\text{fric,deep }}=k_{\text{deep }}{\underset{\_}{v}}_h^2.$$(11)

Simulation tests for the sponge layer are discussed in detail in Carone et al. (2020).

Atmospheric initialization. The atmosphere was initialized using the analytic temperature profile of Parmentier et al. (2015). Therefore, the interior and equilibrium temperatures were used. The equilibrium temperature for zero albedo,

$$T_{\text{eq},0}=T_{\star }[\text{K}]\sqrt{\frac{R_{\star }[\text{m}]}{2a_p[ \text{m}]}},$$(12)

was computed from stellar and orbital parameters (see Table 1). For WASP-18 b, adopting a stellar effective temperature of T_⋆ = 6400 K, stellar radius of R_⋆ = 1.23 R_⊙, and semimajor axis of a_p = 0.0203 AU yields T_eq,0 ≈ 2402 K. T_eq,0 was then applied to estimate the corresponding interior temperature of T_int ≈ 573 K following Thorngren et al. (2019). These temperatures define the initial radiative–convective profile using the analytical model of Parmentier et al. (2015), which yields a potential temperature of the hot adiabat of Θ_ad = 6940 K at 1 bar. Synchronous rotation and a surface gravity of g = 190 m s⁻² are assumed. This setup provides a hot, physically consistent initial state in the deep layers. Sainsbury-Martinez et al. (2019) have shown that the time required for the deep atmosphere to reach equilibrium after cooling from a hot initial condition is less than the time required for heating it from a cold initial state to a hotter one.

Radiative transfer treatment. The ExoRad framework uses the full radiative transfer treatment expeRT/MITgcm (Carone et al. 2020; Schneider et al. 2022b) based on petitRADTRANS (Mollière et al. 2019, 2020). The implementation of radiative transfer in ExoRad is described in detail in Sect. 2.2 of Schneider et al. (2022b). The model solves the radiative transfer equation

$$\underset{\_}{n}·\nabla I_v={\alpha }_v^{\text{tot }}\left(S_v-I_v\right),$$(13)

with the unit vector, n, the intensity, I_ν [Wm⁻²sr⁻¹Hz⁻¹], the frequency, ν [Hz], the source function, S_ν [W m⁻² sr⁻¹ Hz⁻¹], and the inverse mean-free path of the light beam, ${\alpha }_v^{\text{tot}}\left[{\text{m}}^{-1}\right]$. Equation (13) describes how a beam of radiation traveling through a planetary atmosphere loses energy through absorption, gains energy through emission, and redistributes energy through scattering. The equation was solved for photons originating in the planetary atmosphere or scattered from incoming stellar radiation. Direct stellar intensity extinction was modeled separately using exponential decay, and the total intensity was obtained by summing all contributions (intensities of planetary and scattered stellar photons, and attenuated stellar intensity). The mean stellar attenuated intensity in a plane-parallel atmosphere is

$$J_v^{\star }=\frac{I_v^{\star }\left(p_{\text{gas }}=p_{\text{top }}\right)}{4\pi }\exp \left(\frac{-{\tau }_v}{{\mu }_{\star }}\right),$$(14)

with the stellar intensity at the top of the atmosphere $I_v^{\ast }(p_{\text{gas}}=p_{\text{top}})$, calculated using the PHOENIX stellar model spectrum (Husser et al. 2013) via petitRADTRANS. For a tidally locked planet, the cosine of the zenith angle – that is, the angle between the normal vector on top of a given atmospheric column on the planet and the incoming stellar light – is µ_⋆ = cos θ cos ϕ, with latitude θ and longitude ϕ, when the origin of the coordinate system is set to the substellar point. The optical depth is

$${\tau }_v={\displaystyle \int \frac{{\kappa }_v^{\text{tot}}}{g}}d{p}_{\text{gas,}}$$(15)

where ${\kappa }_v^{\text{tot}}={\alpha }_v^{\text{tot}}{\rho }_n\left[{\text{m}}^2{\text{kg}}^{-1}\right]$ is the total gas-phase opacity given from absorption and scattering coefficients. The received bolometric stellar flux was computed via

$$F^{\star }=4\pi {\mu }_{\star }{\displaystyle {\int }_v{J}_v^{\star }}dv\left[ \text{W} {\text{m}}^{-2}\right].$$(16)

The radiative transfer equation (Eq. (13)) was solved using the Feautrier method (Feautrier 1964) with accelerated lambda iteration (Olson et al. 1986) to obtain the planetary flux and iterated to convergence of the planetary intensity field. Convergence is achieved when the relative change in the source function falls below 2% of the reciprocal local lambda operator, providing a good balance between accuracy and performance (Schneider et al. 2022b). The previous radiative time step’s source function serves as the initial guess, allowing later time steps to converge in only a few iterations (fewer than three), compared to 50–100 iterations in the beginning (the first 100 time steps).

Then the bolometric flux, F^pla [Wm⁻²], was computed by integrating the radiation field, $F_v^{\text{pla}}$:

$$F^{\text{pla}}={\displaystyle {\int }_v{F}_v^{\text{pla}}}dv=-4\pi {\displaystyle {\int }_v{H}_v^{\text{pla}}}dv.$$(17)

$H_v^{\text{pla}}$ is the directional flux of radiation at frequency ν emerging from the planet. The total bolometric flux, F^net [Wm⁻²], was then obtained as

$$F^{\text{net }}=F^{\text{pla }}+F^{\star }.$$(18)

Following Showman et al. (2009), F^net was evaluated on vertically staggered cell interfaces, using quadratic Bézier interpolation for temperature profiles (Lee et al. 2022) to enhance accuracy and stability in particular at the terminators, where the zenith angle of incoming radiation is approaching zero. F^net was computed for every second grid column, with linear interpolation between columns and direct computation at tile borders due to the C32 grid geometry.

The total gas-phase opacity, ${\kappa }_v^{\text{tot}}$, is represented by correlated-k tabulated line gas opacities² resolved into 11 wavelength bins from 0.26 µm to 300 µm (Kataria et al. 2013; Schneider et al. 2022b) corresponding to the S1 resolution for the following species: H₂O (Barber et al. 2006), Na, and K including pressure broadening (Allard et al. 2019), CO₂ (Yurchenko et al. 2020), CH₄ (Yurchenko et al. 2017), NH₃ (Coles et al. 2019), CO (Li et al. 2015), H₂S (Azzam et al. 2016), HCN (Barber et al. 2014), SiO (Barton et al. 2013), PH₃ (Khomenko et al. 2014), FeH (Wende et al. 2010), TiO (McKemmish et al. 2019), and VO (McKemmish et al. 2016), as well as H⁻ scattering suitable for an ionized atmosphere (see Table 1 in Schneider et al. 2022b). The gas line opacities were pre-calculated on a grid of gas pressure and temperature (T_gas = 100 K … 4000 K with linear spacing ∆T = 3.9 K and p_gas = 700 bar … 10⁻⁵ bar as set by the vertical atmospheric pressure grid), assuming local chemical equilibrium using easyCHEM (Lei & Mollière 2025). Furthermore, the opacities also include collision-induced absorption from H₂-H₂ (Borysow et al. 2001; Borysow 2002; Richard et al. 2012) and H₂-He (Borysow et al. 1988; Borysow & Frommhold 1989; Borysow et al. 1989), Rayleigh scattering by H₂ (Dalgarno & Williams 1962) and He (Chan & Dalgarno 1965), and opacities for H⁻ free-free and bound-free photon absorption in the presence of free electrons (Gray 2008; Jacobs et al. 2022). Atmospheric chemistry assumes solar metallicity ([Fe/H] = 0) and a solar C/O ratio of 0.55 under local thermodynamic equilibrium (LTE) conditions (Asplund et al. 2009).

As the upper (p_gas < 1 bar) atmospheric temperatures can evolve significantly, radiative fluxes were recalculated every fourth dynamical time step; that is, a radiative time step of τ_rad = 100 s was used.

Planetary system properties. For the equatorial magnetic field strength at the reference radius R_ref [m], B₀ = 5 G was chosen (Coulombe et al. 2023). Coulombe et al. (2023) showed that an internal magnetic field of at least 5 G is required within the framework of a GCM with active magnetic drag (Beltz et al. 2022b) to reproduce the observed white light curve, as such a field effectively suppresses any noticeable longitudinal shift of the hotspot from the substellar point. However, this threshold is likely model-dependent and may vary with the adopted drag formulation and atmospheric conductivity profile.

For R_ref, the transit radius of R_p,obs = 1.19 R_Jup (Table 1) was adopted, which was derived from optical light curves corresponding to the altitude where the slant optical depth reaches unity (e.g., Fortney 2005). This contrasts with Jupiter, where the radius is conventionally defined at p_gas = 1 bar (Seidelmann et al. 2007). The optical transit radius is wavelength-dependent and sensitive to clouds and hazes, affecting atmospheric opacities and possibly leading to an overestimated radius (see, e.g., Burrows et al. 2007). Unlike at Jupiter, there is no fixed physical pressure level for the radius across different exoplanets due to their large diversity. However, the transit radius is the directly observable radius and is consistent with the observed transit depth and the derived planetary mass–radius relation. Using representative parameters for WASP-18 b (R_d = 4343 J kg⁻¹ K⁻¹, g = 190 m s⁻²; see Table 1) and the T_gas(p_gas) profile at ϕ = 0^◦ and θ = 0^◦ from ExoRad, the geometric height difference, ∆z, can be calculated by integrating the hydrostatic equation $dz=-\frac{R_d{T}_{\text{gas}}(p_{\text{gas}})}{g}d\ln p_{\text{gas}}$ between p_obs = 10⁻³ bar and p_target = 1 bar. Mapping R_p,obs at p_obs = 10⁻³ bar to p_target = 1 bar gives ∆z ≈ 467 km, which results in less than 0.6% of R_p,obs. The overview of the simulation parameters and the setup of the different drag treatments in the performed simulations is given in Tables 1 and 2.

Numerical stability and drag implementations. In all simulations with active or anisotropic magnetic drag, it was found that numerical stability imposes a lower limit on the drag timescale, τ_drag [s]. The drag timescale is defined by τ_drag = min(ρ_n/K_H , ρ_n/K_P), where ρ_n [kg m⁻³] is the gas density, and K_p [kg m⁻³ s⁻¹] and K_H[kg m⁻³ s⁻¹] are the total Pedersen and Hall drag coefficients (Eqs. (39), (40)). Specifically, when τ_drag < 100 s, the model becomes unstable and crashes due to the strong damping in the low-pressure upper atmosphere. Such short timescales are often reached on the dayside at p_gas < 10⁻² bar, where ionization is high and diffusivity is low (gray line in Fig. 3). To maintain numerical stability across the simulation runs, a global lower bound of τ_drag ≥ 100 s is imposed, which corresponds to four times the model time step (∆t = 25 s). Similar lower bounds of τ_drag are commonly used in GCMs with active drag (e.g., Rauscher & Menou 2013; Coulombe et al. 2023).

4 Modeling the magnetic drag in the changing ionization environments of exoplanet atmospheres

The differing stellar irradiation between the dayside and night-side of WASP-18 b introduces a strong horizontal pressure gradient in the upper atmosphere. This gradient drives fast zonal winds from the hot to the cold hemispheres, with speeds reaching several kilometers per second. High temperatures in the day-side atmosphere lead to thermal ionization of alkali metals (e.g., Na, K) and other species (e.g., Al, Ti; Helling et al. 2019a), which create a partially ionized plasma. Figure 1b illustrates the thermal ionization fraction, x_e, for WASP-18 b based on the 3D ExoRad results. Rodríguez-Barrera et al. (2015) consider x_e > 10⁻⁷ as a threshold when the gas starts to experience plasma behavior. The calculated x_e in Fig. 1b indicates significant ionization already for upper atmospheric regions where p_gas < 1 bar. Furthermore, the electron plasma frequency exceeds the electron-neutral collision frequency on the dayside for p_gas < 0.1 bar (dotted blue line in Fig. 2), indicating that the ionized medium behaves as a plasma in this area (e.g., Baumjohann & Treumann 1996). In the partially ionized upper atmosphere, frictional heating may arise from collisions between the charged particles (ions, electrons) and the dominant neutrals due to their relative drift velocity.

The importance of the effect of anisotropic drag on the atmospheric flows of (U)HJs has also been recently highlighted by Christie et al. (2025). Their work, in which the anisotropic drag is directly derived from the generalized Ohm’s law (see Appendix B), similarly demonstrates that the Hall effect can alter global atmospheric dynamics. Our work provides an alternative derivation starting from the fundamental equations of motion for the coupled ion, electron, and neutral fluids. This approach, which is based on the framework established for studying Earth’s ionosphere-thermosphere system (e.g., Rees 1989; Song et al. 2001; Zhu et al. 2005; Schunk & Nagy 2009), explicitly parametrizes the drag in terms of the collision frequencies and gyrofrequencies of the existent particle species. This formulation represents an alternative presentation of the same underlying physics that emphasizes the role of anisotropic, collisionally driven transport and provides a physically motivated drag parametrization suitable for the implementation in GCMs.

In this section, different magnetic coupling regimes are discussed through the magnetic Reynolds number and a parametrized form of the anisotropic magnetic drag suitable for implementation in ExoRad is derived. In ExoRad, similar to many other GCMs (SPARC/MITgcm (Showman et al. 2009), THOR (e.g., Mendonça et al. 2016), RM-GCM (Rauscher & Menou 2010), and the Met Office Unified Model (UM; e.g., Wood et al. 2014; Edwards 1996; Drummond et al. 2018)), only the neutral fluid dynamics is evolved and the momentum equations for charged species are not solved explicitly. Instead, an additional source term representing the coupling effect through locally changing ionization is derived.

Before the model is presented, the key plasma parameters that are relevant for the characterization of the magnetic coupling are provided here calculated from the simulation without magnetic drag for the benefit of the reader. The electron plasma frequency, ω_pe [Hz], describes the rate at which electrons oscillate and is given by

$${\omega }_{pe}=\sqrt{\frac{n_e{e}^2}{m_e{ϵ}_0}},$$(19)

with the vacuum permittivity, ϵ₀ [C²kg⁻¹m⁻³s²], the elementary charge, e [C], and the electron number density, n_e [m⁻³]. ω_pe ranges from 5 × 10⁶ … 9 × 10¹¹ Hz on the dayside atmosphere of WASP-18 b. The relation ω_pe/ν_en between ω_pe and the electron–neutral collision frequency, ν_en (Eq. (42)), describes whether the plasma behaves more as a collisionless fluid (ω_pe/ν_en > 1) or a collisional fluid (ω_pe/ν_en < 1). It ranges from 6 × 10⁻⁴ to 50 in the dayside atmosphere. The gyrofrequency, ω_cs [Hz], gives the rate at which charged particles (index s = i for ions and s = e for electrons) gyrate around magnetic fields and is given by

$${\omega }_{cs}=\frac{\left|q_s\right|\left|{\underset{\_}{B}}_{\text{dip}}\right|}{m_s},$$(20)

with q_s [As] the charge of species s, m_s [kg] the mass of species, s, and |B_dip| the magnitude of the dipolar field (Eq. (28)). ω_ce ranges from 9×10⁷ … 2×10⁸ Hz and ω_ci from 1.5×10³ … 3×10³ Hz. The magnetization parameter, k_s (Leake et al. 2014; also called Hall parameter), describes the degree to which a charged species is magnetized and is given by

$$k_s=\frac{{\omega }_{cs}}{v_{sn}},$$(21)

with ν_sn the collision frequency between species s and neutrals (Eqs. (41), (42)). k_i ranges from 2 × 10⁻⁹ to 2 × 10⁻² and k_e from 4 × 10⁻⁷ to 3.8 on the dayside atmosphere of WASP-18 b. The efficiency of magnetic induction relative to diffusion is quantified by the magnetic Reynolds number, R_M, shown in Fig. 1a. R_M compares the timescale of magnetic field advection by the flow to the timescale of magnetic diffusion (see Eq. (A.4)) and is specifically adopted here as

$$R_M=\left|{\underset{\_}{v}}_{\text{h}}\right|H_p{\mu }_0{\sigma }_P,$$(22)

with the pressure scale height, $H_p=\frac{R_d{T}_{\text{gas}}}{g}[\text{m}]$, the magnetic permeability, ${\mu }_0\left[\frac{N}{A^2}\right]$, and the Pedersen conductivity, σ_P [S/m] (Eq. (B.2)).

The ionization fraction, x_e, quantifies the proportion of particles that are ionized compared to the total number of particles in the gas and is defined as

$$x_e=\frac{p_{\text{el}}}{p_{\text{el}}+p_{\text{gas}}}\approx \frac{n_e}{n_n},$$(23)

with the gas and electron pressure p_gas [Pa] and p_el [Pa] and assuming an electron temperature of T_e ≈ T_gas, p_gas ≫ p_el, and quasi-neutrality (n_e ≈ n_i). In the 3D GCM model it is assumed that all species (electrons, ions, and neutrals) have the same temperature. The electron pressure, p_el, is calculated with the code GGchem (Woitke et al. 2018). GGchem calculates the gas-phase chemical equilibrium using elemental abundances, the local gas temperature, and gas pressure to determine the number densities of all atomic and molecular species, including their ions and the electron pressure. The ionization fraction is then interpolated to the gas pressure-temperature profile in ExoRad for further use in computing magnetic drag forces at each grid point in the simulation. This is faster than re-running GGchem.

Fig. 1 (a) Spatial variations in magnetic Reynolds number (RM (Eq. (22)), left) and (b) thermal ionization fraction (xe, right) for WASP-18 b. Large areas of the dayside are in the intermediate and high-RM regimes, while ionization varies over many orders of magnitude, motivating an anisotropic drag treatment that accounts for local conductivity and magnetic geometry. The lines show gas temperature-pressure profiles averaged over all latitudes and over different regions in the atmosphere: dayside (−90◦< ϕ ≤ 90◦), nightside (|ϕ| > 90◦), morning (−97.5◦≤ ϕ ≤ −82.5◦) and evening terminator (82.5◦≥ ϕ ≥ 97.5◦). The results of the simulation with anisotropic drag were used in the calculation of RM, xe, and the averaged temperature-pressure profiles shown here.

Fig. 2 Magnetization parameter (Eq. (21)) profiles for ions (ki) and electrons (ke) averaged over all latitudes and over dayside atmosphere (−90◦< ϕ ≤ 90◦). The dotted blue line (ωpe/νen) shows the relation between the electron plasma frequency and the electron–neutral collision frequency averaged over the dayside atmosphere. The plasma coupling ratios were calculated from the ExoRad simulation for WASP-18 b with anisotropic drag. The dashed gray line shows where the plasma coupling ratio reaches unity.

4.1 Where does magnetic drag affect atmosphere environments?

The efficiency and geometry of the momentum exchange strongly depend on x_e, the local magnetic field orientation, and the relative importance of inductive and diffusive processes. The existing active drag parameterization used in GCM studies of hot Jupiters (e.g., Rauscher & Menou 2013) scales with the local magnetic diffusivity, η (Eq. (A.21)). While this approach captures the braking of the wind flow on the dayside, it treats the drag as isotropic, only affecting the zonal component of the wind velocity but taking into account the variations in the ionization on the day- and nightside. This approach is valid only for R_M ≪ 1 (e.g., Dietrich et al. 2022), where the plasma is only weakly coupled to the field, the magnetic diffusion dominates, and the field lines move freely through the gas. In this case the magnetic field only weakly affects the flow. However, this model does not take into account that the interaction between a conducting fluid and a magnetic field is anisotropic, which becomes important for the intermediate R_M regime (R_M ≈ 1). Charged particles are constrained to move along magnetic field lines, while their movement perpendicular to the field is restricted. This leads to an anisotropic drag force on the neutral gas. In Fig. 1a, a map of the Reynolds number for the modeled atmospheric structure shows that large regions of the dayside extend into the intermediate-to-high R_M regime. Inductive effects become important where R_M > 1 and drag is anisotropic, while the cooler nightside and deeper layers have R_M ≪ 1. The anisotropic magnetic drag is represented by separating the Pedersen and Hall contributions, which are the dissipative and non-dissipative components of the magnetic drag, respectively. The Pedersen component corresponds to drag perpendicular to the magnetic field, leading to flow damping and kinetic energy loss, whereas the Hall component corresponds to drag perpendicular to both the magnetic field and the flow velocity, modifying the flow direction without energy dissipation. As the starting point, the momentum equations for ions, electrons, and neutrals were taken (Eqs. (24)–(26)). Then drift velocities of the charged particles with respect to the neutrals were derived (Eq. (31)) and were used to parametrize the magnetic drag force (Eq. (36)) acting on the neutral gas. This formulation of the anisotropic drag is applicable across low and intermediate R_M regimes and is a significant physical improvement over isotropic drag models, providing a more realistic representation of the directional property of the Lorentz force in the large areas of the atmosphere where the low and intermediate R_M approximation is valid. In the region where R_M ≫ 1, the magnetic field lines are frozen into the plasma and advection dominates over diffusion. Here, the atmospheric gas and magnetic field are fully coupled: the gas and field move together, inducing strong currents and feedback reactions on the flow Rogers (2017). For this, a fully self-consistent MHD treatment would be required for accuracy. However, it should be noted that R_M must be interpreted carefully in the presence of strong magnetic drag. Strong drag can suppress the horizontal flow, resulting in small values of R_M that represent a reduced velocity rather than weak magnetic coupling. Therefore, a small R_M does not necessarily demonstrate if a magnetic drag parametrization is valid in this regime, as the flow might be in the ideal MHD limit. To analyze the system properly, not only R_M but the plasma parameters introduced before need to be considered (Eqs. (19)–(23)).

4.2 Momentum exchange between charged and neutral species and frictional heating

The following derivations are described in a right-handed spherical coordinate system defined by the radial distance, r [m], latitude, θ [^◦], and longitude, ϕ [^◦], with the corresponding unit vectors, ${\underset{\_}{\widehat{e}}}_r$ (radially outward), ${\underset{\_}{\widehat{e}}}_{\theta }$ (positive northward), and ${\underset{\_}{\widehat{e}}}_{\varphi }$ (positive eastward). This coordinate system, which is in height coordinates, differs from the p coordinates in which the HPEs (Eqs. (1)–(4)) are solved primarily because the vertical coordinate itself changes the meaning of vertical derivatives and the different surfaces (gas pressure vs. geometric height) on which the horizontal gradients of fields are defined.

In the planet’s rotating reference frame, the momentum equation for the neutral species is given by (see Song et al. 2001)

$${\rho }_n\frac{d{\underset{\_}{v}}_n}{dt}=-\nabla p_{\text{gas}}+{\underset{\_}{F}}_n+{\rho }_i{v}_{in}\left({\underset{\_}{v}}_i-{\underset{\_}{v}}_n\right)+{\rho }_e{v}_{en}\left({\underset{\_}{v}}_e-{\underset{\_}{v}}_n\right),$$(24)

where ${\underset{\_}{v}}_n=\left(v_{n,r},v_{n,\theta },v_{n,\varphi }\right)$ is the velocity vector of the neutral species with the radial velocity, v_n,r [m s⁻¹], the zonal velocity, v_n,ϕ [m s⁻¹], and the meridional velocity, v_n,θ [m s⁻¹], p_gas [Pa] the gas pressure of the neutral species, F_n the sum of external forces acting on the neutrals (e.g., gravity, Coriolis, centrifugal), ρ_i [kg m⁻³] the ion mass density, ρ_e [kg m⁻³] the electron mass density, and ρ_n [kg m⁻³] the neutral mass density, and ν_in, ν_en [s⁻¹] the ion-neutral and electron-neutral collision frequencies, respectively. The last two terms on the right-hand side represent the frictional drag resulting from the relative motion between charged and neutral species. In this paper, only gas species are considered. The momentum Eq. (24) without the frictional terms in height coordinates is equivalent to the horizontal and vertical momentum (Eqs. (1) and (2)) in p coordinates using the hydrostatic balance assumption in the vertical direction.

The momentum equations for electrons (subscript s = e) and ions (s = i) are given by (see, e.g., Song et al. 2001)

$$n_e{m}_e\frac{d{\underset{\_}{v}}_e}{dt}=-\nabla P_e+{\underset{\_}{F}}_e+e{n}_e\left(\underset{\_}{E}+{\underset{\_}{v}}_e\times \underset{\_}{B}\right)-n_e{m}_e{v}_{en}\left({\underset{\_}{v}}_e-{\underset{\_}{v}}_n\right),$$(25)

$$n_e{m}_i\frac{d{\underset{\_}{v}}_i}{dt}=-\nabla P_i+{\underset{\_}{F}}_i+e{n}_e\left(\underset{\_}{E}+{\underset{\_}{v}}_i\times \underset{\_}{B}\right)-n_e{m}_i{v}_{in}\left({\underset{\_}{v}}_i-{\underset{\_}{v}}_n\right),$$(26)

where e [C] is the elementary charge, F_s the external forces exerting on the charged species, s, and n_s [m⁻³], P_s [Pa], and m_s [kg] the number density, pressure, and mass of each charged species, respectively. E and B are the electric and magnetic fields. Here, the terms with electron-ion collision frequency, ν_ei, and ion-electron collision frequency, ν_ie, were neglected, since our estimates have shown that for our parameter regime ν_en ≫ ν_ei and ν_in ≫ ν_ie. Since ExoRad does not evolve ion or electron dynamics, a steady-state approximation $\left(\frac{d}{dt}=0\right)$ for their drift velocities was derived using a reduced momentum balance equation for electrons and ions, assuming singly charged ions, quasi-neutrality, and neglecting pressure gradients and inertial terms (e.g., Schunk & Nagy 2009):

$$e\left(\underset{\_}{E}+{\underset{\_}{v}}_s\times \underset{\_}{B}\right)-m_s{v}_{sn}\left({\underset{\_}{v}}_s-{\underset{\_}{v}}_n\right)\approx 0.$$(27)

We further assumed that the large-scale electric field in the planetary rest frame is negligible (E ≈ 0). This implies that the electric field in the local frame of the neutral fluid is E′ = E+ v_n × B ≈ v_n× B. This approach does not take into account the feedback from the global current system, which would modify the large scale electric field and currents. The main limitation of the model is that it neglects the polarization electric field and therefore does not ensure current closure. In a partially ionized atmosphere, collisions between neutrals and charged particles introduce charge separation, which creates a polarization electric field. This field is built up until it balances the differential drifts of ions and electrons, though it maintains current continuity by diverting currents within the ionosphere. By neglecting the polarization field, the model only describes the local Lorentz force component of the coupling between charged particles and neutrals, without taking into account any self-consistent electrodynamic feedback in the ionosphere. Therefore, the calculated magnetic drag describes a local forcing approximation rather than the global response of an ionospheric current system. However, this parametrization provides a physically based and computationally efficient method of studying the main large-scale effects of magnetic coupling on the atmospheric dynamics and thermal structure. The model separates the impact of direct collisional drag between charged particles and neutrals, which dominates in the upper atmosphere. Nevertheless, a more complete approach, i.e., solving the Poisson’s equation, would enable global current closure and feedback from large-scale electric fields (see Sect. 6). These effects could modify the magnitude and spatial distribution of the frictional drag and heating. This, however, is outside the scope of the present paper.

We further assumed a dipolar large-scale planetary magnetic field anti-aligned with the planetary rotation axis. The field is purely poloidal, with the magnetic flux density given by

$${\underset{\_}{B}}_{\text{dip}}(r,\theta )=B_0{\left(\frac{R_{\text{ref}}}{r}\right)}^3\left(-2\sin \theta {\underset{\_}{\widehat{e}}}_r+\cos \theta {\underset{\_}{\widehat{e}}}_{\theta }\right).$$(28)

To calculate the dipolar magnetic field, the equatorial magnetic field strength, B₀ [T], was defined at the reference radius, R_ref [m] (Sect. 2). The northward component, B_dip,θ [T], is positive in the direction of increasing θ. The magnitude of the magnetic field, |B_dip|, and its unit vector, $\underset{\_}{\widehat{b}}=\left({\widehat{b}}_r,{\widehat{b}}_{\theta },{\widehat{b}}_{\varphi }\right)$, are

$$\left|{\underset{\_}{B}}_{\text{dip}}\right|=B_0{\left(\frac{R_{\text{ref}}}{r}\right)}^3\sqrt{1+3{\sin }^2\theta },$$(29)

$$\widehat{\underset{\_}{b}}=\frac{{\underset{\_}{B}}_{\text{dip}}}{\left|{\underset{\_}{B}}_{\text{dip}}\right|}=\frac{\left(-2\sin \theta {\underset{\_}{\widehat{e}}}_r-\cos \theta {\underset{\_}{\widehat{e}}}_{\theta }\right)}{\sqrt{1+3{\sin }^2\theta }}.$$(30)

To simplify the implementation, the r dependence of the magnetic field is neglected in the simulations of ExoRad. Treating B_dip as constant with altitude introduces only a small error of less than 4% in |B_dip| in the modeled atmospheric region of WASP-18 b, corresponding to a radial extent of ∆z ≈ 1075 km. This estimate is specific to the geometry and scale height of WASP-18 b. In (U)HJs with more extended atmospheres or smaller radii, the variation in B_dip with altitude could be more pronounced and should be accounted for in the analysis. Additionally, it is assumed for the drag calculation that the neutral wind has no radial component (v_n,r ≈ 0 m/s), as vertical winds are weak compared to horizontal winds. Solving Eq. (27) for the drift velocity Δv_s = v_s −v_n yields

$$\Delta {\underset{\_}{v}}_s=\pm \frac{k_s}{1+k_s^2}\left({\underset{\_}{v}}_n\times \underset{\_}{\widehat{b}}\right)-\frac{k_s^2}{1+k_s^2}{\underset{\_}{v}}_{n,\perp },$$(31)

where the upper sign applies to ions (s = i) and the lower one to electrons (s = e). The velocity of the neutrals perpendicular to B_dip is given by

$${\underset{\_}{v}}_{n,\perp }={\underset{\_}{v}}_n-\left({\underset{\_}{v}}_n·\underset{\_}{\widehat{b}}\right)\underset{\_}{\widehat{b}}.$$(32)

Equation (31) shows that if k_s ≪ 1, collisions dominate and the charged species are strongly coupled to the neutrals. In this regime, magnetic drag is weak. For k_s ≫ 1, the charged particles complete many gyro-orbits before being deflected by collisions. They are magnetized and coupled primarily to the magnetic field, leading to significant relative drift and enhanced magnetic drag on the neutrals. Substituting Eq. (31) into the momentum exchange terms of Eq. (24), the total drag force density exerted by charged particles on the neutral fluid is

$$\begin{array}{l}{\underset{\_}{F}}_{\text{drag }}={\rho }_i{v}_{in}\left(\frac{k_i}{1+k_i^2}\left({\underset{\_}{v}}_n\times \underset{\_}{\widehat{b}}\right)-\frac{k_i^2}{1+k_i^2}{\underset{\_}{v}}_{n,\perp }\right)\\ +{\rho }_e{v}_{en}\left(-\frac{k_e}{1+k_e^2}\left({\underset{\_}{v}}_n\times \underset{\_}{\widehat{b}}\right)-\frac{k_e^2}{1+k_e^2}{\underset{\_}{v}}_{n,\perp }\right).\end{array}$$(33)

To estimate the contributions from ion and electron drag, the Pedersen K_P,s [kg m⁻³s⁻¹] and Hall drag coefficients, K_H,s [kg m⁻³s⁻¹], for each species are introduced:

$$K_{P,s}={\rho }_s{v}_{sn}\frac{k_s^2}{1+k_s^2},$$(34)

$$K_{H,s}={\rho }_s{v}_{sn}\frac{k_s}{1+k_s^2},$$(35)

so that the total drag force density simplifies to

$${\underset{\_}{F}}_{\text{drag }}=(K_{H,i}-K_{H,e})({\underset{\_}{v}}_n\times \underset{\_}{\widehat{b}})-(K_{P,i}+K_{P,e}){\underset{\_}{v}}_{n,\perp }.$$(36)

The first term represents the Hall drag, which acts perpendicular to both v_n and B_dip, while the second term corresponds to the Pedersen drag, which opposes the perpendicular motion of the neutrals to B_dip. Although ion and electron velocities are not evolved dynamically, their influence is included through the Hall and Pedersen coefficients, which depend on local x_e, plasma parameters, and magnetic field strength of the dipolar field. With Eq. (23) the electron and ion mass density (ρ_s = m_sx_en_n) can be calculated.

ExoRad assumes vertical hydrostatic equilibrium, employing gas pressure as the vertical coordinate. It does not explicitly solve the vertical momentum equation. Instead, vertical velocities are derived to ensure mass continuity. While vertical magnetic drag forces can occur (e.g., Christie et al. 2025), their inclusion would violate the vertical hydrostatic equilibrium. However, the resulting vertical acceleration is expected to be negligibly small compared to gravity in the weakly ionized regime considered here. Therefore, the radial drag is neglected and only the meridional and zonal components of the drag force density are computed:

$$F_{\text{drag },\theta }=K_H{v}_{n,\varphi }{\widehat{b}}_r-K_P{v}_{n,\theta }{\widehat{b}}_r^2$$(37)

$$F_{\text{drag },\varphi }=-K_H{v}_{n,\theta }{\widehat{b}}_r-K_P{v}_{n,\varphi },$$(38)

where

$$K_H=K_{H,i}-K_{H,e}\left[\text{kg} {\text{m}}^{-3} {\text{s}}^{-1}\right],$$(39)

$$K_P=K_{P,e}+K_{P,i}\left[\text{kg} {\text{m}}^{-3} {\text{s}}^{-1}\right],$$(40)

are the total Hall and Pedersen coefficients, respectively. The geometric relation ${\widehat{b}}_{\theta }^2=1-{\widehat{b}}_r^2$ was used, assuming the azimuthal component of the magnetic field to be negligible $\left({\widehat{b}}_{\varphi }=0\right)$. The expression in Eqs. (37) and (38) was implemented in ExoRad as an external drag force acting on the horizontal momentum Eq. (1) via ${\underset{\_}{F}}_v=\frac{1}{{\rho }_n}{\underset{\_}{F}}_{\text{drag}}$ and ν_en are

$$v_{in}=4.6\times {10}^{-16}{n}_n,$$(41)

$$v_{en}=6.967\times {10}^{-14}{n}_n{T}_{\text{gas }}^{0.1},$$(42)

where n_n [m⁻³] is the neutral number density. The electron-H collision rate follows Koskinen et al. (2010), based on Danby et al. (1996). The chemical atmosphere gas composition strongly varies between the day- and nightside. While the dayside is dominated by atomic H, the nightside is dominated by H₂ (Fig. 8 of Helling et al. 2019a). In this work, ν_en was computed assuming atomic hydrogen as the dominant neutral species, which is justified on the dayside of WASP-18 b. On the nightside, where molecular hydrogen dominates, ν_en might be underestimated, as collision frequencies differ (Eq. (12) in Koskinen et al. 2010). However, the electron density and therefore the magnetic coupling strongly decrease on the nightside and the resulting impact on the magnetic drag is small. The ion-neutral rate is based on a rigid-sphere model for non-resonant collisions (Chapman 1956; Günzkofer et al. 2023), assuming a mean atomic mass number of 32 that was estimated by summing the relevant atomic masses that make up an ion as given in GGchem. However, the effective ion mass may vary with altitude and longitude as the dominant ion species is determined by the gas temperature and elemental abundances. In the plasma regime considered here, ions are collisionally coupled to the neutrals (k_i ≪ 1; Fig. 2), so the resulting variations in ion mass have a comparatively weak effect on the conductivities relative to the dominant dependence on electron density. Therefore, the assumed constant mass number is sufficient for the present work; however, future studies including altitude-dependent ion composition would be valuable. Resonant collisions (e.g., Na⁺–Na) are negligible due to the dominance of hydrogen in the neutral atmosphere.

The total drag force density (Eq. (36)) provides both a sink of momentum and a source of frictional heating, and is particularly important where the magnetization parameters, k_i and k_e, are large. In Fig. 3, the Hall and Pedersen drag coefficients according to Eqs. (34) and (35) were calculated using the gas pressure (p_gas) and temperature (T_gas) output from the ExoRad simulation of the WASP-18 b atmosphere with anisotropic magnetic drag, as well as x_e output from GGchem. The coefficients shown are averaged over all latitudes and over the dayside atmosphere (−90^◦< ϕ ≤ 90^◦). For p_gas ≥ 10⁻³ bar, K_H,e and K_H,i (orange lines) are of similar magnitude and largely cancel out due to opposite signs (K_H, magenta line), making the electron Pedersen drag (K_P,e, dotted cyan line) the dominant contributor. At higher altitudes where p_gas < 10⁻³ bar, k_e > 1 and k_i < 1 (see Fig. 2) but the small electron mass suppresses K_H,e, so K_H,i dominates the total Hall drag (K_H) and the Hall drag dominates the total Pedersen drag (K_P). To maintain the numerical stability in ExoRad (Sect. 2), the drag timescale is constrained to τ_drag = min(ρ_n/K_H, ρ_n/K_P) ≥ 100 s. Consequently, the imposed lower limit reduces the effective drag strength at low pressures. Therefore, the dayside-averaged K_P deviates from its physically expected value at p_gas < 10⁻² bar (gray line in Fig. 3), and K_H at p_gas ≲ 10⁻³ bar.

In addition to affecting neutral momentum, magnetic drag contributes to energy dissipation via frictional heating. The heating of the neutral atmosphere due to magnetic drag is often referred to as Joule or Ohmic heating, but it is more accurately described as frictional heating resulting from momentum exchange between drifting charged particles and neutrals and not electromagnetic dissipation, which results from electric currents interacting with the electrical resistivity of the gas (Vasyliūnas & Song 2005). For energy to be conserved, the corresponding rate at which energy is dissipated into heat per unit mass of the neutral gas is equal to the rate at which Pedersen drag force works on the fluid, given by

$$q_{\text{fric }}=-\frac{1}{{\rho }_n}{\underset{\_}{F}}_{\text{drag }}·{\underset{\_}{v}}_n=\frac{1}{{\rho }_n}{K}_P{\underset{\_}{v}}_{n,\perp }^2.$$(43)

q_fric is always positive, as it represents a conversion of neutral kinetic energy into thermal energy. The Hall drag force is a deflective force: it changes the direction of the flow but does not remove kinetic energy from the flow. The Pedersen drag component is frictional and removes kinetic energy from the flow, which is converted into heat. To include the frictional heating in the energy equation in ExoRad, q_fric is evaluated at each time step and added to the potential temperature tendency equation.

Fig. 3 Hall (KH,s) and Pedersen (KP,s) drag coefficients for electrons (s = e) and ions (s = i) averaged over all latitudes and over dayside atmosphere (−90◦< ϕ ≤ 90◦). KH = KH,i − KH,e and KP = KP,i + KP,e are the total Hall and Pedersen coefficients. The coefficients were calculated using the gas pressure-temperature profile from the ExoRad simulation run of WASP-18 b with anisotropic drag. The gray line shows KP obtained directly from ExoRad. It diverges from the calculated KP (blue line) in the upper atmosphere because the drag timescale in ExoRad is limited to τdrag ≥ 100 s. The same applies for KH from ExoRad (not shown here).

5 Impact of magnetic drag on winds and temperatures in the atmosphere of WASP-18 b

Magnetic drag is expected to be crucial in shaping UHJs atmospheres, yet its effects remain insufficiently studied and constrained. Here, it is shown how magnetic drag modifies the wind pattern and thermal structure in the atmosphere of WASP-18 b and how different magnetic drag models compare to each other. The extent to which the magnetic drag disrupts the equatorial super-rotation, how it modifies heat redistribution between day and night, and whether it can change the flow pattern will be explored. The results show that magnetic drag generally weakens equatorial super-rotation, and that anisotropic drag can lead to a flow decoupling and deflection between the day- and nightside.

5.1 Atmospheric circulation and different drag treatments

Figure 4 shows the horizontal wind speeds and streamlines for the different drag treatments (no drag, uniform drag, active drag, and anisotropic drag) for WASP-18 b as a test case planet. Figure 5 shows the corresponding mean zonal wind structures. In the absence of magnetic drag (when momentum and energy transfers are unaffected by the magnetic drag), the circulation is dominated by a pronounced eastward (super-rotating) equatorial jet extending across the planet. However, the strength of the jet is weakened in the upper dayside atmosphere (p_gas ≤ 10⁻³ bar: Figs. 4a, 5a), likely as a result of radiative cooling and heating or vertical wind shear. The upper atmosphere is more strongly affected by direct stellar irradiation and less by advection compared to the denser lower atmosphere. The strong irradiation then heats the lower atmosphere and determines the thermal structure via radiative processes. Steep thermal gradients from irradiation create vertical wind shear of the zonal wind. Such shear leads to momentum redistribution and weakening of the zonal jet in the upper atmosphere. Thus, radiative processes play an important role in the upper atmosphere, while advective heat redistribution is more dominant in the lower atmosphere. This circulation pattern is common in hot Jupiter simulations without magnetic drag (e.g., Showman et al. 2009; Showman & Polvani 2010, 2011), driven by a Matsuno–Gill–type response (Matsuno 1966; Gill 1980). Peak horizontal wind speeds exceed 10 km/s.

If a constant magnetic drag acts uniformly on the horizontal momentum and energy transfer (uniform drag), the equatorial super-rotating jet is largely suppressed (Figs. 4b, 4f, and 5j). The circulation is instead dominated by eastward and westward flows of comparable magnitude, which mostly cancel in the zonal mean (Fig. 5j). Similar to Beltz et al. (2022b, see their Fig. 8 for the total atmosphere), a weaker eastward jet is present at p_gas < 1 bar in higher latitudes (|ϕ| > 50^◦) on the day-side atmosphere (Fig. 5b). On the nightside (Fig. 5f), a weak and narrow equatorial jet is visible for p_gas < 10⁻³ bar. The circulation is concentrated near the day–night terminators (longitudes ≈ ± 90^◦), with enhanced horizontal wind speeds. The flow diverges out of the substellar region. Uniform drag removes momentum throughout the atmosphere uniformly, damping the wave–driven feedback that sustains the equatorial super-rotation. Hence, a strong global uniform drag damps the equatorial jet throughout the entire atmosphere. However, such a model does not capture the underlying physics, since it neglects the strong spatial variability of the ionization fraction, particularly on the day- and nightside, resulting in magnetic coupling of both the dayside and nightside.

A drag confined only to the zonal direction (active drag; Figs. 4c and 4g) opposes the zonal flow and decouples the day-to nightside flow completely (at p_gas < 10⁻³ bar), and therefore suppresses the equatorial jet on the dayside (Fig. 5c). At low pressures higher in the atmosphere, where drag timescales are extremely short, the zonal (east–west) flow is efficiently damped and the flow on the dayside is redirected toward the poles, producing a largely meridional (north–south) circulation on the dayside (Figs. 4c and 5c), consistent with previous studies (e.g., Beltz et al. 2022b). A zonal magnetic drag is based on the assumption that only meridional currents affect the atmospheric dynamics. The vanishing of the zonal magnetic drag at the equator due to the sin θ dependence in Eqs. (A.14) and (A.15) has previously been discussed by Christie et al. (2025) and is visible in the streamline in the zonal direction at θ = 0^◦ on the dayside in the upper atmosphere (Fig. 4c. At higher pressures (p_gas = 10⁻¹ bar; Fig. 4g), where drag is weaker, zonal flow partially recovers, but meridional flow remains the dominant component.

A more complete representation of the effect of a planetary dipole magnetic field is achieved with a magnetic drag force acting in both the zonal and meridional directions (anisotropic drag). The drag strongly damps the flow, especially at low gas pressures (p_gas = 10⁻³ bar) on the dayside (Figs. 4d and 5d).

The resulting magnetic drag forces effectively isolate the day-side from the nightside circulation at the equatorial region at the morning terminator (ϕ ≈ −90^◦; Fig. 4d). At higher latitudes (ϕ > |30^◦|) at the morning terminator, the gas flow is restored from the dayside to the nightside. At the evening terminator (ϕ ≈ 90^◦), the streamlines are directed from the dayside to the nightside. Despite this strong damping, an eastward equatorial jet remains on the nightside and extends over a broader latitude range than in the active drag case (Fig. 5h). The wind speed shows dayside stagnation with limited redistribution, coupled to a nightside jet sustained by momentum convergence from deeper levels. This indicates that anisotropic drag firstly acts to dissipate kinetic energy, slow down, and deflect the flow on the dayside, while still allowing the equatorial jet to be maintained on the nightside. The Pedersen drag primarily suppresses zonal winds, producing significantly reduced flow speeds on the dayside, and the Hall drag generates a non-dissipative magnetic coupling between the meridional and zonal flow components (Eqs. (37) and (38)). In the upper atmosphere, where the Hall drag strength becomes comparable to the Pedersen drag, the Hall drag introduces an asymmetry in the flow pattern. This asymmetry is visible in Fig. 4d in the dayside streamlines between −90^◦ < ϕ < 0^◦ and 0^◦ < ϕ < 90^◦. In the longitude range −90^◦ < ϕ < 0^◦, the Hall drag deflects the flow toward more polar-directed motion, whereas in 0^◦ < ϕ < 90^◦ the stream lines show increased curvature and are redirected toward the evening terminator. If the magnetic dipole is aligned with the planetary rotation axis (not shown here), ${\widehat{b}}_r$ (Eq. (30)) changes sign while maintaining the same magnitude. In that case, the Pedersen damping does not change, but the magnetic coupling produced by the Hall drag between v_n,θ and v_n,ϕ reverses sign, leading to a circulation pattern, which is mirrored on the day-side with enhanced curvature occurring in −90^◦ < ϕ < 0^◦ and more meridional flow in 0^◦ < ϕ < 90^◦. This behavior demonstrates the strong dependence of the atmospheric flow on the magnetic field geometry in the upper atmosphere. The resulting flow asymmetry determines the distribution of momentum across the terminator regions and indirectly affects the flow on the nightside. The flow perturbations due to the Hall effect on the nightside do not arise from strong local Hall drag, which is weak due to low ionization fraction, but are instead a nonlocal response from the redirected momentum on the dayside. At the poles on the dayside where the magnetic field lines are nearly vertical, the flow experiences a stronger drag in the meridional direction than at the equator where the field lines are horizontal, meridional magnetic drag force vanishes, and only the zonal velocity is damped by the Pedersen drag component. This is seen in Eqs. (37) and (38), where the polar flow structure is produced by the dominance of the radial magnetic field component, which reaches $\left|{\widehat{b}}_r\right|=2$ at the poles. Both Pedersen (in the meridional direction of the magnetic drag force) and Hall drag terms are increased in the polar regions. The Pedersen drag provides strong damping, particularly in the meridional direction, where the damping scales as ${\widehat{b}}_r^2$, producing strong suppression of polar motions. On the other hand, the Hall drag couples the zonal and meridional velocity components through terms proportional to ${\widehat{b}}_r$, creating a rotation of the horizontal momentum. The Hall term rotates the velocity vector against strong Pedersen damping, producing the flow curvature and directional asymmetry seen in the polar regions (Fig. 4d). At gas pressures of 1 mbar, the overall wind structure is similar to that shown by Christie et al. (2025, their Fig. 13). However, unlike in the results shown here, the wind flow in Christie et al. (2025) for their anisotropic approach does not experience longitudinal asymmetries due to the Hall drag and the partial decoupling of the equatorial jet at the morning terminator. The hottest planet considered in Christie et al. (2025) reaches lower dayside temperatures than WASP-18 b, resulting in a negligible contribution of the Hall term in the drag force.

In both the active drag and anisotropic drag cases, the equatorial jet reestablishes on the dayside for p_gas > 10⁻² bar (active drag, Fig. 5c) and p_gas > 10⁻¹ bar (anisotropic drag, Fig. 5d). While the variation in jet width with pressure is similar in the anisotropic drag and no drag cases (compare Figs. 5e and 5h), the jet in the active drag case is narrower. Plaschzug et al. (2026) introduced the width of the equatorial jet as well as the wind speed as two of four climate characteristics to help characterize the potential diversity of climate regimes in the ensemble of gas giant exoplanets (see their Fig. 5). Table 3 shows how jet width and wind speeds are affected by magnetic drag. The comparison between the parameters of a planet with an F host star and a global temperature of 2400 K from Fig. 5 of Plaschzug et al. (2026) and the calculated parameters from this study for the case in which magnetic effects are absent (no drag; Table 3) shows that the wind jet width is significantly larger, the evening-morning temperature difference is smaller, and the day-nightside temperature difference is larger in this study. This discrepancy might come from different initial conditions for the planet and the exclusion of TiO and VO in the simulations for UHJs in Plaschzug et al. (2026). Including TiO and VO leads to upper dayside atmosphere thermal inversion, as is seen for example in the red line of Fig. 1. Deline et al. (2025) have shown that including magnetic drag and the strong local gas-phase absorber TiO and VO are important to fit the data from optical to IR wavelengths for WASP-18 b. The wind velocity is the highest when the magnetic drag is absent and lowest when the drag force acts uniformly on the horizontal flow, implying that the magnetic drag significantly reduces the maximum wind speeds.

Fig. 4 Maps of horizontal wind speed |vh| and streamlines (white lines) at two pressure levels (pgas =0.1 bar and pgas =0.001 bar) for the four different magnetic drag treatments (no drag, uniform drag, active drag, and anisotropic drag). The substellar point is centered at each map. The vertical dashed black lines mark the area at longitude ±90◦. The output of each velocity component is time-averaged over 100 days of simulation time to eliminate the small-scale fluctuations.

Fig. 5 Mean zonal wind for the four different magnetic drag treatments in different regions of the atmosphere: dayside atmosphere (first row), nightside atmosphere (second row), and total atmosphere (third row). The zonal wind velocity, u, is averaged at a specific latitude, and pressure level across longitudes on the dayside (−90◦≤ ϕ ≤ 90◦), on the nightside (|ϕ| > 90◦), and all longitudes. The black line shows the zero-wind contour, the boundary between super-rotation and counter-rotation. The output of the zonal wind velocity is time-averaged over 100 days of simulation time to eliminate the small-scale fluctuations.

5.2 Gas temperature and different drag treatments

In Fig. 6, horizontal temperature maps for different drag treatments and pressure levels and the hotspot regions are shown. When active drag, anisotropic drag, or no drag is applied to the atmospheric flow, a hotspot splitting in the upper atmosphere (p_gas = 10⁻³ bar) is observed almost symmetric around the equator where two distinct temperature maxima regions instead of one appear. This kind of splitting was also seen in other models (e.g., Carone et al. 2020; Komacek & Showman 2016). In Fig. 6e, a strong eastward equatorial jet is visible (see white arrows), which transports heat efficiently from the substellar point eastward, resulting in an eastward shift of the peak temperature (see cyan contours). The wind direction around the equator is almost purely zonal. When this wind flow slightly changes so that wind also points in the meridional direction, these winds shear the hottest region apart and redistribute it meridionally (Fig. 6a). The hotspots can be seen then as two symmetric lobes around the equator. If magnetic drag acts uniformly in the atmosphere, the hotspot area is almost centered at the substellar point because the equatorial jet is suppressed (Figs. 6b and 6f). Here, the momentum damping is constant.

The active and anisotropic drag introduce a nonlinear feedback that arises because neutral winds adjust the induced electric field $({\underset{\_}{E}}^{\prime }={\underset{\_}{v}}_n\times \underset{\_}{B})$, which drives currents. These currents depend on the locally varying conductivities. Then they act back on the winds via damping and Hall deflection (only in the anisotropic drag case). Therefore, these cases suppress equatorial super-rotation on the dayside at low pressures (Figs. 5c, d) and distort zonal flows, leading to no efficient eastward advection of heat. When strong magnetic drag (Pedersen + Hall) acts on the flow, it damps the wave response. The circulation then becomes less wave-driven and more controlled by pressure-gradient forces. In Table 3 the location of the maximum temperatures for the different drag treatments at p_gas = 10⁻¹ bar is summarized. In the case of active and anisotropic drag, a shift of T_max eastward of the substellar point by 7.5^◦ is seen, which is less compared to the case in which no magnetic drag is applied to the atmosphere (12.5^◦). The reason is that the magnetic drag suppresses or reduces the eastward equatorial jet. This was also proposed by previous studies; for example, Rogers & Komacek (2014) and Kataria et al. (2015). In the active drag and anisotropic drag cases, an asymmetry of the hotspot regions emerges (Figs. 6c, 6d, 6g, 6h). In the active drag model, the primary hotspot region is displaced by ∼20^◦ south of the equator. A secondary hotspot region develops ∼20^◦ north of the equator, but its temperature remains about 10 K cooler than the southern hotspot at p_gas = 10⁻¹ bar. Hence, locally changing drag forces and energy inputs may introduce another level of observable asymmetries; namely, in the form of a primary and a secondary hotspot region. Furthermore, the transition region between the cooler nightside and hot dayside at the morning terminator is shifted eastward at p_gas = 10⁻¹ bar around the equatorial region for the no drag, active drag, and anisotropic drag cases, with the shift in the no drag being the strongest and the one in the anisotropic drag being the weakest (compare the transition region to the position of the dashed black line at ϕ = −90^◦ in Figs. 6e, 6g, and 6h). A similar eastward shift is observed at the evening terminator. This eastward shift of the day–night boundary appears due to the zonal advection of heat by the equatorial jet, which transports the thermal hotspot downstream of the substellar point. When magnetic drag is introduced, Lorentz forces act against the zonal flow and dissipate kinetic energy, thereby weakening the jet and reducing the longitudinal offset of the thermal transition region. Consequently, stronger or more spatially extended drag (anisotropic drag) suppresses the eastward shift of the terminator.

Figure 7 shows the gas temperature-pressure profiles averaged across all latitudes within different longitude regions for the four drag treatments. All models show a strong difference between the day- and nightside in the upper atmosphere, with the dayside reaching almost 3000 K and the nightside remaining below 1100 K. The profiles of the terminator regions lie in between, with the evening terminator being warmer than the morning terminator for the no drag, active drag, and anisotropic drag models, consistent with the eastward-directed flow that advects hot material. When a spatially uniform magnetic drag is applied, the terminator profiles almost overlap for p_gas < 10⁻¹ bar. The temperature differences between the day- and night-side and between the evening and morning terminator, the two remaining climate characteristics (Plaschzug et al. 2026), are displayed in Table 3. It shows that the difference between the averaged temperatures on the morning and evening terminator at p_gas = 10⁻¹ bar is strongest in the anisotropic drag case, with the evening terminator being 20% hotter than the morning terminator, and weakest in the uniform drag treatment. Temperature inversions are present on the dayside and absent on the nightside. On the nightside, the magnetic drag also affects the temperature profiles. The temperatures are lower when magnetic drag is acting on the atmospheric dynamics (active drag and anisotropic drag) compared to the case in which magnetic field does not affect the dynamics (no drag). The nightside temperature of the uniform drag model is larger than for the other cases and almost constant for p_gas < 10⁻¹ bar, which is due to the suppression of zonal heat redistribution. In the uniform drag case, the temperature inversion at the dayside is stronger and lies deeper in the atmosphere (p_gas = 10 bar) than for the other drag treatments. This is also the case for the evening and morning terminator regions. At p_gas > 10 bar, the profiles converge across the regions for all drag treatments, indicating that the deeper atmosphere is much less affected by the drag.

Fig. 6 Maps of gas temperature and velocity field at two gas pressure levels (pgas = 0.1 bar and pgas = 0.001 bar) for the four different magnetic drag treatments. The substellar point is centered at each map and is marked with a white star. The velocity field is shown by the white arrows. The length of the arrow scales with the relative wind velocity. The vertical dashed black lines mark the area at longitude ±90◦. The cyan contours represent the hotspot region (Tgas,max, Tgas,max − 1K) with Tgas,max being the largest temperature for the chosen pressure level. The temperature is time-averaged over 100 days of simulation time to eliminate the small-scale fluctuations.

Fig. 7 Gas temperature-pressure profiles for four different drag treatments averaged over all latitudes and over different regions in the atmosphere: dayside (−90◦< ϕ ≤ 90◦), nightside (|ϕ| > 90◦), morning terminator (−97.5◦≤ ϕ ≤ −82.5◦), and evening terminator (82.5◦≥ ϕ ≥ 97.5◦). All drag treatments show pronounced day-night differences and temperature inversions in the upper atmosphere. Nightside temperatures are cooler when active or anisotropic magnetic drag is applied compared to no drag.

6 Discussion

The implementation of the anisotropic magnetic drag demonstrates that magnetic effects can significantly reshape the atmospheric circulation of UHJs. The main outcome is that the magnetic drag disrupts the equatorial super-rotation on the day-side, thereby decoupling the dayside and nightside circulations at the equatorial region in the upper atmosphere, deflecting the flow on the dayside, and altering the efficiency of heat redistribution. The drag dissipates kinetic energy and dampens the flow, while maintaining an eastward jet of similar width on the nightside compared to an atmosphere that is not affected by the magnetic drag. Additionally, the dayside–nightside temperature difference increases and the hotspot and terminator shift toward the east becomes weaker. The implementation of anisotropic magnetic drag demonstrates an important step toward including plasma-neutral interactions in ExoRad, but several model limitations affect the interpretation of the simulation results.

Large-scale planetary electric field. A main simplification in the present model is the assumption of a negligible large-scale planetary electric potential $({\underset{\_}{E}}_{\text{pol}}=0)$ in the planetary rest frame. In reality, both dynamo-generated and polarization electric fields can arise from charge separation and currents. These fields influence the charged particle velocities, and hence the collisional drag felt by neutrals. By enforcing ${\underset{\_}{E}}_{\text{pol }}=0$, the charged particle’s motion is controlled by ${\underset{\_}{v}}_n\times \underset{\_}{B}$. This includes the direct collisional drag but excludes feedback from electric fields. A self-consistent approach of the electric field would require one to solve a complex elliptic equation (Poisson’s equation) for the global polarization electric field across the atmosphere at every time step (e.g., Koskinen et al. 2014). Therefore, the current continuity equation ($\nabla ·\underset{\_}{j}=0$ with the current density, $\underset{\_}{j}=\underset{\_}{\underset{\_}{\sigma }}·{\underset{\_}{E}}^{\text{'}}$, the conductivity tensor, $\underset{\_}{\underset{\_}{\sigma }}$, and the electric field seen in the frame of the neutrals, $\underset{\_}{{\underset{\_}{E}}^{\prime }={\underset{\_}{E}}_{\text{pol}}+{\underset{\_}{v}}_n\times \underset{\_}{B}}$ has to be solved and the equation to be calculated then becomes $\nabla ·(\underset{\_}{\underset{\_}{\sigma }}·\nabla {\Phi }_{\text{pol}})=\nabla ·(\underset{\_}{\underset{\_}{\sigma }}·({\underset{\_}{v}}_n\times \underset{\_}{B}))$, where ${\underset{\_}{E}}_{\text{pol }}$, which is electrostatic, is expressed as a gradient of a scalar potential ${\underset{\_}{E}}_{\text{pol}}=-\nabla {\Phi }_{\text{pol}}$. To solve the Poisson equation is computationally expensive and beyond the scope of the present work.

To make a rough estimate of the error that is made by neglecting ${\underset{\_}{E}}_{\text{pol }}$, Ohm’s law was considered: $\underset{\_}{j}={\sigma }_P{\underset{\_}{E}}_{\perp }^{\prime }+{\sigma }_H\left(\underset{\_}{\widehat{b}}\times \underset{\_}{E^{\prime }}\right)+{\sigma }_{\Vert }{\underset{\_}{E}}_{\Vert }^{\prime }$, where σ_P is the Pedersen (Eq. (B.2)), σ_H is the Hall (Eq. (B.3)), and σ_∥ is the parallel conductivity (Eq. (A.23)), and the different components of the electric field are parallel $\left({\underset{\_}{E}}_{\Vert }^{\prime }\right)$ and perpendicular $\left({\underset{\_}{E}}_{\perp }^{\prime }\right)$ to the magnetic field. ${\underset{\_}{E}}_{\text{pol}}$ arises from charge separation. The initial motional electric field, ${\underset{\_}{E}}_0={\underset{\_}{v}}_n\times \underset{\_}{B}$, drives two primary currents (Pedersen: ${\underset{\_}{j}}_{\text{P,prim }}={\sigma }_P{\underset{\_}{E}}_0$, Hall: ${\underset{\_}{j}}_{\text{H,prim }}={\sigma }_H\left(\underset{\_}{\widehat{b}}\times {\underset{\_}{E}}_0\right))$ and causes charges to build up at the insulating boundaries (non-conducting nightside at terminator region, non-conducting gas above and below the ionosphere). The charge build-up creates ${\underset{\_}{E}}_{\text{pol }}$ opposing the initial charge flow. This means that ${\underset{\_}{E}}_{\text{pol}}$ drives a secondary Pedersen current $\left({\underset{\_}{j}}_{\text{P,sec }}={\sigma }_P{\underset{\_}{E}}_{\text{pol }}\right)$ that cancels out the primary Hall current $\left({\underset{\_}{j}}_{\text{H},\text{ prim }}={\sigma }_H\left(\underset{\_}{\widehat{b}}\times {\underset{\_}{E}}_0\right)\right)$ to ensure $\nabla ·\underset{\_}{j}=0$. A simplified calculation of the polarization electric field with approximately uniform conductivities would require one to calculate ${\underset{\_}{j}}_{\text{P},\text{sec}}=-{\underset{\_}{j}}_{\text{H},\text{ prim }}$, which gives ${\underset{\_}{E}}_p\approx -\frac{{\sigma }_H}{{\sigma }_P}\underset{\_}{\widehat{b}}\times {\underset{\_}{E}}_0\approx -k_e\underset{\_}{\widehat{b}}\times {\underset{\_}{E}}_0·{\underset{\_}{E}}_{\text{pol}}$ drives a secondary Hall current $\left({\underset{\_}{j}}_{\text{H},\text{sec}}={\sigma }_H\left(\underset{\_}{\widehat{b}}\times {\underset{\_}{E}}_{\text{pol}}\right)\approx \frac{{\sigma }_H^2}{{\sigma }_P}{\underset{\_}{E}}_0\right)$. These currents flow in the same direction as the initial Pedersen drag current. The total drag current is then ${\underset{\_}{j}}_{\text{drag,tot }}={\underset{\_}{j}}_{\text{P,prim }}+{\underset{\_}{j}}_{\text{H,sec }}={\sigma }_P\left(1+\frac{{\sigma }_H^2}{{\sigma }_P^2}\right){\underset{\_}{E}}_0$. The concept described here is analogous to the Earth’s equatorial electojet current system (e.g., Baumjohann & Treumann 1996; Muralikrishna & Kulkarni 2006). In this very simplified calculation, it was shown that neglecting ${\underset{\_}{E}}_{\text{pol}}$, the anisotropic drag model underestimates both the magnetic drag force $({\underset{\_}{F}}_{\text{drag}}=\underset{\_}{j}\times \underset{\_}{B})$ and the frictional heating rate by a factor of $1+k_e^2.k_e$ is shown in Fig. 2. So if k_e ≈ 1 for p_gas < 10⁻³ bar on the dayside atmosphere, then the magnetic drag force and the frictional heating are underestimated by a factor of 2. The model with anisotropic drag includes the initial electric field allowing the calculation of the primary magnetic drag but it neglects the secondary fields to ensure current closure. Therefore, the model underestimates the magnetic drag and frictional heating in the upper dayside atmosphere.

Nonideal MHD effects. Other magnetic effects such as ambipolar diffusion have also been neglected. The nonideal MHD effects can modify $\underset{\_}{E}$ and lead to modifications of the local Joule heating. To evaluate their relative importance, which influences the evolution of the magnetic field, the prefactors of the low-frequency Ohm’s law in the neutral frame, which is derived from the momentum equations of the electrons, ions, and the neutrals, can be approximated by (Leake et al. 2014; Ballester et al. 2018; Khomenko et al. 2014; Forteza et al. 2007)

$$\begin{array}{l}{\underset{\_}{E}}^{v_n}={\eta }_O{\mu }_0\underset{\_}{j}+\frac{{\eta }_H{\mu }_0}{\left|B\right|}\underset{\_}{j}\times \underset{\_}{B}-\frac{{\eta }_A{\mu }_0}{|B{|}^2}(((\underset{\_}{j}\times \underset{\_}{B})\times \underset{\_}{B})\\ -\frac{{\eta }_H{\mu }_0}{\left|B\right|}\nabla p_e+\frac{1}{{\rho }_i{v}_{in}}\left(\nabla \left(2p_e\right)-\frac{{\rho }_i}{{\rho }_n}\nabla p_n\right)\times \underset{\_}{B},\end{array}$$(44)

where ${\underset{\_}{E}}^{v_n}=\underset{\_}{E}+{\underset{\_}{v}}_n\times \underset{\_}{B}$ is the electric field in the neutral rest frame and $\underset{\_}{j}=\nabla \times \underset{\_}{B}/{\mu }_0$ is the electric current density. The terms on the right-hand side correspond to Ohmic diffusion, Hall effect, ambipolar diffusion, Biermann battery, and pressure-gradient contributions. The diffusion coefficients, scaled to units of m² s⁻¹, are given by

$${\eta }_O=\frac{1}{{\mu }_0{\sigma }_{\Vert }}\approx \frac{\left|B\right|}{{\mu }_{0}en{k}_e},{\eta }_H\approx \frac{\left|B\right|}{{\mu }_{0}en},{\eta }_A\approx \frac{|B{|}^2}{{\rho }_i{v}_{in}{\mu }_0}.$$(45)

Here the assumptions k_e ≫ k_i, m_i ≫ m_e, n_e ≈ n_i = n, with the plasma number density, n [m⁻³], and the plasma pressure, p = p_i + p_e ≈ 2p_e [Pa], have been applied. Note that in the estimation shown here, there is no difference made between the magnetic dipole field and the induced field, so B = B_dip. Effective diffusivities were introduced for the battery and pressure-gradient terms to enable a comparison with the diffusive terms: ${\eta }_{pe}^{\text{eff }}=\frac{p_e}{\left|B\right|ne}\left[ {\text{m}}^2 {\text{s}}^{-1}\right]$, ${\eta }_{G,1}^{eff}=\frac{2p_e}{{\rho }_i{v}_{in}}\left[ {\text{m}}^2 {\text{s}}^{-1}\right]$, and ${\eta }_{G,2}^{eff}=\frac{p_n}{{\rho }_n{v}_{in}}\left[ {\text{m}}^2 {\text{s}}^{-1}\right]$. As is shown in Fig. 8, the Hall effect and Ohmic diffusion dominate in the upper dayside atmosphere (p_gas < 10⁻³ bar), with the Hall diffusivity exceeding the Ohmic diffusivity at the upper atmosphere. The other effects remain several orders of magnitude smaller than the Hall effect and Ohmic diffusion in the region of interest. Ambipolar diffusion has only a minor influence in the WASP-18 b case with a magnetic field strength of 5 G but possibly becomes relevant at p_gas < 10⁻⁴ bar compared to Ohmic diffusion and Hall effects due to the decreasing gas number density. Since η_A is dependent on |B|², increasing the magnetic field strength would enhance the effect of the ambipolar diffusion and it could become of a similar magnitude to the Ohmic diffusion at, for example, p_gas < 10⁻³ bar for |B| = 50 G. However, ambipolar diffusion is linked to the feedback of the other MHD effects as well as the configuration of the dipolar and induced magnetic fields, and therefore the effect of the ambipolar diffusion might be different. Savel et al. (2024) proposed a new method of constraining the planetary magnetic field by measuring the ambipolar velocity in the photosphere via comparing heavy ion velocity with the neutral gas velocity using high-resolution spectroscopy. Therefore, this ambipolar diffusion effect can be relevant for studying the potential spectroscopical detectability of the ion-neutral velocity in the photosphere (Soriano-Guerrero et al. 2025). However, the interpretation of the results of this method strongly relies on the underlying model of the drift velocity between ions and neutrals (Christie et al. 2025, Appendix C). Pressure-gradient terms are negligible for the conditions applied here.

Self-excited atmospheric dynamo. Another limitation of the drag parameterization is the neglect of the effect of the atmospheric dynamo. For R_M ≫ 1 (Fig. 4), atmospheric flows can produce toroidal fields comparable in strength to the imposed dipole field (Rogers & Komacek 2014). Several possible scenarios for atmospheric magnetic induction were explored by Dietrich et al. (2022), indicating that both the strength and geometry of the Lorentz force could show large variation from a static dipole configuration. In the present simulations, R_M exceeds unity at p_gas < 2 bar and peaks near R_M ≈ 70 at ∼0.3 bar, suggesting that MHD effects could be dynamically important in this region. An approximate representation of these effects could be added in the model by modifying the magnitude of the magnetic field by a first-order correction; for example, $B=B_{\text{dip }}\sqrt{1+{\left(ϵR_M\right)}^2}$, where ϵ (0 < ϵ ≤ 1) is a parametrization that captures the enhancement of magnetic coupling when atmospheric induction becomes important. Here the balance between induction and diffusion was applied to estimate the induced magnetic field, B_ind ∼ ϵR_MB_dip. ϵ can be calibrated by comparing the magnetic field profiles to MHD results (e.g., Rogers & Komacek 2014; Soriano-Guerrero et al. 2025).

Dipole magnetic field. The assumption of a dipolar planetary field remains reasonable as a first-order approximation for WASP-18 b. The morphology of the magnetic field depends on the balance between Coriolis and inertial forces in the convecting interior, characterized by the Rossby number $Ro=\frac{U}{\Omega L}$ (Christensen & Aubert 2006). Taking characteristic values of U ∼ 0.1 m s⁻¹ for the convective velocity, L ∼ 5 × 10⁷ m for the Jupiter-like shell thickness (Elias-López et al. 2025), and Ω = 7.7 × 10⁻⁵ rad s⁻¹ for the angular velocity of WASP-18 b yields Ro ∼ 3 × 10⁻⁵. This value lies below the critical threshold of Ro ∼ 0.1, consistent with a dipole-dominated dynamo regime. Higher Rossby numbers (Ro > 0.1) indicate weaker, multipolar fields. However, U and L are uncertain and model-dependent, and the dynamo region of WASP-18 b is not directly constrained. Internal heating due to tides or Ohmic dissipation may affect the thickness and structure of the dynamo shell, which potentially favors a multipolar field (see, e.g., Stanley & Bloxham 2006). While the dipole assumption is reasonable for modeling atmospheric drag, the possibility of a multipolar magnetic field cannot be ruled out.

Fig. 8 Diffusion coefficients from the Ohm’s law (Eq. (6)) averaged over all latitudes and over dayside atmosphere (−90◦< ϕ ≤ 90◦). The coefficients were calculated from the output of the ExoRad simulation with anisotropic drag for the atmosphere of WASP-18 b.

7 Summary

In this study, a physically motivated model for the anisotropic magnetic drag is presented to investigate how an anti-aligned planetary dipole magnetic field influences the atmospheric circulation of (U)HJs. This approach addresses the main research questions of (1) how magnetic coupling and the resulting drag modifies the wind circulation and hotspot shift, (2) to what extent the magnetic drag disrupts the equatorial super-rotation, and (3) how it alters heat distribution between the day- and nightside and the climate characteristics. By connecting the anisotropic drag approach to ionospheric conductivities and magnetization, the approach provides a theoretical framework for linking MHD processes to climate characteristics in UHJ atmospheres. Therefore, the magnetic drag parametrization was derived from a multi-fluid description in which electrons and ions are treated as magnetized species in steady-state momentum balance, collisionally coupled to the neutrals. The resulting ion and electron drift velocities relative to the neutrals were used to determine anisotropic drag in terms of Hall and Pedersen components acting on the neutral flow. The parametrization was implemented in the meridional and zonal momentum equations as well as in the temperature evolution equation of the cloud-free 3D GCM ExoRad, and applied to WASP-18 b. The results were compared with previously published drag treatments (uniform drag and active drag) and with simulations with an absent magnetic field influence. They demonstrate that the atmospheric circulation, in particular the presence and width of the equatorial super-rotating jet, strongly depend on the magnetic drag. The main findings are:

• Anisotropic drag primarily dissipates kinetic energy and slows the flow on the dayside, while still allowing an efficient eastward equatorial jet to be maintained on the nightside, in contrast to uniform drag that suppresses overall wind flow. The presence or absence of anisotropic drag does not affect the width of the nightside jet but does affect the horizontal velocity of the jet;

• Active drag completely decouples the dayside and night-side circulation in the upper atmosphere, leading to a slower nightside jet, a polar flow at the dayside, and a weaker eastward shift of the day–night boundary compared to the case in which the magnetic field is absent;

• Anisotropic drag changes the dayside flow in particular for p_gas ≤ 10⁻³ bar. The inclusion of the Hall term and the meridional drag component (Pedersen) leads to a deflection of the flow perpendicular to the dipole field. The result is a direct flow (albeit at a weaker level compared to uniform drag) from the dayside toward the nightside at the evening terminator at all latitudes and at the morning terminator at high latitudes (> |30^◦|);

• Identifying the wind flow direction over the terminators is thus important to identify the coupling regime (uniform: diverging flow from the substellar point, active: strong decoupling at terminator regions, anisotropic: only decoupled at equatorial morning terminator, otherwise dayside-to-nightside flow);

• Because anisotropic drag includes more drag effects, this leads to an overall stronger damping of wind flow with smaller maximum velocities compared to the active drag;

• With anisotropic drag the evening–morning terminator averaged temperature difference at p_gas = 10⁻¹ bar is stronger compared to the case without magnetic drag, with the evening terminator being 20% hotter than the morning terminator;

• Both active and anisotropic drag reshape the circulation such that two hotspots appear north and south of the equator, with smaller eastward shifts compared to the case in which magnetic drag is absent. The results show that local variations in drag forces and energy can create observable asymmetries, producing both primary and secondary hotspot regions;

• The morning-to-evening terminator gas temperature difference, ∆T_gas, at p_gas = 10⁻³ bar (Table 3), is crucial to identify day-to-night side flow patterns arising with different drag treatments. Uniform drag with radial flow exhibits almost no contrast (∆T_gas ≈ 0 K), active drag with polar flow exhibits high contrast (∆T_gas = 145 K), and anisotropic drag due to partial horizontal flow toward the nightside exhibits lower contrast (∆T_gas ≈ 100 K).

While the parametrization of the magnetic drag includes the influence of Pedersen and Hall drag, it neglects several induction processes such as dynamo action, ambipolar diffusion, and ion and electron pressure effects. These processes likely become important in the partially ionized upper atmosphere, where R_M ≳ 1, and may strongly determine both the efficiency and anisotropy of magnetic drag. Observationally, such uncertainties might impact predictions of phase curve amplitudes, hotspot offsets, transit asymmetries, and Doppler wind signatures. Because the model does not include the self-consistent evolution of the magnetic field, and thus neglects inductive feedback, it likely underestimates magnetic stresses. It is therefore not valid in atmospheric regions where R_M ≫ 1. However, for R_M ≲ 1, the parametrization provides reasonable estimates of the magnitude and direction of magnetic drag. Although the approach neglects the polarization electric field which arise to ensure global current closure, it remains a suitable and computationally efficient method of calculating the first-order effects of magnetic fields in a GCM. Understanding these effects and the extent to which they affect observations can provide a constraint on the magnetic field strengths of specific UHJs in future work.

Acknowledgements

A.B. thanks Maxim Khodachenko for useful discussions. The post-processing of GCM data has been performed with gcm-toolkit (Schneider et al. 2022b). The simulations (project id 72245) were performed on the Austrian Scientific Computing (ASC) infrastructure, in particular the Vienna Science Cluster (VSC). We thank the referee for very useful comments that improved this manuscript.

References

Appendix A Active magnetic drag treatment

Here, the active magnetic drag treatment is reviewed, in which the interaction between atmospheric winds and magnetic fields is modeled by calculating the induced currents and resulting Lorentz forces. This approach was suggested by Perna et al. (2010a) and first implemented in hot Jupiter simulations by Rauscher & Menou (2013). It was later extended to UHJs by Beltz et al. (2022b). In contrast to the uniform drag treatment, the active drag approach incorporates spatially and temporally varying diffusivity, providing a more physically realistic representation of magnetic drag in atmospheres that are only partially ionized. The model is based on the following assumptions: (a) The planetary magnetic field is a dipole aligned with the planet’s rotation axis. (b) Zonal winds dominate the atmospheric circulation. (c) The magnetic Reynolds number is small (R_M ≪ 1), so the field is not significantly distorted by magnetic advection and is diffusion-dominated. (d) Since the system is diffusion-dominated other effects such as Hall and ambipolar diffusion are neglected. Under these assumptions, the nonideal MHD induction equation simplifies to

$${\partial }_t\underset{\_}{B}=\nabla \times (\underset{\_}{v}\times \underset{\_}{B}-\eta \nabla \times \underset{\_}{B}),$$(A.1)

where $\underset{\_}{B}[ \text{T}]$ is the magnetic field, $\underset{\_}{v}\left[\frac{m}{s}\right]$ the plasma bulk velocity, and $\eta (r,\theta ,\varphi )=\frac{1}{{\mu }_0\sigma (r,\theta ,\varphi )}\left[\frac{m^2}{s}\right]$ is the magnetic diffusivity with the magnetic permeability ${\mu }_0\left[\frac{N}{A^2}\right]$ and the electrical conductivity $\sigma \left[\frac{S}{m}\right]$. Here, σ corresponds to the parallel conductivity σ_∥ in the model of Perna et al. (2010a, see also Christie et al. 2025). For the model atmosphere used here, this agrees well with the Pedersen conductivity σ_P on the dayside at higher pressure levels (p_gas > 10⁻³ bar), but σ_P starts to deviate from σ_∥ at lower pressures (p_gas ≤ 10⁻³ bar). Therefore, using σ = σ_P is more appropriate for capturing the resistive effect that occurs when the gyrofrequency exceeds the collision frequency. Perna et al. (2010a) uses the velocity of neutrals rather than the plasma bulk velocity $\underset{\_}{v}={\underset{\_}{v}}_n$, because the plasma is weakly ionized and collisions between charged particles and neutrals are frequent. The neutral gas dominates the mass and momentum budget which leads to charged particles being strongly coupled to the neutrals through collisions. The diffusive term introduces spatial gradients of η, which must be carefully treated in planetary atmospheres where conductivity varies both radially and horizontally (Perna et al. 2010a). All vector components and derivatives are expressed within the right-handed spherical coordinate system defined in Sect. 4.2. Note that this approach differs from standard formulations, in which colatitude is often used. The diffusive term is given by

$$\begin{array}{l}\nabla \times (\eta \nabla \times \underset{\_}{B})=-\eta {\nabla }^2\underset{\_}{B}+\frac{\partial \eta }{\partial r}{e}_r\times (\nabla \times \underset{\_}{B})+\frac{1}{r}\frac{\partial \eta }{\partial \theta }{\underset{\_}{e}}_{\theta }\times (\nabla \times \underset{\_}{B})\\ +\frac{1}{r\cos \theta }\frac{\partial \eta }{\partial \varphi }{\underset{\_}{e}}_{\varphi }\times (\nabla \times \underset{\_}{B})\end{array}$$(A.2)

It is useful to introduce the radial, meridional, and zonal magnetic diffusivity scale heights:

$$H_{\eta ,r}=\frac{\eta }{\partial \eta /\partial r},H_{\eta ,\theta }=\frac{r\eta }{\partial \eta /\partial \theta },H_{\eta ,\varphi }=\frac{r\cos \theta \eta }{\partial \eta /\partial \varphi }.$$(A.3)

These scales can be used to evaluate whether horizontal gradients can be neglected relative to the radial one. The calculated scale heights from the model atmosphere showed that the horizontal scale heights dominate the radial scale height (i.e., $H_{\eta ,\theta },H_{\eta ,\varphi }>H_{\eta ,r}$) and therefore horizontal diffusivity gradients in Eq. A.2 can be neglected. In order to estimate the relative importance of the advective (first term on the right-hand side of Eq. A.1) and diffusive terms (second term on the right-hand side of Eq. A.1) the magnetic Reynolds number R_M is used for dimensional analysis. The typical length scale for R_M has to be chosen carefully (see therefore, e.g., Liu et al. 2008 and Perna et al. 2010a). The pressure scale height $H_p=\frac{R_d{T}_{\text{gas }}}{g}[ \text{m}]$ is smaller than H_η,r for most of the upper model atmosphere (p_gas < 0.1 bar) and similar in magnitude in the lower atmosphere, therefore H_p is used to calculate the magnetic Reynolds number:

$$R_M=\frac{v{H}_p}{\eta }.$$(A.4)

Note that R_M may exceed unity in the model atmosphere showing that this approach is not suitable to model such hot atmospheres. When R_M ≪ 1, the system is in quasi-steady-state and time derivatives can be neglected (∂_t = 0). Assuming negligible dynamo action, the magnetic field is decomposed into a dipolar background $\left({\underset{\_}{B}}_{\text{dip }}\right)$ and an induced component (${\underset{\_}{B}}_{\text{ind }}$; see, e.g., Batygin and Stevenson 2010):

$$\underset{\_}{B}={\underset{\_}{B}}_{\text{dip}}+{\underset{\_}{B}}_{\text{ind}},\nabla \times {\underset{\_}{B}}_{\text{dip}}=0.$$(A.5)

Since the induced field is typically small when R_M ≪ 1, the induction equation is linearized around the background field:

$$\nabla \times \left(\underset{\_}{v}\times {\underset{\_}{B}}_{\text{dip }}\right)=\nabla \times \left(\eta \nabla \times {\underset{\_}{B}}_{\text{ind }}\right).$$(A.6)

Using Ampère’s law $\left({\underset{\_}{j}}_{\text{ind }}=\frac{1}{{\mu }_0}\nabla \times {\underset{\_}{B}}_{\text{ind }}\right)$, the induction equation becomes

$$\nabla \times \left(\eta {\mu }_0{\underset{\_}{j}}_{\text{ind }}\right)=\nabla \times \left(\underset{\_}{v}\times {\underset{\_}{B}}_{\text{dip }}\right).$$(A.7)

Focusing on the zonal (ϕ) component and assuming dominant horizontal currents (i.e., j_r ≪ j_θ), the dominant balance becomes

$${\mu }_0{\partial }_r\left(r\eta j_{\text{ind},\theta }\right)={\partial }_r\left(r{v}_{\varphi }{B}_{\text{dip},r}\right)+{\partial }_{\theta }\left(v_{\varphi }{B}_{\text{dip},\theta }\right).$$(A.8)

Integrating over radius from the layer of interest to the model top at radius R, an expression for the meridional current is obtained:

$$\begin{array}{l}{j}_{\text{ind},\theta }(r,\theta ,\varphi )=\frac{R\eta (R,\theta ,\varphi )}{r\eta (r,\theta ,\varphi )}{j}_{\text{ind},\theta }(R,\theta ,\varphi )\\ -\frac{1}{{\mu }_{0}r\eta (r,\theta ,\varphi )}{\displaystyle {\int }_r^{R}d}{r}^{\prime }\left({\partial }_r^{\prime }\left(r^{\prime }{v}_{\varphi }{B}_{\text{dip},r}\right)+{\partial }_{\theta }\left(v_{\varphi }{B}_{\text{dip},\theta }\right)\right).\end{array}$$(A.9)

The first term on the right-hand side is a boundary current in the uppermost modeled level and is set to zero for simplicity (Perna et ail 2010). Furthermore, it is assumed that the radial length scale of diffusivity variations (H_η(r)), is smaller than the characteristic meridional length scales of the zonal wind and magnetic field. As a result, radial gradients dominate over meridional ones, allowing to neglect $\frac{1}{r}\partial /\partial \theta$ terms in comparison to ∂/∂r in the equation above. Assume that ${\partial }_r^{\prime }\left(r^{\prime }{v}_{\varphi }{B}_{\text{dip},r}\right)$ slowly varies with radius and that the atmosphere is very thin compared to the radius of the planet, then it can be approximated by

$$\begin{array}{l}{j}_{\text{ind},\theta }(r,\theta ,\varphi )=-\frac{1}{{\mu }_{0}r\eta (r,\theta ,\varphi )}{\displaystyle {\int }_r^{R}d}{r}^{\prime }\left({\partial }_r^{\prime }\left(r^{\prime }{v}_{\varphi }{B}_{\text{dip},r}\right)\right)\\ \approx -\frac{(R-r)}{{\mu }_{0}r\eta (r,\theta ,\varphi )}{\partial }_r\left(r{v}_{\varphi }{B}_{\text{dip},r}\right)\\ \approx -\frac{v_{\varphi }}{{\mu }_0\eta (r,\theta ,\varphi )}{B}_{\text{dip},r}.\end{array}$$(A.10)

The resulting Lorentz force appears in the zonal momentum equation (where for the mass density ρ the neutral gas mass density ρ_n was assumed) as

$$\rho \frac{d}{dt}{v}_{\varphi }\approx -j_{\text{ind},\theta }{B}_{\text{dip},r}.$$(A.11)

Substituting the idealized expression for the meridional current j_ind,θ = −v_ϕB_dip,r/(µ₀η) gives a drag force proportional to the zonal wind velocity:

$$\rho \frac{d}{dt}{v}_{\varphi }\approx \frac{1}{\eta {\mu }_0}{v}_{\varphi }{B}_{\text{dip},r}^2.$$(A.12)

Assuming a dipolar field with radial component $B_{\text{dip},r}=2B_0\sin \theta \frac{R_p^3}{r^3}$ the equation becomes

$$\rho \frac{d}{dt}{v}_{\varphi }\approx \frac{1}{\eta {\mu }_0}{v}_{\varphi }{\left(2B_0\sin \theta \frac{R_p^3}{r^3}\right)}^2.$$(A.13)

Furthermore, the magnetic field strength is assumed to be independent of radius, simplifying the radial variation of τ_drag. This yields a magnetic drag timescale of the form

$${\tau }_{\text{drag }}=\frac{{\mu }_0\rho \eta }{B_{\text{dip},r}^2}.$$(A.14)

τ_drag³ is used in this study for the "active drag" treatment. Note that this approach differs from the parametrized model used in Rauscher & Menou (2013) and Beltz et al. (2022b), where the drag force is based on local flow geometry, magnetic field strength, and diffusivity via

$${\tau }_{\text{drag, Perna }}=\frac{{\mu }_0\rho \eta }{B_0^{\text{2}}\left|\sin \theta \right|}.$$(A.15)

The latitudinal dependence of the radial magnetic flux density is not fully maintained in their formulation of the drag timescale and interprets drag as being proportional to the component of the wind perpendicular to the magnetic field lines. At the equator (θ = 0^◦), the magnetic drag vanishes (τ_drag = τ_drag,Perna → ∞).

The energy equation is also modified by the kinetic energy which is lost due to the magnetic drag and returned to the atmosphere in form of localized ohmic heating with

$$q_{\text{fric,active }}=v_{\varphi }^2/{\tau }_{\text{drag }}.$$(A.16)

While the active drag approach is more physically motivated than the empirical Rayleigh drag, it is still an approximation. In particular, the assumption of dominant zonal winds, weak induced fields, and R_M ≪ 1 may not hold in all regions of the atmosphere as shown for WASP-18 b. However, this model was still used for simulations of (U)HJs where the condition of R_M ≪ 1 was not given (e.g., Beltz et al. 2022b; Coulombe et al. 2023).

A.1 Thermal ionization

Due to extreme temperatures on the dayside (> 3000 K), thermal collisions are sufficiently energetic to ionize atoms and molecules and form an ionosphere. Thermal ionization rather than photoionization is the dominant process creating strong ionospheric effects in our model, influencing the atmospheric dynamics and magnetic interactions. However, Koskinen et al. (2014) showed that photoionization can dominate over thermal ionization on the dayside atmosphere above the 100 mbar level. In order to determine the thermal ionization fraction x_e,Saha in a gas in LTE, models of Menou (2012), Rauscher & Menou (2013) and others adopted the approximate solution of the Saha equation from Sato (1991). Furthermore, they assumed that the ionization fraction is small (x_e,Saha ≪ 1). Therefore, only single ionization from the ground state is considered (Menou 2012). The Saha equation for a plasma in thermodynamic equilibrium is given by

$$\frac{n_{i,j}^{+}{n}_e}{n_{n,j}}={\left(\frac{m_e{k}_{B}T}{2\pi {\hslash }^2}\right)}^{\frac{3}{\text{2}}}\frac{2Z_{i,j}}{Z_{n,j}}\exp \left(\frac{-I_{i,j}}{k_{B}T}\right),$$(A.17)

with the electron mass m_e [kg], the Boltzmann’s constant $k_B\left[\frac{J}{K}\right]$, the gas temperature T_gas [K], the reduced Plank’s constant ℏ[Js], the electron number density n_e [m⁻³], the ionized number density of the jth chemical species $n_{i,j}^{+}\left[{\text{m}}^{-3}\right]$, the atom number density n_n,j [m⁻³] of the jth species, their partition functions Z_i,j and Z_n,j, and the first ionization potential I_i,j [J] of the jth species. For the calculation of x_e the first 28 elements from the periodic table (hydrogen to nickel) are taken into account (Menou 2012; Rauscher & Menou 2013). The first ionization potential I_i,j ranges from 7×10⁻¹⁹ J (for K) to 4×10⁻¹⁸ J (for He) The partition functions Z_i,j and Z_n,j depend on the temperature in a complex way (Sato 1991). For simplicity, Z_i,j and Z_n,j are set to unity according to Menou (2012) since the specific atmospheric composition is unknown. Menou (2012) and Rauscher & Menou (2013) use a sum of approximate Saha expressions across multiple elements. They evaluate ionization fractions for each element $\left(x_{e,\text{ Saha },j}=n_{i,j}^{+}/n_{\text{tot },\text{j}}\right)$ independently using uncoupled Saha equations. The total number density of the jth element is given by $n_{\text{tot },\text{j}}=n_{n,j}+n_{i,j}^{+}$. They assume that only the jth species contributes to the electrons and therefore use $n_e=n_{i,j}^{+}$, which is not mathematically consistent, but here a useful simplification in a weakly ionized plasma in order to derive x_e,Saha. The Saha equation (Eq. A.17) can be described in terms of the ionization fraction for each ion species: