© The Authors 2024

1 Introduction

Hot Jupiters are massive gaseous planets (0.5–10 Jupiter masses) orbiting very close to their host stars (semimajor axis ≤0.1 AU). Due to their mass and proximity to the central star, hot Jupiters are the perfect natural laboratories to study star-planet interactions (SPIs).

Typically, hot Jupiters’ atmospheres are heated by the radiation arising from the central star (Burrows et al. 2000; Owen & Jackson 2012; Buzasi 2013; Lanza 2013). In some cases, the heated planetary gas has enough thermal energy to escape from the gravitational field of the planet (Lammer et al. 2003). This phenomenon is called photoevaporation. Evidence of these outflows has been reported for example for the systems HD 189733b and HD 209458b (Vidal-Madjar et al. 2004; Ehrenreich et al. 2008; Lecavelier Des Etangs et al. 2010, 2012; Bourrier et al. 2013). The rate of photoevaporation depends on the amount of flux emitted by the star and the density of the planet.

In general, the star generates a wind whose physical characteristics (density and velocity) depend on the star’s age (Matt & Pudritz 2008; Cohen & Drake 2014). The planetary outflow expands into and interacts with the stellar wind. This interaction could produce observable signatures that may result in enhanced stellar activity. Shkolnik et al. (2003, 2005, 2008) and Walker et al. (2008) analyzed some systems with hot Jupiters and reported evidence of chromospheric activity in phase with the planetary period. Furthermore, X-ray observations also suggest a correlation between the planetary phase and the stellar activity in HD 189733A (Pillitteri et al. 2015, 2010, 2011, 2014). HD 189733 is a wide binary system composed of a K1-type star (HD 189733A) and an M star (HD 189733B) (Bakos et al. 2006). HD 189733A hosts a transiting hot Jupiter (HD 189733b) in a 2.2-day orbit (Bouchy et al. 2005). This system represents one of the best targets to study SPI because of its proximity, strong activity, and the presence of a transiting planet, which allow the observation of planetary photoevaporation. HD 189733A has an X-ray luminosity on the order of 10²⁸ erg s⁻¹ and a corona brighter than the solar one, with a mean temperature of around 5 MK. The age of HD 189733 A, derived from its level of activity and X-ray emission, appears to be about 1 Gyr, making it notably younger than the M companion HD 189733B. This discrepancy has been explained by the exchange of angular momentum between the star and the planet (Pillitteri et al. 2010; Poppenhaeger & Wolk 2014).

Pillitteri et al. (2010, 2011, 2014) observed X-ray flares at several epochs in the planetary phase range ϕ ≈ 0.54–0.52, and at 0.64. Moreover, the flares appear to be colder than the stellar corona (kT < 1.0 keV or T < 12 MK), suggesting that they have a different origin than stellar activity (Pillitteri et al. 2015). It was inferred that the planetary magnetosphere may affect the coronal structure of the star (Pillitteri et al. 2010). In particular, Pillitteri et al. (2010) suggested that the planetary magnetosphere may perturb the stellar magnetic field twisting the magnetic field lines and, as a result, producing energetic flares via magnetic reconnection (Pillitteri et al. 2010). This hypothesis is also well supported by theoretical models (Shkolnik et al. 2008; Lanza 2009). Pillitteri et al. (2014) inferred the flare properties from the oscillation timescales and concluded that the flaring loop generated after the transit is long ≈25% of the distance between the star and the planet. Moreover, the density of the emitting region is higher than that of the solar corona as observed from OVII and NeIX spectral lines (Pillitteri et al. 2011, 2014). A parallel interpretation of the origin of the X-ray emission is proposed by Pillitteri et al. (2015). They collected data at far-UV wavelengths with the Hubble Space Telescope Cosmic Origin Spectrograph and observed two flaring episodes, at ϕ ≈ 0.525 and ϕ ≈ 0.588, that appear to coincide with the X-ray activity observed previously. They suggest that part of the evaporating material from the planet is intercepted by the stellar gravitational field and accretes onto the star, forming an accretion column. The impact between the supersonic accretion column with the stellar surface generates a hotspot that is responsible for the enhanced UV and X-ray emission.

