© The Authors 2023

1 Introduction

The evolution of planets is shaped by billions of years of interactions with their host stars and with neighbouring planets. In addition to the central gravitational force of the host star, there are several types of interactions taking place between a star and an exoplanet (Vidotto 2020), some of which are not present, or simply very weak, in the present-day Solar System planets.

In the case of the GJ 436 system, the focus of the current study, the warm Neptune GJ 436b experiences a strong interaction with its host star due to the combined effects of high-energy stellar irradiation that is deposited in the planetary atmosphere and a subsequent interaction of this atmosphere with the stellar wind (Bourrier et al. 2016). Spectroscopic transit observations in Lyman-α (Kulow et al. 2014; Ehrenreich et al. 2015; Lavie et al. 2017; Dos Santos et al. 2019) revealed that GJ 436b is enshrouded by a giant cloud of escaping atmosphere that trails GJ 436b along its orbit. Three-dimensional models suggest that this trailing structure takes the form of a giant comet-like tail (Bourrier et al. 2016; Shaikhislamov et al. 2018; Khodachenko et al. 2019; Villarreal D’Angelo et al. 2021). By modelling the observations of spectroscopic transits of GJ 436b with atmospheric escape models that account for atmospheric interaction with the stellar wind, local properties of the stellar wind at the orbital distance of the planet, such as the wind speed and density, can be estimated; they can then be used to infer global properties of stellar winds, such as their mass-loss rates (Vidotto & Bourrier 2017). All these quantities shape the space weather around a planet.

Indeed, one of the major strengths in observing and modelling star-planet interactions in general is that they can provide alternative ways to physically characterise planetary systems. Stellar wind mass-loss rates, which are difficult to constrain in the case of cool dwarf stars (Vidotto 2021), can also be estimated by investigating star-planet interactions through planet-induced radio emission (Kavanagh et al. 2021; Pineda & Villadsen 2023) or auroral radio emission from the exoplanet (Vidotto et al. 2019). Investigating star-planet interactions through chromospheric hot spots associated with anomalous stellar activity (e.g. Shkolnik et al. 2005; Cauley et al. 2019) or exoplanetary radio emission (e.g. Farrell et al. 1999; Turner et al. 2021) allows one to probe planetary magnetism, a quantity that is believed to be very important in the context of planetary habitability but is still rather poorly known in planets outside the Solar System.

In this work we focus on star–planet interactions that are mediated by the stellar wind and its embedded magnetic field. In the literature, this type of interaction is referred to as ‘star–planet (electro)magnetic interactions’, ‘stellar wind-planet interactions’ or ‘electrodynamic star-planet interactions’, but they are only different names of the same process. In the context of star–planet interactions mediated by the magnetised host star‘s wind, different physical processes take place depending on the wind regime in which the planet orbits.

In the regime of sub-Alfvénic motion of a planet through the wind of its host star, a magnetic ‘connectivity’ between the planet and the star can take place. In this regime, the orbiting planet can trigger magnetohydrodynamic (MHD) waves that can travel towards the star, or alternatively (and even possibly concurrently) it can also trigger magnetic reconnection events in the case of a magnetised planet. Both of these mechanisms release energy that is radiated away and could be detected in the form of, for example, chromospheric hot spots (Lanza 2009; Strugarek et al. 2022), planet-induced coronal radio emission (Saur et al. 2013; Pineda & Villadsen 2023), anomalous flare events (Fischer & Saur 2019), and even potentially the lack of energetic flare events, as recently proposed by Loyd et al. (2023) in the context of GJ 436. We note that all these signatures originate on or near the host star, and not on the planet.

The inverse takes place in the regime of super-Alfvénic motion of a planet through the wind of its host star. In this regime, the interaction between the planet and the stellar wind can lead to radio emission originating from the planet itself (more precisely, from its magnetosphere). This is because, in the super-Alfvénic regime, the magnetic connectivity between the star and planet can no longer form. This is the situation currently experienced by the Solar System planets, where we see that the planetary auroral radio emission is correlated with the solar wind energy that is dissipated in the magnetosphere of Solar System planets (e.g. Farrell et al. 1999). Despite many attempts to detect such signatures in exoplanets, it was only recently that a potential detection of radio emission from a hot Jupiter was reported (Turner et al. 2021, see also Elekes & Saur 2023 for models of these observations).

Fundamental to the nature of stellar wind-planet interactions is thus the regime of the orbital motion of the planet: sub- or super-Alfvénic. To determine the regime, we need to obtain more realistic stellar wind properties, such as densities, velocities, and magnetism. These stellar wind quantities are also what shape the space weather of planets and are thus also relevant for, for example, understanding the flux of energetic particles impacting planetary atmospheres (e.g. Mesquita et al. 2021; Rodgers-Lee et al. 2023). In the case of GJ 436, its wind has been investigated in several works (Vidotto & Bourrier 2017; Mesquita & Vidotto 2020; Mesquita et al. 2021), but it is only recently that the surface magnetic field of the host star GJ 436 was mapped (Bellotti et al. 2023); this allows us to now model, for the first time, the 3D structure of the stellar wind, including the observed surface magnetism. Magnetic fields play an important role in the winds of cool dwarf stars, as they are responsible for the mechanisms that heat and accelerate these winds. Section 2 presents our 3D MHD simulations of the wind of GJ 436. We then use the results of our wind models to compute the space weather conditions (Sect. 3) along the highly misaligned orbit of GJ 436b. GJ 436b is a Neptune-mass planet orbiting very close to its slowly rotating M dwarf host. It has a puzzling observed orbital architecture – a non-circular orbit with an eccentricity of 0.152 (Trifonov et al. 2018) and strong orbital misalignment (Bourrier et al. 2018, 2022) – that is likely the result of a delayed migration following Kozai-Lidov interactions with an outer massive companion (Bourrier et al. 2018; Attia et al. 2021). In Sect. 4 we then compute the strengths of star-planet interactions. In the final section of this paper (Sect. 5), we discuss how our computations compare to the recent observations of Loyd et al. (2023) and present our concluding thoughts.

Fig. 1 Observationally reconstructed radial magnetic field of GJ 436 from Bellotti et al. (2023) using Zeeman-Doppler imaging. This map is incorporated at the boundary of our stellar wind simulations.

2 The wind of GJ 436

The environment surrounding GJ 436b consists of particles and magnetic fields that form the wind of the host star. Here, we characterise this environment by performing data-driven 3D MHD simulations of the stellar wind. Because winds of cool dwarf stars are ultimately driven by stellar magnetism, at the boundary of our model (which is located close to the stellar surface, at the transition region), we include the reconstructed surface field shown in Fig. 1 (Bellotti et al. 2023). This map was reconstructed using the Zeeman-Doppler imaging technique (Donati & Brown 1997; Donati & Landstreet 2009) and, thus, represents the large-scale field of GJ 436 in March-June 2016. For this simulation, we used the Alfvén Wave Solar Model (AWSoM; van der Holst et al. 2014) implemented in the numerical code BATS-R-US (Powell et al. 1999; Toth et al. 2012). We refer the reader to van der Holst et al. (2014) for further details of the model. Below we describe the main input parameters adopted in the simulations of the wind of GJ 436.

In addition to the stellar surface magnetic map, other model inputs are the stellar properties (mass, radius and rotation period; cf. Table 1), the density and temperature at the transition region, and the Alfvén wave properties (boundary wave flux and the correlation length of the Alfvénic wave turbulence). Of these parameters, the least constrained ones are the wave properties, as we discuss below.

