11institutetext: Leiden Observatory, Leiden University, PO Box 9513, 2300 RA,
Leiden, the Netherlands
11email<EMAIL_ADDRESS>22institutetext: Observatoire Astronomique
de l’Université de Genéve, Chemin Pegasi 51b, 1290 Versoix, Switzerland
33institutetext: Department of Physics, College of Science, United Arab
Emirates University, P.O. Box No. 15551, Al Ain, UAE
33email<EMAIL_ADDRESS>44institutetext: Science Division, Directorate
of Science, European Space Research and Technology Centre (ESA/ESTEC),
Keplerlaan 1, 2201 AZ, Noordwijk, The Netherlands 55institutetext: Institut de
Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS,
IRAP/UMR 5277, 14 avenue Edouard Belin, F-31400, Toulouse, France
66institutetext: Laboratoire Univers et Particules de Montpellier, Université
de Montpellier, CNRS, F-34095, Montpellier, France

# The space weather around the exoplanet GJ 436 b.
II. Stellar wind–exoplanet interactions

A. A. Vidotto 11 0000-0001-5371-2675 V. Bourrier 22 0000-0002-9148-034X R.
Fares 33 0000-0002-3301-341X S. Bellotti 114455 0000-0002-2558-6920 J. F.
Donati 55 0000-0001-5541-2887 P. Petit 55 0000-0001-7624-9222 G. A. J. Hussain
44 0000-0003-3547-3783 J. Morin 66 0000-0002-4996-6901

(Received ; accepted )

The M dwarf star GJ 436 hosts a warm-Neptune that is losing substantial amount
of atmosphere, which is then shaped by the interactions with the wind of the
host star. The stellar wind is formed by particles and magnetic fields that
shape the exo-space weather around the exoplanet GJ 436 b. Here, we use the
recently published magnetic map of GJ 436 to model its 3D Alfvén-wave driven
wind. By comparing our results with previous transmission spectroscopic models
and measurements of non-thermal velocities at the transition region of GJ 436,
our models indicate that the wind of GJ 436 is powered by a smaller flux of
Alfvén waves than that powering the wind of the Sun. This suggests that the
canonical flux of Alfvén waves assumed in solar wind models might not be
applicable to the winds of old M dwarf stars. Compared to the solar wind, GJ
436’s wind has a weaker acceleration and an extended sub-Alfvénic region. This
is important because it places the orbit of GJ 436 b inside the region
dominated by the stellar magnetic field (i.e., inside the Alfvén surface). Due
to the sub-Alfvénic motion of the planet through the stellar wind,
magnetohydrodynamic waves and particles released in reconnection events can
travel along the magnetic field lines towards the star, which could power the
anomalous ultraviolet flare distribution recently observed in the system. For
an assumed planetary magnetic field of $B_{p}\simeq 2$ G, we derive the power
released by stellar wind-planet interactions as $\mathcal{P}\sim 10^{22}$ –
$10^{23}$ erg s-1, which is consistent with the upper limit of $10^{26}$ erg
s-1 derived from ultraviolet lines. We further highlight that, because star-
planet interactions depend on stellar wind properties, observations that probe
these interactions and the magnetic map used in 3D stellar wind simulations
should be contemporaneous for deriving realistic results.

###### Key Words.:

Stars: winds, outflows – Stars: individual: GJ 436 – planet-star interactions
– planetary systems

##  1 Introduction