This idea stems from the work of Matsakos & Königl (2015), who developed an SPI magnetohydrodynamic (MHD) model. They did not focus specifically on the HD 189733 system, but developed an idealized star–planet system that they used to describe the different kinds of SPIs. Depending on the physical conditions of the system, various types of interactions may come into play. The planetary outflow might be forced by a strong stellar wind in a cometary tail that trails the planet along the orbit (Matsakos & Königl 2015) and produces observable signatures in phase with the planetary period (Mura et al. 2011; Kulow et al. 2014). Alternatively, the planetary wind might produce a bow shock in the stellar wind during the expansion (Vidotto et al. 2010, 2011a,b; Llama et al. 2011, 2013). In some cases, it might then be captured by the gravitational field of the central star. When this occurs, the stellar wind is forced to spiral down into the star. The stellar magnetic field funnels the falling gas into an accretion column-like structure. The impact of this accretion flow with the stellar surface might produce a hotspot that in turn generates an observable excess of UV and X-ray radiation (Lanza 2013; Matsakos & Königl 2015; Pillitteri et al. 2015). The hotspot orbits in sync with the planetary period but experiences a phase shift of ≈90deg.

Despite the observational evidence supported by theoretical works, there are other explanations for the X-ray activity in HD 189733. Route & Looney (2019) claim that there is no statistical evidence for a bright hotspot synchronized to the planetary orbital period. Using Kolmogorov-Smirnov and Lomb-Scargle periodogram analyses, they found no evidence for persistent hotspots that have locations synchronized to the planetary orbital period. They suggest that the bright regions that persist for a few rotational periods are entirely consistent with the normal evolution of active regions on stars. In addition, they claim that the accretion rate onto HD 189733a is at least two orders of magnitude lower than the typical values observed in classical T Tauri stars (CTTSs). According to Route & Looney (2019), this low accretion rate would produce an undetectable emission.

In this paper, we analyze the dynamic of HD 189733 and the observability of the UV and X-ray band of the SPI. Our questions are: (a) whether the star is accreting material from the planetary winds, and (b) whether this kind of phenomenon produces visible signatures such as those observed (Pillitteri et al. 2010, 2011, 2014, 2015). To address these questions, we developed an MHD model that describes the star–planet system in HD 189733A. We then synthesized high energy emissions from the results of the simulation.

The paper is structured as follows: in Sect. 2, we describe the MHD model used in this work, in Sect. 3 we discuss the results from the simulation, including the results from the synthesis of emission and the implication for the observability of this kind of SPI. Finally, in Sect. 4 we draw our conclusions.

2 Model

Star-planet systems with evaporating hot Jupiters can be fully described by hydrodynamic or MHD models, as discussed in (Matsakos & Königl 2015, see also the references therein). In this work, we adopted a MHD model analogous to the one described in Matsakos & Königl (2015), but, as previously discussed, with a setup mimicking the system HD 189733A.

HD 189733 is a wide binary system located at a distance of ≈19.3 pc from Earth. The two stars of the system are separated by an average distance of D_s ~ 220 AU. HD 189733A is a K1.5V-type star with mass M_s = 0.805 M_⊙ radius R_s = 0.76 R_⊙, and a rotational period of ≈ 11.95 days (Henry & Winn 2008). This star hosts a hot Jupiter with mass M_p = 1.13 M_J that orbits at a distance of D_p ~ 0.031 AU with a period of P_p = 2.219 days (Bouchy et al. 2005). Winds originate from both the star and the planet and propagate through the interplanetary medium.

In the following, we focus on the planetary system HD 189733A. Due to the fact that D_p ≪ D_s we assumed that the effects of the companion star are negligible in the dynamics of the star-planet system.

2.1 Equations

We adopted a spherical coordinate system (R,θ,ϕ) centered at the center of the star. Since M_s ≫ M_p we considered the center of the star as the center of mass of the entire system. The planet orbits the star in the θ = π/2 plane. The system is described in the reference frame corotating with the planet. In this reference frame, the model solves the MHD equations of conservation of mass, momentum, and energy, as well as the equation of induction for the magnetic field and the equation of entropy: $$\frac{\partial }{\partial t}\rho +\nabla \cdot \rho \mathit{u}=0$$(1) $$\frac{\partial }{\partial t}\rho \mathit{u}+\nabla \cdot \left(\rho \mathit{uu}-\mathit{BB}+\mathit{I}{p}_t\right)=\rho \mathit{g}$$(2) $$\frac{\partial }{\partial t}\rho E+\nabla \cdot [\left(\rho E+p_t\right)\mathit{u}-\mathit{B}(\mathit{u}\cdot \mathit{B})]=\rho \mathit{u}\cdot \left(\mathit{g}+{\mathit{F}}_{\text{ext}}\right)$$(3) $$\frac{\partial }{\partial t}\mathit{B}+\nabla \cdot (\mathit{uB}-\mathit{Bu})=0$$(4) $$\frac{\partial \rho \sigma }{\partial t}+\nabla \cdot (\rho \sigma \mathit{u})=0$$(5)