The inner boundary of our model is at the transition region, and values considered for solar wind models (van der Holst et al. 2014; Oran et al. 2017) for ion/electron density (n_tr = 5 × 10¹¹ cm⁻³), temperature (T_tr = 5 × 10⁴ K), and, consequently, gas density (P_tr ≃ 1.4 dyn cm⁻²), were adopted. Our assumed pressure is not too dissimilar to the value assumed in the chromospheric models of the M dwarf star GJ832 by Fontenla et al. (2016), who used P_tr = 0.8 dyn cm⁻². For a temperature of 5 × 10⁴ K, this pressure value translates to a density of 10¹¹ cm⁻³, which is also similar to density values found at this temperature from the chromosphere+wind models of Sakaue & Shibata (2021). However, there are also models in the literature that adopt lower pressure values of P_tr = 0.02−0.04 dyn cm⁻² (Mauas et al. 1997; Peacock et al. 2019), which would imply transition region densities of a few times 10⁹ cm⁻³ at a temperature of 5 × 10⁴ K. Depending on which chromospheric model we use to guide our assumption of n_tr, we could have transition region densities that vary from ~10⁹ to 10¹¹cm⁻³. O Fionnagain et al. (2021) explored the effects that changing this parameter had on the wind of a sub-giant star λ And and they found that varying n_tr by a factor of 100 led to similar wind properties (compare their models A2 and B0, or B1 and D0, from their Table 2). We therefore do not expect that assuming a different value for n_tr within the quoted range would change our wind solutions substantially.

As they propagate through magnetic field lines, Alfvén waves transfer their energy and momentum to the wind plasma, thus causing its heating and acceleration. In AWSoM, the wave dissipation is modelled through a turbulent cascade process. Its correlation length is assumed to increase with flux tube radius; thus $L_{\perp }=L_{\perp ,0}\sqrt{B_{\star }/B}$. Here, we assume the proportionality factor is the same as that used in the AWSoM solar wind models, namely $L_{\perp ,0}\sqrt{B_{\star }}=1.5\times {10}^5$ $\text{m\hspace{0.17em}}\sqrt{T}=1.5\times {10}^4\text{\hspace{0.17em}km\hspace{0.17em}}\sqrt{\text{G}}$. Based on the parametric study of O Fionnagain et al. (2021), we also do not expect that changing the value of this base quantity by a factor of ~ 10 will substantially affect the general properties of our wind solution (see also Shoda et al. 2018, who investigated the effects of the correlation length on wind properties such as its maximum temperature, mass-loss rate, and velocity).

An important quantity of our model is the boundary value of the Alfvén wave flux: $$S_{A,0}={\rho }_{\text{tr}}\langle {\upsilon }_{\perp ,0}^2\rangle {\upsilon }_{A,0},$$(1)

where ${\langle {\upsilon }_{\perp ,0}^2\rangle }^{1/2}$ is the average velocity perturbation of the Alfvén wave, ρ_tr = m_pn_tr is the mass density at the transition region, and the Alfvén velocity is ${\upsilon }_{A,0}=B_{\star }/\sqrt{4\pi {\rho }_{\text{tr}}}$, where the index ‘0’ indicates value computed at the boundary of the simulation. Here, B_⋆ is the local value (i.e. at each longitude and latitude) of the surface magnetic field as derived in the surface magnetic map (Fig. 1). The previous equation can be rewritten as $$\frac{{\langle {\upsilon }_{\perp ,0}^2\rangle }^{1/2}}{65\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}={\left(\frac{5\times {10}^{11}}{n_{\text{tr}}\left[{\text{cm}}^{-3}\right]}\right)}^{1/4}{\left(\frac{S_{A,0}/B_{\star }\left[\text{erg}/\left({\text{cm}}^2\text{\hspace{0.17em}s\hspace{0.17em}G}\right)\right]}{1.1\times {10}^7}\right)}^{1/2}.$$(2)

We note that the ratio S_A,0/B_⋆ is assumed to be uniform at the inner boundary, implying that the average velocity amplitude of the Alfvén wave and the energy density of the wave (ρ_tr$\langle {\upsilon }_{\perp ,0}^2\rangle$) are also uniform at the inner boundary. The choice of solar wave flux adopted in previous models of the solar wind with AWSoM (Oran et al. 2013, 2017; van der Holst et al. 2014) was guided by Hinode observations of velocity fluctuations at the chromospheric region (De Pontieu et al. 2007). Indeed, if one assumes that the non-thermal motions of ions at the transition region are associated with Alfvén waves (e.g. Oran et al. 2017), then one way to constrain the Alfvén wave properties in other stars is by measuring the non-thermal broadening of lines formed in the transition region (Boro Saikia et al. 2023). In this case, the non-thermal velocity ξ is related to the average velocity perturbation of the Alfvén wave as ${\xi }^2=\frac{1}{2}\langle {\upsilon }_{\perp ,0}^2\rangle$.

The ratio S_A,0/B_⋆ is the most important free parameter in our model, and its value affects the wind temperature, its mass-loss rate, and the size of the Alfvén surface, which in turn affects the angular momentum loss of the star (Boro Saikia et al. 2020; O Fionnagain et al. 2021; Kavanagh et al. 2021). Given that, we ran two wind models with two different values of S_A,0/B_⋆: 1.1 × 10⁶ erg (cm² s G)⁻¹ and 1.1 × 10⁷ erg (cm² s G)⁻¹. The latter is the same value as used in AWSoM solar wind simulations (Sokolov et al. 2013; van der Holst et al. 2014). From now on, we refer to these models as ‘I’ and ‘II’. As we show in Sect. 3, our derived mass-loss rates and computed non-thermal velocities of Model I show better agreement with observations, indicating that the adopted lower Alfvén wave flux is more appropriate to describe the wind of GJ 436.

For GJ 436, which has a surface radial field reaching about 28 G, our choice of parameters translates to correlation lengths of L_⊥,0 ≃ 2800 km¹ and surface Alfvén fluxes of up to S_A,0 ≃ [3 × 10⁷,3 × 10⁸] erg cm⁻²s⁻¹ for Models I and II, respectively. The surface amplitudes of the wave perturbation are 2 km s⁻¹ (Model I) and 65 km s⁻¹ (Model II), uniform at the lower boundary. The corresponding surface energy densities of the wave are 11 and 35 erg cm⁻³. With these parameters, the predicted non-thermal velocities at the transition region are ξ ≃ 14 and 45 km s⁻¹ for Models I and II, respectively. The other remaining parameters of the model, namely the stochastic heating parameter, the collisionless heat conduction parameter and the collisionless heat flux parameter have the same values as in Table 1 of van der Holst et al. (2014). The energy equations are solved in their non-conservative form (i.e. solving for the pressure instead of total energy).

Figure 2 shows the 3D outputs of our simulations, after they have reached steady state. The grey streamlines represent the stellar wind magnetic field lines; this geometry results from the self-consistent interaction between field lines and the stellar wind flow. The rotation axis lies along positive z_⋆ − we note that the magnetic field geometry is approximately that of an aligned dipole close to the stellar surface, with elongated closed field lines at large distances and open fields around the stellar rotation poles. The stellar wind mass-loss rate we obtained are 1.1 × 10⁻¹⁵ M_⊙ yr⁻¹ and 2.5 × 10⁻¹⁴ M₀ yr⁻¹ for Models I and II, respectively. For comparison, the solar wind mass-loss rate is 2 × 10⁻¹⁴ M_⊙ yr⁻¹ − given that all of our input values for the Alfvén wave and transition region are solar for Model II, and the surface magnetic flux of GJ 436 is similar to solar², it is not surprising that the mass-loss rate derived in Model II is similar to that of the solar wind. In Sect. 3 we perform a more detailed comparison of our derived mass-loss rates with those derived in other studies of GJ 436. Figure 2 also shows a cut in the y_⋆z_⋆ plane, where colours indicate the stellar wind speed, where we see that along open field lines the wind is accelerated more quickly, with low-speed streams around closed field line regions. As expected, higher velocities are achieved in Model II, with higher Alfvén wave fluxes. The Alfvén surface is shown by the translucent surface. For context, we plot a white circumference at the orbital distance of GJ 436b (for simplicity, we assume a circular orbit). We note most of the time the planet would be in the region of open magnetic field lines and that the orbital distance is nearly always within the Alfvén surface for both models, which can have implications for star-planet interactions (Folsom et al. 2020) and induced radio emission on the star (Kavanagh et al. 2021, Kavanagh et al. 2022), for example. We come back to this in Sect. 4.