The evolution of planets is shaped by billions of years of interactions with
their host stars and also 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 436 b 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 subsequent interaction of this atmosphere with
the stellar wind (Bourrier et al., 2016). Spectroscopic transit observations
in Lyman-$\alpha$ (Kulow et al., 2014; Ehrenreich et al., 2015; Lavie et al.,
2017; Dos Santos et al., 2019) revealed that GJ 436 b is enshrouded by a giant
cloud of escaping atmosphere that trails GJ 436 b 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 436 b 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; these 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 rather 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 to 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 can be found 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 the 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 magnetohydrodynamics
(MHD) waves that can travel towards the star, or alternatively (and even
possibly concurrently) it can also trigger magnetic reconnection events in
case of a magnetised planet. Both of these mechanisms release energy that are
radiated away and could be detected in the form of, e.g., 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.
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 the planet
through the wind of its host star. In this regime, the interaction between a
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 experienced by the solar-system
planets currently, 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). In spite of many
attempts to detect such signature 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, i.e., sub- or super-Alfvénic. To
determine such a 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
relevant also for, e.g., 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), which allows us to now model
for the first time the three-dimensional (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 (Section 3) along the highly
misaligned orbit of GJ 436 b. GJ 436 b is a Neptune-mass planet orbiting very
close to its slowly rotating M dwarf host. It has a puzzling observed orbital
architecture – 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 Section 4, we then compute the strengths of star-planet
interactions. In the final section of this paper (Section 5), we discuss how
our computations compare to the recent observations of Loyd et al. (2023) and
present our concluding thoughts.

##  2 The wind of GJ 436

The environment surrounding GJ 436 b 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 Figure
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 use 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.

Figure 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 addition to the stellar surface magnetic map, other model inputs are the
stellar properties (mass, radius and rotation period; c.f. Table1), 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_{\rm tr}=5\times
10^{11}\leavevmode\nobreak\ \mathrm{cm}^{-3}$), temperature ($T_{\rm
tr}=5\times 10^{4}$ K), and, consequently, gas density ($P_{\rm tr}\simeq
1.4\leavevmode\nobreak\ \mathrm{dyn\leavevmode\nobreak\ cm}^{-2}$), 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_{\rm tr}=0.8\leavevmode\nobreak\ \mathrm{dyn\leavevmode\nobreak\
cm}^{-2}$. For a temperature of $5\times 10^{4}$ K, this pressure value
translates to a density of $10^{11}\leavevmode\nobreak\ \mathrm{cm}^{-3}$,
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_{\rm tr}=0.02$
– $0.04\leavevmode\nobreak\ \mathrm{dyn\leavevmode\nobreak\ cm}^{-2}$ (Mauas
et al., 1997; Peacock et al., 2019), which would imply transition region
densities of a few times $10^{9}\leavevmode\nobreak\ \mathrm{cm}^{-3}$ at a
temperature of $5\times 10^{4}$ K. Depending on which chromospheric model we
use to guide our assumption of $n_{\rm tr}$, we could have transition region
densities that vary from $\sim 10^{9}$ to $10^{11}\mathrm{cm}^{-3}$. O
Fionnagain et al. (2021) explored the effects that changing this parameter had
on the wind of a sub-giant star $\lambda$ And and they found that varying
$n_{\rm 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_{\rm tr}$ within the quoted
range would change our wind solutions substantially.

Table 1: Parameter adopted in our models for GJ 436 and GJ 436 b. Parameter (Star) | value | units | reference
---|---|---|---
Mass ($M_{\star}$) | $\simeq 0.45$ | $M_{\odot}$ | Maxted et al. (2022)
Radius ($R_{\star}$) | $\simeq 0.42$ | $R_{\odot}$ | Maxted et al. (2022)
Rotation period ($P_{\rm rot}$) | $\simeq 44.1$ | days | Bourrier et al. (2018)
Stellar magnetism | c.f. Fig. 1 |  | Bellotti et al. (2023)
Age | $6^{+4}_{-5}$ | Gyr | Torres (2007)
Parameter (Planet) | value | units | reference
Mass ($M_{\rm pl}$) | 21.72 | $M_{\oplus}$ | Maxted et al. (2022)
Radius ($R_{\rm pl}$) | 3.85 | $R_{\oplus}$ | Maxted et al. (2022)
Orbital period ($P_{\rm orb}$) | $\simeq 2.644$ | days | Lanotte et al. (2014)
Semi major axis ($a_{\rm orb}$) | 14.56 | $R_{\star}$ | Maxted et al. (2022)
… | $\simeq 0.028$ | au | …
Spin-orbit angle ($\Psi$) | $103.2$ | degrees | Bourrier et al. (2022)

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 being the same as that used in the AWSoM solar wind
models, namely $L_{\perp,0}\sqrt{B_{\star}}=1.5\times
10^{5}\leavevmode\nobreak\ \mathrm{m}\sqrt{\mathrm{T}}=1.5\times
10^{4}\leavevmode\nobreak\ \mathrm{km}\sqrt{\mathrm{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 $\sim 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_{\rm tr}\langle v_{\perp,0}^{2}\rangle v_{A,0}\,,$ (1)

where $\langle v_{\perp,0}^{2}\rangle^{1/2}$ is the average velocity
perturbation of the Alfvén wave, $\rho_{\rm tr}=m_{p}n_{\rm tr}$ is the mass
density at the transition region, and the Alfvén velocity is
$v_{A,0}=B_{\star}/\sqrt{4\pi\rho_{\rm tr}}$, where the index ‘0’ indicates
value computed at the boundary of the simulation. Here, $B_{\star}$ is the
local value (i.e., at each longitude and latitude) of the surface magnetic
field as derived in the surface magnetic map (Figure 1). The previous equation
can be rewritten as

$\frac{\langle v_{\perp,0}^{2}\rangle^{1/2}}{65\leavevmode\nobreak\
\textrm{km/s}}=\left(\frac{5\times 10^{11}}{n_{\rm
tr}[\textrm{cm}^{-3}]}\right)^{1/4}\left(\frac{S_{A,0}/B_{\star}[\textrm{erg/(cm}^{2}\textrm{\,s\,G})]}{1.1\times
10^{7}}\right)^{1/2}\,.$ (2)

Note that the ratio $S_{A,0}/B_{\star}$ 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 ($\rho_{\rm tr}\langle v_{\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 to 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 $\xi$ is related to the
average velocity perturbation of the Alfvén wave as
$\xi^{2}=\frac{1}{2}\langle v_{\perp,0}^{2}\rangle$.

The ratio $S_{A,0}/B_{\star}$ is the most important free parameter in our
model and its value affects the wind temperature, its mass-loss rate and 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 run two wind models with two different values of
$S_{A,0}/B_{\star}$: $1.1\times 10^{6}\leavevmode\nobreak\ \textrm{erg
(cm}^{2}\textrm{\,s\,G)}^{-1}$ and $1.1\times 10^{7}\leavevmode\nobreak\
\textrm{erg (cm}^{2}\textrm{\,s\,G)}^{-1}$. 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 will
show in Section 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 into correlation lengths of $L_{\perp,0}\simeq
2800\leavevmode\nobreak\ \mathrm{km}$,111The scaling of the correlation length
$L_{\perp,0}\sqrt{B_{\star}}$ has been related to the distance between
magnetic flux tubes on the solar surface (Hollweg, 1986). Hollweg (1986)
empirically estimated a value of
$L_{\perp,0}\sqrt{B_{\odot}}=7520\leavevmode\nobreak\ {\rm km}\sqrt{\rm G}$
for the Sun, or $L_{\perp,0}\sim 4300$ km for $\sim 3$-G field. We do not know
how $L_{\perp,0}$ would change for GJ 436, but because of its smaller radius,
we naively expect that the distance between magnetic flux tubes would be
smaller than solar for GJ 436, consistent to the value of $\sim
2800\leavevmode\nobreak\ \mathrm{km}$ we use. Based on the parametric study of
O Fionnagain et al. (2021), we do not expect a variation in this parameter by
a factor of $\sim 10$ would substantially affect the general properties of our
wind solution. and surface Alfvén fluxes of up to $S_{A,0}\simeq[3\times
10^{7},3\times 10^{8}]\leavevmode\nobreak\ \textrm{erg
cm}^{-2}{\mathrm{s}}^{-1}$, for Models I and II, respectively. The surface
amplitudes of the wave perturbation are 21 km s-1 (Model I) and 65 km s-1
(Model II), uniform at the lower boundary. The corresponding surface energy
densities of the wave are 11 and 35 erg cm-3. With these parameters, the
predicted non-thermal velocities at the transition region are $\xi\simeq 14$
and $45$ km s-1, 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_{\star}$ – 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\times
10^{-15}\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$ and
$2.5\times 10^{-14}\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\
\rm{yr}^{-1}$, for Models I and II, respectively. For comparison, the solar
wind mass-loss rate is $2\times 10^{-14}\leavevmode\nobreak\
M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$ – 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,222From the magnetic
map of GJ 436, we find a surface magnetic flux of $\Phi_{0}=1.4\times 10^{23}$
Mx, which is very similar to solar values (for an average large-scale field
strength of $\sim 3$ G, the solar surface flux is $\Phi_{0}=3\times 4\pi
R_{\odot}^{2}\simeq 1.8\times 10^{23}$ Mx). In stellar wind theory, the mass-
loss rate is linearly proportional to the open magnetic flux ($\Phi_{\rm
open}$), which is also proportional to $\Phi_{0}$ (e.g., Vidotto et al.,
2014). This is also seen in AWSoM models of winds of Sun-like stars
(Evensberget et al., 2023). Because Model II uses the same solar wind
parameters and GJ436 has a similar surface flux to the Sun, this model
produces a mass-loss rate that is thus similar to that of the Sun. it is not
surprising that the mass-loss rate derived in Model II is similar to that of
the solar wind. In Section 3, we will 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_{\star}z_{\star}$ 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 436 b (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,
2022), for example. We will come back to this in Section 4.

Figure 2: Three dimensional 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_{\star}$-axis. The two-dimensional cut is
placed at the $x_{\star}z_{\star}$ plane, and for indication, the orbital
distance of the planet is shown by the translucent sphere at
$14.56\leavevmode\nobreak\ R_{\star}\simeq 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\times
10^{-15}\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$ and
$2.5\times 10^{-14}\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\
\rm{yr}^{-1}$.

##  3 Space weather conditions of GJ 436 b

The space weather of GJ 436 b 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 436 b. Because GJ 436 b has a very misaligned orbit, we
first compute the position of the planet as a function of time, following the
coordinate transformations presented in Section 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}=2458947.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}=2457464.4967$, corresponding to the first
spectropolarimetric observation (Bellotti et al. 2023), and
$\\{x_{\star},y_{\star}\\}=\\{R_{\star},0\\}$ in our simulation coordinates),
the planet was at an orbital phase of $\varphi_{0}=0.1745$ after the preceding
transit midpoint. GJ 436 b has a nearly polar orbit with its orbital axis and
stellar spin misaligned by $\Psi=103.2^{\circ}$ (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.

Figure 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 colour bar),
during which the planet completes about 16 orbital revolutions. The grey
circles indicate the approximate times of mid-transit.

As presented in Section A, in the case of GJ 436 b, which has a highly
misaligned orbit, the time it takes for the planet to be at the same
substellar 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 GJ 436 b. 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 Figure 3. Note
that to compute the ram pressure, we use the relative velocity between the
stellar wind velocity $\@vec{u}$ and the Keplerian motion of the planet
$\@vec{v}_{K}$: $\Delta\@vec{u}=\@vec{u}-\@vec{v}_{K}$. For GJ 436 b, the
Keplerian velocity is
$\@vec{v}_{K}=118(\sin\Psi\hat{\theta}+\cos\Psi\hat{\varphi})\leavevmode\nobreak\
$km s-1, where we adopted a stellar mass of $\simeq 0.45\leavevmode\nobreak\
M_{\odot}$ and an orbital distance of $14.56\leavevmode\nobreak\
R_{\star}\simeq 0.028$ au (Maxted et al., 2022).

Figure 4: The stellar wind proton density, velocity, total pressure (sum of
ram, including planet’s own orbital motion, thermal and magnetic pressures)
and total magnetic field strength at the orbit of GJ 436 b. 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 Figure 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.

The range of wind densities at the orbit of the planet is $[0.12,5.0]\times
10^{3}\leavevmode\nobreak\ {\rm cm}^{-3}$ for Model I and $[2.8,81]\times
10^{3}\leavevmode\nobreak\ {\rm cm}^{-3}$ for Model II. Their time averages
over one stellar rotation period are $510$ and $6200\leavevmode\nobreak\ {\rm
cm}^{-3}$ for Models I and II, respectively (horizontal grey lines shown in
Figure 4). Similarly, for wind velocities, we found that the range of values
that the planet experience are $[30,770]\leavevmode\nobreak\ $km s-1 and
$[220,830]$ km s-1 for Models I and II, respectively, with averages of $510$
and $640$ km s-1. When modelling the transmission spectroscopic transit of GJ
436 b, Bourrier et al. (2016) derived densities of $2^{+2.2}_{-1.2}\times
10^{3}\leavevmode\nobreak\ {\rm cm}^{-3}$ and wind speeds of $85^{+6}_{-16}$
km 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 Figure 4 crosses at the
approximate mid-transit times at the same level of ion density ($2\times
10^{3}\leavevmode\nobreak\ {\rm cm}^{-3}$) and speed ($85$ km s-1) 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^{+1.3}_{-0.75}\times 10^{-15}\leavevmode\nobreak\
M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$ (Vidotto & Bourrier, 2017) better
agree with the values derived in Model I ($1.1\times
10^{-15}\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$), but
not with Model II ($2.5\times 10^{-14}\leavevmode\nobreak\
M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$). Our chosen values of
$S_{A,0}/B_{\star}$ correspond to non-thermal velocities $\xi\simeq 14$ and
$45$ km s-1 at the transition region (see Section 2). Recently, Boro Saikia et
al. (2023) reported $\xi$ ranging from $\sim 5$ to $\sim 30$ km s-1 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_{\star}=1.1\times 10^{7}\leavevmode\nobreak\ \textrm{erg
(cm}^{2}\textrm{\,s\,G)}^{-1}$, derived from solar wind models, might not be
applicable to the winds of (old) M dwarf stars.

The bottom right panel of Figure 4 shows the stellar wind magnetic field
strength at the orbit of GJ 436 b. 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^{2}/8\pi$). We also
compute the ram pressure of the stellar wind ($\rho\Delta u^{2}$). 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 Figure 4. We find average total pressures of $8.2\times 10^{-6}$
dyn/cm-2 and $59\times 10^{-6}$ dyn/cm-2 for Models I and II, respectively.
These values will be used in Section 4, when we compute the strength of star-
planet interactions and the size of the magnetosphere of the planet.

##  4 Stellar-wind–planet interactions

Because of its close proximity to the star, and its sub-Alfvénic motion, GJ
436 b can magnetically interact with its host star. Recently, Loyd et al.
(2023) studied the presence of star-planet interaction signatures in GJ 436 in
UV spectroscopic observations. While they did not find line fluxes modulated
with the orbital period of GJ 436 b (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 GJ 436, compared
to other M dwarfs. They suggested that this could be caused by the presence of
GJ 436 b, 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 $\sim 3\times 10^{24}$ erg s-1. Converting this number
to a bolometric power, they found an upper limit of $\lesssim 10^{26}$ erg s-1
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 now proceed to estimate 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 start by calculating the (unsigned)
Poynting flux of the stellar wind that impacts on the planet

$P=\frac{|(\Delta\@vec{u}\times\@vec{B})\times\@vec{B}|}{4\pi}=\frac{B^{2}\Delta
u\sin\theta}{4\pi}\,.$ (3)

In the previous expression, $B\sin\theta$ is the component of the stellar wind
magnetic field that is perpendicular to $\Delta\@vec{u}$, $\theta$ is the
angle between $\Delta\@vec{u}$ and $\@vec{B}$ – all these quantities are
computed at the position of the planet, which lies at an orbital distance
$a_{\rm orb}=14.56\leavevmode\nobreak\ R_{\star}$. We find that the average
Poynting fluxes computed at the orbit of GJ 436 b are $P\simeq 160$ and $650$
erg cm-2 s-1 (or $0.1$ to $0.6$ W m-2) 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, that can then power electron-cyclotron radio emission (e.g. Vidotto et
al., 2019). In the case of GJ 436 b, though, our wind models indicate that the
motion is sub-Alfvénic (c.f. Section 2) and we use the Poynting flux to
compute the energy released in the Alfvén wing scenario and the reconnection
scenario. 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.

###  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}=\Delta u/v_{A}<1$, the dissipated power is (Saur et al., 2013)

$\mathcal{P}_{\rm wg}\simeq P(2\pi
r_{M}^{2})\alpha^{2}M_{A}\sin\theta=\frac{\Delta{u}^{2}{B}\sin^{2}\theta}{2}r_{M}^{2}\alpha^{2}\sqrt{4\pi\rho}\,.$
(4)

Here, $\alpha$ 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 $\alpha=1$, thus deriving the maximum power
$\mathcal{P}_{\rm 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., sum of the magnetic, thermal and ram pressures, see bottom
left panel of Figure 4), we have (e.g., Vidotto et al., 2014)

$\frac{r_{M}}{R_{p}}\simeq f\left[\frac{(B_{p}/2)^{2}/(8\pi)}{p_{\rm
tot}(a_{\rm orb})}\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\simeq
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 436 b is currently poorly
constrained: Bourrier et al. (2016) derived an upper limit for the equatorial
field strength of $\sim 3$ G, consistent with the upper limit estimates of
$\sim 10$ G from Loyd et al. (2023). Assuming $B_{p}=2$ G leads to an average
magnetosphere size of $5.2\leavevmode\nobreak\ R_{p}$ and
$3.7\leavevmode\nobreak\ R_{p}$, for stellar wind Models I and II,
respectively, as seen on the top panel of Figure 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
$\mathcal{P}_{\rm wg}\simeq 1.2\times 10^{22}$ erg s-1 and $4.8\times 10^{22}$
erg s-1 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}\propto B_{p}^{1/3}$, we have that
$\mathcal{P}_{\rm 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}$, i.e., nearly a factor of 5, resulting in $\mathcal{P}_{\rm
wg}\simeq 5.6\times 10^{22}$ erg s-1 and $22\times 10^{22}$ erg s-1 for Models
I and II, respectively.

Figure 5: Strength of the stellar-wind planet interactions computed in this
work. Top: the 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:
the maximum dissipated power of Alfvén waves travelling towards the star (Saur
et al., 2013). Bottom: the 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 Figure 3. The grey arrows in the upper horizontal axes indicate approximate
times of mid-transits.

###  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$ (Equation 3) over the cross-
section of the planet’s magnetosphere (Lanza, 2009)

$\mathcal{P}_{\rm rec}=\gamma_{\rm rec}P(\pi R_{M}^{2})=\frac{1}{4}\gamma_{\rm
rec}r_{M}^{2}{B^{2}\Delta u\sin\theta}\,,$ (6)

where $\gamma_{\rm 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
$\gamma_{\rm rec}$ would differ at different points in the orbit, for
simplicity here we take its maximum value of $\gamma_{\rm rec}=1$, so that our
power estimates $\mathcal{P}_{\rm rec}$ are upper limits. The results of this
scenario are shown in the bottom panel of Figure 5, where we see that the
average powers released in the reconnection scenario are $\mathcal{P}_{\rm
rec}\simeq 8\times 10^{22}$ erg s-1 and $17\times 10^{22}$ erg s-1 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
$\mathcal{P}_{\rm rec}\propto B_{p}^{2/3}$, thus for a 20-G field, we find an
increase in dissipated power of nearly a factor of 5, reaching average values
of $\mathcal{P}_{\rm rec}\simeq 38\times 10^{22}$ erg s-1 and $80\times
10^{22}$ erg s-1 for Models I and II, respectively.

###  4.3 Interconnecting magnetic loop between star and planet

Here, 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 Section 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^{+5.9}_{-7.6}$ 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 paper, 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) that shows that, in March-June 2016, GJ 436
presented a large-scale field that largely resembles an aligned dipole, with a
maximum (absolute) field strength of $\sim 27\leavevmode\nobreak\ $G (Figure
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 GJ 436 b, a warm-Neptune orbiting at a
distance of $0.028\leavevmode\nobreak\ \rm{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_{\star}=1.1\times 10^{7}\leavevmode\nobreak\ \textrm{erg
(cm}^{2}\textrm{\,s\,G)}^{-1}$, 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 Figure 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\times 10^{-15}\leavevmode\nobreak\
M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$ and $2.5\times
10^{-14}\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$, 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
allow us to derive the conditions of the stellar wind at the orbit of GJ 436
b. 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, GJ 436 b has a nearly polar orbit. This means that during one orbit, the
planet probes regions of open field lines and closed field lines. Taking into
account the spin-orbit angle, we calculated the position of the planet during
consecutive orbits: during one stellar rotation ($\sim 44$ days) the planet
orbits about $\sim 16$ times, experiencing different space weather conditions
at each orbit (see Figures 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 the GJ 436 b. For example, the many peaks and troughs we see in the
local wind properties (Figure 4) are associated to the planet crossing the
stellar equatorial plane.

Regardless of the two wind models we choose, our 3D simulations indicate that
the motion of GJ 436 b 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 on GJ 436 b to be $P\simeq
160$ and $650$ erg cm-2 s-1 (or $0.1$ to $0.6$ W m-2) 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
release energy, Lanza 2009). For both scenarios, we assume that the planet is
magnetised with an equatorial dipolar magnetic field of 1G, leading to a
magnetospheric size of about $5.2\leavevmode\nobreak\ R_{p}$ and
$3.7\leavevmode\nobreak\ R_{p}$, for stellar wind Models I and II,
respectively. Our estimated maximum powers released through star-planet
interactions are of the order of $10^{22}$ to $10^{23}$ erg s-1 (approximately
$10^{-10}$ to $10^{-9}$ of the stellar bolometric luminosity), with the
largest values associated to the reconnection scenario. It is interesting to
note that the powers we derived are not too 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 (Figure 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 increase 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 increase from $\sim 10^{23}$ erg s-1 to
$10^{26}$ erg s-1, the latter of which is the maximum power derived by Loyd et
al. (2023) due to 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 ($\lesssim 10^{3\times{3/2}}\simeq 3\times 10^{4}$ G).
Because of this, unfortunately, we cannot derive any meaningful upper limit to
the planetary magnetic field.

According to the models studied here, GJ 436 b orbits mostly within the Alfven
surface, i.e., it has a sub-Alfvenic orbital motion, except for small portions
of the orbit around the magnetic equator, where the orbit is outside the
Alfven surface. When a planet moves from sub- to super-Alfvenic motion, it is
likely that different types of star-planet interactions can take place. An
immediate effect should be seen by the formation of a bow shock around the
planet, when it moves from sub- to super-Alfvenic motion (e.g. Vidotto et al.,
2011a). Cohen et al. (2014) studied the sub-/super-Alfvenic 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 field of their host stars can experience induction
heating also in their interiors. Additionally, it is also expected that
magnetic connectivity can switch between on and off states – this can be due
to variations along one planetary year (e.g., the planet moving in and out of
the Alfven surface), or variations caused to the Alfven surface because of 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 search for signatures of star-planet interactions.

We also explored 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). This scenario has been recognised as the one that releases the largest
star-planet interaction powers, compared to the other scenarions (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
$\sim 6$ G would reproduce the powers found by Loyd et al. (2023), agreeing
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
436 b 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 Figure 4), we found that our Model I, with the
lowest surface flux of Alfvén waves, better matches the spectroscopic transit
values. Our derived mass-loss rate of $1.1\times 10^{-15}\leavevmode\nobreak\
M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$ is compatible with the one derived
from spectroscopic transits of $1.2^{+1.3}_{-0.75}\times
10^{-15}\leavevmode\nobreak\ M_{\odot}\leavevmode\nobreak\ \rm{yr}^{-1}$
(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 velocity
$\xi$ 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
v_{\perp,0}^{2}\rangle^{1/2}/\sqrt{2}$. If this can be extended to solar-like
stars, this implies that by measuring $\xi$ 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 found an increase in $\xi$ with stellar rotation (or activity). In their
sample, they have about 5 stars with similar masses, radii and effective
temperature to GJ 436, namely: GJ 832, LTT 2050, GJ 3470, GJ 176, GJ 849.
Three of these stars (GJ 832, GJ 176, 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 $\xi$ ranging between $14.7$ and $17.8$ km s-1 for
these objects. These values are very similar to the non-thermal velocity of
$14$ km s-1 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 smaller 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 star $\sigma T_{\rm 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 values of wave
fluxes for slowly rotating M dwarfs. Based on these arguments and the
agreement of the local wind values we derive here with results from
Ly-$\alpha$ transit modelling (Bourrier et al., 2016; Vidotto & Bourrier,
2017), we favour a lower flux of Alfvén waves to be more appropriate to the
case of GJ 436. This might imply that the ‘canonical’ value of
$S_{A,0}/B_{\star}=1.1\times 10^{7}\leavevmode\nobreak\ \textrm{erg
(cm}^{2}\textrm{\,s\,G)}^{-1}$, derived from solar wind models, might not be
applicable to the winds of old M dwarf stars.

Our models indicate that the planet GJ 436 b 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 investigating 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 436 b, 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 Figure 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-$\alpha$ transit of GJ 436
b remains stable over time (Dos Santos et al., 2019).

For this stability to occur in transits observed over 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 GJ 436 has an activity
cycle with a period of about 7 to 8 years (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 released powers in star-planet interactions
(e.g., Klein et al., 2022). Ideally, observations that probe star-planet
interactions mediated by the stellar wind/magnetic field (e.g., planetary
radio emission, chromospheric hot-spots, interaction of evaporated atmospheres
with stellar winds) and the magnetic map used in stellar wind models should be
conducted contemporaneously.

###### Acknowledgements.

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

## References

  * Adams (2011) Adams, F. C. 2011, ApJ, 730, 27
  * Attia et al. (2021) Attia, O., Bourrier, V., Eggenberger, P., et al. 2021, A&A, 647, A40
  * Bellotti et al. (2023) Bellotti, S., Fares, R., Vidotto, A. A., et al. 2023, A&A, 676, A139
  * Boro Saikia et al. (2020) Boro Saikia, S., Jin, M., Johnstone, C. P., et al. 2020, A&A, 635, A178
  * Boro Saikia et al. (2023) Boro Saikia, S., Lueftinger, T., Airapetian, V. S., et al. 2023, ApJ, 950, 124
  * Bourrier et al. (2016) Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., Tanaka, Y. A., & Vidotto, A. A. 2016, A&A, 591, A121
  * Bourrier et al. (2018) Bourrier, V., Lovis, C., Beust, H., et al. 2018, Nature, 553, 477
  * Bourrier et al. (2022) Bourrier, V., Zapatero Osorio, M. R., Allart, R., et al. 2022, A&A, 663, A160
  * Carolan et al. (2021) Carolan, S., Vidotto, A. A., Villarreal D’Angelo, C., & Hazra, G. 2021, MNRAS, 500, 3382
  * Cauley et al. (2019) Cauley, P. W., Shkolnik, E. L., Llama, J., & Lanza, A. F. 2019, Nature Astronomy, 3, 1128
  * Cherenkov et al. (2017) Cherenkov, A., Bisikalo, D., Fossati, L., & Mostl, C. 2017, ApJ, 846, 31
  * Cohen et al. (2014) Cohen, O., Drake, J. J., Glocer, A., et al. 2014, ApJ, 790, 57
  * Cranmer (2017) Cranmer, S. R. 2017, ApJ, 840, 114
  * Cravens (2004) Cravens, T. E. 2004, Physics of Solar System Plasmas
  * De Pontieu et al. (2007) De Pontieu, B., McIntosh, S. W., Carlsson, M., et al. 2007, Science, 318, 1574
  * Donati & Landstreet (2009) Donati, J. & Landstreet, J. D. 2009, ARA&A, 47, 333
  * Donati & Brown (1997) Donati, J.-F. & Brown, S. F. 1997, A&A, 326, 1135
  * Dos Santos et al. (2019) Dos Santos, L. A., Ehrenreich, D., Bourrier, V., et al. 2019, A&A, 629, A47
  * Ehrenreich et al. (2015) Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459
  * Elekes & Saur (2023) Elekes, F. & Saur, J. 2023, A&A, 671, A133
  * Evensberget et al. (2023) Evensberget, D., Marsden, S. C., Carter, B. D., et al. 2023, MNRAS, 524, 2042
  * Fares et al. (2010) Fares, R., Donati, J., Moutou, C., et al. 2010, MNRAS, 406, 409
  * Farrell et al. (1999) Farrell, W. M., Desch, M. D., & Zarka, P. 1999, J. Geophys. Res., 104, 14025
  * Finley et al. (2018) Finley, A. J., Matt, S. P., & See, V. 2018, ApJ, 864, 125
  * Fischer & Saur (2019) Fischer, C. & Saur, J. 2019, ApJ, 872, 113
  * Folsom et al. (2020) Folsom, C. P., O Fionnagain, D., Fossati, L., et al. 2020, A&A, 633, A48
  * Fontenla et al. (2016) Fontenla, J. M., Linsky, J. L., Witbrod, J., et al. 2016, ApJ, 830, 154
  * Hazra et al. (2022) Hazra, G., Vidotto, A. A., Carolan, S., Villarreal D’Angelo, C., & Manchester, W. 2022, MNRAS, 509, 5858
  * Hollweg (1986) Hollweg, J. V. 1986, J. Geophys. Res., 91, 4111
  * Jardine & Cameron (2008) Jardine, M. & Cameron, A. C. 2008, A&A, 490, 843
  * Kavanagh et al. (2021) Kavanagh, R. D., Vidotto, A. A., Klein, B., et al. 2021, MNRAS, 504, 1511
  * Kavanagh et al. (2019) Kavanagh, R. D., Vidotto, A. A., O Fionnagain, D., et al. 2019, MNRAS, 485, 4529
  * Kavanagh et al. (2022) Kavanagh, R. D., Vidotto, A. A., Vedantham, H. K., et al. 2022, MNRAS, 514, 675
  * Khodachenko et al. (2019) Khodachenko, M. L., Shaikhislamov, I. F., Lammer, H., et al. 2019, ApJ, 885, 67
  * Kislyakova & Noack (2020) Kislyakova, K. & Noack, L. 2020, A&A, 636, L10
  * Kislyakova et al. (2018) Kislyakova, K. G., Fossati, L., Johnstone, C. P., et al. 2018, ApJ, 858, 105
  * Klein et al. (2022) Klein, B., Zicher, N., Kavanagh, R. D., et al. 2022, MNRAS, 512, 5067
  * Kulow et al. (2014) Kulow, J. R., France, K., Linsky, J., & Loyd, R. O. P. 2014, ApJ, 786, 132
  * Kumar & Fares (2023) Kumar, M. & Fares, R. 2023, MNRAS, 518, 3147
  * Lanotte et al. (2014) Lanotte, A. A., Gillon, M., Demory, B. O., et al. 2014, A&A, 572, A73
  * Lanza (2009) Lanza, A. F. 2009, A&A, 505, 339
  * Lanza (2013) Lanza, A. F. 2013, A&A, 557, A31
  * Lavie et al. (2017) Lavie, B., Ehrenreich, D., Bourrier, V., et al. 2017, A&A, 605, L7
  * Lothringer et al. (2018) Lothringer, J. D., Benneke, B., Crossfield, I. J. M., et al. 2018, AJ, 155, 66
  * Loyd et al. (2023) Loyd, R. O. P., Schneider, P. C., Jackman, J. A. G., et al. 2023, AJ, 165, 146
  * Mauas et al. (1997) Mauas, P. J. D., Falchi, A., Pasquini, L., & Pallavicini, R. 1997, A&A, 326, 249
  * Maxted et al. (2022) Maxted, P. F. L., Ehrenreich, D., Wilson, T. G., et al. 2022, MNRAS, 514, 77
  * McCann et al. (2019) McCann, J., Murray-Clay, R. A., Kratter, K., & Krumholz, M. R. 2019, ApJ, 873, 89
  * Mesquita et al. (2021) Mesquita, A. L., Rodgers-Lee, D., & Vidotto, A. A. 2021, MNRAS, 505, 1817
  * Mesquita & Vidotto (2020) Mesquita, A. L. & Vidotto, A. A. 2020, MNRAS, 494, 1297
  * Nicholson et al. (2016) Nicholson, B. A., Vidotto, A. A., Mengel, M., et al. 2016, MNRAS, 459, 1907
  * O Fionnagain et al. (2021) O Fionnagain, D., Vidotto, A. A., Petit, P., et al. 2021, MNRAS, 500, 3438
  * Oran et al. (2017) Oran, R., Landi, E., van der Holst, B., Sokolov, I. V., & Gombosi, T. I. 2017, ApJ, 845, 98
  * Oran et al. (2013) Oran, R., van der Holst, B., Landi, E., et al. 2013, ApJ, 778, 176
  * Pagano et al. (2009) Pagano, I., Lanza, A. F., Leto, G., et al. 2009, Earth Moon and Planets, 105, 373
  * Peacock et al. (2019) Peacock, S., Barman, T., Shkolnik, E. L., et al. 2019, ApJ, 886, 77
  * Pineda & Villadsen (2023) Pineda, J. S. & Villadsen, J. 2023, Nature Astronomy
  * Powell et al. (1999) Powell, K. G., Roe, P. L., Linde, T. J., Gombosi, T. I., & de Zeeuw, D. L. 1999, J. Chem. Phys., 154, 284
  * Preusse et al. (2006) Preusse, S., Kopp, A., Buechner, J., & Motschmann, U. 2006, A&A, 460, 317
  * Rodgers-Lee et al. (2023) Rodgers-Lee, D., Rimmer, P. B., Vidotto, A. A., et al. 2023, MNRAS, 521, 5880
  * Rumenskikh et al. (2022) Rumenskikh, M. S., Shaikhislamov, I. F., Khodachenko, M. L., et al. 2022, ApJ, 927, 238
  * Sakaue & Shibata (2021) Sakaue, T. & Shibata, K. 2021, ApJ, 919, 29
  * Saur et al. (2013) Saur, J., Grambusch, T., Duling, S., Neubauer, F. M., & Simon, S. 2013, A&A, 552, A119
  * See et al. (2019) See, V., Matt, S. P., Folsom, C. P., et al. 2019, ApJ, 876, 118
  * Shaikhislamov et al. (2018) Shaikhislamov, I. F., Khodachenko, M. L., Lammer, H., et al. 2018, MNRAS, 481, 5315
  * Shkolnik et al. (2008) Shkolnik, E., Bohlender, D. A., Walker, G. A. H., & Collier Cameron, A. 2008, ApJ, 676, 628
  * Shkolnik et al. (2005) Shkolnik, E., Walker, G. A. H., Bohlender, D. A., Gu, P.-G., & Kuerster, M. 2005, ApJ, 622, 1075
  * Shoda et al. (2018) Shoda, M., Yokoyama, T., & Suzuki, T. K. 2018, ApJ, 853, 190
  * Sokolov et al. (2013) Sokolov, I. V., van der Holst, B., Oran, R., et al. 2013, ApJ, 764, 23
  * Strugarek (2016) Strugarek, A. 2016, ApJ, 833, 140
  * Strugarek et al. (2022) Strugarek, A., Fares, R., Bourrier, V., et al. 2022, MNRAS, 512, 4556
  * Suzuki (2018) Suzuki, T. K. 2018, Publications of the Astronomical Society of Japan, 70, 34
  * Torres (2007) Torres, G. 2007, ApJ, 671, L65
  * Toth et al. (2012) Toth, G., van der Holst, B., Sokolov, I. V., et al. 2012, Journal of Computational Physics, 231, 870
  * Trifonov et al. (2018) Trifonov, T., Kürster, M., Zechmeister, M., et al. 2018, A&A, 609, A117
  * Turner et al. (2021) Turner, J. D., Zarka, P., Grießmeier, J.-M., et al. 2021, A&A, 645, A59
  * van der Holst et al. (2014) van der Holst, B., Sokolov, I. V., Meng, X., et al. 2014, ApJ, 782, 81
  * Vidotto (2020) Vidotto, A. A. 2020, IAU Symposium, 354, 259
  * Vidotto (2021) Vidotto, A. A. 2021, Living Reviews in Solar Physics, 18, 3
  * Vidotto & Bourrier (2017) Vidotto, A. A. & Bourrier, V. 2017, MNRAS, 470, 4026
  * Vidotto et al. (2012) Vidotto, A. A., Fares, R., Jardine, M., et al. 2012, MNRAS, 423, 3285
  * Vidotto et al. (2019) Vidotto, A. A., Feeney, N., & Groh, J. H. 2019, MNRAS, 488, 633
  * Vidotto et al. (2011a) Vidotto, A. A., Jardine, M., & Helling, C. 2011a, MNRAS, 414, 1573
  * Vidotto et al. (2014) Vidotto, A. A., Jardine, M., Morin, J., et al. 2014, MNRAS, 438, 1162
  * Vidotto et al. (2011b) Vidotto, A. A., Jardine, M., Opher, M., Donati, J. F., & Gombosi, T. I. 2011b, MNRAS, 412, 351
  * Vidotto et al. (2010) Vidotto, A. A., Opher, M., Jatenco-Pereira, V., & Gombosi, T. I. 2010, ApJ, 720, 1262
  * Villarreal D’Angelo et al. (2014) Villarreal D’Angelo, C., Schneiter, M., Costa, A., et al. 2014, MNRAS, 438, 1654
  * Villarreal D’Angelo et al. (2021) Villarreal D’Angelo, C., Vidotto, A. A., Esquivel, A., Hazra, G., & Youngblood, A. 2021, MNRAS, 501, 4383
  * Vourlidas et al. (2007) Vourlidas, A., Davis, C. J., Eyles, C. J., et al. 2007, ApJ, 668, L79

##  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 Figure 6.

Figure 6: 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 $\\{x_{\rm
pl},y_{\rm pl},z_{\rm pl}\\}$ in red contains the orbital plane of the planet
$x_{\rm pl}y_{\rm pl}$ and orbital spin axis along $z_{\rm pl}$. The
observer’s inertial reference frame $\\{x_{\rm obs},y_{\rm obs},z_{\rm
obs}\\}$ (in blue) is centred on the star, with the observer located at
$+x_{\rm obs}$, and the stellar rotation axis aligned along $z_{\rm obs}$. The
reference frame $\\{x_{\star},y_{\star},z_{\star}\\}$ (in black) co-rotates
with the star. The spin-orbit angle is $\Psi$.

We start from the orbital reference frame $\\{x_{\rm pl},y_{\rm pl},z_{\rm
pl}\\}$, such that the $x_{\rm pl}y_{\rm pl}$-plane contains the orbit of the
planet and the orbital axis is along $z_{\rm pl}$. The origin of this
reference frame is centred on the star. In this reference frame, the
(circular) planetary motion is simply described as

$\begin{bmatrix}x_{\rm pl}\\\ y_{\rm pl}\\\ z_{\rm
pl}\end{bmatrix}=\begin{bmatrix}R_{\rm orb}\cos\varphi_{\rm pl}\\\ R_{\rm
orb}\sin\varphi_{\rm pl}\\\ 0\end{bmatrix}\,,$ (7)

where the phase of the orbit is $\varphi_{\rm pl}=\Omega_{\rm
orb}t+\varphi_{0}$, with $t$ describing the time and $\Omega_{\rm
orb}=2\pi/P_{\rm orb}$ the orbital rotation rate. Here, we assume a non-null
initial phase $\varphi_{0}$, which we computed later in this section to find
the location of GJ 436 b 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_{\rm obs},y_{\rm
obs},z_{\rm obs}\\}$, whose origin is also centred on the star, and $z_{\rm
obs}$ is aligned with the stellar rotation. The observer is located at
positive $x_{\rm obs}$ and we assume that $x_{\rm obs}\parallel x_{\rm pl}$.
The angle between the orbital spin and the stellar spin axis is the spin-orbit
angle $\Psi$, thus, the angle between $z_{\rm obs}$ and $z_{\rm pl}$ is also
$\Psi$. In this reference frame, the orbital motion of the planet is described
as

$\begin{bmatrix}x_{\rm obs}\\\ y_{\rm obs}\\\ z_{\rm
obs}\end{bmatrix}=\begin{bmatrix}1&0&0\\\ 0&\cos\Psi&\sin\Psi\\\
0&-\sin\Psi&\cos\Psi\end{bmatrix}\begin{bmatrix}x_{\rm pl}\\\ y_{\rm pl}\\\
z_{\rm pl}\end{bmatrix}\,,$ (8)

where we performed an anticlockwise rotation around the axis $x_{\rm
obs}\equiv x_{\rm pl}$ by an angle $\Psi$. According to our reference frame,
the mid-transit occurs when the planet is at $\\{x_{\rm obs},y_{\rm
obs},z_{\rm obs}\\}=\\{R_{\rm orb},0,0\\}$.

The reference frame co-rotating with the star
$\\{x_{\star},y_{\star},z_{\star}\\}$ is also centred on the star, with the
stellar spin axis along $z_{\star}$, which is parallel to $z_{\rm obs}$. To
transform from the observer’s reference frame to the stellar co-rotating
frame, we perform a clockwise rotation around the axis $z_{\star}\equiv z_{\rm
obs}$ by an angle $\Omega_{\star}t$

$\begin{bmatrix}x_{\star}\\\ y_{\star}\\\
z_{\star}\end{bmatrix}=\begin{bmatrix}\cos(\Omega_{\star}t)&-\sin(\Omega_{\star}t)&0\\\
\sin(\Omega_{\star}t)&\cos(\Omega_{\star}t)&0\\\
0&0&1\end{bmatrix}\begin{bmatrix}x_{\rm obs}\\\ y_{\rm obs}\\\ z_{\rm
obs}\end{bmatrix}\,.$ (9)

Therefore, to derive the orbital path of the planet in the reference frame of
the star, we substitute (7) and (8) into (9), thus obtaining

$\frac{x_{\star}}{R_{\rm orb}}=\cos(\Omega_{\rm
orb}t+\varphi_{0})\cos(\Omega_{\star}t)-\sin(\Omega_{\rm
orb}t+\varphi_{0})\sin(\Omega_{\star}t)\cos\Psi$ (10) $\frac{y_{\star}}{R_{\rm
orb}}=\cos(\Omega_{\rm
orb}t+\varphi_{0})\sin(\Omega_{\star}t)+\sin(\Omega_{\rm
orb}t+\varphi_{0})\cos(\Omega_{\star}t)\cos\Psi$ (11) $\frac{z_{\star}}{R_{\rm
orb}}=-\sin(\Omega_{\rm orb}t+\varphi_{0})\sin\Psi\,.$ (12)

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

Note that, in the case of aligned systems, $\Psi=0$ and we have

$\displaystyle{x_{\star}}=R_{\rm orb}\cos(\Omega_{\rm
orb}t+\varphi_{0}+\Omega_{\star}t)$ (13) $\displaystyle{y_{\star}}=R_{\rm
orb}\sin(\Omega_{\rm orb}t+\varphi_{0}+\Omega_{\star}t)$ (14)
$\displaystyle{z_{\star}}=0\,.$ (15)

I.e., in aligned systems, the planet re-encounters the same stellar wind
property once every $2\pi/(\Omega_{\rm orb}+\Omega_{\star})$, 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 showed in Section 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).

To explore the interconnecting loop scenario, we run 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 enforce that the
stellar magnetic field lines are closed up to the planet’s orbit, we choose a
source surface of $\simeq 15R_{\star}$, i.e., beyond the planet’s orbit. Using
the surface magnetic map from Figure 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 Figure 7, where we see that the coronal magnetic field is on average
$7.5\times 10^{-3}$ G.

Figure 7: Top: magnetic field of the star at the orbital distance of the
planet, assuming a potential field extrapolation, with a source surface of
$\simeq 15R_{\star}$ (in this case, no wind model is used). Bottom: the
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 Figure 3. The grey arrows in
the upper horizontal axes indicate approximate times of mid-transits.

With this, we can then compute the power released in the interconnecting loop
model using (Lanza 2013)

$\mathcal{P}_{\rm sb}=(2\pi R_{p}^{2})f_{\rm
open}\left[\frac{B_{p}^{2}v_{K}}{4\pi}\right]\,,$ (16)

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_{\rm open}=1-\left(1-\frac{3\zeta^{1/3}}{2+\zeta}\right)^{1/2}\,$ (17)

with $\zeta=B/B_{p}$. Here, $B$ is the stellar magnetic field extrapolated out
to the orbit of GJ 436 b and $B_{p}$, 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
$B_{p}$ is that at such a height, which is smaller than the surface field.
Therefore, Equation (16) 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\times 10^{25}$ erg s-1. For a
magnetic field of $B_{p}=6$ G, the power released increases to $10\times
10^{25}$ erg s-1 (this value is similar to the maximum power derived by Loyd
et al. 2023). This is shown in the bottom panel of Figure 7.

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 Section 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}\lesssim 6$ G can reproduce the maximum power
observed in star-planet interactions in the GJ 436 system (Loyd et al. 2023).