where $$p_t=P+\frac{\mathit{B}\cdot \mathit{B}}{2},E=ϵ+\frac{\mathit{u}\cdot \mathit{u}}{2}+\frac{\mathit{B}\cdot \mathit{B}}{2\rho }$$(6)

Here, ρ is the plasma density, u the plasma velocity, p_t the total plasma pressure, B the magnetic field, g the gravity acceleration vector, and F_ext an inertial force that appears in our non-inertial rotating frame. F_ext = F_coriolis + F_centrifugal has a Coriolis and centrifugal components given by F_Coriolis = 2(Ω_fr × υ) and F_centrifugal = − [Ω_fr × (Ω_fr × R)]. E is the total plasma energy density, e the thermal energy density, and σ = p/ρ^γ the plasma entropy. We used the ideal gas law P = (γ − 1)/ρ∈, where the polytropic index γ = 3/2.

The calculation was performed using the PLUTO code, a modular Godunov-type code for astrophysical plasmas (Mignone et al. 2012). The code uses parallel computers using the Message Passage Interface (MPI) libraries. The MHD equations are solved using the MHD module available in PLUTO with the Harten-Lax-van Leer Riemann solver. The time evolution is solved using the Hancock method. To follow the magnetic field evolution and to maintain the solenoidal condition, we used the eight waves technique (Powell 1994; Powell et al. 1999).

2.2 Initial and boundary conditions

To impose the initial conditions we first set up the stellar wind which is generated at the surface of the star and expands in a low density static ambient medium (ρ = 10⁻¹⁸ g cm⁻³). The wind is generated at the stellar surface, which is described in the region where R ≤ 1.1 R_⊙. In particular, for R < 0.9 R_⊙ ρ, p and υ are fixed to generate the stellar wind and B is fixed to describe the dipole configuration for the stellar magnetic field. In regions where 0.9_⊙ < R < 1.1_⊙, ρ, p and υ are still fixed to generate the wind, but in order to increase the code stability, we allowed the magnetic field in this region to evolve freely. We let this setup evolve until it reached a stationary condition (t ≈ 1P_p). Then we used this evolved setup and set the planet, its wind, and the unperturbed magnetic field as the initial conditions for our simulation. At t = 0, the planetary wind is generated from the planetary surface. The planetary wind is prescribed in regions inside 1.5R_p where ρ, p, υ and B are fixed. The stellar and planetary interiors are not involved in the simulation and they are prescribed as internal boundaries. As a result, the physical quantities are fixed and the solution is overwritten at each time step. We approximated the star and the planet as rotating solid bodies whose axes rotate parallel to the z-axis. The parameters that characterize the star and planetary winds are summarized in Table 1.

The primary objective of this study is to investigate the intricate interplay between the planetary wind and the stellar wind and its potential to generate accretion events onto the star, which consequently produce discernible X-ray signatures. To achieve this, we opted for an exceptionally high evaporation rate ($\dot{M}\approx {10}^{-9}{M}_{{\text{J\hspace{0.17em}yr}}^{-1}}$) for the planet and a relatively low density stellar wind.

The initial configuration of the global magnetic field results from the combination of the stellar and planetary magnetic fields. The stellar magnetic field is modeled as a Parker spiral, while the planetary magnetic field is represented by a dipole configuration. We explored the influence of the magnetic field strength on shaping the dynamics of the stellar wind. This choice allowed us to assess the impact of the magnetic field strength on the SPI. By incorporating these aspects into our study, we aim to gain deeper insights into the interaction between the planetary wind and the stellar wind, the occurrence of accretion events onto the star, and the resultant observable X-ray signatures. Table 2 shows all the different cases explored in this work.

The inner boundary in R is defined to describe the stellar surface and to generate the stellar wind. The external boundary in R and the boundaries in θ are prescribed as gradient equals 0 (outflow condition) for density, pressure and velocity. The magnetic field at the outer boundary is calculated via parabolic extrapolation using the three closest grid cells. The boundaries for ϕ are located at 0 and 2π; for this reason, they are assumed to be periodic.

2.3 Spatial grid