Fig. 2 3D view of the magnetised stellar wind for Models I and II. The surface magnetic field from March-June 2016 reconstructed by Bellotti et al. (2023) is used at the boundary of the models. The grey streamlines are the magnetic field lines embedded in the wind. The rotation axis of the star is along the z⋆-axis. The 2D cut is placed at the x⋆z⋆ plane, and for reference the orbital distance of the planet is shown with the translucent sphere at 14.56 R⋆ ≃ 0.028 au. The contour shows the wind velocity, and the Alfvén surface is shown by the translucent outer surface. Models I and II predict stellar wind mass-loss rates of 1.1 × 10−15 M⊙ yr−1 and 2.5 × 10−14 M⊙ yr−1.

3 Space weather conditions of GJ 436b

The space weather of GJ 436b is formed by the particles and ambient magnetic field of the stellar outflow. Here, we compute the properties of stellar wind at the position of GJ 436b. Because GJ 436b has a very misaligned orbit, we first compute the position of the planet as a function of time, following the coordinate transformations presented in Appendix A. To calculate the phase of the orbit with respect to the rotational phase of the star, we consider the transit midpoint at HJD_t = 2458 947.26212(12) (Maxted et al. 2022), with an orbital period of 2.64389803(27) d (Lanotte et al. 2014). At stellar rotation phase zero (HJD_s = 2457 464.4967, corresponding to the first spectropolarimetric observation (Bellotti et al. 2023), and {x_⋆, y_⋆} = {R_⋆, 0} in our simulation coordinates), the planet was at an orbital phase of φ₀ = 0.1745 after the preceding transit midpoint. GJ 436b has a nearly polar orbit with its orbital axis and stellar spin misaligned by Ψ = 103.2° (Bourrier et al. 2022). Figure 3 shows the trajectories of the planet as seen in the reference frame of the star during one rotational period of the star (about 44 days), during which the planet completes about 16 orbital revolutions.

As presented in Appendix A, in the case of GJ 436b, which has a highly misaligned orbit, the time it takes for the planet to be at the same sub-stellar point is a complex combination of orbital and rotation periods, as well as the obliquity of the orbit. Because of this, the local condition of the stellar wind is constantly changing. Figure 4 shows the stellar wind proton density, velocity, total pressure (ram + thermal + magnetic) and magnetic field strength, all computed at the orbit of GJ436b. The solid lines are for Model I, and the dashed lines for Model II. Colour represents the time evolution following the same colour scheme as shown in Fig. 3. We note that to compute the ram pressure, we use the relative velocity between the stellar wind velocity u and the Keplerian motion of the planet v_K: Δu = и − v_K. For GJ 436b, the Keplerian velocity is v_K = 118($\sin \Psi \widehat{\theta }+\cos \Psi \widehat{\phi }$) km s⁻¹, where we adopted a stellar mass of ≃0.45 M_⊙ and an orbital distance of 14.56 R_⋆ ≃ 0.028 au (Maxted et al. 2022).

The range of wind densities at the orbit of the planet is [0.12,5.0] × 10³ cm⁻³ for Model I and [2.8,81] × 10³ cm⁻³ for Model II. Their time averages over one stellar rotation period are 510 and 6200 cm⁻³ for Models I and II, respectively (horizontal grey lines shown in Fig. 4). Similarly, for wind velocities, we found that the range of values that the planet experience are [30, 770] km s⁻¹ and [220, 830] km s⁻¹ for Models I and II, respectively, with averages of 510 and 640 km s⁻¹. When modelling the transmission spectroscopic transit of GJ 436b, Bourrier et al. (2016) derived densities of $2_{-1.2}^{+2.2}\times {10}^3{\text{cm}}^{-3}$ and wind speeds of ${85}_{-16}^{+6}\text{km\hspace{0.17em}}{\text{s}}^{-1}$. These observations are taken during transit, so the values they derived should be compared to our model values around mid-transit times. For this comparison, we added to the top panels of Fig. 4 crosses at the approximate mid-transit times at the same level of ion density (2 × 10³ cm⁻³) and speed (85 km s⁻¹) derived in the work of Bourrier et al. (2016). We see that these crosses overlap really well with values of Model I (solid lines), but not with Model II (dashed lines), implying that Model I, with the lowest surface flux of Alfvén waves, better matches the spectroscopic transit values from Bourrier et al. (2016). Indeed, from these observations, the estimated mass-loss rate of ${1.2}_{-0.75}^{+1.3}\times {10}^{-15}{M}_{\odot }{\text{yr}}^{-1}$ (Vidotto & Bourrier 2017) better agree with the values derived in Model I (1.1 × 10⁻¹⁵ M_⊙ yr⁻¹), but not with Model II (2.5 × 10⁻¹⁴ M_⊙ yr⁻¹). Our chosen values of S_A,0 /B_⋆ correspond to non-thermal velocities ξ ≃ 14 and 45 km s⁻¹ at the transition region (see Sect. 2). Recently, Boro Saikia et al. (2023) reported ξ ranging from ∼5 to ∼30 km s⁻¹ for stars with similar rotation periods and effective temperatures as GJ 436. In this case too, the non-thermal velocities from Model I have a better agreement with the range reported in Boro Saikia et al. (2023). Altogether, our derived mass-loss rates and non-thermal velocities suggest that the ‘canonical’ value of S_A,0/B_⋆ = 1.1 × 10⁷ erg (cm² s G)⁻¹, derived from solar wind models, might not be applicable to the winds of (old) M dwarf stars.

The bottom-right panel of Fig. 4 shows the stellar wind magnetic field strength at the orbit of GJ 436b. We find an average stellar magnetic field strength of 0.013 G and 0.026 G for Models I and II, respectively. We use these values to then compute the magnetic pressure (B²/8π). We also compute the ram pressure of the stellar wind (ρΔu²). Adding these two pressures to the thermal pressure, we obtain the total pressure of the stellar wind at the position of the planet, which is shown at the bottom-left panel of Fig. 4. We find average total pressures of 8.2 × 10⁻⁶ dyn/cm⁻² and 59 × 10⁻⁶ dyn/cm⁻² for Models I and II, respectively. These values will be used in Sect. 4, when we compute the strength of star-planet interactions and the size of the magnetosphere of the planet.

Fig. 3 Trajectories of the planet as seen in the reference frame of the star during one rotational period of the star (about 44 days; see the colour bar), during which the planet completes about 16 orbital revolutions. The grey circles indicate the approximate times of mid-transit.

4 Stellar wind–planet interactions

