Magnetohydrodynamic modelling of star–planet interaction and associated auroral radio emission

ABSTRACT

We present calculations of auroral radio powers of magnetized hot Jupiters orbiting Sun-like stars, computed using global magnetohydrodynamic (MHD) modelling of the magnetospheric and ionospheric convection arising from the interaction between the magnetosphere and the stellar wind. Exoplanetary auroral radio powers are traditionally estimated using empirical or analytically derived relations, such as the radiometric Bode’s law (RBL), which relates radio power to the magnetic or kinetic energy dissipated in the stellar wind–planet interaction. Such methods risk an oversimplification of the magnetospheric electrodynamics giving rise to radio emission. As the next step towards a self-consistent picture, we model the stellar wind–magnetosphere–ionosphere coupling currents using a 3D MHD model. We compute electron-cyclotron maser instability-driven emission from the calculated ionospheric field-aligned current density. We show that the auroral radio power is highly sensitive to interplanetary magnetic field (IMF) strength, and that the emission is saturated for plausible hot Jupiter Pedersen conductances, indicating that radio power may be largely independent of ionospheric conductance. We estimate peak radio powers of 10¹⁴ W from a planet exposed to an IMF strength of 10³ nT, implying flux densities at a distance of 15 pc from Earth potentially detectable with current and future radio telescopes. We also find a relation between radio power and planetary orbital distance that is broadly consistent with results from previous analytic models of magnetosphere–ionosphere coupling at hot Jupiters, and indicates that the RBL likely overestimates the radio powers by up to two orders of magnitude in the hot Jupiter regime.

MHD, methods: numerical, planets and satellites: aurorae, planets and satellites: detection, planets and satellites: magnetic fields, radio continuum: planetary systems

1 INTRODUCTION

The direct detection of exoplanets in large parts of the electromagnetic spectrum is hindered by the high-luminosity contrast ratio between the star and planet. Evidence from the Solar system planets, however, indicates that the radio waveband presents a luminosity ratio much more conducive to direct detection, with non-thermal emission from Jupiter of similar intensity to Solar radio bursts (Zarka 2007). Historically, the search for exoplanetary radio emission has focused primarily on Jupiter-like exoplanets in close orbit (3–10 stellar radii) around their parent star. Such planets, given the moniker ‘hot Jupiters’, have been targeted by several studies examining the possibility of detectable auroral radio emission (e.g. Farrell, Desch & Zarka 1999; Zarka et al. 2001; Grießmeier et al. 2004; Lazio et al. 2004; Griessmeier et al. 2005, 2007; Zarka 2007; Jardine & Cameron 2008; Nichols 2011, 2012; Hallinan et al. 2012; Vidotto & Donati 2017). Although the search for radio emission has not been limited exclusively to hot Jupiters, the plausibly strong magnetic fields and intense stellar wind conditions present at these objects is thought to be favourable to the generation of radio emission through star–planet interaction. Auroral radio emission from exoplanets is envisaged to be generated by the electron-cyclotron maser instability (ECMI), the same mechanism responsible for driving radio emission from auroras at magnetized Solar system planets (Wu & Lee 1979; Treumann 2006; Imai et al. 2008; Lamy et al. 2011). This intense, coherent electromagnetic radiation is emitted at the local cyclotron frequency, and therefore Jupiter-like exoplanets are the prime candidates for directly detectable emission, since their potentially high intrinsic magnetic field strengths (∼B_Jup) are required to produce emission above the Earth’s ionospheric cut-off frequency of ∼10 MHz.