Our simulations need a grid that adequately resolves the small-scale planetary wind structures in a huge domain. For this purpose, a uniform grid is too expensive from a numerical point of view. Our strategy follows from Matsakos & Königl (2015). The system is described in a proper corotating frame ofreference in which the star and planet are fixed. In this reference frame the planet is located at (r, θ, ϕ) = (D_p, π/2, π). We also made use of a nonuniform spherical grid centered at the center of the star in order to increase the spatial resolution close to the planet and to the stellar surface. The grid is shown in Fig. 1

The radial coordinate, r, extends from the stellar radius (i.e., r_min = 0.805 R_⊙) to r_max = 15 R_⊙ and it is divided into five parts: the first ranges from r = 0.805 R_⊙ to r = 2 R_⊙ with 36 points, and the resolution of mesh decreases with r, meaning there is higher resolution close to the stellar surface. The second grid ranges from r = 2 R_⊙ to r = 5 R_⊙ with 60 points with a uniform resolution of 0.05 R_⊙. The third part is composed of 14 points and ranges from r = 5 R_⊙ to r = 5.5 R_⊙. The fourth part of the grid that ranges from r = 5.5 R_⊙ to r = 9.0 R_⊙ is uniform with 150 points, which corresponds to a resolution of 0.03 R_⊙. The fifth part of the grid ranges up to r = 15 R_⊙ with 100 points and the resolution decreases with r (see Fig. 1).

The θ coordinate ranges between θ_min = 35° and θ_max = 145°. It is composed of three parts: a central uniform and more resolved grid with 100 points ranging from θ_min = 80^° to θ_max = 100° and two grids with 40 points from θ_min = 35° to θ_max = 80° and θ_min = 100° to θ_max = 145°, respectively (see Fig. 1, bottom panel).

An analogous scheme is adopted for the ϕ coordinate. In particular, the ϕ coordinate is divided into three grids: a central uniform and more resolved grid with 100 points from ϕ_min = 170° to ϕ_min = 190° and two with 150 points each ranging from ϕ_min = 0° to ϕ_min = 170° and ϕ_min = 190° to ϕ_min = 360°. This grid ensures the highest resolution close to the planet and in the planetary plane (see Fig. 1, top panel).

Fig. 1 Structure of the numerical grid used in the model. Top panel: slice of one quadrant of the xz plane. To enhance clarity and due to symmetry, only half of the domain is shown. Bottom panel: slice in the xy (equatorial) plane. The dark red dot represents the planet. To enhance clarity, the black grid is four times coarser than the adopted resolution (i.e., each box contains 4 × 4 grid points).

2.4 Synthesis of emission

The code outputs the density, pressure, and velocity and magnetic field components for each cell in the grid, as well as two passive tracers used to identify the planetary wind and the stellar wind. From the model results we synthesized the emission in the X-ray band [0.1–10 keV]. From the 3D domain, we calculated the emission measure from the j_th cell, with density ρ_j and volume dV_j as ${\text{EM}}_j={\rho }_j^2 \text{d}{V}_j$. For each cell we synthesized the emerging spectrum by multiplying the EM by the spectrum per unit of EM obtained from the CHIANTI atomic lines database (Landi et al. 2013), assuming solar abundances. To obtain the total X-ray emission we integrated the spectra in the band [0.1– 10 keV] obtained for each cell in the whole spatial domain. Here we are assuming that the material in the domain is optically thin.

3 Results

In this section, we present the results of the numerical simulations and the synthesis of the emission for the 4 cases. In the following, we consider the simulation with the lowest value of B_s and B_p as reference case. The other simulations are named as in Table 2. The simulations do have the same duration, in fact each run was stopped as soon as a stationary regime is reached.

3.1 Dynamics

Here, we focus on the dynamics of the planetary wind during the simulations. Each simulation was interrupted once a steady state was reached (the configuration of the system do not change for at least half planetary orbit). Figures 2 and 3 show the evolution of the planetary winds during the orbits for the reference case (Bs5-Bp1), while Fig. 4 shows the evolution of the magnetic field (stellar+planetary). Analogous two movies with the full evolution of the system are present as additional files. In this scenario, the evolution is similar to the type III case in Matsakos & Königl (2015).