Because of its close proximity to the star, and its sub-Alfvénic motion, GJ 436b can magnetically interact with its host star. Recently, Loyd et al. (2023) studied the presence of star-planet interaction signatures in GJ436 in ultraviolet spectroscopic observations. While they did not find line fluxes modulated with the orbital period of GJ 436b (traditionally used as evidence of star-planet interactions; e.g. Pineda & Villadsen 2023), they found an enhancement of low energy flares and a lack of more energetic flares in GJ436, compared to other M dwarfs. They suggested that this could be caused by the presence of GJ 436b, whereby the reduction of the number of high energy flares in GJ 436 would be caused by the planet triggering early flaring events (i.e. before enough energy would have been built up in a larger stellar flare). As a consequence, high-energy flares in GJ 436 would not have time to store energy and would be prematurely released, explaining the excess of low energy flares and lack of high energy flares in the system. In their analysis, they also derived a maximum power for the star-planet interaction, detected in the far ultraviolet lines, of −3 × 10²⁴ erg s⁻¹. Converting this number to a bolometric power, they found an upper limit of ≲10²⁶ erg s⁻¹ for the maximum power released in star-planet interactions in the GJ 436 system.

To compare with the derived power from Loyd et al. (2023), in this section we compute the strength of star-planet interactions. Here, we use two different scenarios for this calculation. One scenario is based on the idea that magnetic reconnection occurring between stellar and planetary magnetic field lines releases magnetic energy (Jardine & Cameron 2008; Lanza 2009; Vidotto et al. 2010). This energy can accelerate electrons, which travel towards the star. This idea has been used to explain anomalous hot spots at stellar chromospheres (Lanza 2009), as reported in observational works (e.g. Shkolnik et al. 2005; Pagano et al. 2009; Cauley et al. 2019). The other scenario we also explore is based on the idea that a planet moving at sub-Alfvénic speed through the magnetised wind of its host star can trigger MHD Alfvén waves, which can travel towards the star along magnetic field lines (Preusse et al. 2006; Saur et al. 2013; Strugarek 2016). As they travel, these waves dissipate their energies, which are radiated away.

We next estimated the power released in the star-planet interaction using the two scenarios mentioned above. To compute the strength of star-planet interactions in both scenarios, we started by calculating the (unsigned) Poynting flux of the stellar wind that impacts on the planet: $$P=\frac{\left|\left(\text{Δ}u\times B\times B\right)\right|}{4\pi }=\frac{B^2\text{Δ}u\text{\hspace{0.17em}sin\hspace{0.17em}}\theta }{4\pi }.$$(3)

In this expression, B sin θ is the component of the stellar wind magnetic field that is perpendicular to ∆u, θ is the angle between ∆u and B − all these quantities are computed at the position of the planet, which lies at an orbital distance a_orb = 14.56 R_⋆. We find that the average Poynting fluxes computed at the orbit of GJ 436b are P ≃ 160 and 650 erg cm⁻² s⁻¹ (or 0.1 to 0.6 W m⁻²) for Models I and II, respectively. In the context of planets moving at super-Alfvénic speeds, the Poynting flux can be used to estimate the amount of power dissipated in the magnetosphere of a planet, which can then power electron-cyclotron radio emission (e.g. Vidotto et al. 2019). In the case of GJ 436b, though, our wind models indicate that the motion is sub-Alfvénic (cf. Sect. 2) and we use the Poynting flux to compute the energy released in the Alfvén wing scenario and the reconnection scenario. We note that the Poynting flux only depends on the characteristics of the wind at a certain orbital distance and does not depend on the properties of the planet, such as its size or magnetic field strength, contrary to the power released in star-planet interactions, as we see below.

Fig. 4 Stellar wind proton density, velocity, total pressure (sum of ram, including the planet‘s own orbital motion and the thermal and magnetic pressures), and total magnetic field strength at the orbit of GJ 436b. The solid lines are for Model I and the dashed lines for Model II. Colour represents the time evolution and follows the colour scheme shown in Fig. 3. The horizontal grey lines are the time average of the stellar wind properties over one stellar rotation period (about 16 orbital periods of the planet). The grey arrows in the upper horizontal axes indicate approximate times of mid-transits. In the upper panels, the crosses indicate the wind values derived in the work of Bourrier et al. (2016), and they are plotted at the mid-transit times. Note how they overlap really well with the values derived in our Model I, indicating that Model I provides a better description of the wind of GJ 436.

4.1 Scenario 1: Alfvén wings

A planet moving at sub-Alfvénic speed through the magnetised wind of its host star can trigger MHD Alfvén waves, which can travel towards the star along magnetic field lines. As they travel, these waves can then dissipate their energies and the maximum power that is radiated away can be estimated by integrating the Poynting flux of the stellar wind over the cross-section of two Alfvén ‘wings’. In the limit of sub-Alfvénic speeds, with Alfvén-Mach number M_A = Δu/υ_A < 1, the dissipated power is (Saur et al. 2013) $${𝒫}_{\text{wg}}\simeq P\left(2\pi r_M^2\right){\alpha }^2{M}_A\text{\hspace{0.17em}sin\hspace{0.17em}}\theta \text{\hspace{0.17em}=}\frac{\text{Δ}{u}^{2}B{\text{sin}}^2\theta }{2}{r}_M^2{\alpha }^2\sqrt{4\pi \rho }.$$(4)

Here, α is a dimensionless parameter related to the conductive properties of the planet and it can be understood as the strength of the interaction; here, we assume α = 1, thus deriving the maximum power Ƥ_wg. In the expression above, the planet is assumed to be magnetised, with a magnetospheric size r_M. Through pressure balance between the magnetic pressure of the planet and the total pressure of the stellar wind (i.e. the sum of the magnetic, thermal, and ram pressures; see the bottom-left panel of Fig. 4), we have (e.g. Vidotto et al. 2014) $$\frac{r_M}{R_p}\simeq f{\left[\frac{{\left(B_p/2\right)}^2/\left(8\pi \right)}{p_{\text{tot}}\left(a_{\text{orb}}\right)}\right]}^{1/6},$$(5)

where B_p is the polar field strength of the planet’s dipolar magnetic field (i.e. twice the value of the equatorial field strength) and f ≃ 2^2/6 is a correction factor used to account for the effects of currents (e.g. Cravens 2004). The value of B_p for GJ 436b is currently poorly constrained: Bourrier et al. (2016) derived an upper limit for the equatorial field strength of ∼3 G, consistent with the upper limit estimates of ∼10 G from Loyd et al. (2023). Assuming B_p = 2 G leads to an average magnetosphere size of 5.2 Rp and 3.7 R_p for stellar wind Models I and II, respectively, as seen on the top panel of Fig. 5. The maximum power that is radiated by Alfvén waves travelling towards the star is shown on the middle panel of this figure, where we see on average that the radiated power is Ƥ_wg ≃ 1.2 × 10²² erg s⁻¹ and 4.8 × 10²² erg s⁻¹ for Models I and II, respectively. It is worth noting that the radiated power depends on the assumed magnetic field strength of the planet indirectly through r_M. As r^M∝ $B_p^{1/3}$, we have that ${𝒫}_{\text{wg}}\propto B_p^{2/3}$. Thus, for a 20 G polar field strength (i.e. 10 G at the equator), the radiated power increases by 10^2/3, which is nearly a factor of 5, resulting in Ƥ_wg ≃ 5.6 × 10²² erg s⁻¹ and 22 × 10²² erg s⁻¹ for Models I and II, respectively.

4.2 Scenario 2: Star-planet magnetic reconnection