Many previous studies estimating expected exoplanetary radio emission employ the radiometric Bode’s law (RBL), an empirical scaling relation extrapolated from Solar system measurements between incident Poynting or kinetic energy flux and emitted radio power (Farrell et al. (Farrell et al. 1999; Zarka et al. 2001; Lazio et al. 2004; Zarka 2007). Estimates of the radio power from hot Jupiter auroras based on the RBL range between 10¹⁴ and 10¹⁶ W, up to five orders of magnitude stronger than Jupiter’s equivalent ECMI-driven emission (Zarka 2007), implying radio fluxes that may be detectable with the next generation of radio telescopes. Two primary processes are assumed to mediate radio emission driven by star–planet interaction: either Alfvén waves, such as mediates emission in the sub-Aflvénic Io–Jupiter interaction; or magnetic reconnection, such as occurs at Earth’s magnetosphere. In considering radio emission propagated by Alfvén waves, Saur et al. (2013) and Turnpenney et al. (2018) estimated total radiated energy fluxes from exoplanets of up to 10¹⁹ W, translating to radio powers of 10¹⁷ W assuming an ECMI efficiency of 1 per cent. Nichols & Milan (2016), using an analytic model, considered an Earth-type Dungey cycle process of magnetic reconnection, computing ionospheric field-aligned currents (FACs) and resulting radio powers arising from ionospheric convection. They found that saturation of the convection-induced electric potential limited the dissipated power, and predicted auroral radio emission from hot Jupiters ∼2 orders of magnitude smaller than RBL-based predictions. In this paper, we use a global 3D magnetohydrodynamic (MHD) model to calculate the magnetosphere–ionosphere coupling currents arising from the interaction between Sun-like stars and hot Jupiters, and hence estimate the resulting auroral radio emission. By using a numerical global MHD model, this study extends the analytic work of Nichols & Milan (2016) to enable a self-consistent calculation of the FACs at exoplanets. The numerical model, as used in this study, takes a set of input boundary conditions and generates a 2D map of the ionospheric FAC density distribution, from which radio power is calculated in post-processing. The overall power is expected to be strongly influenced externally, by the interplanetary magnetic field (IMF) strength, and internally, by the ionospheric Pedersen conductance. Therefore, we run simulations across a broad range of these two parameters, comparing the results with those of Nichols & Milan (2016).

This paper begins with an overview of the theoretical background relevant to the MHD model, along with the formulation used to compute radio powers from the FAC output of the model. Then follows a description of the method employed, before results of the modelling work are presented and discussed. A case intermediate between Earth and hot Jupiter exoplanets is first studied, before cases more appropriate to exoplanets are examined.

2 THEORETICAL BACKGROUND

2.1 Magnetohydrodynamic model

Global 3D MHD simulations, based on first principle physics, are a powerful tool for modelling the dynamics and evolution of magnetic fields and plasma flows in space weather and astrophysical phenomena. In use since the 1980s, early models relied on techniques such as finite-difference methods to solve a system of discretized MHD equations (Van Leer 1979; Leboeuf et al. 1981; Wu, Walker & Dawson 1981; Brecht et al. 1982). Computational solutions of MHD equations require discretization of the MHD equations, which inherently introduces errors into the solution. Modern MHD models use advanced solution techniques which improve the efficiency of the solution and minimize such discretization errors (Gombosi et al. 2004). Such methods rely on approximations to solve a Riemann problem, the form of initial value problem presented in MHD numerical analysis over a finite volume. The 3D MHD solver used in this work is the ‘Block Adaptive Tree Solarwind Roe Upwind Scheme’ (bats-r-us) software first outlined by Powell et al. (1999), and developed at the University of Michigan. The computational scheme of bats-r-us is based on the same elements used in many state-of-the-art MHD codes, and this section describes each of those elements in the scheme.

A governing set of 3D ideal MHD equations is first defined. Various forms of these equations are expressible, and the form chosen is dictated by factors which will ultimately aid in the computational solution of these equations. The set of ideal MHD equations used inbats-r-us are expressed in a gasdynamic conservative form (Gombosi et al. 2004):

$$\begin{eqnarray*} \frac{\partial \rho }{\partial t} + \nabla \cdot (\rho \boldsymbol{u}) = 0 \end{eqnarray*}$$

(1)

$$\begin{eqnarray*} \frac{\partial (\rho \boldsymbol{u})}{\partial t} + \nabla \cdot (\rho \boldsymbol{u u} + p {\bf I})= \frac{1}{\mu _0} \boldsymbol{j} \times \boldsymbol{B} \end{eqnarray*}$$

(2)

$$\begin{eqnarray*} \frac{\partial \boldsymbol{B}}{\partial t} + \nabla \times \boldsymbol{E} = 0 \end{eqnarray*}$$

(3)

$$\begin{eqnarray*} \frac{\partial E_{\mathrm{gd}}}{\partial t} + \nabla \cdot \left[ \boldsymbol{u} (E_{\mathrm{gd}} + p) \right] = \frac{1}{\mu _0} \boldsymbol{u} \cdot (\boldsymbol{j} \times \boldsymbol{B}), \end{eqnarray*}$$

(4)

where I is the identity matrix and the total gasdynamic energy E_gd, is given by

$$\begin{eqnarray*} E_{\mathrm{gd}} = \frac{1}{2} \rho u^2 + \frac{1}{\gamma - 1} p, \end{eqnarray*}$$

(5)

where γ is the ratio of the specific heats. This equation set contains an expression for the conservation of mass (equation 1), conservation of momentum (equation 2), an expression of Faraday’s law (equation 3), and an energy equation (4). These partial differential equations are manipulated into a non-dimensional, symmetrizable form, the full details of which can be found in Powell et al. (1999). A computational domain is divided into cells over which the MHD equations are integrated. From the symmetrizable form Powell (1994) showed that it is possible to derive a Roe-type approximate Riemann solver for the 3D equations. First described by Roe (1981), this is a method for solving partial differential equations by estimating the flux at the interface between two computation cells in some finite-volume discretized domain. Such solvers are required in magnetohydrodynamics, since iterative solutions are costly, and therefore approximations must be made.

The computation domain is divided into a grid of Cartesian cells, and the cells are structured into blocks typically consisting of between 4 × 4 × 4 and 12 × 12 × 12 individual cells. The block-adaptive technology of bats-r-us allows the computational grid to be adapted based on pre-specified physical criteria, such that blocks can be refined in regions where interesting physical features emerge. Adaptive mesh refinement is extremely effective when the problem being treated contains disparate length-scales, and also removes any initial grid-based bias in the solution.

The Space Weather Modelling Framework (swmf) is a software package which integrates several different physics domains extending from the solar surface to the planetary upper atmosphere (Tóth et al. 2005). bats-r-us is used in a number of these components where MHD solutions are required, and physically meaningful combinations of components can be coupled together to study a wide variety of space weather events and phenomena (Tóth et al. 2012). While developed originally to study Sun–Earth events, the swmf has since been adapted and applied to other Solar system planets, satellites, and comets (e.g. Huang et al. 2016; Jia & Kivelson 2016; Tóth et al. 2016; Jia et al. 2019), and may reasonably be used for the study of extrasolar astrophysical systems where the physics domains are appropriate, with some adaptation where may be required.

2.2 Field-aligned current and radio power

This study utilizes the global magnetosphere (GM) component of the swmf, coupled with the ionospheric electrodynamics (IE) component (Ridley, Gombosi & DeZeeuw 2004). The GM domain constructs the magnetic environment and plasma dynamics around the planet, and contains features such as the bow shock, magnetopause, and magnetotail. Upstream boundary condition for the GM component can be obtained from coupling with the inner heliosphere component of the swmf, but in this work are simply input into the model based on reasonable values, as will be discussed below. Currents from the GM component are mapped down along magnetic field lines to provide the FAC boundary conditions for the IE component. The domain of the IE component is a height-integrated spherical surface. Formally, this component is a two-dimensional electric potential solver, which computes conductances and particle precipitation from FACs. The process can be summarized as follows: (1) FACs are calculated by |$\nabla \times \boldsymbol{B}$| at 3.5 R_P, a value also employed by Ridley et al. (2004), where |$\boldsymbol{B}$| is the local magnetic field; (2) the currents are then mapped down along field lines to a nominal ionospheric altitude of ∼110 km using the planetary dipolar field, and are scaled by the ratio B_I/B_3.5, where B_I and B_3.5 are the magnetic field strengths at the ionosphere and 3.5 planetary radii, respectively. (3) Next, a height-integrated ionospheric conductance map is generated and the electric potential is calculated, which is then mapped out along magnetic field lines to the simulation’s inner boundary at 2.5 R_P, where flow velocities and electric fields are prescribed (Ridley et al. 2001, 2004).

Of the several variables output from the IE component, this work is principally concerned with the ionospheric FAC density j_||, and the cross-polar cap potential (CPCP) Φ_CPC outputs. By integrating the total upward or downward FAC density output from the IE component over one hemisphere, the total current, I_tot, flowing into or out of the ionosphere is obtained by

$$\begin{eqnarray*} I_{\mathrm{tot}} = \int ^{2 \pi }_0 \int ^{\pi /2}_0 R_{\mathrm{P}}^2 j_{||} \sin \theta\, \mathrm{d}\theta \,\mathrm{d}\phi , \end{eqnarray*}$$

(6)

where θ and ϕ are the conventional spherical coordinates of colatitude and azimuth, respectively. Total auroral power is also calculated in post-processing by integrating precipitating electron energy flux E_f over one hemisphere:

$$\begin{eqnarray*} P_{\mathrm{tot}} = \int ^{2 \pi }_0 \int ^{\pi /2}_0 R_{\mathrm{P}}^2 E_{\mathrm{f}}(\theta , \phi) \sin \theta \,\mathrm{d} \,\theta \mathrm{d}\phi . \end{eqnarray*}$$

(7)

Qualitatively, the precipitating electron energy flux is the kinetic energy carried by the downward-flowing electrons associated with the upward FAC. The maximum FAC which can be carried by unaccelerated electrons in an isotropic Maxwellian velocity space distribution is given by

$$\begin{eqnarray*} j_{\Vert i0} = e n \left(\frac{W_{\mathrm{th}}}{2 \pi m_{\mathrm{e}}} \right)^{1/2}, \end{eqnarray*}$$

(8)

where e and m_e are the electron charge and mass, respectively, and W_th and n are the thermal energy and number density of the magnetospheric electron source population, respectively, for which we employ canonical Jovian values established from Voyager measurements of W_th = 2.5 keV and n = 0.01 cm⁻³ throughout this work (e.g. Scudder, Sittler & Bridge 1981). We note that these values may differ significantly at hot Jupiters, although no reliable estimates exist at present. The source plasma number density is in relation to the evacuated auroral field lines, and is therefore expected to be much reduced from the ambient plasma density. Qualitatively, the effect of varying these parameters, however, will be to increase the precipitating electron energy flux and auroral power with increasing plasma thermal energy, and with decreasing plasma density. The Jovian values are employed here in the absence of any information for realistic values at hot Jupiters, although future work may investigate a range of these parameter values to determine quantitatively the effect on auroral radio emission. In general, the FACs will be larger than can be carried solely by unaccelerated electrons, and therefore must be driven by a field-aligned electric potential. In common with previous works computing intense exoplanetary and ultracool dwarf auroral radio emissions (Nichols 2011, 2012; Nichols et al. 2012; Turnpenney et al. 2017), we employ Cowley’s (2006) relativistic extension of Knight’s (1973) current–voltage relation for parallel electric fields given by

$$\begin{eqnarray*} \left(\frac{j_{\Vert i}}{j_{\Vert i0}} \right) = 1 + \left(\frac{e \Phi _{\mathrm{{min}}}}{W_{\mathrm{th}}} \right) + \frac{\left(\frac{(e \Phi _{\mathrm{min}}}{W_{\mathrm{th}}}\right)^2}{2 \left[ \left(\frac{m_{\mathrm{e}} c^2}{W_{\mathrm{th}}} \right)+ 1\right]}, \end{eqnarray*}$$

(9)

where Φ_min is the minimum field-aligned voltage required to drive the current j_‖i in the ionosphere, and c is the speed of light in a vacuum. The corresponding precipitating electron energy flux is given by

$$\begin{eqnarray*} \frac{E_{\mathrm{f}}}{E_{\mathrm{f}0}} = 1 +\left(\frac{e\Phi _{\mathrm{min}}}{W_{\mathrm{th}}} \right) + \frac{1}{2} \left(\frac{e\Phi _{\mathrm{min}}}{W_{\mathrm{th}}} \right)^2 + \frac{ \left(\frac{e\Phi _{\mathrm{min}}}{W_{\mathrm{th}}} \right)^3 }{ 2 \left[ 2 \left(\frac{m_{\mathrm{e}} c^2}{W_{\mathrm{th}}} \right) + 3 \right]}, \end{eqnarray*}$$

(10)

where E₀ is the maximum unaccelerated electron energy flux, corresponding to equation (8), given by

$$\begin{eqnarray*} E_{\mathrm{f}0} = 2 n W_{\mathrm{th}} \left(\frac{W_{\mathrm{th}}}{2 \pi m_{\mathrm{e}}} \right)^{1/2}. \end{eqnarray*}$$

(11)

Assuming, in common with observations of Jupiter and Saturn, an ECMI efficiency of 1 per cent (Lamy et al. 2011), the emitted auroral radio power is

$$\begin{eqnarray*} P_\mathrm{r} = \frac{P_{\mathrm{tot}}}{100}. \end{eqnarray*}$$

(12)

Finally, the spectral flux density is calculated by

$$\begin{eqnarray*} F_{\mathrm{r}} = \frac{P_{\mathrm{r}}}{1.6 s^2 \Delta \nu }, \end{eqnarray*}$$

(13)

where s is the distance to the exoplanetary system from the Earth, Δν is the radio emission bandwidth, and a beaming angle of 1.6 sr is assumed on the basis of Jupiter’s observed ECMI emission beaming angle (Zarka, Cecconi & Kurth 2004). Since cyclotron maser emission is generated at the local electron-cyclotron frequency, the bandwidth is the difference between the magnetic field strength at the location of the field-aligned potential and the ionosphere. This formulation assumes that the potential is located high enough up the field line from the ionosphere that the field strength there is much smaller than the field strength in the ionosphere. This assumption is valid beyond a few planetary radii owing to the r⁻³ dependence for a dipole planetary field. We therefore assume that the bandwidth Δν is equal to the electron-cyclotron frequency in the polar ionosphere and hence given by

$$\begin{eqnarray*} \Delta \nu = \frac{e B_\mathrm{ i}}{2 \pi m_{\mathrm{e}}}, \end{eqnarray*}$$

(14)

where B_i is the ionospheric magnetic field strength.

3 METHODOLOGY

Input parameters were chosen for the swmf to simulate a magnetized hot Jupiter interacting with the IMF and stellar wind of a solar-type star. A planet of Jupiter mass (⁠|$1.9 \times 10^{27} \, \mathrm{kg}$|⁠), radius (⁠|$71,492 \, \mathrm{km}$|⁠), and equatorial magnetic field strength (⁠|$426,400\, \mathrm{nT}$|⁠) was specified, with a dipole magnetic field aligned with the planetary spin axis. The orientation of the magnetic field is opposite to that of Jupiter, i.e. the planetary field is pointing northward at the equator. Each simulation was run with a planetary plasma density of |$10^7 \, \mathrm{cm^{-3}}$|⁠, and a temperature of |$8,000\, \mathrm{K}$| at the inner boundary of the simulation, values consistent with those based on modelling of atomic hydrogen, heavy atoms, and ions surrounding hot Jupiters (Muñoz 2007; Koskinen et al. 2013; Shaikhislamov et al. 2016). The incident plasma velocity of the simulations was set to |$250\, \mathrm{km\, s^{-1}}$|⁠. This value represents the impinging plasma flow velocity taking into account both the stellar wind outflow and the Keplerian orbital velocity of the planet in the hot Jupiter regime. We note that although in reality the incident plasma velocity will be predominantly azimuthal due to the high orbital velocity of planet, in this work this incident plasma velocity is prescribed as an input in the model, such that in the results that follow (i.e. Figs 1 and 2) the X-axis is aligned with the incident plasma velocity. This transposition is purely for convenience of modelling, and does not affect the resulting auroral radio power calculations.

Since the stellar wind plasma temperature close to the star in the orbital distance regime of hot Jupiters is approximately the coronal value, a stellar wind plasma temperature input of |$2\, \mathrm{MK}$| was employed throughout this study. The stellar wind density was set to |$10^4\, \mathrm{cm^{-3}}$| in all simulations, a value appropriate for the hot Jupiter regime based on analytic modelling by Nichols & Milan (2016).

In the work presented here, simulations were run for an entirely southward IMF orientation, with B_sw values ranging from 1 to |$10^6\, \mathrm{nT}$|⁠, and with a constant ionospheric Pedersen conductance between 1 and |$10^4 \, \mathrm{mho}$|⁠. The assumption of a constant Pedersen conductance provides a reasonable first-order approximation to the global average. A uniform zero Hall conductance was used throughout, which along with a constant non-zero Pedersen conductance forms the simplest ionospheric conductance model to approximate a realistic magnetospheric configuration. Such a configuration has been used as a standard ionospheric model in many previous MHD simulations (e.g. Fedder & Lyon 1987; Jia, Kivelson & Gombosi 2012a; Jia et al. 2012b).

The simulations were run on a 3D Cartesian computational grid of |$256 \, \times 256 \, \times 256 \, R_{\mathrm{P}}$|⁠. Although bats-r-us may be run in a time-accurate mode, since this work focuses on hypothetical events, we instead used iterative local time-stepping, in which each computation cell takes different time-steps and the simulation progresses for a fixed number of iterations to converge on a steady-state solution. At the inner boundary of each simulation the grid was initially highly resolved near the planet with cells of 1/8 R_P in size, while the remainder of the grid is incrementally more coarsely resolved moving out further from the planet. This initial resolution is entirely geometric, i.e. it is not based on any physical criteria, but rather is determined based on expectations of where interesting regions of the solution will emerge requiring high resolution. Each simulation was allowed to run for 3000 iterations before the grid was refined using the adaptive mesh refinement facility within the MHD code. Refinement added 10 per cent more cells in regions of large ∇P and |$\nabla \times \boldsymbol{B}$| before the simulation resumed running. Refinement based on these criteria was performed every 300 iterations up to a total of 6000 iterations, giving a final grid containing approximately 15 million cells. After the final refinement the simulation was then allowed to run for a total of 50 000 iterations, by which point the solution had reached an approximate steady state.

To validate the approach described above, the model was first tested by replicating earlier results of Ridley et al. (2004) and Ridley (2007), the details of which can be found in Appendix  A. This work builds on those studies with vastly increased IMF strengths and Pedersen conductances.

4 RESULTS

4.1 Magnetospheric structure

A series of simulation results of the magnetospheric field morphology and plasma density is shown in Fig. 1 for IMF B_sw values ranging from 1 to |$10^5\, \mathrm{nT}$|⁠. In each run, the Pedersen conductance was initially |$\Sigma _{\mathrm{P}} = 10^4 \, \mathrm{mho}$|⁠, representing the highly conductive ionospheres expected at hot Jupiters. The remainder of the input parameters were as stated in the previous section. The magnetic field lines are traced in the Y = 0 plane, with the incident plasma flowing from right to left. Note that in these plots the field line spacing does not necessarily represent magnetic field strength. The lower end of the IMF range represents conditions analogous to the IMF experienced at the Earth. In the region typically associated with hot Jupiters, i.e. 3–10 stellar radii, an IMF of approximately 10³–|$10^4 \, \mathrm{nT}$| is expected (Nichols & Milan 2016). Higher IMF strengths have been examined which represent either planets orbiting the star extremely closely (<3 stellar radii), or planets orbiting stars with exceptionally strong magnetic fields. For instance, at a distance of 10 R_*, assuming a predominantly radial field such that B_sw ∝ 1/r², a 10⁶ nT IMF equates to a star with a surface magnetic field strength of 10⁸ nT, approximately three orders of magnitude great than the solar magnetic field.

Figure 1.

Plasma density with traces showing magnetic field lines in the X–Z plane (in GSM co-ordinates). The planet is situated at Z = X = 0 in each plot, and the stellar wind is flowing from right to left. The co-ordinates are given in units of planetary radii. Note that the density scale differs between the individual plots.

Fig. 1(a) shows that when the planet is exposed to |$B_{\mathrm{sw}} = 1 \, \mathrm{nT}$| the magnetosphere formed is similar to that at the Earth, i.e. with a clearly visible magnetotail and an apparent magnetosphere on the sub-stellar wind side of the planet. As the IMF is increased |$B_{\mathrm{sw}} = 10^2\, \mathrm{nT}$| and |$B_{\mathrm{sw}} = 10^3\, \mathrm{nT}$|⁠, Figs 1(c) and (d) show the lobes of the tail opening as the upstream Alfvén Mach number becomes lower, and the formation of Alfvén wings draped across the planet becomes apparent. At |$B_{\mathrm{sw}} = 10^4\, \mathrm{nT}$|⁠, the flow has become sub-Alfvénic, and the Alfvén wings are formed at a large angle from the equatorial plane. As the IMF is increased further to |$B_{\mathrm{sw}} = 10^5 \, \mathrm{nT}$| the planetary field is dwarfed by the stellar field, and essentially presents no perturbation to the overwhelming IMF.

Fig. 2 shows the magnetospheric field line topology and plasma density in close proximity to the planet for the same simulations as Fig. 1. These plots reveal the compression of the sub-stellar wind side magnetosphere due to the IMF pressure as B_sw is increased, as well as escape of the dense planetary plasma along open field lines. Note that in Figs 2(e) and (f) the substellar magnetosphere has collapsed below the |$2.5 \, R_{\mathrm{p}}$| inner boundary of the simulation.

Figure 2.

As for Fig. 1 but for a smaller scale to show the magnetospheric topology in close proximity to the planet.

4.2 Ionospheric electrodynamics

Fig. 3 shows a set of ionospheric plots of FAC density and CPCP at a constant Pedersen conductance value of Σ_P = 10⁴ mho, and with B_sw range of 1–10⁵ nT. In each FAC density plot positive upward current is indicated in red, and negative downward current by blue. The morphology and magnitudes of both quantities within these plots vary as B_sw is increased. At low IMF strengths the auroral FAC density is somewhat irregular and diffuse in shape, and is situated at higher colatitudes of ∼20–30^○, slightly displaced towards the sub-stellar wind side of the planet. The upward FAC density peaks for B_sw = 10³ nT at |$199 \,\, \mu \mathrm{A} \, \mathrm{m^{-2}}$|⁠. At B_sw = 10⁴ nT the auroral FAC structure is a narrow oval at ∼30^○ colatitude, with the peak upward FAC density of |$39.97 \,\, \mu \mathrm{A} \, \mathrm{m^{-2}}$|⁠. As B_sw is increased to 10⁵ the magnitude of the upward current density falls significantly to 8.67 |$\mu \mathrm{A} \, \mathrm{m^{-2}}$|⁠.

Figure 3.

Ionospheric maps of (a) FAC density and (b) CPCP for a fixed ionospheric Pedersen conductance of Σ_P = 10⁴ mho, and B_sw ranging from 1 to 10|$^5 \, \mathrm{nT}$|⁠. The minimum and maximum values of the two parameters are indicated below each plot. In each case, contours of absolute values are shown on a logarithmic scale, with negative values coloured blue and positive values coloured red.

A similar saturation and turnover is seen in the CPCP results in Fig. 3(b). From a CPCP potential of ∼3 keV for B_sw = 1 nT, the potential peaks around 6 keV at B_sw = 10³ nT, before a substantial drop to Φ_CPC = 0.457 kV at B_sw = 10⁵ nT. Saturation of the CPCP, first discussed by Hill, Dessler & Wolf (1976), was subsequently tested by Siscoe et al. (2002) against results from MHD simulations, and was modelled in terms of the stellar wind parameters. Various interpretations have been offered for the phenomenon of CPCP saturation. Siscoe et al. (2002) interpreted the saturation as a weakening of the planetary field at the magnetopause due to the field arising from region 1 currents, thus limiting the rate at which reconnection occurs on the sub-stellar wind side of the planet. Alternatively, Kivelson & Ridley (2008) argued that saturation occurs when the solar or stellar wind impedance is greater than the ionospheric impedance, i.e. when the Pedersen conductance Σ_P dominates the Alfvén conductance Σ_A, causing a partial reflection of Alfvén waves incident on the ionosphere from the solar or stellar wind. The available magnetospheric convection potential is given by (Nichols & Milan 2016)

$$\begin{eqnarray*} \Phi _{\mathrm{m}} = \chi R_{\mathrm{mp}} E_{\mathrm{sw}}, \end{eqnarray*}$$

(15)

where E_sw is the stellar wind motional electric field, and χ is the fraction of the magnetopause standoff distance R_mp which constitutes the stellar wind reconnection channel. Observations for Earth determine a value of χ ≈ 0.5 (Milan et al. 2004), and this value is therefore also employed here. In the reflected Alfvén wave interpretation of Kivelson & Ridley (2008), the electric potential across the ionosphere is given by

$$\begin{eqnarray*} \Phi _{\mathrm{CPC}} = \frac{2 \xi \Phi _{\mathrm{m}} \Sigma _{\mathrm{A}}}{\Sigma _{\mathrm{P}} + \Sigma _{\mathrm{A}}}, \end{eqnarray*}$$

(16)

where the width of the interaction channel, specified by Kivelson & Ridley (2008) as 0.1πR_mp is accounted for by the factor ξ = (0.1π/χ). Hence, when Σ_P ≫ Σ_A the CPCP tends towards

$$\begin{eqnarray*} \Phi _{\mathrm{sat}} = \frac{2 \xi \Phi _{\mathrm{m}} \Sigma _{\mathrm{A}}}{\Sigma _{\mathrm{P}}} \end{eqnarray*}$$

(17)

and saturation occurs. For a fixed Pedersen conductance, increased IMF strength leads to a reduced Alfvén conductance, and thus a decrease in CPCP beyond saturation as observed in the results in Fig. 3(b). The saturation effect is also influenced by the decreasing magnetopause standoff distance, but the dominant contributing factor is decreased Alfvén conductance, since the sub-stellar wind side magnetosphere is completely eroded under high IMF strengths.

A notable feature of the plots in Figs 3(a) is the absence of strong region 2 currents. This is an artefact of the MHD code, with several possible causes suggested by Ridley et al. (2001). Since region 2 currents are generated close to the inner boundary of the model, a high resolution is required to produce currents that are even a fraction of the region 1 currents. Increasing the resolution would increase the time taken for simulations to run, and is therefore a trade-off that must be made in consideration of running the models in a timely fashion. Another option that should achieve the same result is to move the inner boundary to a lower value (e.g. from |$2.5 \, R_{\mathrm{P}}$| to |$1.5\, R_{\mathrm{P}}$|⁠). As with increasing the resolution, this solution also increases the run time of the simulations. Another possible cause of the weak region 2 currents is that gradient and curvature drifts of particles at different energies is not addressed by bats-r-us. The pressure gradient that results from the reconfiguration of the plasma by these drifts may be a source of region 2 currents.

4.3 Responses to variable IMF strength and Pedersen conductance

To fully understand the response of the global simulations to the key driving factors, Fig. 4 shows plots of CPCP, total current, maximum ionospheric FAC density, and radio power, versus both B_sw and Pedersen conductance. In Figs 4(a)–(d), the Pedersen conductance was fixed at Σ_P = 10⁴ mho, since Koskinen et al. (2010) showed that hot Jupiters likely possess highly conductive ionospheres (Σ ≈ 10⁴–10⁵ mho), and the modelled parameters are plotted as a function of B_sw from 1 to 10⁶ nT. In Figs 4(e)–(h), a fixed IMF value of B_sw = 10⁴ nT was used, and Σ_P was varied. The modelled CPCP, shown in Fig. 4(a), initially rises slowly with B_sw to a peak of ∼6.5 keV at B_sw = 10³ nT before falling away sharply to a value of 4.84 × 10⁻² keV at B_sw = 10⁶ nT. The analytic expression of equation (16) for CPCP is also plotted for comparison with the results from the swmf simulations. A discrepancy between the simulation results and the analytic model is apparent: although the general profiles are similar, the swmf CPCP values are approximately an order of magnitude smaller than the corresponding analytic values. A point that may be explained by the fact that the analytic model does not account for viscous interactions at the magnetopause boundary (Nichols & Milan 2016).

Figure 4.

(a) CPCP, (b) total FAC, (c) maximum ionospheric FAC density, and (d) auroral radio power as functions of IMF B_sw for an exoplanet with ionospheric Pedersen conductance of Σ_P = 10⁴ mho. (e) CPCP, (f) total current, (g) maximum FAC density, and (h) auroral radio power as functions of Pedersen conductance for an exoplanet exposed to B_sw = 10⁴ nT. Blue lines represent the swmf results, with diamonds denoting the values are which MHD simulations were conducted. Black lines show the analytic results using the model of Nichols & Milan (2016).

Fig. 4(c) shows the maximum FAC density as a function of B_sw. The magnitude of the FAC density in the analytic model is proportional to the CPCP and Pedersen conductance, given by

$$\begin{eqnarray*} j_{\Vert i} \propto \frac{\Sigma _{\mathrm{P}} \Phi _{\mathrm{CPC}}}{R_{\mathrm{P}}}, \end{eqnarray*}$$

(18)

and the full details of this relation can be found in Nichols & Milan (2016). There is reasonable agreement between the swmf and analytic results, with the saturation and turnover occurring at B_sw = 10³ nT for both, although the simulation results peak at a slighter higher value (⁠|${\sim}197 \,\, \mu \mathrm{A \, m^{-2}}$|⁠) than the analytic results (⁠|${\sim}40 \,\, \mu \mathrm{A \, m^{-2}}$|⁠). Integration of FAC density over a hemisphere (using equation 6) demonstrates that the total ionospheric current also slowly increases with B_sw up to a peak of 2.1 × 10¹⁰ A at B_sw = 10³ nT, with a general profile similar to that of CPCP (Fig. 4a). However, the absolute values of I_tot exceed those from the analytic model of Nichols & Milan (2016). Auroral radio power, calculated using equations (7) and (12), is plotted in Fig. 4(d), and shows a peak value of 4.5 × 10¹⁴ W at B_sw = 10³ nT, a value approximately four orders of magnitude greater than equivalent peak ECMI emission from Jupiter’s aurora.

Figs 4(e)–(h) show the same parameters plotted as a function of Pedersen conductance for a fixed IMF strength of B_sw = 10⁴ nT. The results show that the total maximum FAC density, total current, and radio power are all virtually independent of Pedersen conductance for the swmf simulations, since the low Alfvén conductance implied by the high IMF strength means that the condition for saturation described in Section 4.2, namely when Σ_P ≫ Σ_A, is now satisfied for low Σ_P values.

4.4 Variable orbital distance

Planets at different orbital distances from the host star are subject to not only varying IMF strengths, but also varying stellar wind velocity, density, and Pedersen conductance. In this section, the effects on the CPCP, FACs, and radio power are investigated as functions of orbital distance. For a Sun-like star, Nichols & Milan (2016) calculated analytically how stellar wind parameters vary with orbital distance. Using that work as a reference, spot values of stellar wind velocity, density, IMF strength, and Pedersen conductance were taken at five orbital distances from 2 to 126 R_*, and used as inputs for five swmf simulations. The input parameter values used for each run in the simulation set are summarized in Table 1, and Fig. 5 shows the Pedersen conductance and IMF values used, in relation to the cuts in the parameter space representative of the simulation sets in Section 4.3.

Figure 5.

Values of the Σ_P–B_sw parameter space used in the simulation runs. The dashed line indicates simulation sets at constant Σ_P = 10⁴ mho, i.e. those runs shown in Figs 4(a)–(d), and the solid black line indicates simulations at constant B_sw = 10⁴ nT (Figs 4e–h). Values representative of varying orbital distances are marked by the blue diamonds.

Table 1.

Input parameter values for orbital distance model runs at a range of orbital distances.

d (R*)	|$v$|sw (km s−1)	Bsw (nT)	nsw (cm−3)	ΣP (mho)
2	304.8	3.31 × 104	1.38 × 106	2.46 × 105
5.63	225.2	3.62 × 103	2.77 × 104	3.16 × 104
15.87	279.8	158	1.67 × 103	3.53 × 103
44.72	370.1	4.86	148	405
126.0	453.7	10.4	15	47

Fig. 6(a) shows that the CPCP rises from ∼3 kV in the nominal hot Jupiter orbital distance region, to ∼106 kV at an orbital distance of 126 R_*. Total FAC (Fig. 6b) falls with increasing orbital distance from a value of ∼5 × 10¹⁰ A at d = 2 R_* to ∼2 × 10⁹ A at d = 126 R_*, and a similar trend is observed in the results for maximum FAC density (Fig. 6c). Figs 6(a)–(c) show a good agreement between the swmf results and the analytic results of Nichols & Milan (2016). In Fig. 6(d) auroral radio power is shown along with analytic results of Nichols & Milan (2016), and the RBL. The dashed line shows a least-squares polynomial fit to the swmf results, which yields the relation P_r ∝ d^−1.398. Nichols & Milan (2016), using a Parker spiral IMF, found that the power varies as P_r ∝ d^−5/2 in the inner orbital distance range, i.e. before the notch in Fig. 6(d), and as P_r ∝ d^−5/4 in the outer orbital distance range, where the resultant IMF is dominated by the perpendicular field component. The relation found in this study of P_r ∝ d^−1.398 therefore lies between the two power laws determined by Nichols & Milan (2016). The radio powers in the hot Jupiter orbital range (3–10 R_*) of ∼10¹⁵ W are commensurate with the peak radio powers seen in the results of Fig. 4. Finally, Fig. 6(e) plots the spectral flux density F_r calculated using equation (13) at a distance of s = 15 pc, a value chosen since it is apparent that emission from planets significantly beyond this distance would be below the detection threshold of currently available radio telescopes. Horizontal lines in Fig. 6(e) indicate the sensitivities of the Murchison Widefield Array, the Low-Frequency Array, the Very Large Array (VLA), and the Square Kilometre Array (SKA). The results show that radio flux from hot Jupiters located within 15 pc of the Solar system should be detectable with both the VLA, and in the future with the SKA.

Figure 6.

(a) CPCP, (b) total FAC, (c) maximum ionospheric FAC density, (d) auroral radio power, and (e) spectral flux density as functions of orbital distance d from the central star. Black triangles denote the simulation results, solid blue and red curves indicate comparative analytic functions, and the black dashed line in (d) represents a least-squares fit to the swmf results. Sensitivity thresholds of a number of radio telescopes are marked by horizontal lines in (e).

5 DISCUSSION AND CONCLUSIONS

The model used in the study contains some inherent limitations due to simplifications made to the realistic dynamics of star-hot Jupiter interactions. The upper atmospheres of hot Jupiters undergo intensive escape due to ionization and radiation heating driven by stellar X-ray and EUV radiation, giving rise to the so-called planetary wind (e.g. Yelle 2004; Muñoz 2007; Murray-Clay, Chiang & Murray 2009; Koskinen et al. 2013). This planetary wind is additionally shaped by the gravitational interaction between the star and planet. As can be seen from equations (1)–(5), gravitational effects and ionization due to heating and radiation are not included in the MHD model used in this study. Atmospheric escape in the form of a planetary wind would add an additional pressure outward from the planet which would expand the magnetosphere, and we expect this to increase the strength of the radio emission due to an increase in the size of the stellar wind reconnection channel. This study is a first attempt to establish the M–I coupling dynamics using a 3D MHD model in the framework of existing analytic studies, and as such is primarily intended to investigate the broadbrush effect of IMF strength and Pedersen conductance on the magnetospheric FACs and auroral radio emission. Attempting to incorporate the phenomena of gravitational and ionization effects immediately may introduce additional unconstrained free parameters which could obscure the intention of this study, namely to isolate the response of the model outputs to IMF strength and ionospheric Pedersen conductance. A more complete future study should develop this model to incorporate the gravitational and ionization effects on the star-planet interaction and resulting FACs and radio emission.

One of the most notable features of the results shown in Fig. 4 is the approximately order of magnitude discrepancy between the numerical swmf results and analytic model for the CPCP. A possible explanation for this is suggested by Ridley et al. (2004), who remark that the magnetospheres studied in MHD simulations are, by nature, MHD magnetospheres. Despite accurately depicting the general shape of the pressure distribution, the magnitude is typically underestimated by approximately an order of magnitude at the inner magnetosphere. This underestimation is a result of the lack of energy discretization in the MHD simulation, meaning that the code is unable to model high-energy particles. As the FAC densities in the simulations presented here are substantially higher than those typically encountered in similar swmf modelling of Earth magnetosphere, the associated electrons also have correspondingly higher energies. Hence, an inability to model high-energy electron may have more impact to this study than to studies involving less energetic particles at Earth, and thus provide a plausible explanation for the order of magnitude discrepancy seen between numerical and analytic results.

The saturation of FACs and radio power at relatively low ionospheric Pedersen conductances implies that variations in this parameter may be largely inconsequential at hot Jupiters, where anticipated conductance values are of the order 10⁴–10⁵ mho. Analytic modelling finds that Pedersen conductance affects the point at which CPCP, and therefore radio power, saturates, but this effect is not apparent in the results from swmf simulations, where radio power peaks at a point which appears to be independent of Pedersen conductance.

Fig. 4 shows that as Pedersen conductance is increased, CPCP falls while FAC density remains largely constant for the swmf results. This reinforces previous findings that the magnetosphere acts as neither a voltage nor current generator (Fedder & Lyon 1987; Ridley et al. 2004). By considering Ohm’s law, J = σE, if the magnetosphere is acting as a current generator, where J is constant, then increasing σ would result in a linear increase of E. However, Fig. 4 shows a constant current but non-linear decreasing CPCP. On the other hand, if the magnetosphere acts as a voltage generator then increasing σ would correspond to a linear increase of J for a constant E. Neither of these trends are seen in the swmf results, and in the analytic results both current and potential are found to change simultaneously with Σ_P. The results shown in Fig. 4 represent spot values at horizontal and vertical cuts along a B_sw–Σ_P plane, and exhibit differences between the analytic and numerical MHD values. However, the values plotted in Fig. 6 for a realistic variation with orbital distance actually agree remarkably well.

Various scaling laws exist which approximate the planetary magnetic dynamo performance to provide estimations of the magnetic field strengths at hot Jupiters. For example Griessmeier et al. (2005, 2007) predict magnetic moments of hot Jupiters to be approximately 10 per cent of the Jovian field strength. However, a recent study by Cauley et al. (2019) find evidence for hot Jupiter surface magnetic field strengths approximately 10–100 times larger than those predicted by scaling laws. Therefore, while magnetic field strengths at hot Jupiters may be a fraction of the Jovian value, it is also possible that they exceed the Jovian value. In light of the present high degree of uncertainty in the field regarding hot Jupiter magnetic field strengths determined either from scaling laws or from observations of star-planet interaction, the jovian value is employed in this study as a reasonable benchmark alongside which the Jovian values of other presently unconstrained hot Jupiter parameters (i.e. source plasma population number density and thermal energy) may also be employed for consistency.

Note that we have not considered here hydrodynamic outflow owing to the strong irradiation from the host star or stellar–planetary tidal interaction, and as such our work should be considered to be a first step towards a self-consistent picture of the magnetospheric dynamics of a hot-Jupiter. The escaping planetary wind at hot Jupiters should deform the field lines from their dipolar topology, and under certain conditions may distort the magnetosphere into a magnetodisc configuration (Khodachenko et al. 2015). At Jupiter, the modification of the planetary magnetosphere into a disc-like topology affects the FAC system and the location and size of the main auroral oval (Nichols 2011; Nichols, Achilleos & Cowley 2015). Similar effects may occur at hot Jupiters, although the presence of a magnetodisc may not be typical for each hot Jupiter, and the effect depends on a variety of parameters that are beyond the scope of this paper and would require a dedicated future study.

The use of a constant Pedersen conductance in this work serves as a first approximation to a realistic ionosphere, but future work should examine more plausible conductance models. Exoplanets in close orbit around the host star would likely be tidally locked, meaning that one side of the planet permanently faces the star, and is subject to intense ionizing stellar radiation. This could result in an asymmetric Pedersen conductance pattern, with an ionosphere that is highly conductive on the sub-stellar wind side of the planet and has a low conductance on the opposite side. Particle precipitation from the auroral currents may also further ionize the atmosphere, amplifying the conductance. Such a self-consistent ionospheric model may modify the result presented in this paper, but the general findings would not be expected to alter significantly. This work may also be extended to incorporate planetary rotation. Since the simulations presented here were run in a time-independent mode, planetary rotation was not a factor, but with simple modification simulations could be run in a time-accurate mode to investigate the effects of rotation on the FACs and radio power. In time-accurate mode it would also be possible to examine different stellar wind configurations from the entirely southward B_|$z$| IMF in this work. Hence a more realistic Parker spiral type IMF, or a north–south switching IMF could both be implemented. A future study to explore planetary magnetic field strengths greater the jovian value employed here is also warranted. The trade-off to consider in doing so is the increase in run-time for the swmf.

Ultimately, this work is intended to guide observations regarding the feasibility of detecting auroral radio emission from exoplanets. The maximum predicted radio powers in this study of ∼10¹⁴–10¹⁵ W are consistent with the findings of Nichols & Milan (2016), but, in the hot Jupiter regime, are lower than predicted by studies employing the RBL (e.g. Zarka et al. 2001), potentially explaining the lack of detection to date. A particular finding in this work is that such intense radio emission may only occur for a planet exposed to a narrow range of stellar magnetic field strengths. However, this assumes other stellar wind parameters are unchanging, and the simulations for a range of orbital distances (Fig. 6), taking into account the variation of these other stellar wind parameters, show a relation between radio power and orbital distance similar to that found in the analytic work of Nichols & Milan (2016). The flux density calculations in this study suggest that auroral emission directly from hot Jupiters in the local stellar neighbourhood (≤15 pc) may be presently detectable with telescopes such as the VLA, and in the near future with the SKA. Separating an exoplanetary radio signal from background noise may prove challenging, although modulation of the signal at planetary orbital periods could allow light curves to be folded at that period to aid identification of radio signals.

ACKNOWLEDGEMENTS

ST was supported by an STFC Quota Studentship.

APPENDIX A: MHD MODELLING OF EARTH’S IONOSPHERE

Studies by Ridley et al. (2004) and Ridley (2007) investigated the influence of variable IMF strength and Pedersen conductance on ionospheric CPCP and FACs for Earth. Since this paper investigates similar effects, but at much larger scales and magnitudes, in order to test the method used and to validate the results of this work, the key results from Ridley et al. (2004) and Ridley (2007) were first reproduced and compared with results from the original studies.

Simulations using the swmf were initialized with the terrestrial radius, dipole field strength, and plasma density, as per Ridley et al. (2004) and Ridley (2007). Adaptive mesh refinement of the computational grid was performed in a manner similar to that described by Ridley (2007). Fig. A1 shows the results for the ionospheric FAC density and CPCP as functions of Pedersen conductance, which compare closely with the results of Ridley et al. (2004) (Fig. A2). Fig. A3 plots CPCP as a function of IMF strength, and again the results are in reasonable agreement with those of Ridley (2007) (Fig. A4), thus validating the technique employed in this study.

Figure A1.

Results from this work of (a) maximum ionospheric FAC density and (b) CPCP as functions of Pedersen conductance. (c) CPCP versus the maximum ionospheric FAC density. The diamonds indicate the values at which MHD simulations were run.

Figure A2.

Results from Ridley et al. (2004) equivalent to those from A1.

Figure A3.

Ionospheric CPCP as a function of IMF strength. Diamonds indicate values at which MHD simulations were run.

Figure A4.

Results from Ridley (2007) equivalent to those of A3.