The planetary wind evaporates from the planet but is confined by the magnetic field and its expansion is favored by the magnetic field along the planetary orbit (see Figs. 2a and b). During the expansion, Rayleigh-Taylor instabilities develop at the contact surface between the planetary and stellar winds. The planetary wind is pushed out of equilibrium by the Rayleigh-Taylor instabilities and due to the impact with the stellar wind, loses angular momentum, triggering the formation of an accretion column of planetary material falling into the star. The plasma that loses angular momentum is caught by the stellar gravitational field and funneled by the stellar magnetic field giving rise to a structure akin to an accretion column that links the planet to the star. The stellar magnetic field drives the plasma into the accretion column and shapes the topology of the accretion column. The material impacts the stellar surface ahead of the planet at about 45°–60° (see Fig. 2f) which is similar to what is suggested by Pillitteri et al. (2015). Another part of the planetary wind is pushed away by the stellar wind and spirals away from the system forming a cometary tail structure. In this region, the interaction with the plasma perturbs the stellar magnetic field and generates regions in which reconnection phenomena may occur. Furthermore, due to the collision with the stellar wind a front shock at T ≈ 10⁶ K develops and heats up the wind at temperatures between T ≈ 10⁵ K and 10⁶ K. (Fig. 3).

The planetary wind strongly interacts with the magnetic field, the details of this interaction are shown in Fig. 4. The result of the interaction between the planetary and stellar winds is a complex magnetic field configuration. In particular, there are regions between the star and the planet where the magnetic field lines twist and The magnetic field structure is notably complex in the region between the planet and the star. In this area, magnetic field lines twist, and regions with opposite-polarity magnetic fields occur. Reconnection events may take place in these regions (see Figs. 4c, e and f). However, our models are unable to describe magnetic reconnection events, because they lack resistivity effects and the high spatial resolution required to accurately describe these events. Nevertheless, the study of reconnection events is beyond the scope of the present paper and may be the primary focus of future studies.

The dynamics exhibited by models with more intense magnetic fields compared to the reference case reveal simpler and distinctive features. These two cases provide compelling evidence for the fundamental role of the magnetic field in shaping the dynamics of the planetary wind. Figure 5 illustrates the dynamics observed in the cases of Bs5Bp5 and Bs10Bp1.

Unlike the reference case, both Bs5Bp5 and Bs10Bp1 exhibit the absence of an accretion column. Instead, the planetary wind expands, assuming a droplet-like structure that surrounds the planet. In fact, the expansion of the planetary wind is impeded in these cases by the ram pressure of the stellar wind, presents also in all cases, and the magnetic pressure. Additionally, the intense magnetic field suppresses the Rayleigh-Taylor instabilities, which are responsible for the formation of the accretion column in the reference cases. This suppression prevents the material in front of the planet from accreting onto the star and instead pushes it behind, forming a cometary tail that extends from the planet. The material in this tail eventually escapes from the system. A shock also forms in these cases, due to the impact between the planetary and stellar wind. In particular, we observe a bow shock in front of the planet as a result of the interaction between the planetary wind and the stellar wind. This kind of dynamic has been adopted as an explanation for Lyα observation (e.g., Vidotto et al. 2011a,b)

To uniquely identify the role of the magnetic field in determining the dynamics of the planetary wind, we also performed a purely hydrodynamic simulation (HD in Table 2), which is the case for B_s = B_p = 0. The evolution of the HD case is shown in Fig. 6. In this case, the system’s evolution closely resembles the reference case, but with one notable difference: no accretion is observed on the stellar surface. In fact, as the reference case, the planetary wind initially expands; during the expansions due to the interaction with the stellar wind Rayleigh-Taylor instabilities develop. The number of these instabilities is larger compared to the reference case and as a result, the planet is surrounded by a more extended cloud of planetary wind with respect to the reference case (see for comparison Figs. 6e and 2e). In this particular case, the absence of accretion can be attributed to the material being dispersed more efficiently, primarily driven by a higher occurrence of Rayleigh-Taylor instabilities compared to the reference case. These increased instabilities contribute to the effective scattering and distribution of the material, preventing significant accretion onto the stellar surface.

Fig. 2 Pole-on views of the density evolution for the Bs5-Bp1 case. The different panels show the density of the planetary wind in a log color scale (right color bar) at different times as indicated in the top-left corner of each image. The yellow sphere at the center represents the central star. The planet is located to the left of the star and it embedded in the planetary wind. The blue-to-white lines represent the magnetic field lines (left color bar). A movie is available online.

Fig. 3 Same as Fig. 2 but for temperature. A movie is available online.