Magnetic reconnection occurs when magnetic field lines of different polarities interact with each other. Magnetic reconnection has been considered in the case of star-planet interactions by several authors (e.g. Jardine & Cameron 2008; Lanza 2009; Vidotto et al. 2010). While Jardine & Cameron (2008); Vidotto et al. (2010) considered reconnection events to power radio emission from exoplanets, Lanza (2009) studied how reconnection events could trigger chromospheric hot spots. We can estimate the power released in the reconnection by integrating the Poynting flux, P (Eq. (3)) over the cross-section of the planet’s magnetosphere (Lanza 2009): $${𝒫}_{\text{rec}}={\gamma }_{\text{rec}}P\left(\pi R_M^2\right)=\frac{1}{4}{\gamma }_{\text{rec}}{r}_M^2{B}^2\text{Δ}u\text{\hspace{0.17em}sin\hspace{0.17em}}\theta ,$$(6)

where γ_rec is an efficiency of the reconnection related to the angle between the stellar and planetary field lines. Even though in the GJ 436 system the planet will move through different orientations of the stellar wind magnetic field (e.g. going through closed and open field lines), and thus γ_rec would differ at different points in the orbit, for simplicity here we take its maximum value of γ_rec = 1, so that our power estimates Ƥ_rec are upper limits. The results of this scenario are shown in the bottom panel of Fig. 5, where we see that the average powers released in the reconnection scenario are Ƥ_rec ≃ 8 × 10²² erg s⁻¹ and 17 × 10²² erg s⁻¹ for Models I and II, respectively, where we assumed a planetary polar field strength of 2 G. Similarly to the previous scenario we considered, the power released in the reconnection scenario also depends on the planetary magnetic field as ${𝒫}_{\text{rec}}\propto B_p^{2/3}$. Thus, for a 20 G field we find an increase in dissipated power by nearly a factor of 5; it reaches average values of Ƥ_rec ≃ 38 × 10²² erg s⁻¹ and 80 × 10²² erg s⁻¹ for Models I and II, respectively.

Fig. 5 Strength of the stellar wind-planet interactions computed in this work. Top: magnetospheric size of the planet, assuming a planetary magnetic field strength of 2 G at the pole (i.e., 1 G at the equator). Middle: maximum dissipated power of Alfvén waves travelling towards the star (Saur et al. 2013). Bottom: maximum dissipated power from the star-planet magnetic reconnection scenario (Lanza 2009). The solid lines are for Model I, and the dashed lines for Model II. The colour scheme follows the one adopted in Fig. 3. The grey arrows in the upper horizontal axes indicate approximate times of mid-transits.

4.3 Interconnecting magnetic loop between star and planet

We calculated the power originating from star-planet interactions mediated by the stellar wind using the two scenarios discussed before. There is still a third model that we explore in the appendix, namely that of an interconnecting loop that extends from the stellar surface to the planetary orbit (Lanza 2013). For similar planetary field strengths, the interconnecting loop scenario produces maximum powers that are about three to four orders of magnitude larger than the other scenarios. This scenario assumes that the stellar magnetic field remains closed out to the orbital distance of the planet, which is not seen in our wind models.

For the three scenarios explored here (in Sect. 4 and in the appendix), the predicted powers show peaks and troughs, which are due to the planet crossing the magnetic equator of the star (when the field changes polarity). To compute the visibility of this interaction, one would need to take into account inclination of the rotation axis of the star (${35.7}_{-7.6}^{+5.9}$ deg; Bourrier et al. 2022), the timing when the planet is hidden behind the stellar disk, and possible phase lags.

5 Discussion and conclusions

In this study we performed 3D numerical simulations of the wind of GJ 436. Our model incorporates the recently reconstructed radial magnetic field of GJ 436 (Bellotti et al. 2023), which shows that, in March-June 2016, GJ 436 presented a large-scale field that largely resembled an aligned dipole, with a maximum (absolute) field strength of ~27 G (Fig. 1). Stellar outflows, in the form of winds and coronal mass ejections, shape the space weather environment around planets. From our simulations, we derived the space weather conditions around GJ436b, a warm-Neptune orbiting at a distance of 0.028 au from its host star. The main unknown in our wind models is the energy flux of the Alfvén waves that heat and accelerate the stellar wind. We ran two models with different values of the Alfvén wave fluxes to magnetic field ratio: the canonical value of S_A,0/B_⋆ = 1.1 × 10⁷ erg (cm² s G)⁻¹, commonly used in solar wind simulations, and another value that is one order of magnitude smaller. These models were named Model I (with the lower flux) and Model II (with the higher flux), as can be seen in Fig. 2. With the exception of the Alfvén wave properties, the remaining inputs of our models are well constrained by observations, including the large-scale surface field of the star. With our models, we derive stellar wind mass-loss rates of 1.1 × 10⁻¹⁵ M_⊙ yr⁻¹ and 2.5 × 10⁻¹⁴ M_⊙ yr⁻¹ for Models I and II, respectively. As discussed further below, we find that Model I is a better representation of the wind of GJ 436.

In addition to the global properties of the stellar wind, our models also allowed us to derive the conditions of the stellar wind at the orbit of GJ 436b. Even though at the observed epoch the stellar magnetic field, and thus the stellar wind, was approximately axi-symmetric about the rotation axis of the star, G_j 436b has a nearly polar orbit. This means that during one orbit, the planet probes regions of open field lines and closed field lines. Taking the spin-orbit angle into account, we calculated the position of the planet during consecutive orbits: during one stellar rotation (~44 days) the planet orbits about ~16 times, experiencing different space weather conditions during each orbit (see Figs. 3 and 4). There are two main sources of variation in these conditions: one relates to the choice of wind model (lower versus higher Alfvén wave fluxes), and the other relates to the high obliquity of the orbit of GJ 436b. For example, the many peaks and troughs we see in the local wind properties (Fig. 4) are associated with the planet crossing the stellar equatorial plane.

Regardless of the two wind models we chose, our 3D simulations indicate that the motion of GJ 436b through the stellar wind is sub-Alfvénic, which implies that the planet has a direct magnetic connectivity with its host star. In a second part of our study, we investigated the power released in star-planet interactions mediated by the magnetised stellar wind. We first computed the Poynting fluxes of the stellar wind that impacts GJ 436b to be P ≃ 160 and 650 erg cm⁻² s⁻¹ (or 0.1 to 0.6 W m⁻²) for Models I and II, respectively. Next, we computed the maximum energy that is available to power sub-Alfvénic star-planet interactions through two different scenarios: the Alfvén wing scenario (in which perturbations caused in the stellar magnetic field line by the planet travel towards the star as MHD waves, irradiating their energies; Saur et al. 2013) and the magnetic reconnection scenario (in which magnetic reconnection between stellar and planetary magnetic fields releases energy; Lanza 2009). For both scenarios, we assumed that the planet is magnetised with an equatorial dipolar magnetic field of 1G, leading to a magnetospheric size of about 5.2 Rp and 3.7 Rp for stellar wind Models I and II, respectively. Our estimated maximum powers released through star-planet interactions are of the order of 10²² to 10²³ erg s⁻¹ (approximately 10⁻¹⁰ to 10⁻⁹ of the stellar bolometric luminosity), with the largest values associated with the reconnection scenario. It is interesting to note that the powers we derived are not very sensitive to the choice of Alfvén wave surface flux in the wind model (the largest unknown in our wind models) – a change of one order of magnitude in the surface flux of Alfvén waves leads to average powers that are a factor of 4 different for the Alfvén wing scenario and a factor of 2 in the reconnection scenario only (Fig. 5).

We note here that the planetary magnetic field is an assumption in our calculations and that the power released in the two star-planet interaction scenarios we considered increases with $B_p^{2/3}$. Thus, increasing the magnetic field of the planet by a factor of 10 (i.e. to an equatorial field of 10 G) leads to an increase in the estimated powers by 10^2/3 (i.e. nearly a factor of 5). We could keep increasing the choice of planetary field strength until our estimated powers increased from ~ 10²³ erg s⁻¹ to 10²⁶ erg s⁻¹, the latter of which is the maximum power derived by Loyd et al. (2023) for star-planet interactions in the GJ 436 system. To match the observed upper limit, the maximum planetary magnetic field would need to be unreasonably high (≲10^3×3/2 ≃ 3 × 10⁴ G). Because of this, we unfortunately cannot derive any meaningful upper limit to the planetary magnetic field.

According to the models studied here, GJ 436b orbits within the Alfvén surface (i.e. it has a sub-Alfvénic orbital motion) except for small portions of the orbit around the magnetic equator, where the orbit is outside the Alfvén surface. When a planet moves from sub- to super-Alfvénic motion, different types of star–planet interactions can take place. The formation of a bow shock around the planet should be seen as soon as the planet moves from sub- to super-Alfvénic motion (e.g. Vidotto et al. 2011a). Cohen et al. (2014) studied the sub- to super-Alfvénic transition in terms of joule heating in the upper atmosphere of planets, while Kislyakova et al. (2018) and Kislyakova & Noack (2020) demonstrated that planets orbiting in the varying magnetic fields of their host stars can also experience induction heating in their interiors. Additionally, it is believed that magnetic connectivity can switch between on and off states – this can be due to variations over one planetary year (e.g. the planet moving in and out of the Alfvén surface) or variations in the Alfvén surface due to the evolution of the surface magnetic field (Fares et al. 2010). The latter fact has been attributed as the cause of the disappearance of planet-induced chromospheric modulation in planet-hosting stars (Shkolnik et al. 2008), emphasising the importance of contemporaneous observations of stellar magnetism and the search for signatures of star-planet interactions.

We explore a third star-planet interaction scenario in the appendix of this paper – in this scenario, a magnetic loop connecting the star and the planet would suffer a rupture after being stretched by the planet (Lanza 2013). Of the three scenarios, this has been recognised as the one that releases the largest star-planet interaction powers (Cauley et al. 2019; Strugarek et al. 2022). However, this scenario is not applicable to the case where the star–planet interaction is mediated by a stellar wind that distorts the stellar magnetic field from a potential state. In this case, in the calculations shown in the appendix, a potential field model for the stellar magnetic field is assumed (without a stellar wind). Using this scenario, we demonstrate in the appendix that a planetary field strength of ∼6 G would reproduce the powers found by Loyd et al. (2023), supporting with these authors’ estimates.

In our work we adopted two different values of the Alfvén wave flux that drives the stellar wind. One way to constrain the best choice of surface Alfvén wave flux, and thus decide which of the wind models is more realistic, is to compare the mass-loss rate and wind properties we predict with the values obtained from models of the transmission spectroscopic transit of GJ 436b from Bourrier et al. (2016). Comparing the wind velocities and densities that we derived to those derived in Bourrier et al. (2016, shown as grey crosses in the upper panels of Fig. 4), we find that our Model I, with the lowest surface flux of Alfvén waves, matches the spectroscopic transit values better. Our derived mass-loss rate of 1.1 × 10⁻¹⁵ M_⊙ yr⁻¹ is compatible with that derived from spectroscopic transits of ${1.2}_{-0.75}^{+1.3}\times {10}^{-15}{M}_{\odot }{\text{yr}}^{-}$ (Vidotto & Bourrier 2017). In addition to this, our preference towards Model I is also supported by further empirical and theoretical works, as we explain below.

Empirically, for the Sun, it has been suggested that the non-thermal velocities, ξ, of lines formed at the transition region are associated with the rms amplitude of the Alfvén waves (e.g. Oran et al. 2017) as $\xi ={\langle {\upsilon }_{\perp ,0}^2\rangle }^{1/2}/\sqrt{2}$. If this can be extended to solar-like stars, it implies that by measuring ξ one can constrain the wave properties of other stars. Recently, Boro Saikia et al. (2023) measured the non-thermal velocity amplitudes of a sample of 55 low-mass stars. Overall, they find an increase in ξ with stellar rotation (or activity). In their sample, they have five stars with masses, radii, and effective temperatures similar to those of GJ 436, namely: GJ 832, LTT 2050, GJ 3470, GJ 176, and GJ 849. Three of these stars (GJ 832, GJ 176, and GJ 849) have a rotation period of about 40 days, very similar to GJ 436. Using the Si IV line (1393 Å), Boro Saikia et al. (2023) measured ξ ranging between 14.7 and 17.8 km s⁻¹ for these objects. These values are very similar to the non-thermal velocity of 14 km s⁻¹ from Model I.

Theoretically, it is usually assumed that the wave flux is given by a combination of the available turbulent energy flux inside a magnetic flux tube times the magnetic filling factor (e.g. Cranmer 2017; Suzuki 2018). The magnetic filling factor is smaller for slowly rotating stars (See et al. 2019), which are also usually older. Compared to the Sun, the turbulent energy flux is also lower in M dwarfs, given that the energy flux emerging from convective motions (and thus eventually driving the waves) is proportional to the bolometric energy flux of the star, $\sigma T_{\text{eff}}^4$ (e.g. Suzuki 2018). Sakaue & Shibata (2021) performed a parametric study of wave-driven wind models of M dwarfs, showing that the wave amplitude in the lower corona is smaller for M dwarfs than for the Sun.

In conclusion, both models (Sakaue & Shibata 2021) and empirical measurements (Boro Saikia et al. 2023) of the wave amplitude point to lower wave flux values for slowly rotating M dwarfs. Based on these arguments and the agreement of the local wind values we derive here with results from Ly-α transit modelling (Bourrier et al. 2016; Vidotto & Bourrier 2017), we favour a lower flux of Alfvén waves for the case of GJ 436. This might imply that the canonical value of S_A,0/B_⋆ = 1.1 × 10⁷ erg (cm² s G)⁻¹, derived from solar wind models, might not not applicable to the winds of old M dwarf stars.

Our models indicate that the planet GJ 436b experiences a vast range of stellar wind speeds and densities along its orbit. Because the wind shapes the atmospheric material that is being evaporated from the planet, we expect variations in the interaction region between the escaping atmosphere and the stellar wind, and hence, variations in the large structures formed due to this interaction, such as comet-like tails. We see the effects of different stellar wind properties on the morphology of comet-like tails, for example in 3D models that investigate these interactions (e.g. Villarreal D‘Angelo et al. 2014, 2021; McCann et al. 2019; Carolan et al. 2021; Rumenskikh et al. 2022). Time-dependent effects, such as coronal mass ejections, would also alter the morphology of comet-like tails (Cherenkov et al. 2017; Hazra et al. 2022), similar to what we see when the solar wind interacts with the plasma tail of comets (Vourlidas et al. 2007). Because of the large range of stellar wind properties experienced by GJ 436b, similar temporal effects are expected during its orbit. It is interesting to note, however, the remarkable stability of the wind properties at every consecutive mid-transit time (crosses in the upper panel of Fig. 4). If the local wind properties are relatively constant for a few hours, this would lead to a relatively stable interaction zone, potentially explaining why the Ly-α transit of GJ 436b remains stable over time (Dos Santos et al. 2019).

For this stability to occur in transits observed many years apart, the stellar magnetic field cannot vary substantially. To investigate this, we need further spectropolarimetric campaigns to reconstruct the large-scale field of the star and thus study its evolution. We note that GJ436 has an activity cycle with a period of about 7–8 yr (Lothringer et al. 2018; Dos Santos et al. 2019; Kumar & Fares 2023; Loyd et al. 2023), which implies that the stellar magnetic field could vary (in strength and topology) substantially at different phases of the cycle. This temporal variation would affect the wind properties and space weather conditions around the planet (e.g. Vidotto et al. 2012; Nicholson et al. 2016; Finley et al. 2018; Kavanagh et al. 2019) and, consequently, the power released during star-planet interactions (e.g. Klein et al. 2022). Ideally, observations that probe star-planet interactions mediated by the stellar wind and/or magnetic field (e.g. planetary radio emission, chromospheric hot spots, or the interaction of evaporated atmospheres with stellar winds) and the magnetic map used in stellar wind models should be conducted contemporaneously.