Fig. 4 Magnetic field configuration for the system. Top panel: pole-on view of the system. The equatorial surface shows the magnetic field intensity. Middle panel: view of the magnetic field lines. Bottom panel: close-up view of the planetary region. The red sphere is the planet, and the yellow sphere is the star. A 3D interactive graphic is available at https://skfb.ly/oOzVG/.

3.2 Observability

In this section, we present the results of post processing of the data cubes of the simulations and the observables that can trace back to SPI effects in X-rays. One of the goals of this work is to verify if the interpretation of Pillitteri et al. (2015) about the X-ray flares observed in HD189733 are the results of SPI, and in particular of the planetary wind hitting the stellar surface forming high-temperature X-ray emitting shocks.

3.2.1 Hotspots

The simulations show that the only case in which we observe accretion and onto the stellar surface, and thus hotspots, is the reference case with B_p = 1 G and B_s = 5 G. Figure 7 shows the regions where we expect to observe the impact of the material on the stellar surface and the subsequent generation of hotspots. As expected, the hotspot mirrors the dynamics of the planetary wind. In fact, once the accretion column is formed, the hotspot remains stable throughout the entire simulation. It is worth noting that this model does not describe the shocks due to the impact of the accreting gas on the stellar surface. In fact, the model does not include a proper description of the stellar atmosphere which is assumed to be an external boundary (see Sect. 2). Nonetheless, we can discuss the expected dynamics and the location of the hotspots and their observability. In particular, the hotspots appear to be close to the stellar equator, as a result of the Parker spiral magnetic field configuration assumed.