Acknowledgements

A.A.V. thanks P. Loyd and A. Lanza for useful discussions in the preparation of this paper. The authors thank the anonymous reviewer for their constructive recommendations and suggestions. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 817540, ASTROFLOW; grant agreement No 947634, SPICE DUNE; grant agreement No 740651, NewWorlds). RF acknowledges support from the United Arab Emirates University (UAEU) startup grant number G00003269. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. We acknowledge funding from the French National Research Agency (ANR) under contract number ANR-18-CE31-0019 (SPlaSH). We thank SURF (www.surf.nl) for the support in using the National Supercomputer Snellius. This work used the BATS-R-US tools developed at the University of Michigan Center for Space Environment Modeling and made available through the NASA Community Coordinated Modeling Center.

Appendix A Calculation of orbital trajectories for a misaligned orbit

To calculate the orbital trajectories in the reference frame of the star, we compute two coordinate transformations: one from the orbital reference frame to the observer’s frame, and another from the observer’s frame to the stellar frame. These frames are illustrated in Fig. A.1.

We start from the orbital reference frame {x_pl, y_pl, z_pl}, such that the x_ply_pl-plane contains the orbit of the planet and the orbital axis is along z_pl. The origin of this reference frame is centred on the star. In this reference frame, the (circular) planetary motion is simply described as $$\left[\begin{array}{c}{x}_{\text{pl}}\\ y_{\text{pl}}\\ z_{\text{pl}}\end{array}\right]=\left[\begin{array}{c}{R}_{\text{orb}}\cos {\phi }_{\text{pl}}\\ R_{\text{orb}}\text{\hspace{0.17em}sin\hspace{0.17em}}{\phi }_{\text{pl}}\\ 0\end{array}\right],$$(A.1)

where the phase of the orbit is φ_pl = Ω_orbt + φ₀, with t describing the time and Ω_orb = 2π/P_orb the orbital rotation rate. Here, we assume a non-null initial phase φ₀, which we computed later in this section to find the location of GJ 436b at zero phase of the surface magnetic map, which corresponds to the first spectropolarimetric observation as reported in Bellotti et al. (2023).

Now we take the observer’s (inertial) reference frame {x_obs, y_obs, z_obs}, whose origin is also centred on the star, and z_obs is aligned with the stellar rotation. The observer is located at positive x_obs and we assume that x_obs || x_pl. The angle between the orbital spin and the stellar spin axis is the spin-orbit angle, Ψ; thus, the angle between z_obs and z_pl is also Ψ. In this reference frame, the orbital motion of the planet is described as $$\left[\begin{array}{c}{x}_{\text{obs}}\\ y_{\text{obs}}\\ z_{\text{obs}}\end{array}\right]=\left[\begin{array}{ccc}1& 0& 0\\ 0& \text{cos\hspace{0.17em}Ψ}& \sin \text{\hspace{0.17em}Ψ}\\ 0& -\sin \text{\hspace{0.17em}Ψ}& \cos \text{\hspace{0.17em}Ψ}\end{array}\right]\left[\begin{array}{c}{x}_{\text{pl}}\\ y_{\text{pl}}\\ z_{\text{pl}}\end{array}\right],$$(A.2)

where we performed an anti-clockwise rotation around the axis x_obs = x_pl by an angle Ψ. According to our reference frame, the mid-transit occurs when the planet is at {x_obs, y_obs, z_obs} = {R_orb, 0,0}.

Fig. A.1 Reference frames used in the coordinate transformation, with the star given by the orange circle, the planet indicated by the red circle, and its orbital path shown in red. The planet’s reference frame {xpl, ypl, zpl} (in red) contains the orbital plane of the planet xplypl and orbital spin axis along zpl. The observer’s inertial reference frame {xobs, yobs, zobs} (in blue) is centred on the star, with the observer located at +xobs and the stellar rotation axis aligned along zobs. The reference frame {x⋆, y⋆, z⋆} (in black) co-rotates with the star. The spin-orbit angle is Ψ.

The reference frame co-rotating with the star {x_⋆, y_⋆, z_⋆} is also centred on the star, with the stellar spin axis along z_⋆, which is parallel to z_obs. To transform from the observer’s reference frame to the stellar co-rotating frame, we performed a clockwise rotation around the axis z_⋆ = z_obs by an angle Ω_⋆ t: $$\left[\begin{array}{c}{x}_{\star }\\ y_{\star }\\ z_{\star }\end{array}\right]=\left[\begin{array}{ccc}\text{cos}\left({\text{Ω}}_{\star }t\right)& -\text{sin}\left({\text{Ω}}_{\star }t\right)& 0\\ \text{sin}\left({\text{Ω}}_{\star }t\right)& \text{cos}\left({\text{Ω}}_{\star }t\right)& 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}_{\text{obs}}\\ y_{\text{obs}}\\ z_{\text{obs}}\end{array}\right].$$(A.3)

Therefore, to derive the orbital path of the planet in the reference frame of the star, we substitute (A.1) and (A.2) into (A.3), thus obtaining $$\frac{x_{★}}{R_{\text{orb}}}=\cos \left({\text{Ω}}_{\text{orb}}t+{\phi }_0\right)\cos \left({\text{Ω}}_{★}t\right)-\sin \left({\text{Ω}}_{\text{orb}}t+{\phi }_0\right)\sin \left({\text{Ω}}_{★}t\right)\cos \text{Ψ}$$(A.4) $$\frac{y_{★}}{R_{\text{orb}}}=\cos \left({\text{Ω}}_{\text{orb}}t+{\phi }_0\right)\sin \left({\text{Ω}}_{★}t\right)+\sin \left({\text{Ω}}_{\text{orb}}t+{\phi }_0\right)\cos \left({\text{Ω}}_{★}t\right)\cos \text{Ψ}$$(A.5) $$\frac{z_{★}}{R_{\text{orb}}}=-\sin \left({\text{Ω}}_{\text{orb}}t+{\phi }_0\right)\sin \text{Ψ.}$$(A.6)

Figure 3 shows the trajectories of the planet GJ 436b as seen in the reference frame of the star during one rotational period of the star (about 44 days)

We note that, in the case of aligned systems, Ψ = 0 and we have $$x_{★}=R_{\text{orb}}\cos \left({\text{Ω}}_{\text{orb}}t+{\phi }_0+{\text{Ω}}_{★}t\right)$$(A.7) $$y_{★}=R_{\text{orb}}\sin \left({\text{Ω}}_{\text{orb}}t+{\phi }_0+{\text{Ω}}_{★}t\right)$$(A.8) $$z_{★}=0.$$(A.9)

Said differently, in aligned systems, the planet re-encounters the same stellar wind property once every 2π/(Ω_orb + Ω_⋆), as long as the stellar magnetic field does not evolve significantly.

Appendix B Interconnecting loop model