The amount of X-ray emission emitted by an accretion shock for an optically thin plasma is proportional to the post-shock plasma density squared. Under the assumption of strong shock (Zel'dovich & Raizer 2012), ρ_ps = 4ρ_accr. For these reasons, in order to study the observability, we estimated the accretion rate of material onto the stellar surface and the maximum density of the accretion column. As Fig. 8 shows, after the initial transient phase, which takes about 1 orbital period to establish an accreting column out of the planetary outflow, the accretion rate increases to 10¹² g s⁻¹ with a maximum value of 8 × 10¹² g s⁻¹. This value is dramatically lower than the accretion rates measured in CTTSs, which are typically on the order of 10¹⁷ g s⁻¹ (Hartmann et al. 2016). The lower accretion rate compared to the CTTSs is expected; in this case the source of accretion is the planet, which evidently has a smaller reservoir of mass compared to the stellar disk. For a complete analysis of the observability of the accretion events, we also synthesized the maximum density of the accretion column, which is shown in Fig. 8. The maximum density is about 10⁷ cm⁻³ with a maximum value at 2.5 × 10⁷ cm⁻³. According to Zel’dovich & Raizer (2012) this corresponds to a post-shock density of 10⁸ cm⁻³, which is one order of magnitude lower than typical coronal densities; this means that the emission produced would be indistinguishable from the typical coronal emission. It is important to stress that we are in the most favorable conditions to produce accretion: we assumed an extremely high evaporation rate from the planet (see Sect. 2). For this reason, we can conclude that the accretion shock due to SPI can be excluded as the source of the X-ray flare in phase with the planet observed by Pillitteri et al. (2015).

3.2.2 X-ray emission from diffused wind

From the results of the models, we synthesized the X-ray emission arising from the planetary wind (see Sect. 2.4). In the previous section we demonstrate that flare-like emission observed by Pillitteri et al. (2015) cannot be originated from the impact of the planetary material onto the stellar surface. Here we discuss if the flare-like emission originated from the hot wind surrounding the planet.

Figure 9 illustrates the X-ray light curves obtained from different simulations. In each case, there is a noticeable increase in emission during the transient phase, which spans approximately half of the planetary period, bringing the system into a stable dynamical regime. Subsequently, the emission remains stable, but the HD and Bs5_Bp1 cases exhibit more significant variability compared to the other two cases. This variance is attributed to the weaker magnetic field, which does not effectively suppresses instabilities. In fact, the cases Bs5_Bp5 and BslO_Bpl that show stronger magnetic fields present no further variability once they reach the stationary regime.

In all the cases analyzed, the total emission is on average four orders of magnitude lower than the average X-ray luminosity for HD 189733b (≈10²⁸ erg s⁻¹; e.g. Pillitteri et al. 2022). However, it is worth noting that the case with the strongest magnetic field is the one that shows the highest value of X-ray emission. This is due to the fact that the magnetic confinement of the plasma works at its highest efficiency in the models we explored, and produces the highest density region enshrouding the planet.

Even in this case, we can conclude that this is not the origin of the emission observed by Pillitteri et al. (2015).

Fig. 5 Density evolution for the Bs5-Bp5 case (top three panels) and Bs10-Bp1 (bottom three panels). The images show the density of the planetary wind in a log color scale. The yellow sphere at the center represents the central star.

Fig. 6 Density evolution for the HD case. The images show the density of the planetary wind in a log color scale. The yellow sphere at the center represents the central star.

Fig. 7 Map of the density at the stellar surface in log scale. The spot is located about 45º–6Oº ahead of the planet.

Fig. 8 Maximum value of density on the accretion column (in blue) and accretion rate (in red) vs. planetary period.

Fig. 9 X-ray luminosity vs. planetary period for different cases. In purple is the Bs5_Bpl case, in green the Bs5_Bp5 case, in yellow the BslO_Bpl case and in black the HD case.

4 Conclusions

In this work, we have analyzed whether the X-ray emission that appears in phase with the planetary period observed in HD 189733b by Pillitteri et al. (2015) is a result of an SPI occurring on HD 189733b. In particular, we aimed to verify whether the planetary wind is responsible for the flare-like emission observed. For this reason, we developed a 3D MHD model that describes the system HD 189733A and its hot Jupiter with an extremely strong planetary wind; we explored different magnetic field intensities.

The results of this work are the following:

• The planetary wind expands and interacts with the stellar wind and stellar magnetic field. During the expansion the impact with the stellar wind generates Rayleigh-Taylor instabilities that, if strong enough, push the planetary wind out of equilibrium. The planetary wind is then caught by the stellar gravity and accretes onto the stellar surface, forming an accretion-column-like structure that links the star and the planet;

• Accretion of planetary wind onto a star is not common. In only one of the four cases we explored does the material accrete onto the stellar surface. In fact, in order for it to fall into the star, we need the right combination of magnetic field intensities. A too low magnetic field intensity does not efficiently suppress the Rayleigh-Taylor instabilities, leading to a more diffuse cloud of planetary wind. On the contrary, a too high magnetic field intensity completely suppresses the Rayleigh-Taylor instabilities needed to trigger the accretion column;

• Even in the case that shows accretion, the impact region cannot be responsible for the X-ray variability observed by Pillitteri et al. (2015). The density of the material that hits the stellar surface is on the order of 10⁷ cm⁻³; this is too low to generate an emission distinguishable from the typical coronal emission, which is on the order of 10⁸ cm⁻³ for the solar case;

• The interaction of the planetary wind with the stellar wind generates bow shocks that produce a high-temperature region between the star and the planet. However, the X-ray emission generated from this region is on the order of 10²⁴ erg s⁻¹, which is negligible compared to the average luminosity observed in HD 189733A, about 10²⁸ erg s⁻¹ (e.g., Pillitteri et al. 2022).

We can conclude that this kind of SPI cannot be responsible for the X-ray emission bursts observed in HD 189733A. Other work claims that there is no statistical evidence for a bright hotspot synchronized to the planetary period and that the observed variability in HD 189733A is compatible with the normal evolution of active regions on the star (Route & Looney 2019).

However, an SPI that involves accretion onto the stellar surface and interaction between planetary and stellar wind is not the only scenario that can produce X-ray flares that are in phase with the planetary period. It is important to note that the magnetic configuration in the tail of the planetary wind is highly perturbed. In principle, the interaction of the planetary wind with the stellar magnetic field can generate clumpy regions, as observed for the solar case (Petralia et al. 2016). The perturbed magnetic field configuration may also generate recombination effects that heat the clumpy regions and produce observable X-ray emission. In order to investigate this effect and its observability, we would need a more spatially resolved model that includes a resistive term that takes the heating of the plasma due to magnetic reconnection into account.

Movies

Movie 1 associated with Fig. 2 (Dens_Bs5Bp1) Access here

Movie 2 associated with Fig. 3 (Temp_Bs5Bp1) Access here

Acknowledgements

We acknowledge the financial contribution from the agreement ASI-INAF No. 2018-16-HH.0 (THE StellaR PAth project). We acknowledge support from ASI-INAF agreement 2021-5-HH.0 Partecipazione alla fase B2/C della missione ARIEL (Atmospheric Remote-Sensing Infrared Exoplanet Large survey). PLUTO was developed at the Turin Astronomical Observatory and the Department of Physics of Turin University. We acknowledge the HPC facility (SCAN) of the INAF - Osservatorio Astronomico di Palermo for the availability of high-performance computing resources and support. We acknowledge PRACE for awarding access to the Fenix Infrastructure resources at CINECA, which are partially funded from the European Union’s Horizon 2020 research and innovation programme through the ICEI project under the grant agreement no. 800858.

Appendix A Resolution convergence test

The spatial resolution adopted reflects a compromise between the need of high resolution and the computational cost of each simulation. To achieve this goal, we performed a convergence test on the HD case, considering two additional simulations with numerical meshes four times (HD4x) and two times (HD2x) coarser than the one adopted in the paper, respectively. Indeed, this case is characterized by RT instability, which results in the formation of small-scale structures that may be influenced by the numerical grid and, consequently, by the spatial resolution. The convergence test was performed through the tracer defined for the planetary wind to identify zones whose content is made up of original planetary wind material by more than 10%. Figure A.1 shows the total mass of these zones in the computational domain during the evolution in the three cases analyzed.

We find that the results obtained with different spatial resolutions exhibit the same general trend, with differences not exceeding 15%. After the initial transient phase of about half orbital period, during which the mass linearly increases with time as the planetary wind expands through the spatial domain, a stationary condition is achieved around a value of M_pw ≈ 3 × 10¹⁷ g. This suggests that an equivalent amount of mass enters the domain from the planet’s surface and exits from the external boundary of the domain. On the other hand, we note a striking difference between case HD4x and the other two cases: after the initial transient, the total mass M_pw remains constant in HD4x, whereas it exhibits chaotic and similar variability in the other two cases. This is attributed to the onset of RT instability, leading to substantial perturbations in the planetary wind. Remarkably, these perturbations appear very similar in the two cases HD2x and HD, indicating that our simulations successfully capture the fundamental properties of the dynamics. In light of this, we decided not to further increase the spatial resolution and opted for the HD one.

Fig. A.1 Mass of planetary outflow (at least 10% of which is of planetary origin) as a function of time. The green line is the HD4x case, the yellow line is the HD2x case, and the blue line is the HD case.

References

All Tables

All Figures

Fig. 1 Structure of the numerical grid used in the model. Top panel: slice of one quadrant of the xz plane. To enhance clarity and due to symmetry, only half of the domain is shown. Bottom panel: slice in the xy (equatorial) plane. The dark red dot represents the planet. To enhance clarity, the black grid is four times coarser than the adopted resolution (i.e., each box contains 4 × 4 grid points).

In the text

Fig. 2 Pole-on views of the density evolution for the Bs5-Bp1 case. The different panels show the density of the planetary wind in a log color scale (right color bar) at different times as indicated in the top-left corner of each image. The yellow sphere at the center represents the central star. The planet is located to the left of the star and it embedded in the planetary wind. The blue-to-white lines represent the magnetic field lines (left color bar). A movie is available online.

In the text

	Fig. 3 Same as Fig. 2 but for temperature. A movie is available online.
In the text

	Fig. 4 Magnetic field configuration for the system. Top panel: pole-on view of the system. The equatorial surface shows the magnetic field intensity. Middle panel: view of the magnetic field lines. Bottom panel: close-up view of the planetary region. The red sphere is the planet, and the yellow sphere is the star. A 3D interactive graphic is available at https://skfb.ly/oOzVG/.
In the text

	Fig. 5 Density evolution for the Bs5-Bp5 case (top three panels) and Bs10-Bp1 (bottom three panels). The images show the density of the planetary wind in a log color scale. The yellow sphere at the center represents the central star.
In the text

	Fig. 6 Density evolution for the HD case. The images show the density of the planetary wind in a log color scale. The yellow sphere at the center represents the central star.
In the text

	Fig. 7 Map of the density at the stellar surface in log scale. The spot is located about 45º–6Oº ahead of the planet.
In the text

	Fig. 8 Maximum value of density on the accretion column (in blue) and accretion rate (in red) vs. planetary period.
In the text

	Fig. 9 X-ray luminosity vs. planetary period for different cases. In purple is the Bs5_Bpl case, in green the Bs5_Bp5 case, in yellow the BslO_Bpl case and in black the HD case.
In the text

	Fig. A.1 Mass of planetary outflow (at least 10% of which is of planetary origin) as a function of time. The green line is the HD4x case, the yellow line is the HD2x case, and the blue line is the HD case.
In the text