A third star-planet scenario we also explore in this paper is based on the idea presented by Lanza (2013), in which a magnetic loop connecting the star and the planet would suffer a rupture after being stretched by the planet. The idea is that one footpoint of the loop sits at the surface of the star, while the other footpoint lies at the surface of the planet. As the planet moves through its orbit, the magnetic loop is stretched (stressed) and when it breaks, it can release energy, which travels towards the star, giving rise for example to anomalous hot spots at the stellar surface (Shkolnik et al. 2005; Cauley et al. 2019). This is the scenario that Strugarek et al. (2022) named ’stretch and break’; in their most recent work on star-planet interactions in the GJ 436 system, this model was also investigated by Loyd et al. (2023). In this scenario, the model assumes that the stellar closed corona extends up to the distance of the planet, which is not the case in the stellar wind models we presented here – as shown in Sect. 3, the stellar wind stretches open the magnetic field lines of the star at distances within the planetary orbital distance. Additionally, the interconnecting loop scenario assumes that the stellar magnetic field is close to a potential field, so that it can topologically reconnect with the magnetic field of the planet that is also potential. In the presence of a stellar wind though, the coronal magnetic field departs from a potential state (Vidotto et al. 2011b).

Fig. B.1 Interconnecting loop model. Top: Magnetic field of the star at the orbital distance of the planet, assuming a potential field extrapolation, with a source surface of ≃ 15R⋆ (in this case, no wind model is used). Bottom: Maximum dissipated power in the interconnecting loop model, assuming a polar planetary magnetic field strength (dipole) of 2 G (dashed) and 6 G (solid). The colour scheme follows the colour bar shown in Fig. 3. The grey arrows in the upper horizontal axes indicate approximate times of mid-transits.

To explore the interconnecting loop scenario, we ran a model without the presence of a stellar wind, in which we assume that the stellar magnetic field remains potential and closed up to the orbit of the planet. To compute the field strength at the orbit of the planet, we use a potential field source surface extrapolation. In this model, the source surface represents the distance above which the magnetic field lines become open. To make sure the stellar magnetic field lines are closed up to the planet’s orbit, we chose a source surface of ≃ 15R_⋆ (i.e. beyond the planet’s orbit). Using the surface magnetic map from Fig. 1 and our potential field model, the magnetic field of the star at the orbit of the planet is shown in the top panel of Fig. B.1, where we see that the coronal magnetic field is on average 7.5 × 10⁻³ G.

With this, we can then compute the power released in the interconnecting loop model using (Lanza 2013) $${𝒫}_{\text{sb}}=\left(2\pi R_p^2\right)f_{\text{open}}\left[\frac{B_p^2{\upsilon }_K}{4\pi }\right],$$(B.1)

where the Poynting flux (term within brackets) is now across the planet. The fraction of the planetary surface area that has open magnetic field lines (Adams 2011) and thus have interconnecting magnetic loops is $$f_{\text{open}}=1-{\left(1-\frac{3{\varsigma }^{1/3}}{2+\varsigma }\right)}^{1/2},$$(B.2)

with ζ = B/B_p. Here, B is the stellar magnetic field extrapolated out to the orbit of GJ 436b and Bp, as before, is the assumed polar magnetic field strength of a dipolar planetary magnetic field. In the interconnecting loop scenario, the atmosphere of the planet is assumed to be ionised down to the surface, where the planetary magnetic field is B_p. If the atmosphere becomes neutral at a certain height above the surface, then the value of Bp is that at such a height, which is smaller than the surface field. Therefore, Eq. (B.1) provides an upper limit of the power released in the interaction. Assuming B_p = 2 G, we find that the power released in the interconnecting loop model is on average 1.7 × 10²⁵ erg s⁻¹. For a magnetic field of B_p = 6 G, the power released increases to 10 × 10²⁵ erg s⁻¹ (this value is similar to the maximum power derived by Loyd et al. 2023). This is shown in the bottom panel of Fig. B.1.

We confirm the findings of previous studies (Cauley et al. 2019; Strugarek et al. 2022), which showed that the interconnecting loop scenario provides the largest maximum power among the three scenarios explored in this paper (see the first two scenarios in Sect. 4), and is also in line with the estimates provided in Lanza (2013). According to the interconnecting loop scenario, a planetary magnetic field B_p ≲ 6 G can reproduce the maximum power observed in star-planet interactions in the GJ 436 system (Loyd et al. 2023).

References

All Tables

All Figures

	Fig. 1 Observationally reconstructed radial magnetic field of GJ 436 from Bellotti et al. (2023) using Zeeman-Doppler imaging. This map is incorporated at the boundary of our stellar wind simulations.
In the text

Fig. 2 3D view of the magnetised stellar wind for Models I and II. The surface magnetic field from March-June 2016 reconstructed by Bellotti et al. (2023) is used at the boundary of the models. The grey streamlines are the magnetic field lines embedded in the wind. The rotation axis of the star is along the z⋆-axis. The 2D cut is placed at the x⋆z⋆ plane, and for reference the orbital distance of the planet is shown with the translucent sphere at 14.56 R⋆ ≃ 0.028 au. The contour shows the wind velocity, and the Alfvén surface is shown by the translucent outer surface. Models I and II predict stellar wind mass-loss rates of 1.1 × 10−15 M⊙ yr−1 and 2.5 × 10−14 M⊙ yr−1.

In the text

	Fig. 3 Trajectories of the planet as seen in the reference frame of the star during one rotational period of the star (about 44 days; see the colour bar), during which the planet completes about 16 orbital revolutions. The grey circles indicate the approximate times of mid-transit.
In the text

Fig. 4 Stellar wind proton density, velocity, total pressure (sum of ram, including the planet‘s own orbital motion and the thermal and magnetic pressures), and total magnetic field strength at the orbit of GJ 436b. The solid lines are for Model I and the dashed lines for Model II. Colour represents the time evolution and follows the colour scheme shown in Fig. 3. The horizontal grey lines are the time average of the stellar wind properties over one stellar rotation period (about 16 orbital periods of the planet). The grey arrows in the upper horizontal axes indicate approximate times of mid-transits. In the upper panels, the crosses indicate the wind values derived in the work of Bourrier et al. (2016), and they are plotted at the mid-transit times. Note how they overlap really well with the values derived in our Model I, indicating that Model I provides a better description of the wind of GJ 436.

In the text

Fig. 5 Strength of the stellar wind-planet interactions computed in this work. Top: magnetospheric size of the planet, assuming a planetary magnetic field strength of 2 G at the pole (i.e., 1 G at the equator). Middle: maximum dissipated power of Alfvén waves travelling towards the star (Saur et al. 2013). Bottom: maximum dissipated power from the star-planet magnetic reconnection scenario (Lanza 2009). The solid lines are for Model I, and the dashed lines for Model II. The colour scheme follows the one adopted in Fig. 3. The grey arrows in the upper horizontal axes indicate approximate times of mid-transits.

In the text

Fig. A.1 Reference frames used in the coordinate transformation, with the star given by the orange circle, the planet indicated by the red circle, and its orbital path shown in red. The planet’s reference frame {xpl, ypl, zpl} (in red) contains the orbital plane of the planet xplypl and orbital spin axis along zpl. The observer’s inertial reference frame {xobs, yobs, zobs} (in blue) is centred on the star, with the observer located at +xobs and the stellar rotation axis aligned along zobs. The reference frame {x⋆, y⋆, z⋆} (in black) co-rotates with the star. The spin-orbit angle is Ψ.

In the text

Fig. B.1 Interconnecting loop model. Top: Magnetic field of the star at the orbital distance of the planet, assuming a potential field extrapolation, with a source surface of ≃ 15R⋆ (in this case, no wind model is used). Bottom: Maximum dissipated power in the interconnecting loop model, assuming a polar planetary magnetic field strength (dipole) of 2 G (dashed) and 6 G (solid). The colour scheme follows the colour bar shown in Fig. 3. The grey arrows in the upper horizontal axes indicate approximate times of mid-transits.

In the text